Machine Learning of Bacterial Transcriptomes Reveals Responses Underlying Differential Antibiotic Susceptibility

Prior studies revealed interrelationships between medium composition and bacterial antibiotic susceptibility by contrasting drug activity profiles under multiple medium conditions. To investigate the basis for medium-dependent antibiotic activities, we first defined the baseline effect of medium composition on E. coli K-12 MG1655 gene expression by performing RNA-seq during growth in CA-MHB, M9, or R10LB. We identified 1,075 differentially expressed genes (DEGs) between the physiological medium CA-MHB and nutrient-poor M9 minimal medium, comprising nearly one-quarter of the E. coli genome (Fig. 1a). These DEGs spanned every Cluster of Orthologous Groups (COG) category, with 24% of DEGs of unknown function. R10LB yielded transcriptional responses more similar to CA-MHB than M9, since only 556 genes were differentially expressed between the two media (Fig. 1a).

For context, the RNA-seq data set generated in this study was compared to our previously published high-quality compendium containing 278 E. coli expression profiles, PRECISE (12), which primarily contains samples grown in M9 minimal medium or LB. Principal-component analysis (PCA) of the combined data set showed that CA-MHB and R10LB provoked significantly different responses from those probed in PRECISE, including growth on LB media without RPMI (Fig. 1b). As a control, gene expression in glucose M9 minimal medium was similar to expression profiles on identical medium produced previously (20) (Pearson R = 0.92), highlighting the internal consistency of the RNA-seq analysis. This earlier expression profile grown on glucose M9 minimal medium served as the reference condition, to which all further iModulon activities are compared.

We combined the 278 data sets in PRECISE with the 30 new RNA-seq data sets from this work to form a new compendium of 308 profiles (see the supplemental material). Application of ICA identified 98 iModulons, 88 of which were nearly identical to those extracted from the original PRECISE data set (12) (see Fig. S1). These iModulons are presented in iModulonDB (iModulonDB.org), an interactive visualization tool (21). To dissect the transcriptional differences between the three medium formulations, we computed the differential iModulon activities (DiMAs) between each pair of media. DiMAs are iModulons that have large, statistically significant (false discovery rate [FDR] > 0.1) activity changes between two conditions of interest. DiMAs are akin to coarse-grained DEGs but can be directly traced back to specific transcriptional regulators. In contrast to the hundreds of DEGs between the medium conditions, a maximum of 37 iModulons were perturbed between these conditions, simplifying further analysis by an order of magnitude (Fig. 1c).

On average, iModulons could explain 81% of the difference in expression between each pair of medium conditions. To further coarse grain our analysis, we grouped iModulons together based on regulator function (see Table S1 in the supplemental material). Borrowing from multiple linear regression, we calculated how much of the expression deviation could be directly explained by each iModulon and each iModulon group (Fig. 1d to f). The single largest contributor to expression variation between CA-MHB and M9 was the RpoS iModulon, which accounted for 30% of the variation in expression. Conversely, the top contributors to expression variation between CA-MHB and R10LB were the Fur-1, PurR-2, FlhDC, NarL, and Fnr iModulons. In total, these five iModulons explained 34% of the variance between CA-MHB and R10LB (see Table S2).

The major source of expression variation between rich media (CA-MHB and R10LB) and minimal medium (M9) could be explained by the “fear-greed” trade-off, which balances fast growth (i.e., “greed”) and high stress response readiness (i.e., “fear”) (12, 22) (Fig. 2a). In rich media, bacterial cells experience significantly less environmental stress, leading to a strong decrease in the activity level of stress-related iModulons, such as RpoS and the GadEWX acid-response system (23). This reduced stress in turn frees cellular resources that can be allocated to promote growth, such as increased activity levels of iModulons encoding genes related to protein translation. These two processes are molecularly linked through (p)ppGpp, a signaling molecule that is produced when the cell experiences nutrient limitations or other stresses. ppGpp binds to RNA polymerase, changing its conformation to activate genes under the control of the RpoS stress response sigma factor and decrease transcription of genes related to translation (24, 25).

A mutated RNA polymerase (RpoBE672K) observed during adaptive laboratory evolution on glucose minimal media (22) is hypothesized to interfere with the ppGpp binding site on RNA polymerase. The iModulon activities of the WT strain under rich media are similar to the iModulon activities of the RNA polymerase mutant strain; both exhibit increased expression of the genes in the translation and ppGpp iModulons and decreased expression of genes in the stress-related RpoS and GadEWX iModulons. It is important to note that the direction of iModulon activities changes reflect the expression changes in iModulon genes (i.e., both derepression and binding of activators results in higher iModulon activities). Overall, these activities indicate that, as expected, minimal media induce more (p)ppGpp stress alarmone production than rich media.

Growth rate-dependent changes in gene expression have been observed to potentiate antibiotic resistance and tolerance (26). Using the growth rate information from PRECISE, we found that growth rate effects best explained the RpoS and translation iModulon activities (R² = 0.39 and 0.30, respectively). However, growth effects on biosynthesis and energy metabolism iModulon activities were more muted (see Fig. S2a), indicating that the fear-greed iModulons are the primary targets of growth rate-dependent expression changes.

Many iModulons related to carbon uptake were differentially activated between the three growth media, including those related to the global regulators CRP and Cra (Fig. 2b). CRP represses secondary carbon sources in the presence of glucose and is activated by cAMP, whereas Cra represses glycolysis and is deactivated by fructose derivatives (20). iModulons regulating alternative carbon sources, such as ribose (RbsR), glycerol (GlpR), and pyruvate, were active in both CA-MHB and R10LB, even though glucose is a component of RPMI. The Cra iModulon contains many glycolytic genes and was active in both RPMI and CA-MHB. iModulon activities for M9 media with fructose as a carbon source are shown in Fig. 2b, since many carbon uptake iModulon activities for this condition were similar to CA-MHB and R10LB.

Biosynthetic pathways for vitamins and amino acids contributed to 15 and 19% of the expression deviation from M9 for CA-MHB and R10LB, respectively (Fig. 2c). Since R10LB contains 10% LB, we included the LB activities for comparison. The iModulon activities of amino acid and B-vitamin iModulons are quite similar in CA-MHB, R10LB, and LB, suggesting a similar nutritional content. The main deviations between CA-MHB and R10LB appear to be in histidine and branched-chain amino acid biosynthesis that are regulated by tRNA-mediated transcriptional attenuation.

The PurR-2 iModulon, which controls pyrimidine biosynthesis, exhibited significantly lower activities in CA-MHB than in R10LB. Lower PurR iModulon activities are associated with nucleotide supplementation, indicating higher pyrimidine content in CA-MHB than in R10LB (Fig. 2d). On the other hand, de novo purine biosynthesis pathway was downregulated in both CA-MHB and R10LB relative to M9 medium.

iModulons related to anaerobic respiration were another large contributor to expression deviation in CA-MHB compared to the other two media (Fig. 2e). Three iModulons were enriched with the transcription factor Fnr: one iModulon uniquely enriched with the Fnr regulon, one iModulon enriched with genes coregulated by Fnr and the nitrate-responsive regulator NarL, and a final iModulon containing all genes coregulated by Fnr, NarL, and the nickel regulator NikR. An additional iModulon in this group was related to ArcAB, a two-component system that senses the redox state of quinone electron carriers to regulate the tricarboxylic acid (TCA) cycle among other genes (27, 28). ArcA iModulon activities were partially reduced in both CA-MHB and R10LB.

Iron regulation accounted for a large portion of the expression deviation between CA-MHB and R10LB, as represented by two Fur-enriched iModulons and the FecI iModulon (Fig. 2f). Fur binds to free iron (FeII) and represses genes related to iron siderophore synthesis and transport (29). Low iModulon activity, as seen in R10LB, indicates repression of Fur-regulated genes. FecI senses iron (III) and activates the expression of iron (III) uptake genes.

A few additional iModulons were differentially activated across the multiple medium conditions, and coincidently discriminate the biological processes that are perturbed between the reference M9 expression profile in PRECISE and the expression profiles from M9 media generated in this study. Three iModulons related to the Rcs-phosphorelay were activated in the new M9 and CA-MHB expression profiles (see Fig. S2b). The FliA iModulon, which controls chemotaxis, was activated in all three media, whereas the FlhDC iModulon, which controls biosynthesis and assembly of flagella, was only activated in R10LB (see Fig. S2c).

Although the three media induced vastly different gene expression changes, we used iModulons to simplify the analysis from hundreds of differentially expressed genes to a few dozen key regulators. This interpretation provided a tractable understanding of how medium components affect transcriptional regulation and laid the groundwork for investigating responses to antibiotic treatment.

To examine how bacterial transcriptional responses impact antibiotic susceptibility, we conducted standardized AST for several antibiotics against E. coli MG1655 in CA-MHB, M9, and R10LB (Table 1). Briefly, E. coli was grown overnight in each of the three media, subcultured the following morning, and grown until reaching the logarithmic phase (∼0.4 optical density at 600 nm [OD₆₀₀]). Logarithmic phase cultures were diluted to an OD₆₀₀ of 0.002 and exposed to serial dilutions of the nine selected antibiotics (see Materials and Methods) (6). After 20 h, the cell density was measured to identify the amount of drug required to inhibit ≥90% of the growth of the untreated controls. In three cases, cell growth was uninhibited in all drug concentrations (e.g., ampicillin [AMP] on M9); in one case, the MIC could not be determined (rifampin [RIF] on R10LB).

Based on clear medium-dependent alterations in activity, we selected four antibiotics for further investigation. Two selected antibiotics, ceftriaxone (CTR) and trimethoprim-sulfamethoxazole (T/S), had identical susceptibilities in the rich media (CA-MHB and R10LB) and a lower MIC₉₀ in M9. Ciprofloxacin was selected as it was most effective in R10LB and least active in CA-MHB. Meropenem (MEM) resulted in the opposite trend in susceptibility, where it was most effective in CA-MHB and least active in R10LB.

According to the standardized AST as set forth by the Clinical and Laboratory Standards Institute (CLSI) (30), MIC₉₀ was measured using an inoculum density of ∼1 × 10⁵ CFU/ml (Table 1). However, we required a significantly higher cell density (∼2.5 × 10⁷ CFU/ml) to collect enough mRNA for transcriptomic analysis. Since antibiotic susceptibility can be density dependent, we reassessed antibiotic bactericidal activity at the same culture density used for the mRNA isolations (Fig. 3a). From these data, we determined the minimum bactericidal concentration (MBC₉₉) for each antibiotic in each medium type (Fig. 3b). At the higher culture densities, no differences were observed for ciprofloxacin (CIP) and CTR bactericidal activities in CA-MHB or R10LB, whereas T/S was more effective in CA-MHB than in R10LB. CIP, CTR, and T/S had the lowest bactericidal activity in M9. Despite being classified as bactericidal, no tested dose of MEM was able to eradicate MG1655 at these higher culture densities (see Fig. S3a). Corroborating the MBC₉₉ results, MICs at the higher densities used for the RNA-seq analysis were likewise increased ∼5-fold from baseline MICs (see Fig. S3b).

Since the ICA-based decomposition of transcriptomes provided a clear explanation of the complex responses to different medium compositions, we applied a similar approach to investigate the effects of antibiotic treatment. The MIC₉₀ of each antibiotic in each medium was selected as the experimental treatment dose, as indicated by the black bars in Fig. 3a. Cells were treated with antibiotics during exponential growth and sampled for RNA sequencing 30 min after treatment.

Gene expression after antibiotic treatment primarily clustered by medium, even when controlling for medium-specific responses (see Fig. S4a and b). Although MEM treatment did not kill MG1655 at any concentration (see Fig. S3a), it still induced small changes in the transcriptome (see Fig. S4c). However, these effects were not reflected in the MBC₉₉, and we removed MEM-treated transcriptomes from subsequent analyses.

We first compared medium-specific changes in iModulon activity between antibiotic treatment versus untreated control in each of the three growth media (Fig. 4a and Table S3). Regardless of their mechanisms of action, several trends in iModulon activities were observed across all three antibiotics. Antibiotic treatment increased activity of the translation iModulon in M9 and R10LB (Fig. 4b) but had no consistent effect on the RpoS iModulon across the three antibiotic classes. Although none of the selected antibiotics interact directly with the ribosome, this indicated that antibiotic-mediated stress could indirectly increase ribosome levels in the cell.

We also observed large reductions in both Crp-related iModulon activities upon antibiotic treatment in rich media (Fig. 4c). These iModulons encode genes in secondary carbon source catabolism; the Crp-1 iModulon contains scavenging enzymes such as tryptophanase and acetyl coenzyme A synthetase, whereas the Crp-2 iModulon represents expression of phosphotransferase systems for alternative carbon sources. Lower activity of Crp iModulons indicated a reallocation of resources away from secondary carbon source catabolism. Traditional gene set enrichment analysis of the DEGs induced by antibiotic treatment could not identify Crp as a key regulator of the process, likely due to the size and breadth of the Crp regulon.

Across PRECISE (12), we observed a strong negative correlation between the iModulon activities of two major respiration-related transcription factors, Fnr and ArcA (Pearson R = −0.74, P < 10⁻¹⁰; Fig. 4d). In both CA-MHB and R10LB, antibiotic treatment resulted in significant drops in Fnr iModulon activities, coupled with a simultaneous increase in ArcA-regulated gene expression in the TCA cycle. Fnr contains an iron-sulfur cluster that is destabilized in the presence of molecular oxygen (31), rendering Fnr inactive. We hypothesized that the high expression of these genes in R10LB indicated a high respiration rate that left little intracellular molecular oxygen available to destabilize Fnr’s iron-sulfur cluster. To test this hypothesis, we measured oxygen consumption rates in CA-MHB and R10LB using an assay that is based on the ability of oxygen to quench the excited state of extracellular O₂ consumption reagent (Fig. 4e). The normalized fluorescence is proportional to the depletion of oxygen concentration in the media (32, 33). Here lies one of the major differences between the two media; R10LB promotes a higher respiration rate than CA-MHB before antibiotic treatment.

Together, these two iModulons indicate that antibiotic treatment decelerated the high respiratory rate facilitated by R10LB. This process likely also reduced the redox state, as sensed by the ArcAB response system (27). A prior study used biochemical assays to investigate the effects of bacteriostatic antibiotics on E. coli (34) and observed that antibiotic treatment drastically reduced the respiration rate and moderately reduced the redox state of the cell. The drop in respiration was absent in minimal media, where we inferred a lower initial cellular respiration rate from the Fnr and ArcA iModulon activities.

Treatment of E. coli with three distinct antibiotics induced multiple common transcriptional changes, including (i) an increase in ribosomal gene expression, (ii) downregulation of secondary carbon catabolism, and (iii) a decelerated respiration rate in R10LB and CA-MHB. In the next section, we discuss transcriptional changes unique to each antibiotic.

In addition to the medium-specific iModulon activity trade-offs described above, each antibiotic induced specific changes in iModulon activity that could be correlated with differential antibiotic susceptibility.

T/S is a combination of two drugs, both of which inhibit folate synthesis. E. coli was most susceptible to T/S on M9 media and less sensitive on R10LB and CA-MHB. Folate is used to synthesize purines, thymidylate, and methionine (35). However, the purine biosynthesis (PurR-1) and methionine biosynthesis (MetJ) iModulons were mostly unaffected by T/S treatment in any medium type, indicating that E. coli is not transcriptionally altering folate usage in response to T/S inhibition (see Fig. S5). Instead, we observed significant iModulon activity changes in four other amino acid and nucleotide-related iModulons (Fig. 5a). Sulfate utilization (CysB) and arginine biosynthesis (ArgR) iModulon activities decreased after treatment with T/S in M9 media. The glycine cleavage system (GcvA iModulon), which converts excess glycine to serine, was downregulated after treatment in R10LB, and pyrimidine biosynthesis (PurR-2) was downregulated after treatment in both M9 and R10LB. Although the MIC₉₀ of T/S was identical in CA-MHB and R10LB, the underlying transcriptional responses to the aforementioned iModulons were different, indicating that these responses are not directly linked to T/S susceptibility.

The MIC₉₀ of CTR for E. coli was also identical between CA-MHB and R10LB and higher for M9. Treatment with CTR in M9 media induced a medium-specific decrease in CysB iModulon activity (Fig. 5b). However, the PuuR iModulon increased after CTR treatment in all three media types, indicating a medium-independent response to CTR. PuuR responds to the polyamine putrescine and regulates its utilization and transport (36). Polyamine accumulation from sublethal antibiotic treatment has been observed previously and is hypothesized to reduce porin permeability to antibiotics (37).

The MIC₉₀ of CIP was highest on CA-MHB and lowest on R10LB. CIP blocks DNA gyrase, creating double-stranded DNA breaks (38). Subsequently, after adding CIP to the culture, we observed a strong activity increase in the SOS-response regulator LexA iModulon in all three media indicating a media-independent response to CIP (Fig. 5c). CIP treatment also caused medium-specific effects on the two Fur iModulons (Fig. 5d and e; see also Fig. S6). The Fur-1 iModulon responds to iron starvation, derepressing iron (II) and (III) siderophore synthesis and transport systems, ribonucleotide reductases, superoxide dismutase, and iron-sulfur cluster assembly. The Fur-2 iModulon responds to excess iron, further repressing siderophore transport and repressing the energy-transducing Ton system. Together, these two iModulons trace out a trajectory from iron replete (Fur-repressed) to iron-poor (no Fur repression) conditions. CIP treatment in CA-MHB results in lower Fur-2 iModulon activities, indicating excess free iron availability. On the other hand, CIP treatment in R10LB, results in higher Fur-1 iModulon activities, indicating iron starvation. The difference in Fur iModulon activities may contribute to the differential susceptibility to CIP across the three media.

The difference between the Fur-1 and Fur-2 iModulons are reflected in the DEGs between CIP-treated and untreated cells (Fig. 5f). In CA-MHB, CIP treatment increases the activity of the Fur-1 iModulon, and most of the respective DEGs are specific to the Fur-1 iModulon. On the other hand, CIP treatment in R10LB decreases the Fur-2 iModulon activity, resulting in DEGs specific to the Fur-2 iModulon. Due to the two dimensions of the Fur regulon, traditional gene set enrichment analysis of the DEGs does not identify Fur as a statistically significant enriched regulator (FDR < 0.01) in either medium.

E. coli strain MG1655 was grown in three media: (i) M9 minimal medium (47.8 mM Na₂HPO₄, 22 mM KH₂PO₄, 8.6 mM NaCl, 18.7 mM NH₄Cl, 2 mM MgSO₄, and 0.1 mM CaCl₂) supplemented with 0.2% (wt/vol) glucose (M9); (ii) Roswell Park Memorial Institute 1640 (RPMI; Thermo Fisher Scientific) supplemented with 10% LB (R10LB); and (iii) Mueller-Hinton broth (Sigma-Aldrich) supplemented with 25 mg/liter Ca²⁺ and 12.5 mg/liter Mg²⁺ (CA-MHB). Glycerol stocks of E. coli were streaked onto Luria broth agar (LB agar) and grown at 37°C overnight.

Antibiotic susceptibility was determined as described previously (6). Briefly, E. coli-streaked plates were used to inoculate overnight cultures in the desired media. E. coli overnight cultures were used to inoculate pre-experimental cultures the following day. Pre-experimental cultures were grown to the midlogarithmic phase (OD₆₀₀= ∼0.4) and then diluted to an OD₆₀₀ of 0.002 (∼1 × 10⁵ CFU/ml) to create the experimental cultures. Concentrated antibiotic stocks were serially diluted in 1× DPBS (Corning) to create 10× stocks of each antibiotic at the indicated concentrations. Then, 180-μl portions of the experimental cultures were combined with 20 μl of the desired 10× antibiotic stock in Costar flat-bottom 96-well plates (Corning, catalog no. 3370) and incubated with shaking at 200 rpm at 37°C overnight. Bacterial growth (OD₆₀₀) was determined approximately 20 h later utilizing a Enspire Alpha multimode plate reader (Perkin-Elmer). To calculate the MIC₉₀, defined as the amount of drug required to inhibit ≥90% of the growth of the untreated controls, the density of each drug-treated well was compared to untreated control. MICs for high-density cultures were performed as described above, except with a starting OD₆₀₀of  ∼0.05 after dilution.

Mid-log-phase E. coli MG1655 was used to inoculate experimental cultures of either glucose M9, R10LB, or CA-MHB with approximately 2.5× 10⁷ CFU/ml (OD₆₀₀= ∼0.05). Antibiotics were added, and cultures were incubated as described for the antibiotic susceptibility assays. After 20 h, plates were removed from the incubator and wells were serially diluted 10-fold in their respective media from 10⁰ to 10⁻⁷. Then, 20 μl of each serial dilution was spot plated onto LB agar, followed by incubation at 37°C overnight to enumerate the CFU.

Total RNA was sampled from duplicate cultures. Strains were cultured overnight in the respective media as described above. Mid-log-phase cultures were diluted (starting OD₆₀₀ of ∼0.05) into 30 ml of medium for untreated conditions or 30 ml of medium containing the MIC₉₀ of the respective antibiotic for the respective media (Table 1). Flasks were incubated for 30 min at 37°C with shaking. Cell broth was centrifuged, and the supernatant was removed. RNA extraction and library preparation were performed as described previously (12). Raw sequencing reads were performed as described previously (12). Differential expression was performed using DESeq2 (53), with a log₂-fold change cutoff of 1.5 and q value cutoff of 0.05.

Log-transformed transcripts per million (log-TPM) expression levels were concatenated to the PRECISE data set. ICA and iModulon processing were performed as described previously (12). Briefly, we executed FastICA 100 times with random seeds and a convergence tolerance of 10⁻⁷. We constrained the number of independent components (ICs) in each iteration to the number of components that reconstruct 99% of the variance as calculated by PCA. The resulting ICs were clustered using DBSCAN to identify robust ICs, with the following parameters: an epsilon value of 0.1 and a minimum cluster seed size of 50. This process was repeated 10 times, and only ICs that consistently occurred in all runs were kept.

As described in Sastry et al. (12), iModulons were extracted from ICs by iteratively removing genes with the largest absolute value and computing the D’Agostino K² test statistic of the resulting distribution. Once the test statistic fell below a cutoff, which was identified through a sensitivity analysis (12), we designated the removed genes as the “iModulons.”

To associate iModulons to regulators, we compared the set of significant genes in each iModulon to each regulon (defined as the set of genes regulated by any given regulator) using the two-sided Fisher exact test (FDR < 10⁻⁵), as described by Sastry et al. (12). In addition, combined regulon enrichments were calculated to identify joint regulation of iModulons (such as Fnr+NarL), using both intersection (+) and union (/) of up to three regulons. Final iModulon-regulator associations were determined through manual curation of enriched regulators.

Gene set enrichment analysis for differentially expressed genes was performed as described above, using only single regulator enrichments (FDR < 0.01).

Explained variance between two conditions was calculated as follows:where k is the iModulon of interest. The total explained variance was calculated similarly:

The difference between the CEV and the sum of the explained variance is the confounding variance.

Differentially activated iModulons were computed with a similar process as previously detailed (12). For each iModulon, the average activity of the iModulon between biological replicates, if available, was computed. The absolute value of the difference in iModulon activities between the two conditions was then compared to the fitted log-normal distribution of all differences in activity for the iModulon. iModulons that had an absolute value of activity greater than 7 and an FDR below 0.1 were considered to be significant. The number of DIMAs was computed for each unique pair of conditions in the PRECISE 2.0 compendium.

The oxygen consumption assay was performed using Abcam’s extracellular oxygen consumption assay kit (ab197243). Assay was performed as per the manufacturer’s suggested protocol. Exponentially growing bacterial cultures were diluted to an OD₆₀₀ of ∼0.1, and 100 μl of the diluted cultures was added to a flat and transparent bottom 96-well sterile half-area cell culture plate. Dye was added as suggested in the manufacturer’s protocol, and wells were sealed with mineral oil. The time-resolved fluorescence was recorded in a kinetic cycle with 2 min interval with a 30-μs excitation-emission time lag and a 100-μs integration time using a Tecan infinite 200Pro. Plates were incubated at 37°C in a shaking mode (2-mm orbital amplitude). The OD₆₀₀ was recorded in parallel to normalize for the cell number. Three independent replicates were used for every condition.

All raw data are available on NCBI Gene Expression Omnibus (GEO) under accession number GSE159494. Processed data are available in the supplemental material. The code to compute iModulons is publicly available on Github at https://github.com/SBRG/precise-db.