Research ArticleDrug Discovery

High-throughput metabolomic analysis predicts mode of action of uncharacterized antimicrobial compounds

See allHide authors and affiliations

Science Translational Medicine  21 Feb 2018:
Vol. 10, Issue 429, eaal3973
DOI: 10.1126/scitranslmed.aal3973

Boosting the antibiotics pipeline

Antibiotic discovery requires innovative phenotypic assays to identify the mode of action of compounds during large-scale drug screening. Now, Zampieri et al. have developed a high-throughput method to systematically quantify and interpret dynamic metabolome responses of mycobacteria to new antimicrobial compounds. They demonstrate how one can infer the mode of action of uncharacterized antimycobacterial compounds by comparing drug-induced metabolic responses in Mycobacterium smegmatis. Their combined mass spectrometry–based metabolomics and computational workflow scaled with the size of typical compound libraries and hence could facilitate the selection of antimicrobial compound leads with potentially new modes of action.


Rapidly spreading antibiotic resistance and the low discovery rate of new antimicrobial compounds demand more effective strategies for early drug discovery. One bottleneck in the drug discovery pipeline is the identification of the modes of action (MoAs) of new compounds. We have developed a rapid systematic metabolome profiling strategy to classify the MoAs of bioactive compounds. The method predicted MoA-specific metabolic responses in the nonpathogenic bacterium Mycobacterium smegmatis after treatment with 62 reference compounds with known MoAs and different metabolic and nonmetabolic targets. We then analyzed a library of 212 new antimycobacterial compounds with unknown MoAs from a drug discovery effort by the pharmaceutical company GlaxoSmithKline (GSK). More than 70% of these new compounds induced metabolic responses in M. smegmatis indicative of known MoAs, seven of which were experimentally validated. Only 8% (16) of the compounds appeared to target unconventional cellular processes, illustrating the difficulty in discovering new antibiotics with different MoAs among compounds used as monotherapies. For six of the GSK compounds with potentially new MoAs, the metabolome profiles suggested their ability to interfere with trehalose and lipid metabolism. This was supported by whole-genome sequencing of spontaneous drug-resistant mutants of the pathogen Mycobacterium tuberculosis and in vitro compound-proteome interaction analysis for one of these compounds. Our compendium of drug-metabolome profiles can be used to rapidly query the MoAs of uncharacterized antimicrobial compounds and should be a useful resource for the drug discovery community.


Despite rapid technological progress, the discovery of new antibiotic drugs has been stalled for the past 50 years. To combat the growing burden of antibiotic resistance, innovative drug discovery approaches are required to improve and expedite the antibiotic discovery process (1, 2). To overcome the intrinsic difficulties of target-based screens and to move beyond well-known targets, phenotype-based screening methods are a valid alternative (3). In a classical phenotypic screening approach, promising antimicrobial compounds are selected on the basis of their empirical ability to prevent cell growth in vitro (4). However, the lack of efficient, rapid, and systematic methods to investigate the modes of action (MoAs) of newly discovered cell growth inhibitors has often misguided the selection of the most promising lead compounds (5).

Current systematic inference of MoAs in large chemical libraries is based on large-scale in vitro direct biochemical techniques (2), comparison of chemical structures (6), or functional genomics approaches. These include in vitro resistance mutation mapping (7), cytological profiling (8), monitoring of transcriptional responses (9, 10), and phenotypic assays (1114) to probe interactions with other drugs (15, 16), gene deletion/overexpression (1721), or culture media (22). Although such strategies have been successfully used (2225), their applicability and scalability are often hampered by cost, reliability, and their labor-intensive nature (26, 27).

Emerging applications of metabolomics in drug discovery have focused on biomarker identification and therapeutic monitoring, but only a few are aimed at determining drug MoAs (2835). Large-scale applications of classical metabolomics techniques are hampered by limited coverage and laborious sample preparation. To address these limitations and to enable systematic testing of large libraries of bioactive compounds in vivo, we developed an approach based on nontargeted mass spectrometry in combination with ad hoc data mining to enable high-throughput inference of drug MoAs. We first validated our approach on a set of 62 reference compounds with known MoAs applied to the nonpathogenic bacterium Mycobacterium smegmatis and then used this technique to infer MoAs of a library of 212 uncharacterized antimycobacterial compounds.


Dynamic metabolome profiling after treatment with reference compounds

We built a reference base of metabolic responses from 62 reference compounds, including currently used antibiotics and chemical stress agents, with 17 known MoAs (Fig. 1, A to C, and table S1). These MoAs included inhibition of bacterial DNA replication, protein synthesis, cell wall synthesis, and folate synthesis. Reference compounds were administered at three to four concentrations from subinhibition to full inhibition of exponentially growing isogenic M. smegmatis in culture at an optical density (OD595) of 0.4 in 96-well deep plates (table S2 and Fig. 1A). The dynamic metabolome response was assessed by withdrawing 80 μl of bacterial culture 5, 30, 60, 180, 360, and 600 min after exposure to the compounds. To facilitate dynamic sampling at high throughput and reduce the risk of sample processing artifacts, the whole culture broth was extracted without cell separation using cold solvent extraction and then directly injected into a time-of-flight mass spectrometer (Fig. 1A) (36). This procedure enabled rapid profiling of relative dynamic changes of ~15,000 ions, out of which 1006 could be putatively annotated as deprotonated metabolites. In total, we monitored dynamic metabolite changes across 431 compound-treated conditions and two mock treatments [vehicle controls: H2O or dimethyl sulfoxide (DMSO)] over seven time points in three biological replicates (Fig. 1B).

Fig. 1 Antibiotic-induced metabolome responses in M. smegmatis.

(A) Metabolomics workflow. Cells were grown in 700-μl volumes in 96-well plates to an OD595 of about 0.4, before addition of 10 μl of the antimicrobial compound. Cell culture (80 μl) was withdrawn from each well at each sampling time. Forty microliters was used to determine cell density, and the remaining 40 μl was added to cold extraction buffer. Supernatant was directly injected into a time-of-flight mass spectrometer, and relative changes in metabolite intensities were extrapolated from processing of the metabolome data. (B) Compounds tested. Almost half of the compounds tested included different concentrations of reference antimicrobials (yellow) and chemical stress agents (green) with known MoAs; the remainder were compounds from a GSK library used at 10 μM concentration (blue). (C) Distribution of MoAs for the 62 reference compounds. ATP, adenosine 5′-triphosphate. (D) Schematic representation of the drug-metabolome response data set. For each antimicrobial compound tested, the dynamic profile of 1006 metabolites was interrogated. As an example, the top graph illustrates the response of the folic acid biosynthesis intermediate 4-aminobenzoic acid to the antimicrobial para-aminosalicylic acid (PAS) (red, 104 μM; gray, 41 μM; blue, 25 μM). The bottom graph shows the response of the bacterial metabolite mycobactin to the known antimicrobial isoniazid (red, 1.5 mM; gray, 0.22 mM; blue, 0.11 mM). Thick lines represent the results from the impulse model fitting analysis for the three drug concentrations. Metabolic profiles of 4-aminobenzoic acid and mycobactin across all conditions are shown in light gray. (E) Distribution of metabolic response onset times for antibiotics belonging to the seven main antibiotic categories tested in this study. The onset time is defined as the time at which metabolite changes reached half of their maximum change after treatment of M. smegmatis with the compounds. For each perturbation (treatment with compound), metabolites with a model fitting R2 ≥ 0.6 and a maximum absolute log2 fold change ≥ 2 were retained.

Raw mass spectrometry data were normalized by correcting for instrumental and systematic biases. Before identifying drug-responsive metabolites, we eliminated contaminating ions that were identified in spectra for pure compounds dissolved in water. After removal of drug-associated contaminating ions, a linear model was used to decouple the contribution of the drug treatment to the metabolic changes from plate-to-plate variance, differences in extracted biomass, and background noise. A Z-score normalization was applied before estimating the average and SD over the three biological replicates (see Materials and Methods, table S2, and figs. S1 and S2).

To systematically identify metabolites undergoing significant changes upon treatment with a drug or compound, we used a regression analysis originally described for dynamic transcriptomics data (37). For each compound, dynamic profiles of all metabolites were fitted by a model with five parameters, capable of describing different dynamic adaptive changes of metabolites to new steady-state levels. Fit quality estimated by adjusted R2 values and the magnitude of metabolite changes were used to select responsive metabolites (Fig. 1D). On average, for each individual perturbation, about 5% of the detected metabolites exhibited significant changes (adjusted R2 ≥ 0.6 from fitting analysis and absolute log2 fold change ≥ 2) (figs. S3 and S4), and the majority exhibited a significant response in at least one of the treatments (table S3). The onset time of metabolite responses reflected the distance of drug targets to metabolism. For example, very rapid responses occurred for drugs directly targeting metabolic processes, such as inhibitors of folate synthesis, whereas delays in the timing of the response onset occurred upon treatment with antibiotics targeting bacterial protein synthesis or DNA replication (Fig. 1E). Given that a broad common response to bactericidal antibiotics has been suggested (30), we investigated which metabolites exhibited nonspecific responses common to most perturbations. To this end, for each metabolite, we computed the average of fitting R2 values across all tested conditions and the average of largest absolute Z scores detected along the time course in each treatment. Overall, the commonality of metabolic changes across all tested perturbations was very low (Fig. 2A). Most metabolites affected across multiple conditions were likely associated with general stress responses and exhibited a low but significant correlation with growth inhibition (Fig. 2A and fig. S5). For example, the accumulation of glycerol 1-phosphate, the first degradation product of glycerol, exhibited a mild but significant correlation (P ≤ 1 × 10−10) with growth reduction (Fig. 2A). Overall, our data suggested that common changes independent of drug MoAs were an indirect effect of growth-related processes. Given that this conclusion was based on applying most drugs at subinhibitory concentrations, future work with these and other bactericidal drugs should examine the presently controversial issue of a common mechanism of antibiotic-induced cell death (38, 39).

Fig. 2 Commonalities among metabolite changes in response to antimicrobial treatment.

(A) Correlation between growth rate and metabolite abundance in M. smegmatis after treatment with antimicrobial compounds. Each dot represents a metabolite. The two axes represent the mean R2 across all tested conditions and the mean of maximum Z scores across all tested conditions and time points. Color reflects the degree of Spearman correlation between maximum Z score and growth inhibition across all tested conditions. Metabolites with an average R2 ≥ 0.5 and log2 Z score ≥ 0.5 are shown. UDP, uridine 5′-diphosphate; GDP, guanosine diphosphate; CDP, cytidine 5′-diphosphate. (B) Pathway enrichment for metabolome responses to antibiotics with seven known MoAs: (1) cell wall synthesis inhibitors, (2) DNA cleavage, (3) folic acid biosynthesis inhibitors, (4) quinolones, (5) mycolic acid biosynthesis inhibitors, (6) protein synthesis inhibitors, and (7) RNA synthesis inhibitors. Enrichment was performed with the 50 most frequently identified genes for each antibiotic class. The heat map shows enriched KEGG metabolic pathways with q ≤ 0.01.

Metabolome response similarity and proximity of metabolic changes to drug-target

Next, we asked whether drug-specific responses occurred in the proximity of corresponding targets and might therefore help to reveal the underlying MoAs. We used a genome-scale metabolic model of M. smegmatis (40) to calculate the distance between each enzyme-metabolite pair as the minimum number of reactions connecting the two. To assess whether metabolites exhibiting the largest change in magnitude tended to be located in the proximity of an enzyme, we used a weighted scoring function, in which all metabolic changes were weighted by the respective distance to the tested enzyme. In this analysis, the timing of the response was not taken into account as, for each metabolite, the maximum absolute fold change during the time course was selected. For each treatment and individual enzyme, we systematically assessed the significance of the proximity of metabolic changes by a permutation test (table S5). Enzymes with a significant enrichment for proximal metabolic changes could be direct drug targets or indirect mediators of the immediate cellular response. Both types of metabolic changes defined what we called the mode of metabolic interference of the compound. To identify the most overrepresented metabolic pathways across compounds with similar MoAs, we calculated the probability for each enzyme to exhibit local metabolic changes within groups of compounds with the same MoA. The resulting top 50 enzymes with consistent low P values for multiple drugs of the same class, identified across seven main antibiotic categories with known MoAs, were selected and subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (41) and correction for multiple tests (Fig. 2B and fig. S6) (42).

Identified genes in most antibiotic classes showed strong functional dependencies for the drug MoA. For instance, enzymes in folate biosynthesis were overrepresented in antifolate drug treatments, and enzymes related to nucleotide biosynthesis and DNA replication were overrepresented when M. smegmatis was treated with quinolones. In contrast, enzymes in amino acid and aminoacyl-tRNA (transfer RNA) biosynthesis or RNA degradation were enriched upon treatment with protein or RNA synthesis inhibitors, respectively (Fig. 2B). The link between the mode of metabolic interference and the MoA of inhibitors of cell wall synthesis or mycolic acid biosynthesis was less clear, potentially a consequence of the limited lipid coverage by the metabolomics method used (table S6). Nevertheless, the distinct profiles of metabolites that were often proximal to the target suggested that drugs with common MoAs featured similar metabolic signatures, regardless of whether the target was metabolic or not.

To compare the metabolome responses across multiple drugs, we calculated the similarity of metabolic changes between drug pairs by computing pairwise mutual information values for each pair of drugs at the corresponding time points. Next, the pairwise metabolome similarity between time points within 3 to 10 hours after drug exposure was summed up into one unique similarity score between drug pairs (fig. S1 and table S3). The context likelihood of relatedness algorithm (43) was applied to prune false dependencies due to promiscuity in the metabolic responses across different drugs, most likely caused by common indirect growth-related processes (44). By selecting only those perturbations with known targets, we demonstrated that compounds with similar MoAs induced similar metabolic responses (Fig. 3A). On average, compounds targeting nonmetabolic processes, such as quinolones and DNA cleavage inducers, or compounds targeting the 30S rather than the 50S ribosomal subunits (fig. S7) featured more similar metabolic signatures than did antibiotics with similar MoAs but different enzymatic targets (for example, ampicillin versus cycloserine) (Fig. 3, A and B). Overall, we demonstrated that antibiotics with common MoAs had a strong tendency to elicit similar metabolic responses (fig. S8). We also observed that certain classes of antibiotics, like quinolones, exhibited on average higher metabolome similarity than did others, such as cell wall inhibitors (Fig. 3, A and B).

Fig. 3 Pairwise similarity of antimicrobials with respect to metabolic changes induced in M. smegmatis.

(A) Similarity heat map for 62 reference antimicrobials. Similarity calculated between each drug-perturbed condition is represented as a symmetric heat map. Diagonal values are not taken into account and are set to not-a-number (gray). Highlighted boxes correspond to the seven main MoAs of the 62 reference compounds. (B) Magnification of panel in (A) showing antimicrobial compounds that blocked gyrase and cell wall synthesis in M. smegmatis. (C) Receiver operating characteristic (ROC) curve measuring the ability of metabolome-based predictions using the iterative hypergeometric test (45) to discriminate antibiotics sharing similar MoAs. Notably, we considered only MoAs that applied to more than one antimicrobial reference compound. CPR, ciprofloxacin; LVX, levofloxacin; MFL, moxifloxacin; NAL, nalidixic acid; NFL, norfloxacin; OFL, ofloxacin; AMX, amoxicillin; AMP, ampicillin; CCL, cefaclor; CTX, ceftriaxone; OCI, oxacillin; BAC, bacitracin; CYC, d-cycloserine; EMB, ethambutol; FOS, fosfomycin; AUC, area under the curve.

To avoid any bias and to predict the MoAs of uncharacterized compounds, we developed a computational framework in which the distribution of pairwise similarities between the uncharacterized compound and compounds with known MoAs was first rank–transformed, and q values indicating the statistical significance of a MoA as a top-ranking antibiotic class were calculated with an iterative hypergeometric distribution formula (45). To test the ability of metabolome-based similarity to predict the MoAs of potential new antimicrobial compounds, we used a procedure in which a reference agent, from MoAs with more than one reference compound, was treated as an uncharacterized compound and its MoA predicted de novo. Metabolome-based predictions of MoAs exhibited good agreement with known MoAs (area under the ROC curve = 0.87; Fig. 3C). Hence, we concluded that dynamic metabolome changes induced by treating M. smegmatis with antimicrobial drugs reflected functional drug properties and could be used to infer the underlying MoAs. Given that different concentrations of the same antibiotic typically elicited similar metabolite responses (fig. S9), we simplified the screening of new compounds by using a single drug dosage.

Metabolomics-based prediction of MoAs for uncharacterized compounds

Given the potential capacity to distinguish known MoAs through dynamic metabolome responses, we set out to classify MoAs for 212 uncharacterized antimycobacterial compounds with chemical structures that differed from existing antibiotics, which were identified in a large-scale phenotypic screen by GlaxoSmithKline (GSK) (4, 46). These compounds were able to kill Mycobacterium tuberculosis in vitro with minimum inhibitory concentrations (MICs) below 10 μM and with limited toxicity to human cell lines; structures of these compounds have been made available ( M. smegmatis was challenged with 10 μM of each GSK compound. Comparing the measured dynamic metabolome responses to those of the 62 reference compounds and testing the significance of overrepresented antibiotic classes among the most similar reference compounds allowed us to predict MoAs for these GSK compounds. GSK compounds that did not meet a minimum enrichment q value of 0.05 were not classified in any of the tested MoAs. Most GSK compounds exhibited a significant metabolome similarity to more than one MoA, with the majority associated with two MoAs. However, in most cases, the two predicted MoAs were functionally similar, for example, quinolones and DNA cleavage agents or inhibitors of mycolic acid synthesis or cell wall synthesis (fig. S10). Although our methodology could not fully resolve the molecular target of an uncharacterized compound solely from the metabolome profile, it did provide a robust prediction of the affected cellular process. For the GSK compounds with at least one significant (q ≤ 0.05) prediction, we focused on the most significant associated MoA (Fig. 4A and table S4). Although most GSK compounds were selected to be chemically distinct from known antibiotics (fig. S11 and table S2), more than 70% of the tested GSK compounds were classified as having already known MoAs (table S4). Notably, our data did not enable us to differentiate whether these compounds interfered directly with classical antibiotic targets, whether they bound to different proteins in the same cellular process, or whether they indirectly interfered with transcriptional regulation (47), proteome interactions (48), or enzymatic activation (for example, prodrugs) (49).

Fig. 4 Metabolome-based predictions of MoAs for 212 GSK compounds.

(A) Grouping of metabolome similarity–based predictions for the 212 GSK compounds into known MoAs. (B) Impact of antimicrobial drugs on normalized FolA in vitro activity. The measured dihydrofolate conversion rate was normalized to the activity measured with DMSO vehicle only. All tested compounds were dissolved in DMSO: 40 μM streptomycin (STR), 40 and 1333 μM trimethoprim (TRM, TRM-H), 40 μM PAS, and 40 μM of six GSK compounds. (C) RecA promoter activity in exponentially growing E. coli treated with the following: four GSK compounds predicted to be quinolone-like agents (BRL-7940SA, BRL-10988SA, GSK1066288A, and GSK695914A), norfloxacin, ampicillin, DMSO, and three GSK compounds (GSK2534991A, GSK1826825A, and GSK1518999A) predicted to be a protein synthesis inhibitor, a folic acid biosynthesis inhibitor, or with an unknown MoA, respectively. (D) Gyrase activity of M. tuberculosis measured using an in vitro supercoiling assay at different concentrations of moxifloxacin or GSK1066288A. (E) Gyrase activity of E. coli in the presence of DMSO, moxifloxacin, streptomycin, or GSK1066288A. Activity of denatured gyrase was used as a negative control (ø).

To validate our predictions for known MoAs, we focused on antifolate drugs that were among the first effective antimycobacterial agents and quinolones that exhibit a large spectrum of activities. Reference folate biosynthesis inhibitors have different targets (50). Trimethoprim, for example, directly targets the dihydrofolate reductase FolA in a noncompetitive manner, whereas others, like sulfamethizole, are competitive inhibitors of the dihydropteroate synthase FolP. Yet, others, like PAS, are prodrugs with no clear molecular target (51). Because four of the five predicted antifolate GSK compounds exhibited high metabolome response similarity to trimethoprim, we tested the inhibitory activity of the five compounds on FolA. With a coupled in vitro enzyme assay, we demonstrated direct inhibition of M. smegmatis FolA by four GSK compounds (Fig. 4B). The remaining compound, GSK275628A, elicited a metabolome response most similar to PAS (table S4). In the case of a PAS analog, enzymatic activation would be necessary to achieve FolA inhibitory activity (Fig. 4B) (51). Hence, the result from the in vitro enzyme assay was consistent with our metabolome-based predictions, suggesting a similar MoA to a prodrug for GSK275628A.

Given that quinolones are able to kill mycobacterial strains resistant to first-line antimycobacterial drugs (52), we attempted to validate the predictions for four GSK compounds with similar metabolic signatures to quinolones. We tested the induction of DNA damage by each individual GSK compound using a green fluorescent protein (GFP) reporter plasmid of RecA promoter activity (53) in exponentially growing Escherichia coli. RecA is a key regulatory protein that induces the SOS response upon DNA damage (39, 54). Consistently, with our metabolome-based predictions, we demonstrated that three of the four GSK compounds predicted to be quinolone-like agents increased RecA promoter activity (Fig. 4C). As a control, we tested ampicillin, DMSO, and three randomly picked GSK compounds (GSK2534991A, GSK1826825A, and GSK1518999A) predicted to be a protein synthesis inhibitor, folic acid biosynthesis inhibitor, and with an unknown MoA, respectively. All negative controls failed to elicit a RecA response, reinforcing the specificity of the assay and the validity of our predictions (Fig. 4C). The quinolone-predicted GSK compound that failed to elicit RecA overexpression (BRL-7940SA) also did not exhibit any inhibitory activity, potentially because it could not penetrate the cell envelope of E. coli. Next, we asked whether any of the four predicted quinolone compounds could directly inhibit the gyrase complex. We tested direct gyrase binding with an in vitro DNA supercoiling assay (fig. S12). One of the four predicted quinolone-like antibiotics, GSK1066288A, directly inhibited activity of M. tuberculosis and E. coli gyrases, comparable to the fourth-generation quinolone antibiotic moxifloxacin (Fig. 4, D and E, and fig. S13).

Given that the majority (that is, 77%) of the 212 tested GSK compounds exhibited responses similar to known MoAs (Fig. 4A), our results suggested that the GSK library contained only a relatively modest number of compounds that could potentially affect nonconventional cellular processes. Among the about 25% of GSK compounds without a metabolome response similarity to known MoAs, 33 did not show significant response differences (q ≥ 0.05) to the control conditions. These compounds were not active under the tested conditions, possibly hinting at a low penetration, or the metabolic response was beyond our detection or coverage limits. The MoAs of these compounds need to be further characterized directly in M. tuberculosis (fig. S14). For the remaining 16 compounds with a metabolic signature different from any of the reference compounds, we identified enzymes in close proximity to the detected metabolic changes and searched for enriched metabolic pathways (fig. S15 and table S2). These metabolome-based analyses suggested unique abilities of these compounds to interfere with nitrogen metabolism, oxidative phosphorylation, and the metabolism of fatty acids, vitamin B12, or glyoxylate (fig. S15). The highly similar metabolome responses for 6 of 16 compounds suggested a common MoA (Fig. 5A).

Fig. 5 Analysis of metabolite changes in M. smegmatis after treatment with the GSK compound GSK2623870A.

(A) Pairwise similarity between M. smegmatis metabolome response profiles to 16 GSK compounds with no similarity to known MoAs. (B) Schematic representation of trehalose monomycolate exporter protein MmpL3. Transmembrane segments are represented in violet. Circles indicate the locations of amino acid changes associated with resistance of M. tuberculosis to antimycobacterial lead compounds previously found to select for resistance mutations in MmpL3. These compounds include SQ109 (blue), THPP (orange), and SPIROS (purple). The black star indicates the amino acid change associated with resistance to the GSK compound GSK2623870A. The genomes of the M. tuberculosis H37Rv and Beijing GC1237 mutant strains contain an A to G SNP at position 755, which resulted in a tyrosine to cysteine missense mutation at position 252 of the MmpL3 protein. (C) Results from limited proteolysis analysis. Each dot in the volcano plot represents the relative difference in peptide abundance between the treated and untreated proteome extracts. Proteins highlighted in red are known to physically interact with fatty acid synthase FAS-I (figs. S17 and S18). For each protein, the size of the dot reflects the number of interacting partners (80) with significant conformational changes. (D) Rapid metabolic changes induced by the GSK antimicrobial compound GSK2623870A. Each dot corresponds to the R2 and Z-score values of the metabolite 5 min after exposure of M. smegmatis to the antimicrobial compound. Metabolites highlighted in red are involved in fatty acid metabolism.

To investigate the potentially new MoAs for these six GSK compounds, we analyzed their similar responses more deeply. The most common and prominent metabolic changes occurred in α,α-trehalose-6-phosphate and neighboring metabolites in lipid metabolism, carbon metabolism, and the cobalamin/heme biosynthesis pathway (fig. S16, blue dots). GSK1829729A and GSK1829728A are chemically similar to tetrahydropyrazolo[1,5-a]pyrimidine-3-carboxamide (THPP). Recent pull-down experiments with Mycobacterium bovis demonstrated that THPP binds to EchA6, an enzyme homolog of enoyl–coenzyme A (CoA) hydratase that is involved in fatty acid synthesis (55). Consistently, our network-based analysis of metabolic changes induced by GSK1829729A and GSK1829728A identified local changes in the proximity of several enoyl-CoA hydratase reactions (table S5). To test whether the six GSK compounds acted through a MoA similar to that of THPP, we selected M. tuberculosis H37Rv and Beijing GC1237 mutants that were spontaneously resistant to GSK2623870A.

Spontaneous drug-resistant M. tuberculosis mutants have been extensively used to resolve the MoAs of bioactive chemicals, because the mutations can occur directly in the drug target (7) or may impact biological processes functionally related to the drug MoA (17). Mutants occurred at a frequency of 5 × 10−7 on solid media containing GSK2623870A (29.76 μg/ml) (16× MIC). Despite the difference in chemical structures between THPP and GSK2623870A, whole-genome sequencing (WGS) revealed a single-nucleotide polymorphism (SNP) in the mmpL3 gene (Rv0206c) of both H37Rv and Beijing GC1237 M. tuberculosis mutants (Fig. 5B). Notably, similar resistance mechanisms were observed for THPP (56), N-benzyl-69,79-dihydrospiro[piperidine-4,49-thieno[3,2-c]pyran (SPIROS) (56), and N-(2-adamantyl)-N-[(2E)-3,7-dimethylocta-2,6-dienyl]ethane-1,2-diamine (SQ109) (57). MmpL3 is an essential inner membrane transporter in M. tuberculosis, responsible for translocating mycolic acids in the form of trehalose monomycolates across the cell membrane (58, 59). It is also involved in heme uptake (60), potentially explaining the metabolome changes in heme biosynthesis seen in M. smegmatis in response to all six GSK compounds (fig. S16).

MmpL3 has recently been proposed as a promising target for developing drugs against M. tuberculosis (61). MmpL3 has been found to be the target of various antimicrobial drugs with different chemical scaffolds (56, 59, 62), and different mutations in mmpL3 have been associated with drug resistance (56, 57, 59, 63). Moreover, MmpL3 inhibitors have shown synergistic interactions with first- and second-line antimycobacterial drugs (57, 64). Despite the apparent attractiveness of MmpL3 as a drug target (65), the vulnerability of this protein to the action of compounds with different scaffolds has generated controversy regarding the importance of the differences in compound chemical structures and their MoAs (55, 65).

To better characterize the MoA of GSK2623870A, we used a chemical proteomics strategy that combined limited proteolysis with mass spectrometry to probe protein conformational changes on a proteome scale (6668). Protein extracts of M. tuberculosis harvested during exponential growth were incubated with the GSK2623870A compound for 5 min before conducting limited proteolysis. The analysis identified 30,532 peptides mapping to 2564 unique proteins in M. tuberculosis (table S8). For each peptide, we quantified the difference in abundance between the treated and untreated proteome samples. Differences in peptide abundance reflected protein conformational changes induced by the action of the compound. Such conformational changes can reflect direct binding events or protein network rearrangements by disruption of protein-protein interactions. For 69 peptides mapping to 60 unique proteins, we detected significant changes between the treated and untreated conditions (absolute log2 fold change ≥ 2 and P ≤ 0.01) (Fig. 5C). Whereas no significant changes were associated with MmpL3, we found several conformational changes in proteins involved in fatty acid metabolism. In particular, the fatty acid synthase FAS-I protein (Rv2524c) exhibited multiple conformational changes together with many of its interacting proteins, such as FadA, FadA3, and FabD (Fig. 5C and figs. S17 and S18). This finding was consistent with the most rapid metabolic changes (absolute log2 fold change ≥ 1 and R2 ≥ 0.6; 5 min after treatment) induced by GSK2623870A. These metabolic changes occurred for several metabolites in close proximity to FAS-I, exemplified by the depletion of the enzyme substrate malonyl-CoA (Fig. 5D). Proteins with a conformational change included the known penicillin target PbpA, as well as proteins involved in the tricarboxylic acid cycle. Collective evidence suggests a complex MoA for GSK2623870A due to potentially multiple targets or an indirect role for MmpL3 as a drug importer (55).


We describe here a general approach to systematically quantify and interpret dynamic metabolome responses at high throughput to predict the MoAs of small molecules, independent of whether their targets were metabolic or not. The large compendium of drug-metabolome profiles that we generated could be used as a reference data set to assess similarity in the responses of mycobacteria to new antimicrobial compounds. We predicted and experimentally validated compounds in the GSK library and identified new inhibitors of DNA replication and folate biosynthesis. Network-based analysis of dynamic drug-metabolome changes resolved distinct metabolic responses that were specific to drug MoAs and were used here to identify compounds with potential new MoAs against mycobacterial targets. Metabolic changes could be detected in the absence of growth inhibitory conditions (fig. S19), representing a key advantage over classical bacterial growth assays. This feature was crucial to enable identification of a group of six antimycobacterial compounds with similar yet distinguishable metabolic signatures. WGS revealed that MmpL3 mediated resistance to one of the GSK compounds in M. tuberculosis, and limited proteolysis analysis suggested that fatty acid metabolism was the main cellular process in M. tuberculosis targeted by this compound.

Overall, our results suggested that the compound GSK2623870A mainly affected fatty acid metabolism, similar to THPP, expanding the diversity of chemical compounds that potentially could be used to target metabolic processes upstream of mycolic acid biosynthesis. Although the exact functions and roles of MmpL3 in mediating the MoA of this GSK compound still remain to be elucidated, the mutant and proteome results confirmed the initial metabolome-based prediction of an unconventional MoA for this group of compounds. Whereas the 16 antimycobacterial compounds with predicted unconventional MoAs exhibited very modest inhibitory activity in M. smegmatis, they elicited subtle but unique metabolic patterns reflecting the underlying specificity of their MoAs. Although these compounds are promising candidates for discovering new targets, future studies should investigate the direct metabolic impact of these compounds on M. tuberculosis to better characterize their MoAs. Nevertheless, the relevance of these compounds tested in M. smegmatis for M. tuberculosis illustrates that our approach was able to classify compound MoAs in a manner largely independent of the organism. In contrast to traditional chemogenomic methods that are often restricted to model organisms, our metabolome-based approach did not require libraries of mutant bacteria and so can be applied to any bacterial species. For reasons of technical simplicity and cost, it might be preferred to perform the first large-scale profiling of compounds in a related, nonpathogenic bacterium given that many MoAs will be similar and drug efficacy can be driven by factors such as drug penetration.

Current state-of-the-art methods to monitor and interpret dynamic drug-induced cellular responses at the genome level do not scale with the size of typical compound libraries used in drug discovery (2830, 32, 38, 69). Although traditional omics methods such as proteomics, transcriptomics, or genome sequencing (70) provide a large information output, our nontargeted metabolomics workflow enabled rapid measurements at a 10- to 100-fold higher throughput. No single technology provides a general solution for identification of a compound’s MoA or specific molecular target. Our approach, however, could play a crucial role in the initial characterization of new compounds and could provide complementary data for the interpretation of cytological and phenotypic assays (12, 13). One major limitation of our approach is the difficulty in decoupling metabolic changes that are direct consequences of drug-target interactions from indirect adaptive metabolic adaptations. At the cost of a lower throughput, time-resolved metabolic profiles with seconds resolution could better resolve direct drug effects and facilitate interpretation of drug-induced metabolome changes. However, it is worth noting that if drug targets are inactive or have low expression under the conditions tested, then drug effects will be more subtle or impossible to detect. Mimicking in vivo conditions in an in vitro setting can improve the efficacy of phenotypic-based assay screens and maximize the output of a metabolomics approach.

Our metabolome-based screening approach can be applied directly to extract multiple quantitative signatures indicative of functional properties of MoAs in large compound libraries. A major advantage of our systematic and rapid methodology for predicting MoAs in large chemical screens is the discrimination of compounds with previously known MoAs from those that have new MoAs. Knowing the MoAs of compounds at an early stage of drug discovery can guide the selection of the most promising leads, even in cases where the drugs do not yet exhibit strong bactericidal or bacteriostatic effects. For drug development, MoA identification is important for follow-up studies, biomarker development, and drug target identification. Even after a relevant drug target is established, additional functional information from our metabolomics drug profiling strategy could help to identify unwanted off-target effects or establish new roles for the target protein within its biological network. Hence, combining our metabolomics and computational-based methodology with other omics strategies and more targeted validation technologies should help to overcome the current challenge of predicting the biological effects of drugs before clinical trials begin.

Although the set of reference compounds used in this study was restricted to antibiotics with the most common MoAs, future studies may benefit from a larger set of characterized compounds, regardless of their antimicrobial activities. We envisage that by increasing the number of tested organisms and reference compounds, alternative computational approaches, such as machine learning (71), could open up new opportunities in antibiotic discovery. Beyond inference of the MoAs of bioactive compounds, our approach could be used to identify direct and indirect enzymatic targets of nonlethal compounds, paving the way to predictive models and rational design of effective combinatorial treatments (71, 72).


Study design

This study was designed to evaluate the ability of a mass spectrometry approach to predict the MoA of small molecules by monitoring dynamic metabolic responses of mycobacteria after drug treatment. To this end, we established a new combined experimental-computational framework for large-scale profiling of bacterial metabolic responses to treatment with bioactive compounds. This approach was first tested in M. smegmatis on a reference set of 62 compounds applied at different dosages. Rigorous analysis of metabolomics data demonstrated that antibiotics with a similar MoA elicited similar metabolic changes. This strategy was used to characterize the MoAs of a set of 212 new antimycobacterial compounds, the result of an antibiotic discovery effort of the pharmaceutical company GSK. This library was used under the terms and conditions of a material transfer agreement (MTA) between GSK and ETH Institution.

Strains and media

M. smegmatis MC2-155 was used throughout this study. Standard 7H9 minimal medium was used to perform the metabolome drug screening experiment. One liter of medium contained 2.5 g of Na2HPO4, 1 g of KH2PO4, 0.5 g of (NH4)2SO4, 0.5 g of l-glutamic acid, 0.1 g of NaCl, 0.1 g of MgSO4·7H2O, 40 mg of ferric ammonium citrate, 1.8 mg of ZnSO4·7H2O, 1 mg of CuSO4, 1 mg of pyridoxine, 0.5 mg of CaCl2, 0.5 mg of biotin, 2 g of glycerol as carbon source, and 0.05% (v/v) tyloxapol. M. tuberculosis: H37Rv (ATCC 27294), Beijing GC1237 (73), M. smegmatis MC2-155 (74), Mycobacterium marinum (Aranson 1926), Mycobacterium aurum (75), and their derivatives were cultured at 37° or 32°C, in the case of M. marinum, in Middlebrook 7H9 medium supplemented with ADC [0.5% (w/v) bovine serum albumin, 0.2% (w/v) dextrose, 0.085% (w/v) NaCl, 0.0003% (w/v) beef catalase; Difco], 0.5% (w/v) glycerol, and 0.05% (w/v) Tween 80 or in solid Middlebrook 7H11 medium supplemented with OADC [0.05% (w/v) oleic acid, 0.5% (w/v) bovine serum albumin, 0.2% (w/v) dextrose, 0.085% (w/v) NaCl, 0.0003% (w/v) beef catalase; Difco].

Metabolomics drug screening

M. smegmatis was precultivated at 37°C and 300 rpm for 24 hours, diluted 1:80 in 125 ml of fresh medium in 1-liter shake flasks, and grown for ~15 hours. Cell cultures (700 μl) were distributed in 96-well deep plates at an OD of about 0.4 and incubated for 1 hour at 37°C and 300 rpm before exposure to the drug treatment by the addition of 10 μl of drug solution. Sixty-two reference compounds were chosen from well-characterized antibiotics and stress factors, covering most of the currently available antibacterial MoAs (see Fig. 1C and table S1 for more details), together with 212 antituberculosis compounds, obtained from GSK. All stocks were prepared in DMSO solution, except the following ones, which were prepared in H2O due to the low solubility of the compounds in DMSO: CPR, GEN, GET, NEO, H2O2, Parq, NaOH, CuCl, EDTA, CTX, OCI, and PMB. The GSK compounds were used at a single concentration of 10 μM, whereas the reference compounds were administered in at least three different dosages, aiming for ~100, ~50, and ~20% inhibition of growth rate (table S2). All treatments were measured in three biological replicates. DMSO and H2O treatments were included on every plate as control, resulting in 433 different treatments.

Samples were taken immediately before treatment and after 5, 30, 60, 180, 360, and 600 min from drug exposure. At each time point, 40 μl of whole-cell broth was transferred to 120-μl medium on microtiter plates, and OD was measured at 595 nm. In parallel, 40 μl of the cultures was directly transferred to 120 μl of extraction liquid solution containing 50% (v/v) methanol and 50% (v/v) acetonitrile at −20°C. The extraction was carried out by incubating the samples for 1 hour at −20°C. Samples were centrifuged for 5 min at 4000 rpm, and 80 μl of the supernatant was transferred to 96-well storage plates and stored at −80°C.


Growth rates of M. smegmatis upon drug treatment were calculated as the slopes of linear fits to log-transformed OD curves experimentally determined from 30 to 360 min after drug exposure (fig. S5). In few cases, apparent negative growth rates caused by cell death were set to zero. The inhibition ratios reported in table S4 equal one minus the ratio between the growth rate of the treatment and the average growth rate of the controls on the same plate.

Mass spectrometry measurement, spectral data processing, and annotation

The analysis was performed on a platform consisting of an Agilent Series 1100 LC pump coupled to a Gerstel MPS2 autosampler and an Agilent 6550 Series Quadrupole Time of Flight mass spectrometer as described previously (36). Mass spectra were recorded from mass/charge ratio (m/z) 50 to 1000 using the highest resolving power (4 GHz HiRes).

All steps of mass spectrometry data processing and analysis were performed with Matlab (MathWorks). For each acquired spectra, ion peaks were identified with the findpeaks function of the Signal toolbox. Identified peaks with less than 5000 ion counts were excluded from further analysis. Detected ions were matched to a list of metabolites based on the corresponding molar mass. Because there is no available manually curated genome-scale model of M. smegmatis metabolism, we compiled a comprehensive list of metabolites from the following models: M. tuberculosis H37Rv model from the Model SEED database [Seed83332.1 (76);] and the manually curated M. tuberculosis model from (77), E. coli K12 model iJO1366 from (78), and M. smegmatis MC2-155 model from the Biomodels database [BMID000000141548 (79);]. For the full list of metabolites used for annotation, see table S6. The chemical formula of each metabolite was used to calculate the deprotonated monoisotopic molecular mass. Detected ions within a mass tolerance difference of less than 0.003 Da were associated to the nearest reference metabolites. In total, 1006 different reference masses could be detected (see tables S2 to S6 for the reference masses and the associated metabolites).

Supercoiling assay

To test gyrase activity, the M. tuberculosis and E. coli DNA gyrase assay kit from Inspiralis and TopoGEN, respectively, were used. The procedure adopted to monitor DNA gyrase–dependent supercoiling followed kit instructions. Reaction was performed in 30-μl volume by incubating about 250 ng of pHOT1 Relaxed DNA together with the purified gyrase protein at 37°C for 45 min. Moxifloxacin or the compound GSK1066288A was added to the reaction mixtures at different concentrations. Termination, detection, and readout were carried out as described in the kit protocol using agarose gel electrophoresis. Bands were quantified using ImageJ. Raw gel pictures can be found in the Supplementary Materials.

Folate assay

The capability of compounds to inhibit the folate biosynthesis was tested with an enzyme assay. Full protein extract was prepared from M. smegmatis the following way. Cell culture (70 ml) in exponential phase (at about OD of 1) was harvested and centrifuged for 10 min at 4°C and 5000 rpm. The precipitate was resuspended in 8-ml lysis buffer [100 mM tris buffer (pH, 7.5), 4 mM phenylmethylsulfonyl fluoride, 2 mM dithiothreitol, 5 mM MgCl2]. French press was used to extract the cells. The samples were centrifuged for 10 min at 4°C and 5000 rpm, and the supernatant was filtered to concentrate the protein extract to about 2 ml (centrifuged at 4°C and 5000 rpm in an ultrafiltration tube). The reactions were performed in 150 μl on a microtiter plate. The ratio of the buffer/protein extract was 1:2 in the solution. The reaction was started by adding 0.1 and 0.4 mM of dihydrofolate and NADPH (reduced form of nicotinamide adenine dinucleotide phosphate), respectively. A plate reader was used to measure the decrease in the NADPH level at 340 nm for 1 hour, sampling the absorbance every 26 s.

Selection of spontaneous GSK2623870A-resistant mutants of M. tuberculosis H37Rv and Beijing GC1237 and WGS

GSK2623870A spontaneous resistant mutants were selected on 7H11 plates supplemented with OADC and 16× MIC of antibiotic at 37°C after 3 weeks of incubation. Genomic DNA was isolated from saturated cultures of the different mycobacteria and resistant isolates. Sequencing was carried out by GATC, using Illumina MySeq 300 and HiSeq 4000 instruments. Identification of SNPs, insertions, and deletions was done using different bioinformatics tools as GATK modules, PICARD, and SnpEff. Raw sequence reads were deposited in the National Center for Biotechnology Information sequence read archive database under accession number SRP116576.

Statistical analyses

Statistical analyses were performed using Matlab R2015b (MathWorks). In total, we monitored dynamic metabolic changes across 431 compound treatment conditions and two mock treatments (vehicle control H2O or DMSO) over seven time points in three biological replicates. To assess differences in metabolite abundance before and after treatment, we used a multiple linear regression scheme. Similarity of metabolic profiles between different treatments was assessed using an entropy-based measure of similarity (Mutual Information), and prediction of compound MoAs was based on an iterative hypergeometric test described in (45). When necessary, P values were corrected for multiple tests by q value estimation. Metabolomics data and analysis results are available as supplementary tables.


Materials and Methods

Fig. S1. Schematic description of MS data normalization and analysis.

Fig. S2. Analysis of correlations across biological replicates.

Fig. S3. Distribution of responsive metabolites.

Fig. S4. Number of affected metabolites per MoA.

Fig. S5. Distribution of growth inhibitory activities.

Fig. S6. Pathway enrichment for metabolome responses to antibiotics with seven major MoA.

Fig. S7. Pairwise drug similarity.

Fig. S8. Metabolome-based similarity.

Fig. S9. Similarity between compounds with equivalent MoAs as a function of the difference in growth inhibition.

Fig. S10. Distribution of predicted MoAs.

Fig. S11. Pairwise compound chemical distance.

Fig. S12. M. tuberculosis gyrase assay.

Fig. S13. M. tuberculosis and E. coli gyrase assay.

Fig. S14. Distribution of growth inhibitory activities for GSK compounds with classified and unclassified MoAs.

Fig. S15. Pathway enrichment analysis for compounds with potential unconventional MoAs.

Fig. S16. Common metabolic responses across GSK compounds with potential new MoAs.

Fig. S17. Protein-protein interactions among proteins with significant conformational changes detected by limited proteolysis analysis.

Fig. S18. Robustness of results from limited proteolysis.

Fig. S19. Similarity between compounds with equivalent MoAs as a function of growth inhibition.

Fig. S20. Data normalization.

Fig. S21. Schematic representation of the procedure used to estimate pairwise similarity among tested compounds.

Table S1. Antibiotic perturbation list.

Table S2. Metabolome data (perturbation name versus metabolites data matrix).

Table S3. Impulse model fitting results (maximum fold change and R2 matrices).

Table S4. MoA predictions (list of top predicted MoAs and complete matrix of enrichment q values).

Table S5. Gene-drug assignments with P values coming from the network locality analysis.

Table S6. List of metabolites used to annotate peaks in the mass spectra.

Table S7. Comparison of similarity metrics.

Table S8. Results from analysis of limited proteolysis data.

References (8186)


Acknowledgments: We thank M. Zimmerman and D. Sevin for technical help and GSK for providing the compound library. Funding: This work was supported by a Sciex Fellowship to B.S. and an ETH Zürich Postdoctoral Fellowship to M.Z. B.P. was supported by the Wellcome Trust, the Lendület Programme of the Hungarian Academy of Sciences, and grant GINOP-2.3.2-15-2016-00014 (EVOMER). M.V.B. and B.G. acknowledge funding from the European Seventh Framework Program (NAREB project no. 604237). Author contributions: M.Z. designed the project and the experimental profiling of metabolome responses. B.S. performed the metabolomics experiments, supercoiling assay, and in vitro enzyme assays. M.Z. performed the validation experiments in E. coli and designed the analysis. M.Z. and B.S. analyzed the data. B.P. advised on bioinformatics analysis. M.V.B. and B.G. performed and analyzed all experimental data in M. tuberculosis. S.B. and I.P. performed proteome extraction and limited proteolysis analysis. M.V.B., B.G., A.T., P.P., J.L., and S.G. advised on data interpretation. M.Z. and U.S. wrote the manuscript. All authors contributed to the preparation of the manuscript. Competing interests: P.P. is an advisor for Biognosys AG and is an inventor on a patent licensed by Biognosys AG that covers the limited proteolysis–mass spectrometry method used in this study: “Method and tools for the determination of conformation and conformational changes of proteins and of derivatives thereof” (patent no. EP12008011.4). J.L. is an employee of GSK. The other authors declare that they have no competing interests. Data and materials availability: The compounds used in this study are available from GSK under an MTA. The compendium of drug-metabolome profiles is available in the supplementary tables. Metabolome data can be downloaded from with accession number S-BSST113.
View Abstract

Stay Connected to Science Translational Medicine

Navigate This Article