Abstract

   Many drugs can perturb the gut microbiome, potentially leading to
   negative health consequences. However, mechanisms of most microbe-drug
   responses have not been elucidated at the genetic level. Using
   high-throughput bacterial transcriptomics, we systematically
   characterized the gene expression profiles of prevalent human gut
   bacteria exposed to the most frequently prescribed orally administered
   pharmaceuticals. Across >400 drug-microbe pairs, significant and
   reproducible transcriptional responses were observed, including
   pathways involved in multidrug resistance, metabolite transport,
   tartrate metabolism, and riboflavin biosynthesis. Importantly, we
   discovered that statin-mediated upregulation of the AcrAB-TolC efflux
   pump in Bacteroidales species enhances microbial sensitivity to vitamin
   A and secondary bile acids. Moreover, gut bacteria carrying acrAB-tolC
   genes are depleted in patients taking simvastatin, suggesting that
   drug-efflux interactions generate collateral toxicity that depletes
   pump-containing microbes from patient microbiomes. This work provides a
   resource to further understand the drivers of drug-mediated microbiota
   shifts for better informed clinical interventions.
     __________________________________________________________________

   Today, half of all Americans take at least one prescription drug, with
   national spending predicted to reach $400 billion by 2025^[34]1,[35]2.
   The prevalent use of pharmaceuticals is a major contributor to the
   alarming shifts in the gut microbiome, especially in industrialized
   countries^[36]3,[37]4. A recent screen of >1,000 orally administered
   drugs revealed a high frequency of antibacterial activity among
   human-targeted medications (24%), especially antipsychotics^[38]5.
   Clinically, proton pump inhibitors and atypical antipsychotics often
   trigger microbiota changes, with side-effects resembling those of
   antibiotic use (e.g., diarrhea, fungal infection)^[39]5,[40]6. Even
   though the top three prescribed drug classes (i.e., antihyperlipidemic
   agents, antidepressants, analgesics^[41]7) are all linked to gut
   microbiota disturbances, the mechanisms driving these shifts are poorly
   understood^[42]5. Physiologically, altered microbiome composition can
   negatively impact epithelial integrity^[43]8, gut motility^[44]9,
   nutrient availability^[45]10, immune homeostasis^[46]11,[47]12, and
   even treatment response. For instance, depletions in gut bacterial
   diversity have been linked to worse outcomes in immunotherapy drug
   trials^[48]13 and higher susceptibility to infection by pathogens such
   as Clostridiodes difficile^[49]14.

   Gut microbes interact with pharmaceutical compounds in a variety of
   ways. Bacteria can biotransform drugs to impact efficacy^[50]15,[51]16,
   bioavailability^[52]17, and toxicity^[53]18. For example, commensal
   Eggerthella lenta inactivates digoxin, an orally delivered cardiac
   glycoside, by expressing the cardiac glycoside reductase (cgr)
   operon^[54]15. Gut bacteria can also deplete xenobiotics from their
   local environment via bioaccumulation^[55]19 and mitigate the effects
   of toxic medications by employing conserved antibiotic resistance
   mechanisms^[56]5. Within minutes of xenobiotic exposure, microbes often
   exhibit highly specific transcriptional responses^[57]20,[58]21.
   Therefore, transcriptomic measurements can help quickly dissect
   specific drug-microbe responses (e.g., digoxin^[59]15,
   5-fluorouracil^[60]22, and levodopa^[61]23). However, current cost and
   scale of bacterial transcriptomics has precluded high-throughput
   studies using this approach^[62]24.

   Here, we describe the systematic transcriptomic analysis of 409
   drug-microbe pairs to dissect gene-level gut bacterial responses to top
   pharmaceuticals. We developed a high-throughput transcriptomics
   pipeline for optimization in non-model gut microorganisms that was then
   applied to generate a total of 978 individual drug-microbe samples
   (including replicates). By analyzing the transcriptomes of prevalent
   human gut bacteria exposed to top-prescribed orally administered drugs,
   we observed shared and strain-specific responses, including in pathways
   for drug degradation, vitamin biosynthesis, and multidrug efflux.
   Further analysis of a human microbiome cohort dataset confirmed the
   clinical significance of our findings. This study highlights the
   utility of large-scale transcriptomics for functional discovery of gut
   microbiota-xenobiotic interactions.

RESULTS

Transcriptomic map of microbial responses to top medicines

   Of the top 200 drugs prescribed today, 83% are orally delivered and
   expected to interact directly with the gut
   microbiome^[63]5,[64]25,[65]26. We therefore sought to measure the
   transcriptomic responses of prevalent gut bacteria exposed to
   top-prescribed oral pharmaceuticals. We implemented and optimized a
   multiplexed RNA-seq technique^[66]27 for high-throughput
   transcriptomics of non-model gut bacteria and incorporated
   cost-efficient ribosomal RNA (rRNA) depletion for diverse non-model gut
   bacteria that we previously developed^[67]24 ([68]Fig. 1a). This
   modified pipeline can generate high-quality transcriptomes for diverse
   gut bacterial phyla at a cost of <$16 for ~1M reads per sample
   ([69]Supplementary Table 1).

Figure 1. Top prescribed drugs elicit rich transcriptomic responses from
prevalent gut bacterial species.

   Figure 1.
   [70]Open in a new tab

   (a) Schematic of our high throughput transcriptomic pipeline for
   non-model gut microorganisms. (b) Heatmap color indicates delta in
   average optical density (OD) at 48 hours of bacterial strains grown in
   mGAM supplemented with different pharmaceuticals at 100 µM as compared
   to a same-volume solvent control (n=3 per condition). Drug condition is
   indicated on the x-axis and strain identity is indicated on the y-axis.
   Stars indicate >0.3 average absolute change in growth with p[adj]<0.05
   (*), calculated by two-sided independent t-test with Benjamini &
   Hochberg correction (c) Bacterial strains were exposed to drugs (n=2
   per condition) or vehicle controls (n=4 per condition) and
   differentially expressed genes (p[adj]<0.05, calculated as above) were
   identified. Heatmap color represents the log[10] value of the
   differentially expressed gene ratio (DEGR), defined as number of
   differentially expressed genes in a drug condition divided by total
   number of genes within a strain genome, with drug indicated on the
   x-axis (classes grouped by color) and strain identity indicated on the
   y-axis. Bar plot inserts show log[10] values of average DEGRs across
   strains (top bar insert) or drugs (left bar insert). Maximum likelihood
   phylogenetic tree is shown on the left with bacterial family identity
   indicated by color in the figure legend.

   We assembled a panel of 14 representative and highly
   prevalent^[71]28–[72]30 human gut isolates spanning the Bacteroidetes,
   Firmicutes, Actinobacteria, and Proteobacteria phyla (see
   [73]Supplementary Fig. 1, [74]Supplementary Tables 2 & [75]3) and 19
   oral drugs from the top prescribed drugs in 2017 according to the
   Agency for Healthcare Research and Quality ([76]Supplementary Table 4,
   [77]Methods)^[78]31. Our drug list included the top 8 most prescribed
   drugs as well as 10 additional neurotransmitter modulators, which were
   included due to established associations between the psychotherapeutic
   drug class and microbiota compositional changes^[79]5. Lenalidomide, a
   chemotherapeutic with the largest market cap in the small molecule drug
   class^[80]31, was also added.

   We first assessed the antimicrobial activity of the chosen drugs
   against our bacterial panel ([81]Fig. 1b, [82]Supplementary Fig. 2;
   [83]Methods). Drug concentrations were chosen to span median drug
   concentrations within the small and large intestine, which have been
   estimated to approach 20 and 100 µM, respectively^[84]5. We did not
   observe growth inhibition for most strains at 2 or 20 µM concentrations
   ([85]Supplementary Fig. 2). At 100 µM, growth inhibition was observed
   for at least one drug in 18 of 19 strains within 24 hours ([86]Fig.
   1b). Notably, B. longum did not exhibit growth defects in any drug
   condition. Broad spectrum antimicrobial activity was seen for
   simvastatin, amlodipine, and a subset of psychotherapeutics (i.e.,
   sertraline, paroxetine, duloxetine, fluoxetine), in agreement with
   previous reports^[87]5. Interestingly, we observed the greatest
   magnitudes of drug-induced growth defects within Bacteroidetes strains,
   suggesting that these species are more susceptible to drug toxicities
   ([88]Fig. 1b, [89]Supplementary Fig. 2). We further validated that the
   observed toxicity profiles extended to complex bacterial communities by
   exposing the drug panel to a fresh fecal sample from a healthy
   volunteer. We observed similar growth inhibition profiles of 13 panel
   strains grown in fecal community ([90]Supplementary Fig. 3), supporting
   the relevance of our individual strain-level growth measurements in the
   context of a complex consortium.

   Since substantial growth inhibition can confound transcriptional
   measurements, we performed all subsequent transcriptomic assays at 20
   µM, which reasonably approximated intestinal drug concentrations for
   detecting drug responses while minimizing drug toxicities^[91]5,[92]32.
   To further avoid impact of antimicrobial toxicity, cells were harvested
   at 1.5 hours post-exposure, which is shorter than the doubling time of
   several gut bacteria^[93]33. The transcriptomic pipeline was first
   validated using Eggerthella lenta exposed to our drug panel, as well as
   digoxin as a positive control ([94]Supplementary Figs. 4-[95]5,
   [96]Supplementary Table 5). Differential E. lenta gene expression was
   observed in 17 of 20 drug conditions, including the expected
   upregulation of the cgr operon by digoxin^[97]15 ([98]Supplementary
   Fig. 5). The transcriptomic pipeline was thus applied to all
   drug-microbe pairs in biological duplicates ([99]Supplementary Fig. 4b,
   [100]Supplementary Tables 6-[101]8).

   Overall, substantial and consistent transcriptional responses were
   observed across drug classes and bacterial phyla ([102]Supplementary
   Table 6). We used the magnitude of global transcriptional response as a
   common proxy^[103]34 by calculating the ratio of differentially
   expressed genes (DEGs) to the total gene count per genome, which we
   refer to as the DEG ratio or DEGR ([104]Fig. 1c). All drugs produced
   differential expression in at least one strain. The largest aggregate
   DEGR was produced by simvastatin (0.014), followed by sertraline
   (0.010), levothyroxine (0.006), and paroxetine (0.004) ([105]Fig. 1c,
   top bars). Notably, in many cases drugs eliciting large global
   transcriptional changes also exhibited broad-spectrum toxicity in the
   growth screen ([106]Fig. 1b, [107]c; [108]Supplementary Fig. 6).

   Large transcriptional changes are often associated with expression
   changes of global regulators and transcription factors (TFs)^[109]35.
   Indeed, the drug-microbe pairings associated with the highest DEGRs
   (Simvastatin-B. stercoris, Sertraline-F. saccharivorans, Simvastatin-A.
   shahii) also induced the highest numbers of TFs (14, 13, and 12,
   respectively) ([110]Supplementary Fig. 7). However, in several cases,
   bacterial-drug exposures generated differential expression without TF
   modulation, such as D. longicatena exposed to metoprolol tartrate and
   P. dorei exposed to various selective serotonin reuptake inhibitors
   (SSRIs) ([111]Supplementary Fig. 7). In these cases, differential
   expression was not consistently correlated with antimicrobial activity
   ([112]Fig. 1b, [113]c). These results indicate that drug-microbe
   exposures producing large and broad gene expression changes often
   correlate with drug toxicity, while those eliciting specific
   transcriptional responses may not.

   To understand the functional impact of different exposures on the gut
   microbiota, we performed pathway enrichment analysis using the KEGG
   (Kyoto Encyclopedia of Genes and Genomes) database to classify DEGs
   across the drug panel, agnostic of strain identity ([114]Fig. 2). Most
   differentially regulated modules (p[adj]<0.01, FC>4) were associated
   with mechanisms of antibiotic resistance ([115]Fig. 2a). Specifically,
   pathways related to multidrug resistance, transport, and two-component
   systems were enriched. Simvastatin, sertraline, and amlodipine – the
   drugs exhibiting the broadest toxicity in growth screens – strongly
   upregulated multidrug resistance pathways associated with efflux
   transporters. Further, trazodone and levothyroxine, which triggered
   differential expression in Bacteroidetes, Firmicutes, and
   Actinobacteria without impacting growth ([116]Fig. 1b, [117]c), showed
   similarly enriched pathways to more toxic screened compounds ([118]Fig.
   2a). While growth deficits were not detected in trazodone and
   levothyroxine-treated cultures at the maximum concentrations, this
   result suggests that these compounds could exhibit toxicity in vivo at
   concentrations exceeding 100 µM. Notably, post-treatment microbiota
   changes have been detected in hypothyroid patients taking
   levothyroxine^[119]36 and psychiatric patients taking
   trazodone^[120]37.

Figure 2. Pathway enrichment analysis reveals modulation of conserved efflux
pathways by top pharmaceutical compounds.

   Figure 2.
   [121]Open in a new tab

   (a) KEGG modules with significant enrichment (p[adj] <0.05 by two-sided
   independent t-test with Benjamini & Hochberg correction, e-value<10^−5)
   based on pathway analysis of DEG datasets are shown. Drug condition is
   labeled on the x-axis and corresponds to background color of the bubble
   plot. IDs of KEGG modules with significant enrichment are labeled on
   the left y-axis, and KEGG pathway groups are labeled on the right
   y-axis. Bubble color indicates direction of regulation
   (red=up-regulation, blue=down-regulation), bubble intensity indicates
   log[2](FC) enrichment, and bubble size indicates level of significance
   (-log[10] of p[adj]). (b) Top 19 upregulated and top 10 downregulated
   functional orthologs within the DEG dataset across drugs, as annotated
   within the KEGG orthology (KO) database. Bar color indicates direction
   of transcriptional regulation as in (a). -Log[10] of p[adj] (calculated
   as in (a) is shown on the x-axis. K-number assigned by the KO database
   is labeled on the left y-axis, and putative ortholog function is
   labeled on the right y-axis.

   We next examined functional orthologs within pathways that exhibited
   the greatest magnitudes of differential regulation across drug
   conditions ([122]Fig. 2b). Among the top 7 most upregulated orthologs,
   6 corresponded to conserved multidrug efflux pumps ([123]Fig. 2b).
   Among these orthologs, the top 4 (HAE1, AcrA, mexK, oprM) belonged to
   the Resistance-Nodulation-Division (RND) permease superfamily, a drug
   and heavy metal efflux system whose upregulation is associated with
   gram-negative bacterial antibiotic resistance^[124]38,[125]39. The
   fifth and sixth orthologs (bcrB and ABCB-BAC) belonged to the
   ATP-binding cassette (ABC) superfamily, which was also highly
   represented among the top downregulated orthologs across drug
   conditions ([126]Fig. 2b). Interestingly, the ABC superfamily is not
   considered to contribute substantially to bacterial multidrug
   resistance; however, these pumps are highly associated with
   chemotherapy resistance in eukaryotic cells^[127]40. Together, these
   observations suggest that gut bacterial strains utilize conserved
   multidrug efflux pumps to ameliorate toxicities of human-targeted drug
   compounds^[128]5.

   Given the link between the identified RND-type efflux pumps and
   antibiotic resistance in gram-negative species^[129]38,[130]39, we
   wondered whether gut bacteria employ the same resistance pathways to
   ameliorate toxicities in commonly prescribed antibiotic classes. To
   explore this question, we generated transcriptomes for 4 species (A.
   rectalis H1, D. longicatena H1, B. longum H1, P. vulgatus H1) exposed
   to 7 commonly prescribed antibiotics that target bacterial synthesis of
   DNA (ciprofloxacin), proteins (tetracycline, streptomycin,
   erythromycin), or the cell membrane (cefotaxime, ampicillin,
   vancomycin) ([131]Supplementary Figs. 4 & [132]8, [133]Supplementary
   Table 9). Principal coordinates analysis (PCoA) revealed that
   transcriptomic responses to traditional antibiotics clustered by drug
   target across sensitive strains (A. rectalis, D. longicatena, P.
   vulgatus), suggesting that bacteria utilize conserved pathways to
   mitigate toxicities of drugs with similar bacterial targets
   ([134]Supplementary Fig. 8a). Notably, minimal transcriptomic responses
   to streptomycin exposures were observed, consistent with inherent
   streptomycin resistance of anaerobically grown bacterial
   isolates^[135]41. Further, B. longum, the only bacterial strain in our
   panel that did not exhibit growth sensitivity to at least one screened
   pharmaceutical, demonstrated minimal transcriptional responses to
   antibiotic exposures as compared to a vehicle control, a behavior
   consistent with resistance^[136]32. To compare bacterial pathways
   induced by traditional antibiotic exposures, we next performed a KEGG
   pathway enrichment analysis of DEGs in the 4 tested bacteria exposed to
   orally delivered drugs or traditional antibiotics ([137]Supplementary
   Fig. 8b). Membrane-directed antibiotics upregulated pathways associated
   with cell membrane synthesis and modification (vancomycin resistance,
   fatty acid synthesis pathways), while ribosomal-targeted antibiotics
   induced upregulation of ribosomal pathways and downregulation of
   ATP-synthesis. Remarkably, our analysis revealed no overlap in
   bacterial pathways induced by non-antibiotic pharmaceuticals and
   traditional antibiotic compounds, suggesting that non-antibiotic orally
   delivered compounds impact unique mechanisms of prokaryotic multidrug
   resistance.

Drugs impact gut bioavailability, biosynthesis and toxicity

   To obtain a deeper functional understanding of bacterial-drug
   interactions at the operon level, we searched for clusters of DEGs
   transcribed within the same operon. Numerous differentially expressed
   operons were identified across strains ([138]Supplementary Table 8).
   Existing mechanistic studies of gut microbiota drug interactions have
   identified a range of bacterial drug responses that can cause
   differential treatment outcomes, including drug metabolism, toxicity
   mitigation and alteration of the prokaryotic
   metabolome^[139]5,[140]19,[141]42. In order to explore the diversity of
   physiologically relevant bacterial xenobiotic interactions, we selected
   3 representative operons associated with drug metabolism, vitamin
   biosynthesis, and toxicity mitigation for further examination
   ([142]Fig. 3).

Figure 3. Top pharmaceutical compounds impact gut bacterial metabolism,
vitamin production, and mitigation of toxic metabolites.

   Figure 3.
   [143]Open in a new tab

   (a). Top schematic shows the components of the tartrate dehydratase
   (ttd) operon in E. coli. Genes within the operon are labeled on the
   x-axis, and isolates containing ttd analogs are labeled on the y-axis.
   Heatmap color indicates log[2](FC) in transcripts per million (TPM)
   expression. Stars indicate p[adj] <0.05 (*), p[adj] <0.01 (**),
   p[adj]<0.001(***) calculated by two-sided independent t-test with
   Benjamini & Hochberg correction. Gray panels indicate absence of a gene
   analog. (b) Relative abundance of bacterial isolates containing the ttd
   operon was calculated within a dataset of T2D patients^[144]48 (n=323
   total patients). The x-axis indicates patient cohort, grouped by
   healthy controls not taking metoprolol (red dots), T2D patients taking
   metoprolol (green dots), and T2D patients not taking metoprolol (blue
   dots). Relative abundance of isolates is shown on the y-axis. P-values
   calculated by two-sided independent t-test are annotated with brackets.
   (c) Top schematic shows the components of the riboflavin biosynthesis
   (rib) operon in D. longicatena H1. Drugs are labeled on the x-axis
   (grouped and colored by class), and gene components of the rib operon
   are labeled on the y-axis. Heatmap color indicates log[2](FC) in TPM
   expression. Stars indicate significance as calculated and annotated in
   (a). (d) Regulation of acrAB-tolC within Bacteroidales isolates by
   atorvastatin (left panel) and simvastatin (right panel). Bottom
   schematic shows the components of the acrAB-tolC operon in E. coli.
   Genes within acrAB-tolC are labeled on the x-axes, and isolates
   containing AcrAB-TolC analogs are labeled on the y-axis. Heatmap color
   indicates log[2](FC) in TPM expression, as indicated by the bottom
   right red-to-blue color bar. Stars indicate significance as calculated
   and annotated in (a). Dots to the right of strain rows indicate percent
   operon identity as compared to the P. vulgatus AcrAB reference^[145]57,
   as indicated by the top right red-to-yellow color bar.

   First, we observed upregulation of the tartrate dehydratase (ttd)
   operon in E. coli, A. rectalis, and D. longicatena exposed to
   metoprolol tartrate, a beta blocker used to treat hypertension^[146]43
   ([147]Fig. 3a). Metoprolol is chemically formulated either as
   metoprolol tartrate or metoprolol succinate for immediate release or
   extended release, respectively. Which form of metoprolol is
   administered can impact bioavailability, with metoprolol tartrate
   producing higher peak-to-trough variation among patients^[148]44. For
   bacteria, dietary tartrate metabolism through the ttd operon provides
   another carbon source in the gastrointestinal tract^[149]45. Our
   results suggest that certain gut microbes can metabolize the tartrate
   salt of immediate-release metoprolol, which could contribute to
   fluctuations in drug bioavailability within patient cohorts. These
   findings could also have broader implications for other
   tartrate-conjugated drugs^[150]46. Furthermore, upregulation of
   tartrate metabolism by metoprolol could have unintended negative
   consequences in the context of cardiovascular disease^[151]47.
   Increases in bacterial tartrate metabolism have been linked to
   metabolic disorders, including atherosclerotic cardiovascular disease
   (ACVD), type 2 diabetes (T2D) and obesity^[152]47. Enrichment in
   tartrate metabolism has also been associated with higher abundance of
   E. coli, a metoprolol tartrate metabolizer identified in our screen.
   Using a published clinical metagenomic microbiome dataset of T2D
   patients^[153]48, we performed an additional analysis exploring the
   distribution of ttd genes among T2D patients consuming metoprolol. We
   found that T2D patients taking metoprolol had an increased abundance of
   the ttd operon in their gut metagenome compared to patients not taking
   metoprolol (p=0.046) and healthy controls (p=0.008, [154]Fig. 3b).
   Therefore, our findings combined with prior studies suggest that
   treatment with metoprolol tartrate for hypertension, which is the
   strongest predictor for ACVD^[155]49, might inadvertently influence the
   pathophysiology of a metabolic disorder via an increase in microbiota
   tartrate metabolism.

   Next, we identified differential regulation of the riboflavin
   biosynthesis (rib) operon in D. longicatena by several drugs within our
   panel ([156]Fig. 3c). The gut microbiome is an important source of
   riboflavin (vitamin B2) in humans^[157]50. Riboflavin production is an
   essential pathway in bacteria, as downstream metabolites flavin
   mononucleotide (FMN) and flavin adenine dinucleotide (FAD) are
   co-enzymes for oxidases, reductases, and dehydrogenases^[158]51. In
   gram-positive bacteria, riboflavin biosynthesis is downregulated by an
   FMN-responsive riboswitch, which also responds to roseoflavin^[159]51.
   In our transcriptomics dataset, several SSRIs (sertraline, paroxetine,
   fluoxetine, duloxetine), a cardiovascular medication (amlodipine),
   simvastatin, and levothyroxine all downregulate riboflavin production
   in D. longicatena. We also observed upregulation of the rib operon by
   D. longicatena in response to trazodone and atorvastatin exposure.
   Clinically, vitamin B2 deficiency is associated with higher risk of
   psychiatric disorders for which SSRIs are indicated (e.g., depression),
   though whether these depletions are a cause or a result of disease is
   not understood^[160]52,[161]53. Reduced riboflavin concentrations can
   also contribute to hyperhomocysteinemia, a well-studied independent
   risk factor for atherosclerosis^[162]54,[163]55. Our data suggest that
   the administration of certain SSRIs and cardiovascular medications
   could modulate vitamin B2 levels in the setting of mood disorders or
   heart disease, respectively, which could have unintended detrimental
   consequences on the disease state.

   Finally, we observed a strong statin-mediated upregulation of genes
   encoding the AcrAB-TolC efflux pump in all Bacteroidales strains tested
   (B. fragilis, B. stercoris, B. uniformis, P. dorei, P. vulgatus, P.
   distasonis, A. shahii) ([164]Fig 3d). The AcrAB-TolC pump has been
   linked to resistance against several classes of antibiotics as well as
   non-antibiotic orally delivered pharmaceuticals including methotrexate
   and tamoxifen^[165]5,[166]56. The pump is also known to mediate
   bacterial sensitivity to retinol and secondary bile acids^[167]57.
   Given that previous studies have linked statin use with changes in gut
   microbiome composition^[168]25, we thus explored the mechanism by which
   upregulation of AcrAB-TolC by statins could impact microbial physiology
   in the gut.

Statins modulate AcrAB-TolC activity in gut microbes

   Statins have been widely studied for their potential as non-traditional
   antibiotics, with simvastatin garnering particular attention for its
   activity against multidrug resistant pathogens such as Staphylococcus
   aureus^[169]58. Simvastatin alters the composition of gut microbiome in
   patients, but the mechanism driving this microbial shift is not
   understood^[170]25. In our transcriptomic screening, we found that
   simvastatin strongly upregulated acrAB-tolC genes across multiple
   Bacteroidales strains, while atorvastatin significantly upregulated
   acrAB-tolC genes only in P. distasonis ([171]Fig. 3d). Interestingly,
   these transcriptional profiles correlated with strong and modest
   toxicities exhibited in our growth screens by simvastatin and
   atorvastatin, respectively ([172]Fig. 1b), motivating further
   exploration into the role of AcrAB-TolC in the interplay between gut
   microbes and antimicrobial compounds.

   In Bacteroides species, the AcrAB-TolC efflux pump ameliorates toxicity
   of antibiotics and common dietary metabolites such as retinol (vitamin
   A)^[173]57. We thus hypothesized that drugs altering AcrAB-TolC
   expression in Bacteroidales species could result in targeted changes in
   the toxicity of other compounds. Using P. distasonis, we first
   determined the MIC of retinol (>24 µg/mL) and several common
   antibiotics with different mechanisms of action whose resistance is
   mediated by AcrAB-TolC^[174]56 ([175]Supplementary Table 10). In the
   presence of simvastatin, the MICs of several of these compounds were
   significantly shifted, with retinol showing the most pronounced effect
   (>2-fold decrease, [176]Supplementary Fig. 9). We then exposed four
   Bacteroidales strains (P. distasonis H1, P. dorei H1, P. vulgatus H1,
   P. vulgatus ATCC 8482) to retinol in the presence or absence of statins
   and additional drugs (sertraline, trazodone, amlodipine) observed to
   upregulate the acrAB-tolC operon ([177]Fig. 4a; [178]Supplementary Fig.
   10). As a negative control, we co-incubated retinol-exposed
   non-Bacteroidales strains (E. coli, A. rectalis, D. longicatena) with
   simvastatin ([179]Supplementary Fig. 11). At 20 µM, simvastatin
   enhanced the sensitivity to retinol in all Bacteroidales strains, with
   the strongest shift seen in P. distasonis and P. vulgatus (4-fold
   reduction in MIC, [180]Fig. 4a). In non-Bacteroidales strains (E. coli,
   D. longicatena, A. rectalis), retinol MICs were not significantly
   altered by simvastatin co-incubation ([181]Supplementary Fig. 11).
   Notably, we observed a modest shift approaching significance in
   simvastatin-treated A. rectalis, which we attribute to combined global
   toxicity burdens of simvastatin and retinol, unrelated to the
   AcrAB-TolC pump. Sertraline modestly increased retinol sensitivity in
   all Bacteroidales strains, amlodipine showed similar effects in two of
   four strains (P. vulgatus H1, P. dorei H1), atorvastatin modestly
   decreased the MIC of retinol in P. distasonis only, and trazodone did
   not significantly influence MIC of retinol in any strains ([182]Fig.
   4a, [183]Supplementary Fig. 10).

Figure 4. Statin exposure alters Bacteroides sensitivity to common dietary
metabolites via the AcrAB-TolC efflux pump.

   Figure 4.
   [184]Open in a new tab

   (a) Bacteriodales isolates exposed to retinol at various concentrations
   in the presence of 20 µM simvastatin (orange lines), atorvastatin
   (green lines), or vehicle control (blue lines). Graph titles indicate
   strain identity, with retinol concentration shown on the x-axis and
   growth (OD^600, optical density at 600 nm wavelength) relative to a
   vehicle control shown on the y-axis. (b) P. distasonis H1 exposed to
   various deoxycholic acid (DCA) concentrations in the presence of 20 µm
   simvastatin (orange lines) or vehicle control (blue lines). Bile acid
   concentration is shown on the x-axis, and growth relative to a vehicle
   control is shown on the y-axis. Standard deviation bars (n=3
   biologically independent culture replicates per condition) are shown,
   and stars indicate p[adj] <0.05 (*), p[adj] <0.01 (**), and p[adj]
   <0.001 (***) calculated by two-sided independent t-test with Benjamini
   & Hochberg correction. (c) P. vulgatus H1 overexpressing different
   copies of AcrAB-TolC (see schematic) or a plasmid control exposed to
   different concentrations of retinol in the presence of 20 µm
   simvastatin (n=3 biologically independent culture replicates, orange
   lines) or vehicle control (n=3 biologically independent culture
   replicates, blue lines). Retinol concentration is shown on the x-axis,
   and average OD^600 growth across replicates relative to average growth
   of vehicle controls is shown on the y-axis. (d) Structures of
   atorvastatin, simvastatin, and tenivastatin (simvastatin-hydroxy acid).
   The lactone ring in simvastatin is shown in red. (e) P. distasonis H1
   and P. vulgatus ATCC8482 grown in liquid mGAM supplemented with
   simvastatin (orange line), tenivastatin (raspberry line), or
   atorvastatin (green line) at 100 µM concentrations. Standard error bars
   (n=4 biologically independent culture replicates per condition) for
   average OD^600 growth across replicates are shown. (f) P. distasonis H1
   retinol MIC curves in the presence of 20 µM simvastatin (n=2, orange
   line), tenivastatin (n=3 biologically independent culture replicates,
   raspberry line), or vehicle control (n=3 biologically independent
   culture replicates, blue line). p[adj] values and standard deviations
   of average OD^600 growth at each concentration are annotated as in (b).

   Deoxycholic acid (DCA) is a secondary bile acid whose toxicity is also
   modulated by AcrAB-TolC^[185]57. We thus tested whether simvastatin
   enhanced DCA toxicity in P. distasonis. Simvastatin at 20 µM
   significantly limited P. distasonis growth when co-incubated with DCA
   (MIC = 64 µg/ml, 66% reduction) ([186]Fig. 4b). Simvastatin did not
   lower the MIC of other secondary bile acids (i.e., hydroxydeoxycholic
   acid, ursodeoxycholic acid) whose toxicity is not mediated by
   AcrAB-TolC^[187]57 ([188]Supplementary Fig. 12). Importantly, the
   simvastatin-altered MICs of retinol and DCA each fall within estimated
   colon concentration ranges of these metabolites in the human
   gut^[189]57,[190]59.

   Having established that simvastatin significantly upregulates
   AcrAB-TolC in Bacteroidales species, and that this exposure increases
   Bacteroidales sensitivity to retinol and DCA, we next sought to confirm
   that AcrAB-TolC upregulation directly increases sensitivity to dietary
   metabolites. In order to establish causality of AcrAB-TolC upregulation
   in increasing collateral sensitivity of Bacteroidales species to
   retinol, we engineered a P. vulgatus H1 isolate to overexpress
   different copies of AcrAB-TolC or a control cargo, and measured retinol
   MIC for engineered isolates exposed to simvastatin or DMSO ([191]Fig.
   4c, [192]Methods). As compared to plasmid overexpression of a cargo
   control, overexpression of either copy of AcrAB-TolC increased
   collateral sensitivity to retinol, establishing that efflux pump
   overexpression directly increases collateral sensitivity to vitamin A.
   Notably, while minimal toxicity was associated with plasmid
   overexpression in the control P. vulgatus strain, acrAB-tolC expression
   also increased toxicity of simvastatin at 20 µM concentrations, as
   exhibited by P. vulgatus growth at retinol concentrations of zero.
   Further, overexpression of acrAB-tolC copy 2 caused greater enhancement
   of retinol and simvastatin toxicity in P. vulgatus as compared to copy
   1, suggesting distinct efficacy levels of these genomic efflux pump
   copies. Taken together, our transcriptomic dataset revealed a
   statin-mediated interaction with the AcrAB-TolC efflux pump that could
   enhance the antimicrobial effects of common dietary and host
   metabolites in vivo, potentially accounting for the gut microbiota
   alterations seen in statin-treated patient populations^[193]25.

   Simvastatin is generally administered as a lactone prodrug, while
   atorvastatin is administered in the active compound form^[194]58
   ([195]Fig. 4d). Given the role of lactone moieties in many
   antibiotics^[196]60, we wondered whether the transcriptional regulation
   of the AcrAB-TolC efflux pump by simvastatin might be linked to the
   toxicity of its lactone moiety. To test whether the lactone ring in
   simvastatin induces toxicity, we incubated simvastatin-sensitive P.
   distasonis and P. vulgatus strains with the lactone prodrug
   (simvastatin), the non-lactone activated compound (tenivastatin), and
   atorvastatin ([197]Fig. 4e). Interestingly, tenivastatin, which does
   not contain a lactone ring, did not exhibit toxicity against either
   tested strain. Notably, AcrAB-TolC preferentially binds to hydrophobic
   substrates (e.g., lipophilic lactone side chains)^[198]61,[199]62. To
   determine whether AcrAB-TolC modulation of bacterial inhibition by
   dietary metabolites depends on the presence of the toxic simvastatin
   lactone moiety, we exposed a P. distasonis strain to various
   concentrations of retinol in the presence or absence of 20 µM
   simvastatin or tenivastatin ([200]Fig. 4f). Tenivastatin did not lower
   the MIC of retinol, suggesting that the intact lactone ring of the
   unmetabolized prodrug simvastatin is necessary to enhance bacterial
   retinol sensitivity.

   Next, to better understand the prevalence of AcrAB-TolC in the human
   gut microbiome, we assessed its frequency in 4,930 representative
   metagenome-assembled genomes (MAGs) sourced from healthy human
   microbiomes ([201]Fig 5a, [202]Methods)^[203]63. 106 MAGs contained a
   complete AcrAB-TolC homolog (e <10^−5, coverage >90%, identity >40%),
   of which 59.4% were linked to the Bacteroidaceae family. The high
   prevalence of AcrAB-TolC in Bacteroidaceae suggests that
   Bacteroidaceae-dominated individuals may be particularly susceptible to
   simvastatin-mediated microbiota changes.

Figure 5. AcrAB-TolC is linked to gut microbiota shifts in statin-treated
patient populations.

   Figure 5.
   [204]Open in a new tab

   (a) The number of genomes containing acrAB-tolC analogs based on
   Ref^[205]63. Bacterial family is indicated on the x-axis and by bar
   color. The number of acrAB-tolC-containing genomes within each family
   is indicated on the y-axis. (b) The relative abundance of acrAB-tolC
   within treatment groups of the BMIS cohort is shown. Statin or control
   cohorts are indicated on the x-axis. Simvastatin (red) and atorvastatin
   (blue) cohorts are highlighted, and corresponding p-values (calculated
   by two-sided independent t-test) are annotated with colored brackets.
   Each dot represents a different patient metagenome. (c) Relative
   abundance of bacterial species within the BMIS dataset, annotated using
   Ref^[206]63. For each genus, the log[10] relative abundance within
   simvastatin-treated (red) and non-treated (black) individuals is shown
   as a box and whisker plot. Box hinges correspond to the 25^th and 75^th
   percentiles, and whiskers extend to values within 1.5 the interquartile
   range (outliers omitted). Stars indicate p<0.05 (*) as calculated by
   two-sided independent t-test. (d) Relative abundance of Bacteroides
   species (left two panels) or all gut bacteria (right two panels) in
   simvastatin-treated and untreated individuals within the BMIS cohort.
   Within each panel, the x-axis indicates simvastatin (red) or control
   (black) treatment cohort, and relative abundances of species with or
   without the acrAB-tolC operon are labeled on the y-axis. P-values are
   annotated as in (b).

   To explore the clinical relevance of simvastatin-induced upregulation
   of the acrAB-tolC operon, we next quantified the metagenomic abundance
   of this operon’s homologs in stool samples from the Body Mass Index
   Spectrum (BMIS) cohort^[207]25 ([208]Fig. 5b-[209]d). Notably, the
   AcrAB-TolC pump was highly prevalent in fecal metagenomes of the BMIS
   cohort, with 91.6% of study participants showing >20% median relative
   abundance of AcrAB-TolC across treatment groups ([210]Fig. 5b).
   Further, a significant decrease in relative abundance of the acrAB-tolC
   operon (p[adj] = 0.0102) in the gut microbiome of patients receiving
   simvastatin was observed ([211]Fig. 5b). In contrast, acrAB-tolC gene
   abundance was not changed in patients who received atorvastatin,
   fluvastatin, pravastatin, or rosuvastatin. Accordingly,
   simvastatin-treated patients also showed a significant depletion in
   Bacteroides species (p = 0.0030), as well as increases in Alistipes and
   unclassified Firmicutes ([212]Fig. 5c). Moreover, Bacteroides and other
   bacteria containing a complete AcrAB-TolC homolog were significantly
   depleted in simvastatin-treated cohorts as compared to untreated
   controls (p = 0.0013), whereas bacteria missing this operon were not
   depleted (p = 0.699) ([213]Fig. 5d). Notably, a relative abundance
   analysis for the separate gene components of the acrAB-tolC operon in
   statin-treated cohorts identified the most significant depletion (p =
   0.003) in the AcrB2 gene, which encodes a homotrimer docked to the
   intracellular bacterial membrane that binds to pump substrates
   ([214]Supplementary Fig. 13)^[215]61.

   Taken all together, we posit a potential model of statin-mediated
   depletion of gut microbes containing the AcrAB-TolC efflux pump
   ([216]Supplementary Fig. 14). Under normal conditions, gut microbial
   residents containing AcrAB-TolC homologs (e.g., Bacteroidales species)
   utilize this pump to prevent periplasmic and cytoplasmic accumulation
   of metabolites, such as retinol and deoxycholic acid, to toxic levels.
   In the context of simvastatin treatment, acrAB-tolC is upregulated,
   leading to increased collateral sensitivity to dietary metabolites and
   reduced cell viability. Efflux pump upregulation leading to decreased
   bacterial viability has been observed in clinical cohorts^[217]64,
   although the mechanism of this phenomenon has not been elucidated to
   date. Notably, the prototypical acrAB-tolC operon is regulated by the
   self-contained acrR repressor ([218]Fig. 3d), which in E. coli is
   strongly upregulated under conditions of global stress^[219]65. Thus,
   it is possible that by inducing AcrR expression, engineered AcrAB-TolC
   pump overexpression could indirectly impact co-regulation of global
   stress response and lead to reduced cell viability. Regardless, with
   simvastatin treatment, increased expression of AcrAB-TolC causes
   decreased viability of gut bacteria and consequential alterations in
   gut microbiota composition.

Strain differences in drug-mediated transcriptomic responses

   Finally, we sought to explore how drug-microbe transcriptional
   responses differed among different strains in our panel ([220]Fig. 6,
   [221]Supplementary Fig. 15). First, we analyzed the DEGRs across all
   drugs by principal coordinates analysis (PCoA) to quantify similarities
   in global transcriptomic profiles of strains. Interestingly, we found
   that microbial responses to drugs clustered significantly at the family
   level (i.e., Bacteroidaceae, Lachnospiraceae), (p = .003 by PERMANOVA
   test) ([222]Fig. 6a), suggesting that closely related bacteria may
   respond similarly to different drugs. However, a detailed comparison of
   conspecific strains revealed a more variable pattern of response. For
   example, while strains of A. rectalis (H1 vs. ATCC 33656) exhibited
   highly consistent drug responses, strains of B. uniformis (H1 vs. ATCC
   8492) responded more heterogeneously to pharmaceuticals ([223]Fig. 6a).

Figure 6. Diverse patterns of transcriptional response among conspecific
bacterial isolates.

   Figure 6.
   [224]Open in a new tab

   (a) Plot of PCoA results comparing the transcriptomic responses of
   tested gut bacterial isolates to top pharmaceutical compounds, using
   DEGRs within each drug condition as features. Variance contribution of
   first and second principal coordinates are shown on the x- and y- axes,
   respectively. Bacterial family of each isolate is indicated by color.
   Conspecific strains are connected by dotted lines. (b) For each
   possible pair of Bacteroidaceae strains, a comparison of average
   nucleotide identity (ANI, shown on x-axis) and Spearman correlation of
   DEGRs across all drug conditions (shown on y-axis) is shown. Each dot
   indicates a pairwise comparison, with color indicating whether the
   comparison was performed on conspecific (red) or allospecific (green)
   strains. (c) Differential regulation of selected operons by simvastatin
   (purple dots), sertraline (green dots), paroxetine (blue dots), or
   levothyroxine (red) dots are depicted for different conspecific strains
   (purple boxes depict operon expression in strain H1, yellow boxes
   depict operon expression in a publicly available strain). Grey dots
   indicate p[adj] >0.05. Species and operon identities are indicated by
   panel titles, and strain identity is shown on the right y-axis.
   Log[2](FC) in gene expression as compared to a vehicle control is shown
   on the left y-axis and aligned gene IDs are shown on the x-axis.
   Absence of an operon is indicated by a beige panel with red annotation.
   Gray dashed lines and box indicate significance threshold of FC>2.

   To further delineate the strain specificity of drug responses, we next
   quantified the average nucleotide similarity (ANI) between different
   Bacteroidales strain pairs, which we compared with the correlation of
   global transcriptomic response (DEGR values) between strains ([225]Fig.
   6b). We observed higher correlations (R > 0.7) in overall drug response
   between pairs with the greatest ANI values (ANI > 0.95), suggesting
   that the more similar the bacterial genome, the more similar the global
   transcriptomic response across drug conditions ([226]Fig. 6b). In order
   to compare the transcriptomic responses of conspecific strains more
   closely, we next aligned the genomes for all conspecific strains
   ([227]Methods). We observed high correlations (R^2 > 0.79) in
   expression of shared genes between all conspecific strain pairs,
   suggesting shared transcriptional responses between isolates
   ([228]Supplementary Fig. 15).

   We next sought to investigate strain-level drug response at the operon
   level. Remarkably, we identified two common types of transcriptional
   differences between conspecific strains: 1) absence of an operon in one
   isolate, and 2) differential regulation of a shared operon between
   isolates, with both differences possibly having functional impacts. For
   example, sertraline upregulated a mobile element containing Type IV
   secretion system machinery in A. rectalis H1 that was not present in A.
   rectalis ATCC 33656 ([229]Fig. 6c). Similarly, sertraline and
   paroxetine upregulated BepE efflux machinery in B. uniformis H1 that
   was not present in ATCC 8492 ([230]Fig. 6c). On the other hand, several
   drugs (sertraline, paroxetine, levothyroxine, simvastatin)
   downregulated a shared putative L-arabinose utilization operon in A.
   rectalis H1, but this same operon (98.31% average conservation in
   coding sequence) was not differentially expressed in A. rectalis 33656
   ([231]Fig. 6c). Conversely, these compounds downregulated a shared
   putative transketolase operon (99.89% average conservation in coding
   sequence) in A. rectalis 33656 only ([232]Fig. 6c). Interestingly, in
   addition to differences in shared operon expression between isolates,
   we also observed instances of differential regulation of different gene
   copies by the same drug perturbation ([233]Fig. 6c). For example,
   within F. saccharivorans DSMZ 26062, levothyroxine, sertraline,
   simvastatin, and paroxetine exposures triggered simultaneous up- and
   down-regulation of different copies of araQ, a permease protein
   associated with L-arabinose transport (99.25% average conservation in
   coding sequence; [234]Fig. 6c). A high degree of strain-level
   functional diversity could explain the substantial microbiome
   compositional variation often observed in clinical drug studies that
   only rely on 16S taxonomic analysis rather than whole genome
   sequencing. Together, these results underline the importance of
   assessing bacterial drug responses at the strain level with genomic
   information, as well as extending traditional microbiota drug screens
   to multiple conspecific strains to capture the full intra-strain
   diversity of drug responses.

DISCUSSION

   While numerous clinical studies suggest prevalent drug-microbiota
   interactions^[235]5,[236]13,[237]23,[238]25, genetic-level
   understanding of xenobiotic-induced microbiota shifts remains limited.
   In this study, we generated to our knowledge the largest transcriptome
   dataset for exploring gut microbiota-drug responses to common orally
   delivered drugs, totaling >400 bacterial-drug pairwise interactions.
   With this dataset, we uncovered that simvastatin-mediated upregulation
   of the AcrAB-TolC efflux pump generates collateral toxicity to dietary
   metabolites ex vivo. Using a clinical metagenomic dataset, we confirmed
   that AcrAB-TolC-containing gut bacteria are depleted in patients taking
   simvastatin, suggesting that upregulation of acrAB-tolC under statin
   exposure increases cellular toxicity.

   Beyond statins, our data established clinically relevant links between
   cardiovascular medications and SSRIs and the gut microbiome. Our ex
   vivo data and computational analysis of T2D patient data suggest the
   upregulation of tartrate metabolism in gut microbes exposed to
   metoprolol. This result motivates the need for further clinical studies
   to determine whether metoprolol tartrate places patients at higher risk
   for developing T2D and other metabolic disorders. The clinical dataset
   used in our analysis does not distinguish metoprolol tartrate versus
   metoprolol succinate, which are used short- or long-form treatment. As
   metoprolol tartrate exhibits a greater inter-individual variation in
   bioavailability^[239]44, future work exploring the impact of microbiome
   metabolism on the in vivo bioavailability of metoprolol tartrate and
   other tartrate-conjugated pharmaceuticals is warranted. Separately, our
   findings provide the first evidence that drugs can modulate riboflavin
   production by gut bacteria. Mechanistic delineation of how
   pharmaceuticals induce or repress vitamin B2 production could enable
   the development of tools to promote microbial production of this
   essential vitamin^[240]51. As vitamin B2 depletions have been linked to
   depression and atherosclerosis^[241]52,[242]53,[243]55, our data also
   motivate clinical studies of whether rib-modulating medications impact
   key vitamin reservoirs in the microbiomes of psychiatric and
   cardiovascular patients, and whether these interactions impact disease
   pathology.

   From an ecological standpoint, shared transcriptional responses could
   provide mechanistic insights into how pharmaceuticals drive shifts in
   microbiome populations over time. Our pathway meta-analysis
   demonstrated that many drugs upregulate highly conserved bacterial
   multidrug efflux pumps. Efflux-associated pathways were distinct from
   resistance pathways upregulated by traditional antibiotics, suggesting
   novel mechanisms of drug-driven microbiota responses. Our results
   showing efflux-mediated collateral toxicity induced by
   simvastatin-retinol coincubation underline the importance of exploring
   the molecular efflux networks to identify drug-efflux pump interactions
   that can result in gut microbiota shifts. Further, recent data has
   linked chronic prescription drug use and polypharmacy with
   antimicrobial resistance in gut microbiota^[244]4. One work even linked
   the upregulation of efflux pumps by four antidepressants used in our
   study (sertraline, duloxetine, escitalopram, bupropion) to the
   promotion of antimicrobial resistance among gut community
   members^[245]66. These data all suggest that chronic medication-induced
   increases in the activity of efflux pumps could contribute to higher
   rates of antimicrobial resistance among gut microbiota within patient
   populations.

   Our study has several limitations. Ex vivo screens cannot capture the
   full environmental complexity that may exist between microbiota
   community members and the host in vivo, and transcriptomics does not
   directly probe host-mediated drug-microbiota interactions. While many
   of our observed inhibitory drug activities aligned with published
   clinical data^[246]5,[247]25, we detected neither growth nor
   transcriptional changes in gut bacteria exposed to metformin. The role
   of host-derived factors mediating metformin-induced microbiota shifts
   is well documented (i.e., secondary bile acids cause depletion of
   Bacteroides and enrichment of Escherichia species that trigger
   downstream regulation of the farnesoid receptor)^[248]67. This and
   other instances of low transcriptomic signals generated in our study by
   drugs linked to microbiota changes in vivo (e.g., omeprazole) could
   suggest missing host- or community- derived factors, motivating future
   in vivo studies that integrate these other factors (e.g., bile acids
   and dietary vitamins). While we leveraged published clinical microbiome
   data to further validate our ex vivo drug-microbe interaction results,
   our simvastatin findings could benefit from prospective clinical trials
   designed to dissect the impacts of simvastatin and disease status on
   gut microbiota dynamics, and to determine whether AcrAB-TolC is simply
   a biomarker or a functional mediator of simvastatin-associated gut
   microbiota shifts.

   Altogether, this work has demonstrated the utility of high-throughput
   transcriptomics to delineate microbiota-drug interactions. We envision
   that this low-cost and scalable pipeline can be easily applied to any
   microbiota-treatment pairing including xenobiotic, prebiotic, and
   probiotic treatments. These efforts will lead to a greater mechanistic
   understanding of how different environmental exposures impact the gut
   microbiome that in turn can affect host health and response to medical
   interventions.

ONLINE METHODS

Materials and Culture Conditions

   All bacteria used in this study have fully sequenced genomes, the
   information for which can be found in ([249]Supplementary Table 3).
   Natural bacterial isolates used in this study were derived from a
   single fecal sample taken from a healthy individual for an unrelated
   study. Sex and gender of this individual, while reported in the
   unrelated study, were deemed irrelevant in this study and were not
   reported. This work was approved and conducted under Columbia
   University Medical Center Institutional Review Board protocol AAAR0753,
   and written informed consent was obtained from the subject. Additional
   conspecific strains of A. rectalis, F. saccharivorans, B. uniformis, P.
   vulgatus, and B. fragilis, as well as the Eggerthella lenta strain
   described in Ref^[250]15, were obtained from public strain catalogs.
   All publicly available strains used in this study were either purchased
   from the American Type Culture Collection (Manassas, VA) or the Leibniz
   Institute DSMZ-German Collection of Microorganisms and Cell Cultures
   GmbH. All monoclonal isolates were Sanger sequenced (Azenta Life
   Sciences) at the 16S-V4 region pre- and post-experimentation to confirm
   strain identity. Natural isolates used in this study are available from
   the corresponding author upon request.

   For our drug panel, we selected 19 orally administered drugs from the
   top prescribed pharmaceuticals in 2017 according to the Agency for
   Healthcare Research and Quality ([251]Supplementary Table 4)^[252]31.
   This list included the top eight most prescribed drugs, included
   lisinopril, atorvastatin, levothyroxine, metformin, amlodipine,
   metoprolol, omeprazole, and simvastatin. Additional neurotransmitter
   modulators were selected from the top 25 drugs prescribed, as follows:
   sertraline, fluoxetine, citalopram, escitalopram, paroxetine,
   duloxetine, venlafaxine, amitriptyline, bupropion, and trazodone. We
   finally included lenalidomide due to its large market cap^[253]31. All
   chemicals used in this study were purchased from Sigma-Aldrich
   (Burlington, MA), Thermo Scientific Chemicals (Waltham, MA), or Avantor
   (Phillipsburg, NJ) ([254]Supplementary Tables 4, [255]9-[256]10).
   Probes for rRNA depletions were purchased from Integrated DNA
   Technologies (Coralville, IA). Unless otherwise noted, bacterial
   cultures were grown in ½ diluted Gifu Anaerobic Medium Broth, modified
   (mGAM, HyServe 05433) media prepared according to manufacturers’
   instructions. Prior to experimentation, all mGAM media was reduced for
   24 hours under anaerobic conditions (5% H2, 10% CO2, 85% N2) within a
   Coy Laboratory products anaerobic chamber. Chemical plates were
   prepared in 96-well format under aerobic sterile conditions and reduced
   for 3–6 hours prior to experimentation.

Bacterial Transcriptome Preparation and Extraction

   For all transcriptomic experiments, bacterial cultures at exponential
   phase were added to 96- deep-well plates (Axygen P-DW-20-C) containing
   pre-reduced chemicals to reach a final concentration of 20 µM (500 µM
   for metformin). After 90 minutes of exposure in shaken media at 37°C,
   cultures were centrifuged, supernatant removed, and cell pellets were
   transferred to −80°C prior to bacterial RNA extraction.

   Bacterial RNA was extracted using RNAsnap methods^[257]68. Briefly,
   frozen bacterial pellets were suspended in 500 µl RNAsnap mix (95%
   formamide, 18 mM ethylenediaminetetraacetic acid, 0.025% sodium dodecyl
   sulphate, 1% B-mercaptoethanol) before addition of ∼200 µl of 0.1mm
   Zirconia Silica beads (Biospec 11079101Z). Cells were then lysed by
   bead beating for 3× 2.5 minutes in a Biospec Mini Bead Beater (Biospec
   1001) with 5-minute intervals and then subjected to centrifugation at
   4,300 g for 5 minutes. The clean supernatant was then transferred to a
   new deep well plate and RNA was purified using a Zymo ZR-96 RNA Clean &
   Concentrator kit (Zymo R1080) per manufacturer’s instructions.

RNA-seq Library Preparation and Sequencing

   RNA-seq libraries were constructed following a modified RNAtag-seq
   protocol as detailed in Ref^[258]24. Briefly, 400 ng of total RNA
   lysate was subjected to fragmentation in 2X FastAP buffer (ThermoFisher
   EF0651), genomic DNA removal (TURBO^™ DNase, ThermoFisher AM2239), and
   dephosphorylation (FastAP, ThermoFisher EF0651). Fragmented RNA was
   purified using SeraPure SPRI bead cleanup^[259]69 and ligated with
   barcoded first adapter ligation by T4 RNA ligase I (NEB M0437M). After
   pooling and purification with the Zymo RNA Clean and Concentrator-5 kit
   (Zymo R1015), we quantified barcoded RNA using a Qubit RNA HS Assay Kit
   (ThermoFisher [260]Q32855) and then performed RNase H-based ribosomal
   RNA depletion on 400 ng of barcoded RNA sample using a 10:1
   probe-to-RNA ratio, as previously outlined^[261]24.

   rRNA-depleted RNA was subjected to downstream library preparation
   following standard RNAtag-seq protocols^[262]27, including reverse
   transcription (ThermoFisher 18090010) and second adapter ligation (NEB
   M0437M). Ligation products were further amplified with primers
   containing Illumina P5 and P7 adapters and sample indexes, and PCR
   reactions were stopped during exponential amplification. PCR products
   were subjected to gel electrophoresis on E-Gel^™ EX Agarose Gels, 2%
   (ThermoFisher G402002) and expected DNA smears (300–600 bp) were
   excised and extracted using the Monarch^™ DNA gel Extraction Kit (NEB
   T1020L). Resulting libraries were sequenced to a minimum of 4.9X
   ([263]Supplementary Table 6). Sequences of all adapters and primers
   used in library preparation are provided in Ref^[264]24.

RNA-seq Data Analyses

   RNA libraries were analyzed for differential expression analysis as
   outlined in Ref^[265]24, which lists all adapter and primer sequences
   used. Briefly, raw sequencing reads were demultiplexed using
   Sabre^[266]70 and bcl2fastq^[267]71 prior to adapter removal with
   Cutadapt v2.1^[268]72, using parameters ‘-a file:[RNAtagSeq adapter.fa]
   -u 5 –minimum-length 20 –max-n 0 -q 20’ to remove low-quality bases and
   adapters. To mitigate the effect of rRNA reads, we performed alignments
   against the 16S and 23S rRNAs of corresponding strains with
   Bowtie2^[269]73, using the versions and parameters outlined in
   ^[270]24. Genomes of natural isolates used for sequencing alignments
   were de novo assembled as described in^[271]24. Genome assemblies for
   all publicly available isolates were downloaded from the NCBI database,
   with all accession numbers listed in ([272]Supplementary Table 3). The
   number of reads uniquely mapped to each coding sequence was calculated
   using featureCounts v1.6.2^[273]74 without restraint on strandness (-s
   0), and the expression level of each CDS by transcripts per million
   (TPM) was quantified using an in-house script.

   Finally, we used DESeq2^[274]75 to perform differential expression
   analysis on the number of reads uniquely mapped to each coding
   sequence, calculating the FC of gene expression after perturbations.
   Differential gene expression was defined as >2 fold-change (FC) in gene
   expression relative to a vehicle control, with p[adj] <0.05, and
   average fragments per kilobase per million reads aligning to annotated
   open reading frame (FPKMO) > 0.10. P-values of differential expression
   were all adjusted by the Benjamini–Hochberg procedure in DESeq2 using
   default settings. Identified DEGs were used to identify drug-enriched
   modules and pathways using the Kyoto Encyclopedia of Genes and Genomes
   (KEGG) database^[275]76. KEGG enrichment with p[adj]<0.05 and e<.00001
   was considered statistically significant. All analyses were visualized
   in R^[276]77.

Growth Assays

   For all growth experiments, overnight bacterial cultures were
   back-diluted 1:100 before addition to 96-well plates containing 5–25 µl
   of pre-reduced chemicals or bile acids to reach 2, 20, or 100 µM
   concentrations in 1 ml of solution (all concentrations for
   bacterial-drug pairings listed in [277]Supplementary Fig. 2). As an
   exception, higher concentrations of 50, 500 and 2500 µM were used for
   metformin ([278]Methods, [279]Supplementary Fig. 2), which has much
   higher predicted drug concentrations in the intestinal tract^[280]5.
   Further, A. shahii was not included in growth experiments it exhibits
   characteristically poor growth in liquid that prevents optical density
   measurements^[281]78. Cultures were then incubated at 37°C in shaken
   media for 24 hours. For 24-hour temporal growth screens, 190 µl of
   culture was harvested at 6, 10, and 24 hours for measurement of optical
   density at 600 nm (OD^600) using an Epoch2 Microplate Spectrophotometer
   (Agilent Technologies). For all MIC tests, cultures were harvested for
   visual inspection and OD^600 measurement at 24 hours only. Growth
   inhibition was defined in all growth experiments as relative growth
   depletion >30%, with false discovery rate (FDR)-adjusted p value
   (p[adj]) <0.05. To determine statistical differences in relative growth
   between conditions, two-sided independent t-tests with Benjamini &
   Hochberg correction were performed in R to determine FDR-adjusted
   p-values unless otherwise stated.

Analysis of Public Datasets

   To identify species harboring the acrAB-tolC operon, we generated an
   amino acid reference from the P. vulgatus acrAB-tolC operon detailed
   in^[282]57. We then performed a homolog search for this reference
   operon in a public dataset of 4930 representative metagenome-assembled
   genomes (MAGs) characterized in Ref^[283]63. Briefly, we annotated all
   putative protein sequences in the MAG dataset using Prokka^[284]79
   before performing a homolog search with blastP^[285]80. Gene targets
   with an e-value < 10^−5, coverage > 90% and identity > 40% were
   considered hits. We then estimated abundances of all species containing
   the acrAB-tolC operon using metadata from Ref^[286]63. All associated
   p-values were calculated by two-tailed independent t-test in R.

   To examine whether metoprolol administration could be linked to the
   upregulation of the ttd operon in vivo, raw amplicon sequencing data
   from a public dataset (n=145) were downloaded from Sequence Read
   Archive (SRA) under accession code ERP002469. Single-ended raw reads
   were processed by Cutadapt v2.1^[287]72 with the following parameters
   “--minimum-length 24 -u 10 --trim-n -q 15” to remove low-quality bases
   and Illumina adapters. Reads passing quality filtering were then
   aligned against our ttd operon reference by Bowtie2 v2.3.4^[288]73 in
   --very-sensitive mode. Read counts per gene were normalized by gene
   length and sequencing depth (i.e., reads per kilobase per million
   mapped reads) for expression-level quantification. All associated
   p-values were calculated using the Mann-Whitney test in R.

   To examine whether simvastatin-modulated microbiota shifts could be
   linked to acrAB-tolC operon prevalence in vivo, raw amplicon sequencing
   data from the cross-sectional Meta-Cardis Body Mass Index Spectrum
   (BMIS) cohort (n = 888) were downloaded from the EMBL-EBI European
   Nucleotide Archive (ENA) under accession number PRJEB37249.
   Single-ended raw reads were processed by Cutadapt v2.1^[289]72 with the
   following parameters “--minimum-length 24 -u 10 --trim-n -q 15” to
   remove low-quality bases and Illumina adapters. Reads passing quality
   filtering were then aligned against our acrAB-tolC operon reference by
   Bowtie2 v2.3.4^[290]73 in --very-sensitive mode. Read counts per gene
   were normalized by gene length and sequencing depth (i.e., reads per
   kilobase per million mapped reads) for expression level quantification.
   All associated p-values were calculated by two-tailed independent
   t-test in R.

AcrAB-TolC Engineering Studies

   All Bacteroides expression vectors were generated using Gibson assembly
   (NEBuilder 2x HiFi DNA assembly master mix), and polymerase chain
   reaction (PCR) fragments for cloning were generated using Q5 DNA
   Polymerase (NEB). We first generated two PCR fragments from a plasmid
   designed and constructed in our laboratory containing a constitutive
   Bacteroidales promoters described in Ref^[291]81 (additional features
   of this custom vector backbone are described in a separate manuscript
   currently under consideration) to drive expression of the AcrAB-TolC
   gene copies. We then cloned AcrAB-TolC copy 1 (primers
   TCTCGTCAAACAAATATAAATAATATAAACATGAAAATGACAGTAAATAGTATGAAATGT and
   AGAAGGGCACCAATAACTGCCTTAAAAAAATTAATTATTCACGTCCACCGC) and copy 2
   (primers TCTCGTCAAACAAATATAAATAATATAAACATGAAATTTTATTGCAAACCTACGT and
   AGAAGGGCACCAATAACTGCCTTAAAAAAATTATTCTTTTTTGCCTTTGGTCATC) from the P.
   vulgatus H1 genome using PCR. A 3 fragment Gibson assembly was
   incubated at 50 °C for 1 hour to generate our plasmid construct. As a
   control, we cloned a tetracycline (Tet) resistance gene (described in
   Ref^[292]82) into our plasmid construct in place of AcrAB-TolC cargo.
   Finally, AcrAB-TolC-containing plasmids (or Tet plasmid controls) were
   transformed into chemically competent NEB turbo cells. Transformed
   colonies were screened by PCR for the correct insert size length and
   whole plasmids were sequenced on an Illumina Nextseq 500/550 platform
   or using Plasmidsaurus (Eugene, OR). Vector constructs with the correct
   payload were used in subsequent conjugation experiments.

   Before conjugation, donor strains harboring conjugative Bacteroides
   expression vectors were grown from a single colony in 5 ml of LB-Lennox
   media (BD) supplemented with 50 µg ml−1 of carbenicillin and 50 µM DAP
   at 37 °C overnight (~10–16 h). We prepared donor and recipient P.
   vulgatus H1 cells and carried out conjugations as previously described
   in Ref^[293]83 under anaerobic conditions. Transconjugant colonies were
   Sanger sequenced at the 16Sv4 region to confirm strain identity, and
   stable plasmid maintenance was confirmed by colony PCR using the
   primers listed above. Positive P. vulgatus transconjugants were then
   picked into 5 mL ½ mGAM supplemented with 20 µg/ml erythromycin and
   grown overnight. These strains were banked in glycerol stocks (25%
   glycerol final concentration) and stored at −80 °C. Subsequent growth
   curve experiments were done using overnight cultures grown from these
   glycerol stocks in ½ mGAM.

Conspecific Strain Gene Mapping

   To compare the transcriptomic response for strains of the same species,
   we performed gene mapping for five pairs of conspecific strains used in
   this study. Briefly, all protein sequences of the strains were first
   annotated using Prokka^[294]79. Protein alignment was then performed
   for the protein sequences of conspecific strains by blastP^[295]80
   using the parameters “-max_target_seqs 20 -evalue 0.001”. Hits with
   identity > 0.75 and coverage > 0.75 were considered as mapped genes for
   conspecific strains.

Statistics and Reproducibility

   For all isolates generated in this study, individual of origins for
   isolated gut strains were assigned based on the defined identity of
   original feces ([296]Supplementary Table 2), for which covariate
   analysis is not applicable. The specific number of bacterial strains
   used in this study (19) was chosen to allow for inclusion of major
   bacterial phyla within an individual donor microbiome, with replicates
   included as described in [297]Methods to allow for statistical
   significance calculations. For clinical datasets analyzed, cohort
   assignment for individuals was taken from the metadata of the original
   study. Blinding was not possible during experiments as we were
   comparing transcriptional responses of different bacterial isolates to
   different drug conditions. Transcriptomic and growth processing within
   species was blinded as different drug conditions were processed
   together using pooled methods. Data exclusion was based on sequencing
   coverage or genome quality to remove technical artifacts as described
   in the [298]Methods section. The sample sizes in this study are either
   number of bacterial strains tested (with replicates, see below), or
   number of experimental or control cohort individuals taken from a
   publicly available dataset (for which sample size calculations were
   already performed). All sample sizes are listed in figure legends and
   [299]methods sections where applicable. All in vitro assays on
   bacterial strains were performed with multiple (2–4) technical
   replicates, as noted in the [300]methods section and figure legends
   where applicable, in order to confirm replicability and enable
   statistical significance calculations. All analyses of associated data
   were performed with the same parameters and criteria described in
   [301]Methods section above. Data distribution was assumed to be normal,
   but this was not formally tested.

Supplementary Material

   Suppl Matter
   [302]NIHMS1985061-supplement-Suppl_Matter.pdf^ (4.9MB, pdf)
   Supplementary Tables
   [303]NIHMS1985061-supplement-Supplementary_Tables.xlsx^ (21.3MB, xlsx)

ACKNOWLEDGEMENTS