Article Consistency across multi-omics layers in a drug-perturbed gut microbial community Sander Wuyts1,†, Renato Alves1,† , Maria Zimmermann-Kogadeeva1,† , Suguru Nishijima1,† , Sonja Blasche1,2, Marja Driessen1, Philipp E Geyer3 , Rajna Hercog1, Ece Kartal1 , Lisa Maier1 , Johannes B M€uller3, Sarela Garcia Santamarina1,‡ , Thomas Sebastian B Schmidt1 , Daniel C Sevin4 , Anja Telzerow1, Peter V Treit3 , Tobias Wenzel1,§ , Athanasios Typas1 , Kiran R Patil1,2 , Matthias Mann3,5 , Michael Kuhn1,* & Peer Bork1,6,7,8,** Abstract Multi-omics analyses are used in microbiome studies to under- stand molecular changes in microbial communities exposed to dif- ferent conditions. However, it is not always clear how much each omics data type contributes to our understanding and whether they are concordant with each other. Here, we map the molecular response of a synthetic community of 32 human gut bacteria to three non-antibiotic drugs by using five omics layers (16S rRNA gene profiling, metagenomics, metatranscriptomics, metaproteo- mics and metabolomics). We find that all the omics methods with species resolution are highly consistent in estimating relative spe- cies abundances. Furthermore, different omics methods comple- ment each other for capturing functional changes. For example, while nearly all the omics data types captured that the antipsy- chotic drug chlorpromazine selectively inhibits Bacteroidota repre- sentatives in the community, the metatranscriptome and metaproteome suggested that the drug induces stress responses related to protein quality control. Metabolomics revealed a decrease in oligosaccharide uptake, likely caused by Bacteroidota depletion. Our study highlights how multi-omics datasets can be utilized to reveal complex molecular responses to external pertur- bations in microbial communities. Keywords metabolomics; metagenomics; metaproteomics; metatranscriptomics; microbiology Subject Categories Microbiology, Virology & Host Pathogen Interaction; Proteomics DOI 10.15252/msb.202311525 | Received 9 January 2023 | Revised 4 July 2023 | Accepted 6 July 2023 | Published online 24 July 2023 Mol Syst Biol. (2023) 19: e11525 Introduction The human gut microbiota is a complex community of microorgan- isms, which is affected by endogenous and environmental factors such as host genotype, diet, drug treatment and disease status, and in turn, influences host health and disease progression (Kau et al, 2011; Cho & Blaser, 2012; Cani, 2018; Durack & Lynch, 2018; Schmidt et al, 2018; Lindell et al, 2022). Currently, insights into the structure and function of the microbiota community mainly come from 16S rRNA gene profiling and shotgun metagenomics. While 16S rRNA amplicon sequencing offers a cost-efficient way to assess bacterial abundance at a higher taxonomic level, whole-genome shotgun metagenomics resolves the abundance of species and strains, together with the functional potential they encode (Quince et al, 2017; Almeida et al, 2019; Pasolli et al, 2019). In addition, gene and protein expression and metabolite abundance in the com- munity can be quantified with metatranscriptomics (Bashiardes et al, 2016), metaproteomics (Zhang & Figeys, 2019) and metabolo- mics (Zierer et al, 2018; Han et al, 2021), respectively. Ultimately, the combination of these methods should enable the integration of 1 European Molecular Biology Laboratory, Heidelberg, Germany 2 Medical Research Council Toxicology Unit, Cambridge, UK 3 Department of Proteomics and Signal Transduction, Max Planck Institute of Biochemistry, Martinsried, Germany 4 Cellzome, GlaxoSmithKline R&D, Heidelberg, Germany 5 Proteomics Program, NNF Center for Protein Research, Faculty of Health Sciences, University of Copenhagen, Copenhagen, Denmark 6 Max Delbr€uck Centre for Molecular Medicine, Berlin, Germany 7 Yonsei Frontier Lab (YFL), Yonsei University, Seoul, South Korea 8 Department of Bioinformatics, Biocenter, University of W€urzburg, W€urzburg, Germany *Corresponding author. Tel: +49 6221 387 8361; E-mail: mkuhn@embl.de **Corresponding author. Tel: +49 6221 387 8361; E-mail: peer.bork@embl.org †These authors contributed equally to this work ‡Present address: MOSTMICRO Unit, Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa, Oeiras, Portugal §Present address: Institute for Biological and Medical Engineering, Schools of Engineering, Medicine and Biological Sciences, Pontificia Universidad Catolica de Chile, Santiago, Chile � 2023 The Authors. Published under the terms of the CC BY 4.0 license Molecular Systems Biology 19: e11525 | 2023 1 of 17 https://orcid.org/0000-0002-7212-0234 https://orcid.org/0000-0002-7212-0234 https://orcid.org/0000-0002-7212-0234 https://orcid.org/0000-0001-6021-1246 https://orcid.org/0000-0001-6021-1246 https://orcid.org/0000-0001-6021-1246 https://orcid.org/0000-0002-8444-9272 https://orcid.org/0000-0002-8444-9272 https://orcid.org/0000-0002-8444-9272 https://orcid.org/0000-0001-7980-4826 https://orcid.org/0000-0001-7980-4826 https://orcid.org/0000-0001-7980-4826 https://orcid.org/0000-0002-7720-455X https://orcid.org/0000-0002-7720-455X https://orcid.org/0000-0002-7720-455X https://orcid.org/0000-0002-6473-4762 https://orcid.org/0000-0002-6473-4762 https://orcid.org/0000-0002-6473-4762 https://orcid.org/0000-0002-0611-7918 https://orcid.org/0000-0002-0611-7918 https://orcid.org/0000-0002-0611-7918 https://orcid.org/0000-0001-8587-4177 https://orcid.org/0000-0001-8587-4177 https://orcid.org/0000-0001-8587-4177 https://orcid.org/0000-0003-1374-5311 https://orcid.org/0000-0003-1374-5311 https://orcid.org/0000-0003-1374-5311 https://orcid.org/0000-0002-5264-5877 https://orcid.org/0000-0002-5264-5877 https://orcid.org/0000-0002-5264-5877 https://orcid.org/0000-0001-8443-1315 https://orcid.org/0000-0001-8443-1315 https://orcid.org/0000-0001-8443-1315 https://orcid.org/0000-0002-0797-9018 https://orcid.org/0000-0002-0797-9018 https://orcid.org/0000-0002-0797-9018 https://orcid.org/0000-0002-6166-8640 https://orcid.org/0000-0002-6166-8640 https://orcid.org/0000-0002-6166-8640 https://orcid.org/0000-0003-1292-4799 https://orcid.org/0000-0003-1292-4799 https://orcid.org/0000-0003-1292-4799 https://orcid.org/0000-0002-2841-872X https://orcid.org/0000-0002-2841-872X https://orcid.org/0000-0002-2841-872X https://orcid.org/0000-0002-2627-833X https://orcid.org/0000-0002-2627-833X https://orcid.org/0000-0002-2627-833X the major molecular layers of the cell, resulting in a more complete picture of the microbiome (Jansson & Baker, 2016; Heintz-Buschart & Wilmes, 2018). Several studies have shown how a combination of two or more of these omics methods could lead to novel insights regarding the dynamics and inner workings of a microbial commu- nity (Heintz-Buschart et al, 2016; Lloyd-Price et al, 2017; Salazar et al, 2019; Taylor et al, 2020). While multi-omics measurements provide information across molecular layers, their comprehensive integration remains challenging. One challenge is the limited knowl- edge about the concordance of different measurements in complex in natura settings in the absence of ground truth. Another challenge in comparing and integrating multi-omics datasets is the difference in their dynamics in response to perturbations. Although metabolite changes occur on a time scale of seconds, transcriptional changes usually occur on a time scale of minutes, while protein abundance changes take the longest to respond to a perturbation (Gerosa & Sauer, 2011; Choi et al, 2020). Synthetic microbial communities have been increasingly used to obtain a better understanding of the dynamics and species–species interactions (Goldford et al, 2018; preprint: Cheng et al, 2021). Com- pared with a natural gut microbiota, these synthetic communities have lower complexity, higher controllability and reproducibility and a well-defined composition at the strain level, at the cost of being sim- plified representations of natural ecosystems (Roy et al, 2014; Aranda- D�ıaz et al, 2022; Weiss et al, 2022). Yet, they do offer advantages over single species studies, as single species’ behaviour can significantly differ in mono-culture compared with co-culture (D’hoe et al, 2018). The complex interactions between the gut microbiota and non- antibiotic drugs have been elucidated from large-scale human stud- ies and high-throughput laboratory experiments (Rizkallah et al, 2010; Forslund et al, 2015, 2021; Spanogiannopoulos et al, 2016; Wilson & Nicholson, 2017; Zimmermann et al, 2021). This relation- ship is bidirectional, as drugs can influence microbiome composi- tion (Jackson et al, 2018; Maier et al, 2018; Vich Vila et al, 2020; Vieira-Silva et al, 2020), while the gut microbiota can have an impact on a drug’s efficacy and toxicity by altering its chemical structure (Zimmermann et al, 2019a,b; Javdan et al, 2020; Kl€une- mann et al, 2021). The emerging knowledge on drug–microbiota interactions has the potential to influence the future of drug devel- opment and personalised medicine (Doestzada et al, 2018; Weersma et al, 2020; Maier et al, 2021; Zimmermann et al, 2021). We therefore set out to answer the following three questions: How do the different omics methods perform in capturing dynamic changes in microbial communities in response to perturbations? Can we identify the drug’s mechanism of action on the bacteria, or the bacteria’s defensive responses? In which time frame do drugs cause perturbations on the bacteria, as visible in genomes, transcripts, proteins and metabolites? To this end, we designed a controlled time-course experiment with a synthetic community of 32 human gut representatives (Tramontano et al, 2018) in response to three drugs from diverse indication areas: chlorpromazine (antipsy- chotic), metformin (antidiabetic) and niclosamide (anthelmintic), which were previously reported to impair growth of several gut bac- teria (Maier et al, 2018). We followed the response of the defined community to the three non-antibiotic drugs over 4 days on the structural and functional levels across multi-omics layers, based on 16S rRNA gene, metagenome, metatranscriptome, metaproteome and untargeted metabolome profiling. Results Establishment of a synthetic community for drug perturbations To investigate microbial community response to drug perturbations in a controlled system across five omics layers, we combined 32 human gut microbiome representatives (Tramontano et al, 2018) and exposed this community to three different non-antibiotic drugs (Fig 1A; Appen- dix Table S1). The complete experiment was performed twice (run A and run B) as biological replicates, starting from the initial community assembly step from single bacterial cultures. More specifically, seven slow-growing species (inoculated on day 1) were combined with 25 fast-growing species (inoculated on day 3) on day 5 to form a synthetic community (Fig 1A and B). In order to ensure stable community com- position, we performed three culture passages by growing the mixed culture for 48 h and transferring 1% of total volume to a fresh culture medium. Samples for 16S rRNA amplicon sequencing were taken immediately after combining the strains (Inoculummix) and after each passage (Transfers 1–3) to evaluate the stabilisation of the community (Fig 1A; top row). We found that in both runs of the experiment the community reached a stable composition with four highly abundant species after three transfers (relative abundance > 10% for Escherichia coli, Clostridium perfringens, Veillonella parvula and Bacteroides thetaiotaomicron, Fig EV1A). The Bray–Curtis dissimilarity index showed that both runs were highly similar after the third transfer (Fig EV1C). After stabilisation, in each run the community perturbation was performed in duplicate during exponential growth (i.e., 5 h after passaging, as determined by optical density [OD] measurements on the previous transfer; Fig EV1D) by addition of one of the following drugs: (i) 5 mM metformin, a type 2 diabetes drug, (ii) 20 lM chlor- promazine, an antipsychotic drug, or (iii) 20 lM niclosamide, an anthelmintic drug (Fig 1A), while DMSO was used as a control. These are in the range of the estimated colon concentrations, which are available for metformin (1.5 mM) and chlorpromazine (25 lM) (Maier et al, 2018). We chose a higher concentration for metformin based on the reported intestinal concentration and previous data on metformin concentrations sufficient to impair growth of gut micro- biota members in vitro (Bailey et al, 2008b; Maier et al, 2018). The communities were sampled right before the addition of the drugs and 15 min, 30 min, 1 h and 3 h following the drug perturbation (Fig 1A). These time points were chosen to elucidate the early response of the bacterial community to the drug treatment. After 43 h, an additional sample was taken, and the communities were transferred into a fresh culture medium containing the drugs at ini- tial concentrations. A final sample was taken 48 h after this passage (91 h after the initial drug addition). In general, high correlation was evident between technical replicates within the same omics dataset (Appendix Fig S1). Consistency of community composition across omics measurements We first evaluated similarities and differences between the omics measurements in their ability to estimate species abundance. For sequencing-based omics methods, we performed both na€ıve ana- lyses with commonly used computational pipelines that do not use the information about synthetic community composition (DADA2 2 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al for 16S rRNA amplicon sequencing (Callahan et al, 2016), mOTUS v2.5 for metagenomics and metatranscriptomics (Milanese et al, 2019)), and targeted analyses based on mapping to the 32 reference genomes of species comprising our community (Materials and Methods). Within each omics method, both computational approaches produced highly similar results (Appendix Fig S2). As the composition-na€ıve approach only yields genus-level resolution for 16S rRNA sequencing data (Knight et al, 2018), we used the refer- ence genome mapping approach that yields higher resolution for all methods for comparison of community composition across omics types. For consistency, the same methodology (reference genome mapping) was used for metagenomics and metatranscriptomics. For metaproteomics data, we estimated species abundance by summing protein intensities for all proteins assigned to each species and dividing these values by the total protein intensity in each sample, as suggested previously (Kleiner et al, 2017). We compared relative species abundances between all pairs of omics methods except for metabolomics, which by nature repre- sents total metabolite measurements in the community and does not allow to separate compounds by species. Based on the correlation analysis, we found the abundance estimates to be highly similar (minimum Spearman correlation coefficient q = 0.78). Congruence was more pronounced for highly abundant species (Fig 2A). Specifi- cally, metagenomics and metatranscriptomics were the most similar of all pairwise comparisons (q = 0.92). Further, 16S rRNA amplicon sequencing showed high similarity with metagenomics for species with relative abundances higher than 0.001% (q = 0.89). However, for several species with low relative abundances, 16S rRNA sequencing provided higher relative abundance estimates compared to metagenomics, while other species, detected by metagenomics, were not detected with 16S rRNA sequencing. For this observation, no clear taxon-specific or condition-specific effect was found (Fig EV2A–C), indicating that the differences at these low relative abundances are most likely a result of differences in sequencing depth per sample, as has been previously reported (Pereira-Marques et al, 2019; Durazzi et al, 2021). Although metaproteomics is not yet widely used for species abundance estimation, we found the corre- sponding estimates in good agreement with the other omics methods, but only for species with relative abundance above 1% (q = 0.78–0.84; 16 out of 29 species detected across all the samples). This indicates that metaproteomics is less sensitive than sequencing-based methodologies for species abundance estimation, as has also been observed for in natura metaproteomics studies (Zhang & Figeys, 2019). Our results show generally high consistency between omics data types in relative species abundance estimations, and underline that metaproteomics can, in principle, provide robust species abundance estimates, at least for synthetic microbial com- munities, albeit with lower sensitivity. Consistency of functional profiles across omics measurements For each protein-coding gene of each species, we can compare rela- tive abundances across the three molecular layers: gene (metage- nomics), transcript (metatranscriptomics), and protein (metaproteomics). We performed such pairwise comparisons both for individual genes across all species (Appendix Fig S3) and for genes grouped based on the KEGG orthology (Kanehisa et al, 2017; Fig 2B). The correlation between metagenomic and metaproteomic estimates of gene and protein abun- dances was moderate (q = 0.5 for KEGG grouped features and q = 0.48 for all non-zero genes and proteins). Metatranscriptomics and metapro- teomics were the most similar (q = 0.73 for KEGG orthologs and q = 0.60 for transcripts and proteins), followed by metagenomics and metatranscriptomics (q = 0.7 for KEGG orthologs and q = 0.61 for genes and transcripts). To systematically assess how much information on the func- tional level is captured by metagenomics, metatranscriptomics and metaproteomics for different species, we estimated gene and path- way coverage by calculating the proportion of genes or pathways Figure 1. Experimental design and species used in this study. A Schematic overview of the experimental design. B Species cladogram constructed by pruning the relevant species from the GTDB species cladogram (release 95). � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 3 of 17 Sander Wuyts et al Molecular Systems Biology Figure 2. Comparison of species and feature abundances and functional coverage across omics methods. A Scatter plots representing species abundance defined as relative abundance of corresponding omics measurement in each sample. Each dot represents single species abundance in one sample. q, Spearman correlation coefficient. B Scatter plots representing gene, transcript or protein abundance linked through KEGG orthology. Each dot represents a single KEGG ortholog in one sample. C Left: Genome coverage of each of the omics datasets across samples for each species, right: relative species abundance estimated by metagenomics. The fraction of coverage is defined as the number of genes to which at least one read was mapped (for metagenomics and metatranscriptomics), or the number of detected proteins for metaproteomics divided by the total number of genes in the corresponding genome. (Metagenomics n = 75 samples, metatranscriptomics n = 101 samples and metaproteomics n = 112 samples). D KEGG pathway coverage. For the metabolomics dataset, pathway coverage is defined as the number of unique pathway metabolites tentatively detected in at least one sample, divided by the total number of metabolites in the pathway. For metagenomics, metatranscriptomics and metaproteomics, KEGG orthologs are used instead of pathway metabolites. E Heatmap of Mantel correlations across omics methods and Spearman correlation between replicates within each omics method. 4 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al that were detected by each method (Fig 2C and D). We found that 18 out of 32 species had an almost complete coverage (> 90%) in metagenomics, indicating that for these species most of the genes were recovered in all samples measured in this experiment (Fig 2C; in total 101,559 out of 103,921 possible protein-coding genes were detected at least once in the metagenomics dataset). This was not the case for 14 low-abundant species, for which the average gene content coverage was < 20%. For metatranscriptomics, the coverage was generally lower than for metagenomics (91,094 out of 103,921 possible transcripts detected at least once). This is however expected as not all genes are expressed in any given condition. Metaproteomics coverage was found to be much lower than metage- nomics and metatranscriptomics (9,144 out of 103,921 predicted proteins). This may be due to the limited dynamic range: In contrast to mass-spectrometry-based measurements, sequencing-based methods include an amplification step that increases the amount of material and makes it possible to cover rare transcripts and genes. For Escherichia coli, the most abundant species in our synthetic community, the maximum coverage of proteins across all samples did not exceed 30% (1,428 proteins out of 4,978 [29%] predicted proteins compared to 4,978 genes out of 4,978 predicted genes [100%] for metagenomics and 4,962 transcripts out of 4,978 tran- scripts [99%] for metatranscriptomics). This result is lower than state-of-the-art single species proteomics experiments, where around ~ 62% (2,586 detected proteins out of 4,189 predicted pro- teins) of bacterial proteins are captured (Mateus et al, 2020), likely due to the increased sample complexity in the community context, the increased search space of proteins and the presence of highly similar sequences in homologous proteins (where peptides cannot be unambiguously mapped to one protein). Since metabolomics data reflect the total pools of metabolites in the sample and cannot be analysed at the species level, we assessed the coverage of metabolic pathways defined in the KEGG database and compared it to pathway coverages by other omics methods (Fig 2D). For our analysis, we used 1,142 detected ions tentatively annotated as 3,488 possible metabolites by matching their accurate masses against the HMDB database (Wishart et al, 2018). This approach generally provides only low confidence in individual annotations and is unable to distinguish between isomers, yet ensures very broad tentative metabolome coverage. For metabolic pathways annotated in bacterial genomes, we observed an average pathway coverage of 35% for metabolomics, as compared to 44% for metaproteomics and 86% for metatranscriptomics. Even though direct comparison of omics methods is challenging, we believe that the lower coverage for metabolomics has several explanations. First, we measured metabolites in supernatant samples, to capture the drug and its metabolites and the secreted metabolites that play important roles in microbial communities in cross-feeding and sig- nalling (Yu et al, 2022). A more in-depth study could also addition- ally use the cell pellet for metabolomics, for example, to detect bioaccumulation (Kl€unemann et al, 2021). This means that compo- nents of the rich medium masked part of the signal (e.g., amino acids, peptides and polysaccharides), and extracellular products of bacterial metabolism, especially produced by only one or few spe- cies, may therefore be too dilute to be detected. Second, only a sub- set of all metabolites present in the bacterial cell will be secreted outside of the cell. Third, to calculate metabolic pathway coverage, we assumed that each pathway consists of metabolites that are produced or consumed by metabolic enzymes annotated in bacterial genomes, which is likely an overestimation of pathway sizes, since presence of an enzyme-coding gene in the genome does not neces- sarily imply that this enzyme was expressed or that its reactants were present in our experimental conditions. To further compare the samples measured with different omics methods, we performed a Mantel test, which measures a correlation coefficient between sample similarity matrices calculated based on each omics data type individually (Fig 2E). For example, while it is not possible to directly compare matrices of species and protein abundances, it is possible to calculate sample similarity matrices for these two methods that can then be compared with each other. Notably, transcript abundance as measured by metatranscriptomics showed a high correlation (≥ 0.57) with sample distance matrices of all other omics measurements, underlining that this method cap- tured both species abundance and functional information in our experiment. Hierarchical clustering of Mantel correlation coefficients revealed two groups, which shared transcript abundance data from metatranscriptomics as a common member: one group with species abundance data (from metagenomics, metatranscriptomics and metaproteomics) and gene abundance (metagenomics); and the sec- ond group with protein abundance data (metaproteomics) and metabolite abundances. The emergence of these groups can be explained by the nature of the data used to calculate sample dis- tance matrices: species and gene abundances in one group, and functional feature abundances in the other group. Altogether, meta- transcriptomics was found to be the most universal and versatile readout, as it can both provide robust and sensitive estimates of spe- cies abundance, and at the same time reflects functional changes, which are in concordance with protein changes detected by metaproteomics. Chlorpromazine treatment strongly affects community composition After testing the technical consistency between omics measurements in a synthetic microbial community, we explored the impact of drug perturbations on the community composition and the respective responses at species, gene, transcript, protein and metabolite levels. For the control condition and all perturbations (chlorpromazine, metformin and niclosamide), similar dynamic changes in alpha diversity were observed over time. In general, the alpha diversity (inverse Simpson index) increased as the community grew over time after inoculation, however, this increase was lower for chlor- promazine compared with the other drugs and the control condition (Fig EV3A and B). We observed different community dynamics between runs A and B during the exponential phase: E. coli and C. perfringens were the most abundant species in all conditions in run A (Figs 3A and EV3C), while E. coli dominated community com- position during exponential phase in run B. However, community compositions became more similar between the runs at 43 h after drug treatment (Appendix Fig S4). These analyses revealed that the addition of metformin and niclosamide had negligible effects on the community composition, while chlorpromazine treatment shifted the community composition in both runs. To identify differentially abundant species after drug perturba- tion, we analysed the composition of microbiomes by comparing species abundances in drug-treated samples against control samples � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 5 of 17 Sander Wuyts et al Molecular Systems Biology estimated by each omics type (Fig 3B; ANCOM (Mandal et al, 2015)). This analysis revealed that most members of the Bacteroidota phylum (Odoribacter splanchnicus, Parabacteroides distasonis, Phocaei- cola vulgatus, Bacteroides fragilis, Bacteroides thetaiotaomicron and Bacteroides uniformis) were less abundant in chlorpromazine-treated samples. This reduction in Bacteroidota abundance was detected across all four omics methods capturing community composition, indi- cating that each of these methods is capable of detecting strong signals of species abundance change. In addition to Bacteroidota, Fusobac- terium nucleatum was found to be less abundant in chlorpromazine- treated samples. In contrast, the other two drugs did not cause major shifts in relative abundances: although ANCOM test identified signifi- cant changes of abundance of several species, their relative abundance was not changing more than two-fold (Fig 3B). In summary, we found a consistent and substantial depletion of species belonging to the phy- lum Bacteroidota upon chlorpromazine treatment. Multi-omics measurements capture functional response of the community to all three drugs As compositional shifts do not provide information on the mecha- nisms of response of each community member, we investigated these functional responses in more detail by performing differential analysis of metatranscriptomic, metaproteomic and metabolomic datasets after a normalisation step wherein taxonomic abundance effects were reduced (see ‘Gene, transcript and protein counting’ in the Materials and Methods section). The highest number of differen- tially abundant transcripts, proteins and metabolites were found in samples treated with chlorpromazine (adjusted P-value < 0.001 and absolute fold change > 4 compared with the control for metatran- scriptomics, adjusted P-value < 0.05 and absolute fold change > 1.5 for metaproteomics and metabolomics; Fig 4A), which is in line with our findings that chlorpromazine caused the largest disruption to bacterial community (Fig 3B). Transcriptional response to chlor- promazine is detected already after 15 min of treatment across spe- cies belonging to different phyla, suggesting that, although Bacteroidota show the strongest response, other species also adapt their gene expression. In order to evaluate similarities between functional responses across omics data types, we performed pathway enrichment analy- sis of differentially abundant features between drug treatment and controls across all time points using the KEGG pathway annotations (Fig 4B). In general, we detected less overlap between omics layers on the functional level compared to species abundance analysis, as no single pathway was statistically significant in the enrichment analysis of all three functional omics datasets. Across all conditions, five pathways were found to be significantly enriched upon drug treatment compared to the control condition in two omics data types, while 35 pathways were statistically significantly enriched in only one omics dataset. The largest number of significantly enriched pathways was found in chlorpromazine-treated samples for metatranscriptomics data. For our metformin-treated samples, we did not observe substan- tial effects of metformin neither on community composition nor on transcript or protein abundance in our study, at least at the concen- trations used. Only a small number of pathways were significantly overrepresented (pFDR < 0.001 for metatranscriptomics and pFDR < 0.05 for metaproteomics and metabolomics) within the set of up- and downregulated features (transcripts/proteins/metabo- lites) in metformin-treated samples (Fig 4B). Further inspection of putative metabolites involved in these pathways showed that their abundance also decreased upon addition of metformin in the non- bacterial control samples (Appendix Fig S5). This indicates that met- formin primarily interferes with the measurement of these putative metabolites, probably due to their chemical similarity, underlining the importance of including non-bacterial control samples to study drug response. However, we cannot exclude that metformin also interacts with lysine and arginine metabolism pathways in bacteria, as reported before (Forslund et al, 2015; Pryor et al, 2019). In the previous single-strain experiments, metformin at the same concen- tration was shown to have an effect on several of the tested species (Maier et al, 2018). We attribute this discrepancy to the possibility that these species show a different behaviour in a community set- ting compared to a single culture setting, as already shown for other species (D’hoe et al, 2018). Unfortunately, our datasets do not pro- vide any further hypotheses as to what the underlying cause of this protective community effect could be. For niclosamide-treated samples, 10 pathways were significantly enriched (pFDR < 0.001) among regulated transcripts, including amino acid and nitrogen metabolism. Transcripts of nitrogen metab- olism pathway upregulated in the early time points (15 min, 30 min, 1 h, 3 h) were annotated as NAD-specific glutamate dehy- drogenase (belonging to the Cluster of Orthologous Groups COG0334 from the EggNOG database present in B. thetaiotaomicron, P. vulgatus, B. fragilis), hydroxylamine reductase (COG1151 in C. perfringens, B. uniformis) and carbamate kinase (COG0549 in Eggerthella lenta) (Appendix Fig S6). Previously, NAD-specific gluta- mate dehydrogenase was found to be upregulated in response to nitrogen availability in Mycobacterium smegmatis, where it is assumed to have a de-aminating activity (Harper et al, 2010). Fur- thermore, hydroxylamine reductase and carbamate kinase are enzymes belonging to the family of oxidoreductases which both act on nitrogenous compounds. Therefore, the upregulated pathway and its transcripts suggest increased nitrogen metabolism in niclosamide-treated samples. Further examination of our metabo- lomic dataset revealed that niclosamide gets degraded in both runs of the experiment (Fig EV4), which could explain the observed absence of perturbations of the community composition. Nitrore- ductases are known to detoxify niclosamide (Copp et al, 2020). While members of the nitroreductases family (COG0778) are expressed, we did not observe significant changes in their expres- sion levels upon treatment with niclosamide. Additional follow-up experiments are needed to elucidate the mechanisms underlying the microbial degradation of niclosamide and the roles of individual community members. Chlorpromazine induces stress response and metabolic changes in the community Since the number of differentially abundant features and pathways was high in chlorpromazine-treated samples (Fig 4A and B), we tested whether there are features that change concordantly across omics layers. We first compared transcript and protein fold changes upon perturbation, which revealed general agreement between rela- tive changes in gene expression and protein abundance, with tran- script fold changes at each time point correlating more strongly with 6 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al protein changes at later time points (Fig 4C; Appendix Fig S7), likely reflecting the delay between transcription and translation processes. Based on this analysis, we assessed the most prominent and concordant changes between metatranscriptomics and metapro- teomics 15 min and 1 h after chlorpromazine addition, respectively (Fig 4D). The most concordantly downregulated features were pro- teins and genes of Bacteroidota species and F. nucleatum, including ribosomal proteins, elongation factors, and central carbon metabo- lism enzymes gldA (glycerol dehydrogenase), gapdh (glyceralde- hyde 3-phosphate dehydrogenase), and pta (phosphate acetyltransferase), the latter two being downregulated in several species (Fig 4D). Furthermore, the most upregulated features found both in metatranscriptomics and metaproteomics were stress response genes in E. coli, such as the small heat shock proteins IbpA Figure 3. Changes in community composition upon drug perturbation. A Relative species abundance changes over time in the three drug conditions and control. Time 0 indicates timepoint of the drug addition 5 h after the passage in the fresh medium. Relative abundance measured from metagenomics data. B Left, distribution of relative species abundance for each species across all samples (all conditions and timepoints). Right, heatmap of species abundance fold changes measured by different omics methods for each drug condition versus control. Significance of changes estimated by the ANCOM test is indicated by asterisks: *changes detected at 0.7 threshold of W statistic; **changes detected at 0.8 threshold; ***changes detected at 0.9 threshold; nd, not detected. � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 7 of 17 Sander Wuyts et al Molecular Systems Biology and IbpB (Inclusion body-associated proteins A and B), other chap- erones, and ABC transporters. IbpA and IbpB serve as a first line of defence against protein aggregation (Miwa et al, 2021). In addition to ibpA and ibpB, we found upregulation of the transcriptional regu- lator rpoH and the chaperones dnaK and groEL, which are also involved in heat shock response (Yura, 2019; Appendix Fig S8). Figure 4. Functional analysis of transcript, protein and metabolite response after niclosamide, metformin or chlorpromazine treatment. A Number of differentially abundant transcripts, proteins and metabolites. B Pathway enrichment analysis across all conditions and time points. P-values are indicated by asterisks: *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001. C Heatmap representing Spearman correlation between fold changes (relative to control) detected by metatranscriptomics and metaproteomics across all drug perturbations. D Scatterplot depicting protein fold changes (relative to control) detected after 15 min of chlorpromazine exposure by metatranscriptomics versus after 1 h of exposure by metaproteomics. E COG enrichment analysis differentiating between species susceptible to chlorpromazine treatment (Bacteroidota) and non-susceptible species (non-Bacteroidota). COGs that are enriched in upregulated genes are coloured in red, while COGs that are enriched in downregulated genes are coloured in blue. Only COGs that were found to be significantly enriched in at least three out of four early time points are shown. COG names that are coloured are discussed in more detail in the main text. F Di-, tri- and tetra-saccharide abundances as measured by untargeted metabolomics (tentative metabolite annotation is based on m/z values indicated in the panel titles). The lines are coloured according to the experimental conditions (chlorpromazine, metformin, niclosamide and control), and the line type represents whether these are community culture or non-bacterial controls. 8 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al Together, these results show that chlorpromazine causes the activa- tion of a stress response in E. coli, probably due to induction of pro- tein aggregation either directly or indirectly. We then tested whether genes associated with stress response were differently regulated between chlorpromazine-susceptible and non-susceptible species. Two COGs related to the stress response were enriched in upregulated genes in at least two of the four early time points in the depleted (susceptible) species (Fig 4E, annotated in green, Fig EV5). One of them, COG0265, is upregulated by both susceptible and non-susceptible species and encompasses serine proteases (e.g., HtrA proteins such as DegP and DegQ), which repre- sent an important class of chaperones and heat-shock-induced ser- ine proteases, protecting periplasmic proteins. Furthermore, two COGs enriched in upregulated genes were related to (multidrug) transporter activity. COG1538, which contains genes annotated as membrane protein OprM, was the only COG enriched in upregulated genes by Bacteroidota in all four early time points (Fig 4F, annotated in purple). In Pseudomonas aeruginosa, OprM is part of MexAB– OprM, a multidrug efflux pump of the resistance-nodulation-cell division (RND) superfamily, where it plays a central role in multi- drug resistance by transporting drugs from the cytoplasm across the inner and outer membranes outside the cell envelope (Alekshun & Levy, 2007; Tsutsumi et al, 2019). RND-efflux pumps are found in a number of Gram-negative bacteria, for example, AcrAB–TolC is found in E. coli (Du et al, 2018) while Bacteroides fragilis harbours multiple copies of RND pumps BmeABC (Ghotaslou et al, 2018). Further, in addition to COG1538 (OprM homologues), also COG0841 containing homologues of the MexB/AcrB/BmeB protein (Fig 4E, also annotated in purple) was found to be enriched in upregulated genes, both in Bacteroidota and non-Bacteroidota species. These observations suggest an important role of the AcrAB–TolC/MexAB- OprM/BmeABC efflux pumps in determining chlorpromazine sus- ceptibility. Indeed, a recent study showed that chlorpromazine is both a substrate and an inhibitor of the AcrB multidrug efflux pump in Salmonella enterica and E. coli (Grimsey et al, 2020). Together, our results suggest that chlorpromazine could also be an inhibitor of BmeB, the AcrB/MexB homologue in Bacteroidota species and that this, potentially in combination with protein aggregation, could be one of the reasons explaining why Bacteroidota are affected by chlorpromazine treatment. Finally, depletion of Bacteroidota and downregulation of their genes involved in saccharide uptake might explain the enrichment of ‘Starch and sucrose’ and ‘Fructose and mannose metabolism’ path- ways among metabolites more abundant upon chlorpromazine treat- ment compared to the control samples (Fig 4B). As Bacteroides species are known to be capable of metabolising a wide variety of polysaccharides (Schwalm & Groisman, 2017), we believe that the higher abundance of ions tentatively annotated as oligosaccharides after chlorpromazine treatment (Fig 4F) measured by metabolomics is a result of their reduced consumption by these species. At the same time, the abundance of Bacteroidota is decreased compared to control. With the present data, we are not able to quantify the rela- tive contributions of reduced Bacteroidota biomass and lower meta- bolic activity per cell to the observed differences in metabolite levels. Taken together, by integrating multi-omics measurements, we propose that a series of events happens upon treatment with chlor- promazine: (i) a stress response is induced across several bacterial species with overexpression of ibpA and ipbB chaperones being the most pronounced response in E. coli; (ii) this stress response involves upregulation of AcrB/BmeB type of RND pumps, which may be bound and blocked by chlorpromazine in a species-specific manner; (iii) Bacteroidota species are more susceptible to chlor- promazine and are quickly depleted from the community while also downregulating genes involved in saccharide uptake, both of which result in (iv) higher levels of oligosaccharides in the culture medium due to the reduced ability of the perturbed community to utilise them. Discussion In this study, we evaluated the impact of drug perturbations on a synthetic gut microbial community by analysing five different omics data types in a highly controlled in vitro experiment. In general, we found concordance between all omics data types regarding the esti- mation of community composition (taxonomic profiling). To our knowledge, this is the first study to systematically compare taxo- nomic profiles obtained by four omics data types and can thus serve as a baseline for integrating different data types in ‘in natura’ set- tings. Using the synthetic community, we could show a high correla- tion between metagenomics and metatranscriptomics (q = 0.92), similar to a previous study that used only these two omics methods (q = 0.81; Heintz-Buschart & Wilmes, 2018). The taxonomic profiles obtained from our metaproteomics dataset, which is increasingly used in microbiome studies (e.g. Kleiner et al, 2017; Kleikamp et al, 2021), showed correlations between q = 0.78 and q = 0.84 with all other omics for species with a relative abundance higher than 1%. Although the number of detected proteins and the detec- tion limits remain to be improved, we showed that species abun- dance estimates can be derived from metaproteomics in a relatively simple, defined microbial community. For a defined community as used here, 16S amplicon sequencing would be sufficient to capture species abundance. However, when closely related species are used or natural communities are cultured ex vivo, a method based on shotgun sequencing is necessary. If the available resources make it necessary to prioritise, then metatran- scriptomics sequencing would be preferred, as its estimation of spe- cies abundance is highly correlated with metagenomics, but it additionally allows for functional profiling and offers insights into alterations in gene expression. Metabolomics measurements offer a complementary readout of the community functions. In our experiment, metabolomics revealed differences in the degradation of the studied drugs and suggested the links between the abundance of oligosaccharide com- pounds in the culture medium and changes in Bacteroidota abun- dance in response to chlorpromazine treatment. However, since our flow injection approach enabled acquisition of only exact ion masses (MS1), our dataset contains many ambiguous ion annota- tions. Follow up studies collecting tandem mass spectra (MS1 and MS2) and using chemical standards are required to confirm the identities of the changing metabolites (Schymanski et al, 2014). Fur- thermore, our dataset contains only extracellular metabolite mea- surements, and thus only provides indirect information on intracellular functional changes in the bacterial community. Of the three drugs used for perturbation, only chlorpromazine caused a large disturbance in the community composition. Surprisingly, met- formin, which has been shown to alter the gut microbiome in � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 9 of 17 Sander Wuyts et al Molecular Systems Biology patients (Forslund et al, 2015; Wu et al, 2017), did not perturb the community in our study, even though our earlier study suggested that the growth of at least four different species is inhibited by met- formin at the concentration used in monocultures (F. nucleatum, B. longum, P. copri and P. merdae; Maier et al, 2018). This observa- tion hints at a protective effect from the community, although this protective effect is not caused by drug degradation, as metformin concentrations remained high during the course of experiment (Fig EV4), but could be due to other interactions between members of the bacterial community, for example, by altered growth rates (D’hoe et al, 2018). Similarly, niclosamide was expected to cause a depletion of most of the members of the synthetic community, except for E. coli and B. wadsworthia (Maier et al, 2018), which was not observed in this study, also pointing to community-related pro- tection effects. Our metatranscriptomic data revealed an upregula- tion of genes related to nitrogen metabolism, while niclosamide concentration decreased during incubation, which was not observed in the non-bacterial controls. Therefore, we believe that certain spe- cies are capable of degrading niclosamide, which ultimately protected the whole community against possible inhibitory effects of niclosamide treatment. For chlorpromazine, the observed depletion of Bacteroidota spe- cies was in concordance with single species experiments (Maier et al, 2018). The antibiotic activity of chlorpromazine was reported relatively soon after its first usage in the 1950s (Kristiansen & Verg- mann, 1986; Dinan & Cryan, 2018). Its antibiotic mechanism of action is described to be multifold and includes effects on the cell membrane, energy generation and interference with cell replication due to DNA intercalation in E. coli (Grimsey et al, 2020). In our study, several genes and proteins related to protein aggregation were found to be upregulated in the metatranscriptomic and meta- proteomic data in E. coli and other community members. One study already reported protein aggregation of bovine insulin after chlor- promazine treatment (Bhattacharyya & Das, 2001). However, it remains unclear whether chlorpromazine can cause protein aggrega- tion in microbes either directly or indirectly, a hypothesis that should be followed-up in future experiments. Finally, we identified upregulation of RND-type efflux pumps in the Gram-negative bacteria, even in the Bacteroidota species that were severely depleted. It was recently shown that in S. enterica and E. coli, chlorpromazine is both a substrate and an inhibitor of AcrB, the inner membrane transporter of the tripartite system AcrAB–TolC, which is an RND-type efflux pump (Bailey et al, 2008a; Grimsey et al, 2020). Based on our data, we hypothesise that BmeB, the AcrB homologue in Bacteroidota, is also susceptible to chlorpromazine inhibition as we found upregulation of this and related genes, similar to what has been described by others in single species experiments (Grimsey et al, 2020). The suggested mechanism could be of signifi- cance in the battle against the rising multidrug resistance of Bacter- oides fragilis, a commensal bacterium that can act as a virulent pathogen when it escapes its normal niche (Wexler, 2007, 2012; Niestezpski et al, 2019). However, chlorpromazine’s antimicrobial activity generally occurs at concentrations higher than those clini- cally achievable (Grimsey & Piddock, 2019). Therefore, it is possible that, similarly as suggested for S. enterica, chlorpromazine could act as an antimicrobial adjuvant for Bacteroidota where its inhibition of RND-type efflux pumps prevents the extrusion of administered anti- biotics (Grimsey et al, 2020). From the perspective of human health, these results underline the detrimental effect of antipsychotics on the gut microbiome reported before (Dinan & Cryan, 2018). However, the revealed phylum-specific differences provide an opportunity to explore whether complementation of antipsychotic therapy with Bacteroidota-promoting dietary interventions could improve mental health and increase patients’ quality-of-life by restoring a healthy microbiota (Patnode et al, 2019). In conclusion, we directly compared data from multiple omics methods and showed that they agree on species abundance estima- tion of a defined and drug-perturbed microbial community in vitro. Those methods that are able to detect functional information also correlate with each other, albeit to a lower degree. We could also confirm expected time delays between transcriptional and trans- lational responses to perturbations, underlining that these methods reveal biological insights that happen at different time scales. While we were not able to detect the induction of metabolising enzymes in response to drug perturbation, we could detect the upregulation of other resistance mechanisms such as transporters just 15 min after the perturbation (Fig 4). Future studies could therefore investigate a broader panel of drugs with the timepoints established in this study. Although multi-omics analysis of natural communities is hampered by their increasing complexity, combining multiple omics measure- ments allows to measure the response of the community to perturba- tions across molecular layers and provides information that is not achievable by any method alone. Materials and Methods Species and drug selection The species used in this study represent a subset of abundant and prevalent species from the human gut. In total, 32 species were selected based on our previous work (Maier et al, 2018; Tramontano et al, 2018). The bacterial isolates were received from DSMZ, BEI Resources or ATCC and Dupont Health & Nutrition. The drugs were chosen because of their antimicrobial activity (Maier et al, 2018) and diversity in therapeutic usage. Reference genomes Reference genomes were downloaded from RefSeq in March 2019 (release 92) and reannotated using Prokka v1.14.0 (Seemann, 2014). Taxonomic classification was based on GTDB taxonomy release 95 (Parks et al, 2018) and inferred using GTDB-Tk v1.3.0 (Hyatt et al, 2010; Matsen et al, 2010; Price et al, 2010; Eddy, 2011; Ondov et al, 2016; Jain et al, 2018; Chaumeil et al, 2020). Further func- tional annotations (e.g., the KEGG orthology and eggNOG ortholo- gous group) were retrieved using eggNOG-mapper v2.0.1 which is based on eggNOG v5.0 (Huerta-Cepas et al, 2017). A cladogram was built by pruning the species cladogram from GTDB (bac120.tree, release 95) using the ETE toolkit (Huerta-Cepas et al, 2016). Medium and drug preparation mGAM medium was prepared according to manufacturer’s instruc- tions (HyServe GmbH & Co.KG, Germany, produced by Nissui Phar- maceuticals) and all the single species were grown in this medium 10 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al except V. parvula (Todd-Hewitt Broth (Sigma-Aldrich) + 0.6% sodium lactate) and B. wadsworthia (mGAM + 60 mM sodium formate + 10 mM taurine). All media were placed in anaerobic chamber 1 day before use under anoxic conditions (Coy Laboratory Products Inc.) (2% H2, 12% CO2, rest N2). Chlorpromazine (TCI Chemicals) and niclosamide (Santa Cruz Biotechnology) were added from DMSO stock solution. Metformin (Sigma) was added as pow- der directly into the medium after which the medium was filter- sterilised. Final concentrations of each drug were chosen based on previous work (Maier et al, 2018) with concentrations of 5 mM for metformin and 20 lM for chlorpromazine and niclosamide. The higher concentration for metformin is motivated by previously published data, which showed that a concentration of 20 lM was not sufficient to impair growth of gut microbiome members in vitro (Maier et al, 2018; Fig EV4). Experimental set-up and sample collection Species were pre-inoculated in isolation on liquid mGAM medium from pure stocks and incubated at 37°C under anaerobic condi- tions for a period of 3 or 5 days, depending on the growth rate of each species (see Fig 1). The monocultures were subsequently mixed in equal proportions based on their OD and then inoculated in 100 ml of mGAM liquid medium. To allow species to reach a stable state (stabilisation phase), the mixed culture was grown for 48 h after which 1 ml was transferred to fresh medium. In total, three passages were performed and after the second transfer OD measurements were taken to determine the start of the exponen- tial phase. Following the stabilisation phase, the mixed community was inoculated in medium prepared with one single drug or DMSO (con- trol) as soon as the community reached the exponential phase (OD roughly equal to 2–3). The cultures were subsequently sampled (3 mL) at fixed time intervals (0 min, 15 min, 30 min, 1 h, 3 h, 48 h), transferred to fresh medium (with drugs or DMSO) after 48 h and then sampled again 48 h later (or 96 h after the start of the experiment). The whole experiment was performed twice (labelled as run A and run B). 1.5 ml of each collected sample was centrifuged (30 s at max speed) after which the supernatant was removed and the cell pellet was stored at �80°C until further processing for DNA and RNA extractions. For protein and metabolite extraction, again 1 ml of each collected sample was centrifuged (30 s at max speed) and 450 ll of supernatant was used for metabolite extraction while the cell pellet was used for protein extraction (proteins in the cells). The remainder of the samples was frozen at �80°C as backup. DNA and RNA extraction Genomic DNA and total RNA were extracted from the same flash- frozen samples using Allprep Powerfecal DNA/RNA kit (Qiagen, Hilden Germany) following the manufacturer’s protocol but an addi- tional phenol–chloroform extraction step of 700 ll was performed after lysis. DNA yield was measured by using QubitTM dsDNA HS Assay Kit (Qubit, Waltham, Massachusetts, USA), split into two aliquots for ribosomal 16S rRNA amplicon sequencing and metage- nomic shotgun sequencing and was stored at �20°C. RNA yield was measured via Bioanalyzer (Agilent, Santa Clara, California, USA) with Pico and Nano chips depending on the sample concentration and stored at �80°C for further analysis. 16S rRNA amplicon, metagenomic and metatranscriptomic sequencing For 16S rRNA amplicon sequencing, extracted DNA was amplified using primers targeting the V4 region of the 16S rRNA gene on the F515 and R806 primer pair (Caporaso et al, 2011). PCR was performed according to the manufacturer’s instructions of the KAPA HiFi HotStart PCR Kits (Roche, Basel Switzerland) using barcoded primers and a two-step PCR protocol (NEXTflexTM 16S V4 Amplicon- Seq Kit, Bioo Scientific, Austin, Texas, USA). PCR products were pooled and purified using size-selective SPRIselect magnetic beads (0.8 left-sized, Beckman Coulter, Brea, CA, USA). The library was then diluted to 6 pM for sequencing. The library was sequenced on an Illumina (San Diego, USA) MiSeq platform using 2 × 250 bp paired-end reads at Genomics Core Facility (European Molecular Biology Laboratory [EMBL], Heidelberg, Germany). Metagenomic libraries for all samples were prepared using the NEB Ultra II and SPRI HD kits with a targeted insert size of 350, and sequenced on an Illumina HiSeq 4000 platform (Illumina, San Diego, CA, USA) in 2 × 150 bp paired-end with the aim of 1.5 Gbp average setup at the Genomics Core Facility (EMBL, Heidelberg, Germany). RNA samples were depleted for ribosomal RNA using the NEBNext Bacteria rRNA Depletion Kit (New England Biolabs, Ips- wich, Massachusetts, USA). Samples were pooled into a library using the NEBNext Ultra II Directional RNA Library Prep Kit (New England Biolabs) and subsequently sequenced on Illumina NextSeq500 platform (75 bp; single end) at Genomics Core Facility (EMBL, Heidelberg, Germany). Quality control of raw reads was performed using NGLess (Coelho et al, 2019). For metagenomics, reads were trimmed to the longest subread where each base had a Phred score of at least 25. For metatranscriptomics, a sliding window approach was used and reads were trimmed to the longest subread with an average Phred score of 20 (window size: 4 bp). Resulting reads shorter than 45 bp were discarded. To remove possible human contamination, all reads were mapped against a human reference database (release GRCh38.p10, Ensembl; Zerbino et al, 2018) using NGLess and samtools (Li et al, 2009). Reads with an identity threshold ≥ 90% were discarded. For metatranscriptomics specifically, rRNA reads were also removed from the dataset using SortMeRNA (Kopylova et al, 2012) with default parameters. Protein extraction Sample preparation, including protein extraction, digestion and pep- tide purification was performed according to the in-StageTip proto- col (Kulak et al, 2014, 20) with automation on an Agilent Bravo liquid handling platform according to (Geyer et al, 2016). In brief, samples were incubated in the PreOmics lysis buffer (P.O. 00001, PreOmics GmbH) for reduction of disulfide bridges, cysteine alkyl- ation and protein denaturation at 95°C for 10 min. Samples were sonicated using a Bioruptor Plus from Diagenode (15 cycles of 30 s). The protein concentration was measured using a tryptophan assay. In total, 200 lg protein of each organism were further � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 11 of 17 Sander Wuyts et al Molecular Systems Biology processed on the Agilent Bravo liquid handling platform by adding trypsin and LysC (1:100 ratio—lg of enzyme to lg of sample pro- tein). Digestion was performed at 37°C for 4 h. The peptides were purified in consecutive steps according to the PreOmics iST protocol (www.preomics.com). After elution from the solid phase extraction material, the peptides were completely dried using a SpeedVac centrifuge at 60°C (Eppendorf, Concentrator plus). Peptides were suspended in buffer A* (2% acetonitrile [v/v], 0.1% trifluoroacetic acid [v/v]) and sonicated for 30 min (Branson Ultrasonics, Ultrasonic Cleaner Model 2510). Metaproteomics Samples were analysed using a liquid chromatography (LC) system coupled to a mass spectrometer (MS). The LC was an EASY-nLC 1200 ultra-high pressure system (Thermo Fisher Scientific) and was coupled to a Q Exactive HFX Orbitrap MS (Thermo Fisher Scientific) using a nano-electrospray ion source (Thermo Fisher Scientific). Purified peptides were separated on 50 cm HPLC-columns (ID: 75 lm; in-house packed into the tip with ReproSil-Pur C18-AQ 1.9 lm resin [Dr. Maisch GmbH]). For each LC–MS/MS analysis about 500 ng peptides were separated on 100 min gradients. Peptides were separated with a two-buffer-system consisting of buffer A (0.1% [v/v] formic acid) and buffer B (0.1% [v/v] formic acid, 80% [v/v] acetonitrile). Peptides were eluted with a linear 70 min gradient of 2–24% buffer B, followed stepwise by a 21 min increase to 40% buffer B, a 4 min increase to 98% buffer B and a 5 min wash of 98% buffer B. The flow rate was constant at 350 nl/ min. The temperature of the column was kept at 60°C by an in- house-developed oven containing a Peltier element, and parameters were monitored in real time by the SprayQC software (Scheltema & Mann, 2012). First, data-dependent acquisition (DDA) was performed of each single organism to establish a library for the data independent acquisition (DIA) of the community culture samples. The DDA scans consisted of a Top15 MS/MS scan method. Target values for the full scan MS spectra were 3e6 charges in the 300–1,650 m/z range with a maximum injection time of 25 ms and a resolution of 60,000 at m/z 200. Fragmentation of precursor ions was performed by higher- energy C-trap dissociation (HCD) with a normalised collision energy of 27 eV. MS/MS scans were performed at a resolution of 15,000 at m/z 200 with an ion target value of 5e4 and a maximum injection time of 120 ms. Dynamic exclusion was set to 30 s to avoid repeated sequencing of identical peptides. MS data for the community culture samples were acquired with the DIA scan mode. Full MS scans were acquired in the range of m/z 300–1,650 at a resolution of 60,000 at m/z 200 and the auto- matic gain control (AGC) set to 3e6. The full MS scan was followed by 32 MS/MS windows per cycle in the range of m/z 300–1,650 at a resolution of 30,000 at m/z 200. A higher-energy collisional dis- sociation MS/MS scans was acquired with a stepped normalised collision energy of 25/27.5/30 eV and ions were accumulated to reach an AGC target value of 3e6 or for a maximum of 54 ms. The MS data of the single organisms and of the community cultures were used to generate a DDA-library and the direct-DIA- library, respectively, which were computationally merged into a hybrid library using the Spectronaut software (Biognosys AG). All searches were performed against a merged protein FASTA file of the reference genomes annotated using Prokka (see above). Searches used carbamidomethylation as fixed modification and acetylation of the protein N-terminus and oxidation of methionines as variable modifications. Trypsin/P proteolytic cleavage rule was used, permit- ting a maximum of 2 missed cleavages and a minimum peptide length of 7 amino acids. The Q-value cut-offs for both library gener- ation and DIA analyses were set to 0.01. Metabolomics measurements Untargeted metabolomics analysis of cell-free supernatants by flow injection-mass spectrometry was performed as described previously (Fuhrer et al, 2011). Briefly, samples were analysed on a LC/MS platform consisting of a Thermo Scientific Ultimate 3000 LC system with autosampler temperature set to 10°C coupled to a Thermo Sci- entific Q-Exactive Plus Fourier transform MS equipped with a heated electrospray ion source and operated in negative or positive ionisa- tion mode. The isocratic flow rate was 150 ll/min of mobile phase consisting of 60:40% (v/v) isopropanol:water buffered with 1 mM ammonium fluoride at pH 9 for negative ionisation mode or 60:40% (v/v) methanol:water buffered with 0.1% formic acid at pH 2 for positive ionisation mode, in both cases containing 10 nM tauro- cholic acid and 20 nM homotaurine as lock masses. Of note, the LC system was only used to transfer samples from the autosampler to the MS, but did not include a chromatographic column for analyte separation (‘flow injection’). Mass spectra were recorded in profile mode from 50 to 1,000 m/z with the following instrument settings: sheath gas, 35 a.u.; aux gas, 10 a.u.; aux gas heater, 200°C; sweep gas, 1 a.u.; spray voltage, �3 kV (negative mode) or 4 kV (positive mode); capillary temperature, 250°C; S-lens RF level, 50 a.u; resolu- tion, 70 k @ 200 m/z; AGC target, 3 × 106 ions, max. inject time, 120 ms; acquisition duration, 60 s. Spectral data processing includ- ing peak detection and alignment was performed using an auto- mated pipeline in R analogous to previously published pipelines (Fuhrer et al, 2011). To evaluate the impact of the drugs on mea- surements of other metabolites, we also prepared non-bacterial con- trols (i.e., each drug incubated in mGAM culture medium) and analysed them with the same procedure. Detected ions were tenta- tively annotated as metabolites based on accurate mass within a dynamic tolerance depending on local instrument resolving power ranging from 1 mDa at m/z = 50 to 5 mDa at m/z = 1,000 using the Human Metabolome Database (Wishart et al, 2018) as reference considering [M-H] and [M-2H] ions in negative mode or [M+], [M + H], [M + Na] and [M + K] ions in positive mode and up to two 12C to 13C substitutions. We additionally provide mappings to the Microbial Metabolites Database (Wishart et al, 2023) as part of the associated data repository. Of note, this approach precludes the resolution of isomers, of metabolites mapping to the same ion using different adduct assumptions, of unaccounted neutral gains or losses, or of metabolites with slightly distinct masses that neverthe- less map to the same ion within the respective local matching tolerance. Metabolomics data analysis Raw intensity values were quantile-normalised separately for ions acquired in positive and negative modes. For further analysis, the data from the two acquisition polarity modes were combined in one 12 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al http://www.preomics.com table and filtered as follows: only annotated ions were retained; ions annotated to 13C-compounds only were removed; for each metabo- lite, only the ion with the annotation considered most likely was retained (either the ion with the highest correlation with the total ion current, or the ion with the largest mean intensity across sam- ples). We provide both the filtered and the original unfiltered table with metabolite annotations in the associated data repository. Gene, transcript and protein counting Metagenomic and metatranscriptomic reads were mapped against a database of reference genomes containing only the species used in this study, using NGLess and samtools, with a minimum match size of 45 and minimum identity of 97%. Abundance estimates were produced by counting the number of reads mapping to each genome included in the study. If a read mapped to multiple genes, the count was distributed to each of the genes (e.g., if a read maps to gene X and gene Y, gene X and gene Y each get a count of 0.5). Proteins quantification and filtering. Proteins were filtered based on the information from the DDA experiment on which peptides are detected in which single species. Metaproteomics report with pro- tein and peptide quantification obtained from Spectronaut software applied to DIA samples was used as input. For each peptide in the community peptide report file, number of exact protein and species matches was calculated. For each protein, only unique peptides that match to one species were left for quantification. For each protein, the peptides were sorted according to the number of samples in which they were detected. Protein abundance was calculated as the mean of three most commonly measured peptides as suggested before (Ludwig et al, 2018). If the number of peptides was < 3, the protein was discarded. To reduce taxonomic abundance effects in downstream analyses, taxon-specific scaling was performed on metagenomics, metatran- scriptomics and metaproteomics as described by (Klingenberg & Meinicke, 2017). These measurements are all relative, and therefore changes in cell counts or biomass are not taken into account. Species abundance estimation Multiple computational strategies were used to estimate species abundance. Unless stated otherwise, for all analyses the species abundances resulting from read mapping were used. For this approach, first a database of 16S rRNA regions was constructed by manually querying the SILVA rRNA database (Quast et al, 2013) and extracting the representative sequence from each of our 32 spe- cies. Amplicon sequencing reads were then mapped against this database using MAPseq v1.2.4 (Matias Rodrigues et al, 2017). Paired reads were mapped independently and assignments were only considered upon agreement. Abundance estimates were then produced by counting the number of reads mapping to each genome included in the study. For metagenome derived estimates, total counts were normalised by the size of the genome (number of base- pairs). For metatranscriptome derived estimates, additional steps were required. Gene predictions by Prokka/Prodigal were used to calculate the total number of coding bases per genome, after exclu- sion of rRNA regions. Finally, total read counts were normalised by the number of coding bases on each genome. Species abundance was estimated from metaproteomic data by summing up all filtered protein intensities detected per each species, and dividing the sum by the total summed protein intensity in a given sample. In addition, to the approaches based on read mapping, several popular tools were used to estimate species abundance. For amplicon sequencing, DADA2 v1.10 (Callahan et al, 2016) was used with the GTDB database release 86 (Parks et al, 2018) for sequence classification which was limited to genus level classification. Meta- genomic and metatranscriptomic species abundances were esti- mated using mOTUs v2.5 (Milanese et al, 2019) and MetaPhlAn v3 (Beghini et al, 2021). Coverage analyses Gene, transcript and protein coverage were defined as the number of genes/transcripts/proteins that showed a count higher than 0, divided by the total number of predicted genes per species. For path- way coverage, the same approach was used, but genes/transcripts/ proteins were grouped by the KEGG pathways instead and thus divided by the number of KEGG orthologs in one single pathway. The same procedure was repeated for metabolites, but using the number of metabolites per pathway as predicted by KEGG instead of the number of KEGG orthologs. Mantel test Mantel tests were performed to compare each pair of omics data- sets and evaluate the similarity between them. Abundance tables of each omics were transformed into distance matrices using 1�Spearman’s correlation coefficient, and the matrices were com- pared using the mantel function in the vegan package (version 2.5.5) with the default option. For gene (metagenomic), transcript (metatranscriptomic) and protein (metaproteomic) level profiles, features with mean abundances below 1E-7, 1E-7 and 1E-5, respec- tively, were excluded, and only features above those thresholds were included in the analysis. All the features were included in the species-level profiles for each omics. Sixty-one samples that were common among all the omics datasets were used in this analysis. Differential species abundance analysis Differential analysis of species abundance across conditions was performed with ANCOM v. 2.1. Tables of species abundances calcu- lated from each omics measurements were preprocessed with feature_table_pre_process with sample names used as sample vari- ables, condition used as group variable, and parameters out_cut = 0.05; zero_cut = 0.90; lib_cut = 0; neg_lb = TRUE. The ANCOM function was applied to each pre-processed table with condition used as the main variable and time used as the formula for adjustment (with parameters: main_var = “condition’; p_adj_method = ‘BH’; alpha = 0.05, adj_formula = “time”; rand_formula = NULL). P-values were adjusted with Benjamini– Hochberg method (p_adj_method = ‘BH’). The cut-off of 0.7 for the W statistic was used to identify significantly differentially abundant species (detected_0.7 = TRUE). � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 13 of 17 Sander Wuyts et al Molecular Systems Biology Differential transcript, protein and metabolite abundance analysis Differential transcript analysis was performed using DESeq2 v1.26.0 (Love et al, 2014) after taxon-specific scaling (see above). The design formula included the factors run, drug, time point and the interaction term drug:timepoint. Statistical testing was performed with the Wald-test and IHW (Ignatiadis et al, 2016) to control the false discovery rate. Differential protein and metabolite analysis were performed using repeated measures analysis of variance using the lmer func- tion in the ade4 package. The same formula used in the differential transcript analysis was also used in the analysis. To exclude low- abundant features, those that have 0 or NA in at least half of the samples were removed prior to the analysis. P-values were adjusted by the IHW method. Fold changes of proteins and metabolites com- pared to those of controls were calculated based on raw values. Pathway and COG enrichment analysis Pathway enrichment was performed on differentially abundant fea- tures (cut-off for metatranscriptomics abs(log2(fold change)) > 2, pFDR < 0.001, cut-off for metabolomics and metaproteomics abs (log2(fold change)) > log2(1.5), pFDR < 0.05) with Fisher exact test using stats.fisher_exact in Python 3.7.7. P-values were adjusted with Benjamini–Hochberg procedure with multipletests function from statsmodels. For metabolomics, pathway enrichment analysis was performed for ion features rather than metabolite features (e.g. if one ion is annotated to two or more metabolites from the same pathway, it is counted only once in pathway enrichment analysis). For each feature, only one measurement corresponding to the maxi- mum absolute fold change over time was used for pathway enrich- ment analysis. COG enrichment was performed in the R environment using ClusterProfiler (Wu et al, 2021). Data availability The MS-based proteomics data have been deposited to the Proteo- meXchange Consortium via the PRIDE partner repository and are available via ProteomeXchange with identifier PXD036445. Metabo- lomic data has been submitted to MetaboLights under accession number MTBLS3129. Sequencing data is deposited at the European Nucleotide Archive (ENA): PRJEB46619. Preproccessed data files and tables are available on Figshare at https://doi.org/10.6084/m9. figshare.21667763. Code to generate all figures is available at https://github.com/grp-bork/multiomics_Wuyts_2022. Expanded View for this article is available online. Acknowledgements We acknowledge Vladimir Benes, Matthew Hayward, Melanie Tramontano, Thea Van Rossum, Camille Goemans, Carlos Voogdt and Michael Zimmermann for helpful discussions. We gratefully acknowledge support by the EMBL’s Genomics Core facility. The work was supported by the European Molecular Biology Laboratory and has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement number 668031 (to RA, SN, SW). PB: German Federal Ministry of Education and Research (LAMarCK, no. 031L0181A). KRP: UK Medical Research Council (project number MC_UU_00025/11). MZ-K: Postdoc Mobility Fellowship from the Swiss National Science Foundation (P400PB_186795) and a postdoctoral fellowship from the AXA Research Fund. LM, SGS and TW were supported by the EMBL Interdisciplinary Postdoc (EIPOD) program under Marie Sklodowska- Curie Actions COFUND (grant numbers 291772 and 664726). Open Access funding enabled and organized by Projekt DEAL. Author contributions Sander Wuyts: Data curation; software; formal analysis; visualization; methodology; writing – original draft; writing – review and editing. Renato Alves: Conceptualization; data curation; software; formal analysis; visualization; methodology; writing – original draft; writing – review and editing.Maria Zimmermann-Kogadeeva: Data curation; software; formal analysis; visualization; methodology; writing – original draft; writing – review and editing. Suguru Nishijima: Data curation; software; formal analysis; visualization; methodology; writing – original draft; writing – review and editing. Sonja Blasche: Investigation; writing – review and editing.Marja Driessen: Conceptualization; investigation. Philipp E Geyer: Investigation; writing – review and editing. Rajna Hercog: Investigation. Ece Kartal: Investigation; writing – review and editing. Lisa Maier: Conceptualization. Johannes B M€uller: Investigation. Sarela Garcia Santamarina: Investigation; writing – review and editing. Thomas Sebastian B Schmidt: Software; formal analysis; writing – review and editing. Daniel C Sevin: Formal analysis; investigation; writing – review and editing. Anja Telzerow: Investigation. Peter V Treit: Investigation; writing – review and editing. Tobias Wenzel: Investigation. Athanasios Typas: Conceptualization; supervision; writing – original draft; writing – review and editing. Kiran R Patil: Conceptualization; supervision; writing – original draft; writing – review and editing.Matthias Mann: Supervision.Michael Kuhn: Conceptualization; formal analysis; supervision; visualization; writing – original draft; writing – review and editing. Peer Bork: Conceptualization; supervision; funding acquisition; writing – original draft; writing – review and editing. Disclosure and competing interests statement The authors declare no competing interests. PB, AT and MM are members of the Editorial Advisory Board of Molecular Systems Biology. This has no bearing on the editorial consideration of this article for publication. References Alekshun MN, Levy SB (2007) Molecular mechanisms of antibacterial multidrug resistance. Cell 128: 1037 – 1050 Almeida A, Mitchell AL, Boland M, Forster SC, Gloor GB, Tarkowska A, Lawley TD, Finn RD (2019) A new genomic blueprint of the human gut microbiota. Nature 568: 499 – 504 Aranda-D�ıaz A, Ng KM, Thomsen T, Real-Ram�ırez I, Dahan D, Dittmar S, Gonzalez CG, Chavez T, Vasquez KS, Nguyen TH et al (2022) Establishment and characterization of stable, diverse, fecal-derived in vitro microbial communities that model the intestinal microbiota. Cell Host Microbe 30: 260 – 272.e5 Bailey AM, Paulsen IT, Piddock LJV (2008a) RamA confers multidrug resistance in Salmonella enterica via increased expression of acrB, which is inhibited by chlorpromazine. Antimicrob Agents Chemother 52: 3604 – 3611 Bailey CJ, Wilcock C, Scarpello JHB (2008b) Metformin and the intestine. Diabetologia 51: 1552 – 1553 14 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al https://doi.org/10.6084/m9.figshare.21667763 https://doi.org/10.6084/m9.figshare.21667763 https://github.com/grp-bork/multiomics_Wuyts_2022 https://doi.org/10.15252/msb.202311525 Bashiardes S, Zilberman-Schapira G, Elinav E (2016) Use of metatranscriptomics in microbiome research. Bioinform Biol Insights 10: 19– 25 Beghini F, McIver LJ, Blanco-M�ıguez A, Dubois L, Asnicar F, Maharjan S, Mailyan A, Manghi P, Scholz M, Thomas AM et al (2021) Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. eLife 10: e65088 Bhattacharyya J, Das KP (2001) Aggregation of insulin by chlorpromazine. Biochem Pharmacol 62: 1293 – 1297 Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP (2016) DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13: 581 – 583 Cani PD (2018) Human gut microbiome: hopes, threats and promises. Gut 67: 1716 – 1725 Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, Fierer N, Knight R (2011) Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci USA 108: 4516 – 4522 Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH (2020) GTDB-Tk: a toolkit to classify genomes with the genome taxonomy database. Bioinformatics 36: 1925 – 1927 Cheng AG, Aranda-D�ıaz A, Jain S, Yu F, Iakiviak M, Meng X, Weakley A, Patil A, Shiver AL, Deutschbauer A et al (2021) Systematic dissection of a complex gut bacterial community. bioRxiv https://doi.org/10.1101/2021.06.15.448618 [PREPRINT] Cho I, Blaser MJ (2012) The human microbiome: at the interface of health and disease. Nat Rev Genet 13: 260 – 270 Choi B, Cheng Y-Y, Cinar S, Ott W, Bennett MR, Josi�c K, Kim JK (2020) Bayesian inference of distributed time delay in transcriptional and translational regulation. Bioinformatics 36: 586 – 593 Coelho LP, Alves R, Monteiro P, Huerta-Cepas J, Freitas AT, Bork P (2019) NG- meta-profiler: fast processing of metagenomes using NGLess, a domain- specific language. Microbiome 7: 84 Copp JN, Pletzer D, Brown AS, Van der Heijden J, Miton CM, Edgar RJ, Rich MH, Little RF, Williams EM, Hancock REW et al (2020) Mechanistic understanding enables the rational design of salicylanilide combination therapies for gram-negative infections. MBio 11: e02068-20 D’hoe K, Vet S, Faust K, Moens F, Falony G, Gonze D, Llor�ens-Rico V, Gelens L, Danckaert J, De Vuyst L et al (2018) Integrated culturing, modeling and transcriptomics uncovers complex interactions and emergent behavior in a three-species synthetic gut community. Elife 7: e37090 Dinan TG, Cryan JF (2018) Schizophrenia and the microbiome: time to focus on the impact of antipsychotic treatment on the gut microbiota. World J Biol Psychiatry 19: 568 – 570 Doestzada M, Vila AV, Zhernakova A, Koonen DPY, Weersma RK, Touw DJ, Kuipers F, Wijmenga C, Fu J (2018) Pharmacomicrobiomics: a novel route towards personalized medicine? Protein Cell 9: 432 – 445 Du D, Wang-Kan X, Neuberger A, van Veen HW, Pos KM, Piddock LJV, Luisi BF (2018) Multidrug efflux pumps: structure, function and regulation. Nat Rev Microbiol 16: 523 – 539 Durack J, Lynch SV (2018) The gut microbiome: relationships with disease and opportunities for therapy. J Exp Med 216: 20 – 40 Durazzi F, Sala C, Castellani G, Manfreda G, Remondini D, De Cesare A (2021) Comparison between 16S rRNA and shotgun sequencing data for the taxonomic characterization of the gut microbiota. Sci Rep 11: 3030 Eddy SR (2011) Accelerated profile HMM searches. PLoS Comput Biol 7: e1002195 Forslund K, Hildebrand F, Nielsen T, Falony G, Le Chatelier E, Sunagawa S, Prifti E, Vieira-Silva S, Gudmundsdottir V, Krogh Pedersen H et al (2015) Disentangling type 2 diabetes and metformin treatment signatures in the human gut microbiota. Nature 528: 262 – 266 Forslund SK, Chakaroun R, Zimmermann-Kogadeeva M, Mark�o L, Aron-Wisnewsky J, Nielsen T, Moitinho-Silva L, Schmidt TSB, Falony G, Vieira-Silva S et al (2021) Combinatorial, additive and dose-dependent drug-microbiome associations. Nature 600: 500 – 505 Fuhrer T, Heer D, Begemann B, Zamboni N (2011) High-throughput, accurate mass metabolome profiling of cellular extracts by flow injection-time-of- flight mass spectrometry. Anal Chem 83: 7074 – 7080 Gerosa L, Sauer U (2011) Regulation and control of metabolic fluxes in microbes. Curr Opin Biotechnol 22: 566 – 575 Geyer PE, Kulak NA, Pichler G, Holdt LM, Teupser D, Mann M (2016) Plasma proteome profiling to assess human health and disease. Cell Syst 2: 185 – 195 Ghotaslou R, Yekani M, Memar MY (2018) The role of efflux pumps in Bacteroides fragilis resistance to antibiotics. Microbiol Res 210: 1 – 5 Goldford JE, Lu N, Baji�c D, Estrela S, Tikhonov M, Sanchez-Gorostiaga A, Segr�e D, Mehta P, Sanchez A (2018) Emergent simplicity in microbial community assembly. Science 361: 469 – 474 Grimsey EM, Piddock LJV (2019) Do phenothiazines possess antimicrobial and efflux inhibitory properties? FEMS Microbiol Rev 43: 577 – 590 Grimsey EM, Fais C, Marshall RL, Ricci V, Ciusa ML, Stone JW, Ivens A, Malloci G, Ruggerone P, Vargiu AV et al (2020) Chlorpromazine and amitriptyline are substrates and inhibitors of the AcrB multidrug efflux pump. mBio 11: e00465-20 Han S, Van Treuren W, Fischer CR, Merrill BD, DeFelice BC, Sanchez JM, Higginbottom SK, Guthrie L, Fall LA, Dodd D et al (2021) A metabolomics pipeline for the mechanistic interrogation of the gut microbiome. Nature 595: 415 – 420 Harper CJ, Hayward D, Kidd M, Wiid I, van Helden P (2010) Glutamate dehydrogenase and glutamine synthetase are regulated in response to nitrogen availability in Myocbacterium smegmatis. BMC Microbiol 10: 138 Heintz-Buschart A, Wilmes P (2018) Human gut microbiome: function matters. Trends Microbiol 26: 563 – 574 Heintz-Buschart A, May P, Laczny CC, Lebrun LA, Bellora C, Krishna A, Wampach L, Schneider JG, Hogan A, De Beaufort C et al (2016) Integrated multi-omics of the human gut microbiome in a case study of familial type 1 diabetes. Nat Microbiol 2: 16180 Huerta-Cepas J, Serra F, Bork P (2016) ETE 3: reconstruction, analysis, and visualization of phylogenomic data. Mol Biol Evol 33: 1635 – 1638 Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, Bork P (2017) Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol Biol Evol 34: 2115 – 2122 Hyatt D, Chen G-L, Locascio PF, Land ML, Larimer FW, Hauser LJ (2010) Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics 11: 119 Ignatiadis N, Klaus B, Zaugg J, Huber W (2016) Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nat Methods 13: 577 – 580 Jackson MA, Verdi S, Maxan M-E, Shin CM, Zierer J, Bowyer RCE, Martin T, Williams FMK, Menni C, Bell JT et al (2018) Gut microbiota associations with common diseases and prescription medications in a population- based cohort. Nat Commun 9: 2655 Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S (2018) High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun 9: 5114 Jansson JK, Baker ES (2016) A multi-omic future for microbiome studies. Nat Microbiol 1: 16049 � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 15 of 17 Sander Wuyts et al Molecular Systems Biology https://doi.org/10.1101/2021.06.15.448618 Javdan B, Lopez JG, Chankhamjon P, Lee Y-CJ, Hull R, Wu Q, Wang X, Chatterjee S, Donia MS (2020) Personalized mapping of drug metabolism by the human gut microbiome. Cell 181: 1661 – 1679.e22 Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K (2017) KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 45: D353 –D361 Kau AL, Ahern PP, Griffin NW, Goodman AL, Gordon JI (2011) Human nutrition, the gut microbiome and the immune system. Nature 474: 327 – 336 Kleikamp HBC, Pronk M, Tugui C, Guedes da Silva L, Abbas B, Lin YM, van Loosdrecht MCM, Pabst M (2021) Database-independent de novo metaproteomics of complex microbial communities. Cell Syst 12: 375 – 383.e5 Kleiner M, Thorson E, Sharp CE, Dong X, Liu D, Li C, Strous M (2017) Assessing species biomass contributions in microbial communities via metaproteomics. Nat Commun 8: 1558 Klingenberg H, Meinicke P (2017) How to normalize metatranscriptomic count data for differential expression analysis. PeerJ 2017: e3859 Kl€unemann M, Andrejev S, Blasche S, Mateus A, Phapale P, Devendran S, Vappiani J, Simon B, Scott TA, Kafkia E et al (2021) Bioaccumulation of therapeutic drugs by human gut bacteria. Nature 597: 533 – 538 Knight R, Vrbanac A, Taylor BC, Aksenov A, Callewaert C, Debelius J, Gonzalez A, Kosciolek T, McCall L-I, McDonald D et al (2018) Best practices for analysing microbiomes. Nat Rev Microbiol 16: 410 – 422 Kopylova E, No�e L, Touzet H (2012) SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28: 3211 – 3217 Kristiansen JE, Vergmann B (1986) The antibacterial effect of selected phenothiazines and thioxanthenes on slow-growing mycobacteria. Acta Pathol Microbiol Immunol Scand B 94: 393 – 398 Kulak NA, Pichler G, Paron I, Nagaraj N, Mann M (2014) Minimal, encapsulated proteomic-sample processing applied to copy-number estimation in eukaryotic cells. Nat Methods 11: 319 – 324 Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R (2009) The sequence alignment/map format and SAMtools. Bioinformatics 25: 2078 – 2079 Lindell AE, Zimmermann-Kogadeeva M, Patil KR (2022) Multimodal interactions of drugs, natural compounds and pollutants with the gut microbiota. Nat Rev Microbiol 20: 1 – 13 Lloyd-Price J, Mahurkar A, Rahnavard G, Crabtree J, Orvis J, Hall AB, Brady A, Creasy HH, McCracken C, Giglio MG et al (2017) Strains, functions and dynamics in the expanded human microbiome project. Nature 550: 61 – 66 Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550 Ludwig C, Gillet L, Rosenberger G, Amon S, Collins BC, Aebersold R (2018) Data-independent acquisition-based SWATH-MS for quantitative proteomics: a tutorial. Mol Syst Biol 14: e8126 Maier L, Pruteanu M, Kuhn M, Zeller G, Telzerow A, Anderson EE, Brochado AR, Fernandez KC, Dose H, Mori H et al (2018) Extensive impact of non- antibiotic drugs on human gut bacteria. Nature 555: 623 – 628 Maier L, Goemans CV, Wirbel J, Kuhn M, Eberl C, Pruteanu M, M€uller P, Garcia-Santamarina S, Cacace E, Zhang B et al (2021) Unravelling the collateral damage of antibiotics on gut bacteria. Nature 599: 120 – 124 Mandal S, Treuren WV, White RA, Eggesbø M, Knight R, Peddada SD (2015) Analysis of composition of microbiomes: a novel method for studying microbial composition. Microb Ecol Health Dis 26: 27663 Mateus A, Hevler J, Bobonis J, Kurzawa N, Shah M, Mitosch K, Goemans CV, Helm D, Stein F, Typas A et al (2020) The functional proteome landscape of Escherichia coli. Nature 588: 473 – 478 Matias Rodrigues JF, Schmidt TSB, Tackmann J, von Mering C (2017) MAPseq: highly efficient k-mer search with confidence estimates, for rRNA sequence analysis. Bioinformatics 33: 3808 – 3810 Matsen FA, Kodner RB, Armbrust EV (2010) Pplacer: linear time maximum- likelihood and Bayesian phylogenetic placement of sequences onto a fixed reference tree. BMC Bioinformatics 11: 538 Milanese A, Mende DR, Paoli L, Salazar G, Ruscheweyh H-J, Cuenca M, Hingamp P, Alves R, Costea PI, Coelho LP et al (2019) Microbial abundance, activity and population genomic profiling with mOTUs2. Nat Commun 10: 1014 Miwa T, Chadani Y, Taguchi H (2021) Escherichia coli small heat shock protein IbpA is an aggregation-sensor that self-regulates its own expression at posttranscriptional levels. Mol Microbiol 115: 142 – 156 Niestezpski S, Harnisz M, Korzeniewska E, Aguilera-Arreola MG, Contreras- Rodr�ıguez A, Filipkowska Z, Osi�nska A (2019) The emergence of antimicrobial resistance in environmental strains of the Bacteroides fragilis group. Environ Int 124: 408 – 419 Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, Phillippy AM (2016) Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol 17: 132 Parks DH, Chuvochina M, Waite DW, Rinke C, Skarshewski A, Chaumeil P- A, Hugenholtz P (2018) A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol 36: 1 – 35 Pasolli E, Asnicar F, Manara S, Quince C, Huttenhower C, Correspondence NS, Zolfo M, Karcher N, Armanini F, Beghini F et al (2019) Extensive unexplored human microbiome diversity revealed by over 150,000 genomes from metagenomes spanning age, geography, and lifestyle resource. Cell 176: 1 – 14 Patnode ML, Beller ZW, Han ND, Cheng J, Peters SL, Terrapon N, Henrissat B, Gall SL, Saulnier L, Hayashi DK et al (2019) Interspecies competition impacts targeted manipulation of human gut bacteria by fiber-derived glycans. Cell 179: 59 – 73.e13 Pereira-Marques J, Hout A, Ferreira RM, Weber M, Pinto-Ribeiro I, van Doorn L-J, Knetsch CW, Figueiredo C (2019) Impact of host DNA and sequencing depth on the taxonomic resolution of whole metagenome sequencing for microbiome analysis. Front Microbiol 10: 1277 Price MN, Dehal PS, Arkin AP (2010) FastTree 2 – approximately maximum- likelihood trees for large alignments. PLoS One 5: e9490 Pryor R, Norvaisas P, Marinos G, Best L, Thingholm LB, Quintaneiro LM, De Haes W, Esser D, Waschina S, Lujan C et al (2019) Host-microbe-drug- nutrient screen identifies bacterial effectors of metformin therapy. Cell 178: 1299 – 1312.e29 Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO (2013) The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41: 590 – 596 Quince C, Walker AW, Simpson JT, Loman NJ, Segata N (2017) Shotgun metagenomics, from sampling to analysis. Nat Biotechnol 35: 833 – 844 Rizkallah MR, Saad R, Aziz RK (2010) The human microbiome project, personalized medicine and the birth of pharmacomicrobiomics. Curr Pharmacogenomics Person Med 8: 182 – 193 Roy KD, Marzorati M, den Abbeele PV, de Wiele TV, Boon N (2014) Synthetic microbial ecosystems: an exciting tool to understand and apply microbial communities. Environ Microbiol 16: 1472 – 1481 Salazar G, Paoli L, Alberti A, Huerta-Cepas J, Ruscheweyh H-J, Cuenca M, Field CM, Coelho LP, Cruaud C, Engelen S et al (2019) Gene expression changes and community turnover differentially shape the Global Ocean Metatranscriptome. Cell 179: 1068 – 1083.e21 16 of 17 Molecular Systems Biology 19: e11525 | 2023 � 2023 The Authors Molecular Systems Biology Sander Wuyts et al Scheltema RA, Mann M (2012) SprayQc: a real-time LC-MS/MS quality monitoring system to maximize uptime using off the shelf components. J Proteome Res 11: 3458 – 3466 Schmidt TSB, Raes J, Bork P (2018) The human gut microbiome: from association to modulation. Cell 172: 1198 – 1215 Schwalm ND, Groisman EA (2017) Navigating the gut buffet: control of polysaccharide utilization in Bacteroides spp. Trends Microbiol 25: 1005 – 1015 Schymanski EL, Jeon J, Gulde R, Fenner K, Ruff M, Singer HP, Hollender J (2014) Identifying Small molecules via high resolution mass spectrometry: communicating confidence. Environ Sci Technol 48: 2097 – 2098 Seemann T (2014) Prokka: rapid prokaryotic genome annotation. Bioinformatics 30: 2068 – 2069 Spanogiannopoulos P, Bess EN, Carmody RN, Turnbaugh PJ (2016) The microbial pharmacists within us: a metagenomic view of xenobiotic metabolism. Nat Rev Microbiol 14: 273 – 287 Taylor BC, Lejzerowicz F, Poirel M, Shaffer JP, Jiang L, Aksenov A, Litwin N, Humphrey G, Martino C, Miller-Montgomery S et al (2020) Consumption of fermented foods is associated with systematic differences in the gut microbiome and metabolome. mSystems 5: e00901-19 Tramontano M, Andrejev S, Pruteanu M, Kl€unemann M, Kuhn M, Galardini M, Jouhten P, Zelezniak A, Zeller G, Bork P et al (2018) Nutritional preferences of human gut bacteria reveal their metabolic idiosyncrasies. Nat Microbiol 3: 514 – 522 Tsutsumi K, Yonehara R, Ishizaka-Ikeda E, Miyazaki N, Maeda S, Iwasaki K, Nakagawa A, Yamashita E (2019) Structures of the wild-type MexAB– OprM tripartite pump reveal its complex formation and drug efflux mechanism. Nat Commun 10: 1520 Vich Vila A, Collij V, Sanna S, Sinha T, Imhann F, Bourgonje AR, Mujagic Z, Jonkers DMAE, Masclee AAM, Fu J et al (2020) Impact of commonly used drugs on the composition and metabolic function of the gut microbiota. Nat Commun 11: 362 Vieira-Silva S, Falony G, Belda E, Nielsen T, Aron-Wisnewsky J, Chakaroun R, Forslund SK, Assmann K, Valles-Colomer M, Nguyen TTD et al (2020) Statin therapy is associated with lower prevalence of gut microbiota dysbiosis. Nature 581: 310 – 315 Weersma RK, Zhernakova A, Fu J (2020) Interaction between drugs and the gut microbiome. Gut 69: 1510 – 1519 Weiss AS, Burrichter AG, Durai Raj AC, von Strempel A, Meng C, Kleigrewe K, M€unch PC, Rössler L, Huber C, Eisenreich W et al (2022) In vitro interaction network of a synthetic gut bacterial community. ISME J 16: 1095 – 1109 Wexler HM (2007) Bacteroides: the good, the bad, and the nitty-gritty. Clin Microbiol Rev 20: 593 – 621 Wexler HM (2012) Pump it up: occurrence and regulation of multi-drug efflux pumps in Bacteroides fragilis. Anaerobe 18: 200 – 208 Wilson ID, Nicholson JK (2017) Gut microbiome interactions with drug metabolism, efficacy, and toxicity. Transl Res 179: 204 – 222 Wishart DS, Feunang YD, Marcu A, Guo AC, Liang K, V�azquez-Fresno R, Sajed T, Johnson D, Li C, Karu N et al (2018) HMDB 4.0: the human metabolome database for 2018. Nucleic Acids Res 46: D608 –D617 Wishart DS, Oler E, Peters H, Guo A, Girod S, Han S, Saha S, Lui VW, LeVatte M, Gautam V et al (2023) MiMeDB: the human microbial metabolome database. Nucleic Acids Res 51: D611 –D620 Wu H, Esteve E, Tremaroli V, Khan MT, Caesar R, Manner�as-Holm L, St�ahlman M, Olsson LM, Serino M, Planas-F�elix M et al (2017) Metformin alters the gut microbiome of individuals with treatment-naive type 2 diabetes, contributing to the therapeutic effects of the drug. Nat Med 23: 850 – 858 Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L et al (2021) clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb) 2: 100141 Yu JSL, Correia-Melo C, Zorrilla F, Herrera-Dominguez L, Wu MY, Hartl J, Campbell K, Blasche S, Kreidl M, Egger A-S et al (2022) Microbial communities form rich extracellular metabolomes that foster metabolic interactions and promote drug tolerance. Nat Microbiol 7: 542 – 555 Yura T (2019) Regulation of the heat shock response in Escherichia coli: history and perspectives. Genes Genet Syst 94: 103 – 108 Zerbino DR, Achuthan P, Akanni W, Amode MR, Barrell D, Bhai J, Billis K, Cummins C, Gall A, Gir�on CG et al (2018) Ensembl 2018. Nucleic Acids Res 46: D754 –D761 Zhang X, Figeys D (2019) Perspective and guidelines for Metaproteomics in microbiome studies. J Proteome Res 18: 2370 – 2380 Zierer J, Jackson MA, Kastenm€uller G, Mangino M, Long T, Telenti A, Mohney RP, Small KS, Bell JT, Steves CJ et al (2018) The fecal metabolome as a functional readout of the gut microbiome. Nat Genet 50: 790 – 795 Zimmermann M, Zimmermann-Kogadeeva M, Wegmann R, Goodman AL (2019a) Separating host and microbiome contributions to drug pharmacokinetics and toxicity. Science 363: eaat9931 Zimmermann M, Zimmermann-Kogadeeva M, Wegmann R, Goodman AL (2019b) Mapping human microbiome drug metabolism by gut bacteria and their genes. Nature 570: 462 – 467 Zimmermann M, Patil KR, Typas A, Maier L (2021) Towards a mechanistic understanding of reciprocal drug-microbiome interactions. Mol Syst Biol 17: e10116 License: This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. � 2023 The Authors Molecular Systems Biology 19: e11525 | 2023 17 of 17 Sander Wuyts et al Molecular Systems Biology http://creativecommons.org/licenses/by/4.0/ Abstract Introduction Results Establishment of a synthetic community for drug perturbations Consistency of community composition across omics measurements Consistency of functional profiles across omics measurements msb202311525-fig-0001 msb202311525-fig-0002 Chlorpromazine treatment strongly affects community composition �Multi-omics� measurements capture functional response of the community to all three�drugs Chlorpromazine induces stress response and metabolic changes in the community msb202311525-fig-0003 msb202311525-fig-0004 Discussion Materials and Methods Species and drug selection Reference genomes Medium and drug preparation Experimental �set-up and sample collection DNA and RNA extraction 16S rRNA amplicon, metagenomic and metatranscriptomic sequencing Protein extraction Metaproteomics Metabolomics measurements Metabolomics data analysis Gene, transcript and protein counting Species abundance estimation Coverage analyses Mantel�test Differential species abundance analysis Differential transcript, protein and metabolite abundance analysis Pathway and COG enrichment analysis Data availability msb202311525-supitem Acknowledgements Author contributions Disclosure and competing interests statement References msb202311525-bib-0001 msb202311525-bib-0002 msb202311525-bib-0003 msb202311525-bib-0004 msb202311525-bib-0005 msb202311525-bib-0006 msb202311525-bib-0007 msb202311525-bib-0008 msb202311525-bib-0009 msb202311525-bib-0010 msb202311525-bib-0011 msb202311525-bib-0012 msb202311525-bib-0013 msb202311525-bib-0014 msb202311525-bib-0015 msb202311525-bib-0016 msb202311525-bib-0017 msb202311525-bib-0018 msb202311525-bib-0019 msb202311525-bib-0020 msb202311525-bib-0021 msb202311525-bib-0022 msb202311525-bib-0023 msb202311525-bib-0024 msb202311525-bib-0025 msb202311525-bib-0026 msb202311525-bib-0027 msb202311525-bib-0028 msb202311525-bib-0029 msb202311525-bib-0030 msb202311525-bib-0031 msb202311525-bib-0032 msb202311525-bib-0033 msb202311525-bib-0034 msb202311525-bib-0035 msb202311525-bib-0036 msb202311525-bib-0037 msb202311525-bib-0038 msb202311525-bib-0039 msb202311525-bib-0040 msb202311525-bib-0041 msb202311525-bib-0042 msb202311525-bib-0043 msb202311525-bib-0044 msb202311525-bib-0045 msb202311525-bib-0046 msb202311525-bib-0047 msb202311525-bib-0048 msb202311525-bib-0049 msb202311525-bib-0050 msb202311525-bib-0051 msb202311525-bib-0052 msb202311525-bib-0053 msb202311525-bib-0054 msb202311525-bib-0055 msb202311525-bib-0056 msb202311525-bib-0057 msb202311525-bib-0058 msb202311525-bib-0059 msb202311525-bib-0060 msb202311525-bib-0061 msb202311525-bib-0062 msb202311525-bib-0063 msb202311525-bib-0064 msb202311525-bib-0065 msb202311525-bib-0066 msb202311525-bib-0067 msb202311525-bib-0068 msb202311525-bib-0069 msb202311525-bib-0070 msb202311525-bib-0071 msb202311525-bib-0072 msb202311525-bib-0073 msb202311525-bib-0074 msb202311525-bib-0075 msb202311525-bib-0076 msb202311525-bib-0077 msb202311525-bib-0078 msb202311525-bib-0079 msb202311525-bib-0080 msb202311525-bib-0081 msb202311525-bib-0082 msb202311525-bib-0083 msb202311525-bib-0084 msb202311525-bib-0085 msb202311525-bib-0086 msb202311525-bib-0087 msb202311525-bib-0088 msb202311525-bib-0089 msb202311525-bib-0090 msb202311525-bib-0091 msb202311525-bib-0092 msb202311525-bib-0093 msb202311525-bib-0094 msb202311525-bib-0095 msb202311525-bib-0096 msb202311525-bib-0097 msb202311525-bib-0098 msb202311525-bib-0099 msb202311525-bib-0100 msb202311525-bib-0101 msb202311525-bib-0102 msb202311525-bib-0103 msb202311525-bib-0104 msb202311525-bib-0105 msb202311525-bib-0106 msb202311525-bib-0107 msb202311525-bib-0108 msb202311525-bib-0109