Genome-wide mutant profiling predicts the mechanism of a Lipid II binding antibiotic.

Identifying targets of antibacterial compounds remains a challenging step in the development of antibiotics. We have developed a two-pronged functional genomics approach to predict mechanism of action that uses mutant fitness data from antibiotic-treated transposon libraries containing both upregulation and inactivation mutants. We treated a Staphylococcus aureus transposon library containing 690,000 unique insertions with 32 antibiotics. Upregulation signatures identified from directional biases in insertions revealed known molecular targets and resistance mechanisms for the majority of these. Because single-gene upregulation does not always confer resistance, we used a complementary machine-learning approach to predict the mechanism from inactivation mutant fitness profiles. This approach suggested the cell wall precursor Lipid II as the molecular target of the lysocins, a mechanism we have confirmed. We conclude that docking to membrane-anchored Lipid II precedes the selective bacteriolysis that distinguishes these lytic natural products, showing the utility of our approach for nominating the antibiotic mechanism of action.

The need for new antibiotics to treat hospital- and community-acquired bacterial infections has been widely publicized[1]. Nevertheless, antibacterial development has struggled to keep pace with emerging resistance. Multi-drug resistance in Gram-negative and Gram-positive pathogens has severely limited the effectiveness of major antibiotic classes, including fluoroquinolones, β-lactams, and glycopeptides[2,3]. The dwindling number of efficacious drugs to treat bacterial infections necessitates the development of better approaches to produce the next generation of antibacterials.

Target identification is a major bottleneck to advancing antibacterials through clinical development. It is crucial to identify the molecular target of a compound to rule out non-specific mechanisms of action and guide structure-activity studies. Whole genome sequencing can provide the molecular target if resistant mutants can be raised to a compound. When this approach is not feasible due to compound limitations or fails due to a compound’s mechanism, other approaches must be used. Numerous strategies to characterize the mechanism of action of antibacterial compounds have been developed. These include biochemical approaches that compare how a compound affects incorporation of radiolabeled precursors into macromolecules (MMS assays, for macromolecular synthesis)[4], imaging approaches that examine how compound treatment affects cytological profiles (BCP, for bacterial cytological profiling)[5], and functional genomics strategies that systematically evaluate compound activity against arrayed over- and underexpression mutant libraries[6,7]. Functional genomics strategies can nominate individual molecular targets and resistance mechanisms, whereas the other biochemical approaches typically provide information about pathways only; however, arrayed libraries are time-consuming to make, expensive to maintain, and laborious to interrogate with new compounds as each library member is assayed individually.

We thought it might be possible to predict antibiotic mechanism of action using mutant fitness data from pooled Staphylococcus aureus transposon libraries. Next-generation transposon sequencing methods such as Tn-seq can map the locations of all transposon insertions in a pooled mutant library, and it is possible to assess the fitness of each gene knockout under a given condition by comparing sequence reads for that gene in treated and untreated samples[8-11]. Transposon libraries for Tn-seq analysis are typically prepared using a single transposon cassette that generates only inactivation mutants. However, we have developed a S. aureus transposon mutagenesis platform that includes a suite of bar-coded transposon cassettes with outward facing promoters[11]. Depending on the orientation of insertion, a transposon with an outward facing promoter that inserts proximal to a gene may upregulate it. Target upregulation is known to shift the minimum inhibitory concentrations (MIC) of many antibiotics and has been exploited previously to identify targets of antibacterial compounds by either: 1) testing upregulation mutants for a shift in MIC in an arrayed library format, or 2) selecting upregulation mutants from a pooled library by plating on antibiotic, a strategy that achieves spatial separation of transposants[12]. While the latter strategy is efficient, it is very compound-intensive. We thought that if upregulation signatures could be clearly discerned in Tn-seq data, then direct analysis of antibiotic-treated library cultures would have substantial advantages over other methods in terms of efficiency and compound usage; moreover, Tn-seq data provides information concerning mutations that decrease fitness as well, which collectively could provide additional insights into mechanism and intrinsic resistance factors[9,11,13]. Because antibiotics having similar mechanisms can be clustered based on their inactivation mutant fitness profiles[14], we thought it might be possible to use these profiles to predict the mechanism of action for unknown compounds in cases where upregulation signatures are insufficient.

We previously reported a S. aureus transposon library containing 690,000 unique transposon mutants, which was made by combining six sublibraries, each prepared using a different transposon cassette[11] (Supplementary Table 1). The transposon cassettes contained outward-facing promoters of different strengths to induce a gradient of gene expression levels when inserted in the correct orientation upstream of a gene[12]. To develop computational methods that use Tn-seq data to predict mechanism of action, we have treated the transposon library with 32 antibiotics with known mechanisms. We show that the known targets for many of the compounds, as well as other known and new resistance mechanisms, can be identified from directional biases in transposon insertions. We also show that inactivation mutant fitness data can be used to predict mechanism of action. In a striking validation of the latter approach, we have reassigned the mechanism of action of a family of bacteriolytic antibiotics called the lysocins.

RESULTS

Mutant fitness profiles were obtained for 32 antibiotics

We grew the high-density S. aureus transposon library in the presence of 32 different antibiotics (Fig. 1 and Supplementary Table 2) to obtain Tn-seq data (Fig. 2a). All major targets for broad spectrum antibiotics were represented. The cell wall biosynthetic pathway was particularly well represented and included inhibitors that prevent Lipid II biosynthesis by inhibiting intracellular enzymes in the peptidoglycan biosynthetic pathway (fosfomycin and D-cycloserine)[15], by blocking lipid carrier recycling (bacitracin and amphomycin)[16], and by sequestering the lipid carrier in other intermediates (wall teichoic acid export inhibitors)[17,18]. We also tested cell wall inhibitors that allow Lipid II synthesis, but block its assembly into crosslinked peptidoglycan. CDFI and DMPI act by preventing Lipid II export to the cell surface[19]; moenomycin, cefaclor, and oxacillin inhibit late-stage enzymes that polymerize or crosslink peptidoglycan[15]; vancomycin, ramoplanin, and lysobactin also inhibit late steps in peptidoglycan assembly, but do so by binding to Lipid II[15]. For other pathways, we included antibiotics representing several distinct mechanisms (Supplementary Table 2). For example, we included four mechanistically different protein translational inhibitors[20]: chloramphenicol and linezolid, which bind to the peptidyl transferase site in the 50S subunit; tetracycline, which binds to the 30S subunit and prevents aminoacyl-t-RNA binding to the A site; and mupirocin, which prevents protein translation by inhibiting isoleucyl t-RNA synthetase. The three fatty acid biosynthesis (Fab) inhibitors we tested (triclosan, cerulenin, and platensimycin) also inhibit different steps in the Fab pathway[21]. Treatment concentrations were chosen empirically by identifying at least one concentration that resulted in a moderate delay (three to five hours) in growth to stationary phase compared to the untreated control, as well as one concentration that resulted in a substantial delay (~20 hours). DNA was prepared for sequencing as previously described[11]. In general, we found that lower concentrations of antibiotic were useful for revealing compound-susceptible mutants (negative selection) whereas higher concentrations selected for compound-resistant mutants (positive selection).

Figure 1

Schematic showing known sites of action of antibiotics used to treat the S. aureus transposon library

In addition to the antibiotics shown, the library was also treated with daptomycin, CCCP, polymyxin, and gramicidin A. Pathways/processes affected and enzymes inhibited by antibiotics used in this study are summarized in Supplementary Table 2.

Figure 2

Genes upregulated by transposon insertion reveal mechanisms of resistance for the tested antibiotics

(a) Schematic illustrating that gene upregulation, which can increase fitness in the presence of an antibiotic, can be detected by an increase in the number of reads due to insertions of transposons with outward-facing promoters in a single orientation ahead of a gene. (b) and (c) Examples of data showing insertions proximal to known targets or resistance determinants in untreated versus treated samples for DMPI and daptomycin. The data shown is for a single transposon cassette. The same findings were found using at least one other transposon cassette in the same pooled library experiment.

Insertion bias identifies upregulated genes

We reasoned that if gene upregulation conferred a fitness advantage in the presence of an antibiotic, we should observe both a strong bias in the orientation of insertions and an increase in the sequence reads ahead of the upregulated gene compared with the untreated control. We developed an automated method to identify upregulation signatures based on these expectations (Supplementary Fig. 1). The “upregulation signatures” potentially reveal molecular targets as well as other resistance mechanisms (Fig. 2 and Supplementary Fig. 2). While we recognize that a transposon promoter may upregulate a gene several kilobase pairs away from the insertion, especially if the genes are in the same operon, we assumed for this analysis that the genes proximal to strongly biased insertions were most likely to be upregulated. For non-essential candidate genes, if gene upregulation increases resistance to a compound, then inactivating the same gene should increase susceptibility to the compound, and we would expect very few reads in the treated sample compared to the control.

Gene upregulation signatures can reveal known targets

Several classes of upregulated genes previously implicated in antibiotic resistance were identified from upregulation signatures, including genes for known targets. For example, for fosfomycin we identified both alleles of murA[22]; for DMPI and CDFI we identified murJ, which encodes the Lipid II flippase[19,23] (Fig. 2b); and for triclosan we identified fabI, which encodes enoyl-ACP reductase [24] (Supplementary Table 3). We also identified genes for other known resistance mechanisms such as efflux pumps (norA for ciprofloxacin[25]; lmrB2/3, emrA, and yhgE for platensimycin; lmrB2 for sorangicin), antibiotic modification genes (fosB, which encodes an enzyme that inactivates fosfomycin[26]), target modification genes (uppP, which encodes an enzyme that dephosphorylates undecaprenyl pyrophosphate, the target of bacitracin[27]), and genes encoding stress response systems (vraRS for cefaclor[28]). For daptomycin, a calcium-dependent lipopeptide antibiotic that disrupts membrane integrity[29,30], we identified mprF (Fig. 2c), which encodes a polytopic membrane protein that attaches lysine to phosphatidylglycerol to modify cell membrane charge[31]. Consistent with this, other groups have identified genetic changes in both clinical and laboratory daptomycin-resistant isolates that increase MprF activity[29,32]. In the daptomycin-treated sample, we also identified other strong upregulation signatures for genes putatively encoding small proteins of unknown function (SAOUHSC_00969, SAOUHSC_02149, and SAOUHSC_02164).

We confirmed that upregulation signatures identify mutants for which increased fitness is due to increased expression of downstream genes in two ways. First, we raised resistant mutants to four of the antibiotics on plates and identified colonies with transposon insertions in the promoter region of known target genes. Sequencing identified the specific location of the insertion as well as the promoter found in the transposon construct, and semi-quantitative RT-PCR confirmed increased expression of the target genes compared to the wildtype control (Supplementary Fig. 3). For the genes of unknown function identified as conferring a survival benefit in the presence of daptomycin, we showed that overpression from a plasmid increased resistance to daptomycin, but not to several other antibiotics tested, consistent with the Tn-seq data that showed an upregulation signature only for daptomycin (Supplementary Table 3; Supplementary Fig. 4 and 5). Inactivation mutants of two of the genes were also tested and found to have increased susceptibility to daptomycin. We conclude that upregulation signatures in Tn-seq data can identify known targets as well as other resistance mechanisms for many antibiotics even when the mutant pools are large and complex.

Genome-wide fitness profiles can predict mechanism

We sought alternative ways to use the Tn-seq data to predict antibiotic mechanism of action because upregulation of a gene cannot always identify a target. For example, some antibiotics do not have protein targets. Others may have protein targets, but if these are part of a larger complex, and if the target is not the limiting component, then upregulation of its gene may not confer a fitness advantage. For other targets, overexpression may be toxic and upregulation may only confer a fitness advantage under a narrow range of antibiotic concentrations, which may not have been tested. We wondered whether it would be possible to predict mechanism of action of an antibiotic by comparing its mutant fitness profile to the profiles of a panel of antibiotics with known mechanisms (Fig. 3a). To test this possibility, we assigned a normalized fitness value to each gene based on the change in number of reads mapping to it upon antibiotic treatment compared to an untreated control. Reads for all six transposons were combined to obtain these fitness values, which ranged from 0 (depleted compared to the control) to 1 (enriched). Once fitness values were assigned, we explored a supervised machine learning approach that would allow us to categorize the antibiotics used in our training set into different classes based on their known mechanisms (Supplementary Table 4). We found that including mainly the data from library samples grown in antibiotic concentrations that did not severely inhibit growth gave the most consistent results. We compiled a list of the genes with the highest and lowest fitness scores from each Tn-seq file, resulting in a set of 342 genes that were used with the K-nearest neighbors’ algorithm (k=3) to predict mechanism. To validate the method, we used an exhaustive leave-one-out cross-validation approach in which all samples but one were used as the training dataset while the remaining sample was used as the test dataset. The algorithm predicted the correct pathway/mechanism 72% of the time, a respectable outcome (Fig. 3b and Supplementary Table 5); however, there were distinct differences in prediction accuracies depending on the pathway. For some classes of inhibitors, e.g., fatty acid synthesis inhibitors and RNA polymerase inhibitors, all predictions were correct at all concentrations, while for others, e.g., the DNA replication inhibitors and folate pathway inhibitors, the predictions were mostly incorrect. This likely reflects the multitude of cellular processes affected during folate depletion and replication inhibition.

Figure 3

Supervised machine learning predicts mechanism of action of antibiotics

(a) Inactivation mutant fitness values for an antibiotic of unknown mechanism are compared to fitness values obtained for known antibiotics to predict mechanism of action. Boxes represent genes and shading represents fitness of the corresponding inactivation mutants. Each column represents a data set for an antibiotic, with color-coding depicting different antibiotic classes. (b) Summary of prediction accuracy for different antibiotic mechanisms using the leave-one-out cross-validation method, with number of samples in depicted categories shown. FAS, fatty acid synthesis inhibitors; RNA, RNA polymerase inhibitors; PROTEIN, Protein synthesis inhibitors; Cell Wall, peptidoglycan biosynthetic enzyme inhibitors and other lipid II flux perturbing compounds; DNA, DNA gyrase/topoisomerase inhibitors and folate pathway inhibitors (See Supplementary Table 5 for predictions).

The lysocins were predicted to perturb Lipid II flux

The best test of a method for predicting antibiotic mechanism of action is whether it provides useful information on a compound having an unknown mechanism. In the process of isolating gliding bacteria to discover novel anti-infectives, we isolated a strain of Lysobacter enzymogenes and found strong antibacterial activity in crude extracts active against several Gram-positive bacteria. Following activity-guided fractionation, we isolated three natural products with m/z 809.4416 [M+2H]2+, m/z 816.4492 [M+2H]2+, and m/z 823.4570 [M+2H]2+. We did not initially recognize that one of these compounds corresponded to the recently reported lysocin E[33] and so we proceeded with purification, structure elucidation, and evaluation of the antibacterial profiles of all three compounds (Fig. 4a and Supplementary Note). After structure elucidation, the compound with m/z 809.4416 [M+2H]2+ (calculated for C75H118N20O20:809.4410) turned out to be lysocin E, while the compound with m/z 823.4570 [M+2H]2+ (calculated for C77H122N20O20: 823.4567) is lysocin I. The third natural product displaying m/z 816.4492 [M+2H]2+ (calculated for C76H120N20O20:816.4488) is a new lysocin derivative, lysocin J[33], that has so far only been described in a patent[34] with no evidence of its isolation. Lysocins are peptolides that share a common skeletal backbone consisting of 12 amino acids (2 x L-Thr, 2 x D-Arg, L-Ile, L-Leu, L-Ser, L-Glu, D-Gln, D-Trp, Gly and N-Me-D-Phe) where the absolute configuration was determined by Advanced Marfey’s Analysis[35]. Structurally, these three compounds differ in the length of the 3-hydroxy fatty acid chain linked to the amino group of L-Thr1.

Figure 4

The lysocins

(a) Structures of lysocin E, I, and J. (b) Mid-log cultures treated with lysocins at 2x MIC for four hours are clear, indicating lysis. Vancomycin (Van) is not bacteriolytic at the same concentration. (c) and (d) Mutants defective in menaquinone and heme biosynthesis are resistant to lysocin E, with the level of resistance inversely correlated with mutant fitness (Supplementary Fig 10). Panel (c) shows spot dilutions of the indicated strains in the absence (upper panel) and presence (lower panel; 2x MIC) of lysocin E. The small, colorless colonies display a small colony variant phenotype. (d) Summary of liquid MIC results with growth in each well normalized to the untreated control and presented as color density, with dark purple representing full growth. Each column represents two-fold dilutions from 40 μg/mL to 0 of lysocin E. The spot dilution assay and MIC tests are representative of three biologically independent experiments. (e) Antibiotics that affect the carrier lipid cycle of cell wall precursor biosynthesis can lead either to Lipid II depletion by sequestering Und-PP/Und-P or to Lipid II accumulation by preventing utilization of Lipid II. Examples of antibiotics that act at different steps are given.

Lysocin E was previously reported to be active against Gram-positive bacteria, including Mycobacterium spp., in vitro and in silkworm infection models[33]. We found in vitro activity in the same concentration range as reported and confirmed lysocin J and lysocin I as having similar in vitro potency (Supplementary Table 6). While active against M. smegmatis, S. aureus, S. epidermidis, and some other Gram-positive bacteria, the compounds were virtually inactive against Enterococcus spp., Streptococcus pneumoniae, Gram-negative bacteria and yeast/fungi (Supplementary Table 6). Due to their potent activity against laboratory S. aureus strains, including MRSA and VISA (vancomycin-intermediate S. aureus), we tested lysocins against clinical isolates of S. aureus and confirmed their potency (MIC50 1–2 μg/mL; n =15) (Supplementary Table 7). The already reported bactericidal and bacteriolytic effects of lysocin E[33] were also observed for lysocin J and lysocin I. Lysis was extremely rapid, with the density of an S. aureus culture (in cfu/mL) dropping by more than five logs in two hours (Fig. 4b and Supplementary Fig. 6). Strong bactericidal activity was also observed against mycobacteria (Supplementary Fig. 7). Despite these bacteriolytic properties, the compounds were not hemolytic against human red blood cells at concentrations well above the MIC (25x) against S. aureus (Supplementary Fig. 8). Therefore, the compounds are not non-specific membrane disruptors, but recognize features of bacterial cell surfaces as a prelude to lysis.

Based on high level resistance of menA null mutants with some supporting biochemical data,[33] it was proposed that the membrane target of lysocin E is menaquinone, an electron carrier important for aerobic synthesis of ATP. Mutants in menaquinone biosynthetic genes have a reduced membrane potential, low ATP levels, and other changes in metabolism that result in very slow growth; they also display reduced susceptibility to aminoglycosides and cell wall-active antibiotics, including the Lipid II binder nisin (Supplementary Fig. 9)[36-38]. We tested several other electron transport system mutants to determine if the resistant phenotype was specific to men pathway nulls. Five of nine mutants tested, menA, menB, menH, hemX, and hemB, showed lysocin E resistance in spot dilution assays (Fig. 4c and Supplementary Fig. 10a). In liquid culture, the menA, menB, and hemB mutants (Fig. 4d and Supplementary Fig. 10b), which have a small colony variant phenotype[37], showed high level resistance while the other two mutants, which formed larger (i.e., faster growing) colonies, showed lower level resistance. Therefore, the resistant phenotype was not specific to menaquinone depletion, but was also observed for mutants defective in the biosynthesis of heme, the cofactor used for the cytochromes that act as terminal oxidases in the electron transport chain. Therefore, we questioned whether the primary membrane target of the lysocins is, in fact, menaquinone. We tried to use macromolecular synthesis inhibition to identify the pathway affected by the lysocins, but lysis occurred so rapidly that all pathways were affected and the results were inconclusive, as reported [33].

We treated the transposon library with two concentrations of the three lysocin compounds to obtain mutant fitness profiles that might provide additional information about the mechanism of action. Profiles were obtained for two concentrations of each compound (0.25x and 0.5x MIC). Application of the machine learning algorithm predicted an “uncertain” mechanism for two of the samples, but the other four were predicted to perturb Lipid II pools. Compounds in the training set with this mechanism include bacitracin, which depletes Lipid II, and ramoplanin, which accumulates Lipid II. Although the lysocins are rapidly lytic, we observed a modest accumulation of Lipid II at sub-MIC concentrations (Supplementary Fig. 11a)[39]. Lipid II accumulation can only occur if MurJ, the Lipid II flippase is blocked, in which case accumulation is intracellular, or if Lipid II utilization on the cell surface is prevented because the compound in question somehow blocks the action of the polymerases that make peptidoglycan from Lipid II (Fig. 4e). We decided to test whether lysocin E prevents Lipid II polymerization using a prototypical peptidoglycan polymerase, S. aureus SgtB[40].

The lysocins bind to Lipid II

To obtain the necessary substrate to assess mechanism, we isolated S. aureus Lipid II using a two-step extraction procedure[41,42]. We then monitored its polymerization by SgtB at 0, 4 and 8 μM lysocin E as the Lipid II concentration was varied (Fig. 5a). The peptidoglycan polymer was labeled with biotin to enable detection[42]. There was no reaction in the presence of lysocin E until the Lipid II concentration exceeded half the inhibitor concentration (Fig. 5b, 5c and Supplementary Fig. 12). The distinctive shape of the inhibition curves is characteristic of compounds that bind Lipid II tightly: only after free substrate becomes available, which occurs at an inhibitor:substrate ratio that reports on the stoichiometry of binding, does peptidoglycan polymer form[43,44]. In the case of lysocin E, the inhibition curves are consistent with an antibiotic-to-substrate binding stoichiometry of 2:1.

Figure 5

Lysocin E binds to Lipid II

(a) Schematic of the enzymatic method to detect polymerized Lipid II. Products polymerized by SgtB are detected via Western blot after labeling with biotin-D-Lysine (BDL) using S. aureus PBP4. Keys are same as shown in Figure 1. (b) Blot showing that 8 μM lysocin E inhibits formation of peptidoglycan until the Lipid II concentration exceeds 4 μM. (c) Curves showing percent conversion to polymer as a function of Lipid II concentration for 0, 4, and 8 μM lysocin E. See Methods for details and Supplementary Fig. 12 for control and 4 μM blots. These curves are characteristic of a substrate binder with high affinity for Lipid II and a stoichiometry of 2:1 (inhibitor:Lipid II). (d) Schematic of the affinity capture assay developed to confirm binding of lysocin E to Lipid II. Lipid II biotinylated by PBP4 was bound to streptavidin beads, and then lysocin E was applied. (e) Total ion chromatograms showing that lysocin E bound to Lipid II beads could be eluted with 6 M GTC (top panel) or by including non-biotinylated Lipid II in the eluent (bottom panel). Menaquinone did not elute the antibiotic from the Lipid II affinity resin. These data are representative of at least three independent experiments.

We developed a new affinity capture assay to confirm direct binding of lysocin E to Lipid II. To do so, we prepared a Lipid II affinity resin by first exchanging the terminal D-Ala in S. aureus Lipid II for biotin-D-Lys using S. aureus PBP4[39,45] and then attaching biotinylated Lipid II to streptavidin-derivatized magnetic beads (Fig. 5d). Lysocin E was added to the beads, which were washed extensively. We found that we could elute intact lysocin E from the Lipid II affinity resin using 6 M guanidinium thiocyanate (GTC; Fig. 5e). In contrast, lysocin E added to control beds that lacked Lipid II remained in the flow-through (Supplementary Fig. 13). Lysocin E bound to the Lipid II affinity resin could be competed off the resin by adding Lipid II to the wash buffer, but not by adding menaquinone (Fig. 5e). As it was previously shown that lysocin E can interact with menaquinone[33], these results imply either that Lipid II and menaquinone can bind simultaneously to different parts of lysocin E in a non-competitive manner or, if the binding sites overlap, that menaquinone binds much more weakly than Lipid II. We note that the kinetic data and the observation that lysocin E is retained on the affinity column through repeated washings are consistent with a high affinity interaction with Lipid II. Taken together, our studies indicate that the molecular target of the lysocins in bacterial membranes is Lipid II, consistent with the predicted mechanism. The high level resistance of small colony variants lacking menaquinone or heme may in part reflect low Lipid II levels because we found this precursor was almost undetectable in these mutants (Supplementary Fig. 11b).

DISCUSSION

We have reported a functional genomics platform to predict compound mechanism of action that makes use of upregulation as well as inactivation mutant fitness data. Upregulation signatures were found more often in antibiotic samples treated at relatively high concentrations. In part, this is because we used a stringent method requiring strong positive selection pressure to identify upregulated genes. Work is underway to establish new computational methods that include statistically valid approaches to set cutoffs so that less obvious upregulation signatures can be confidently identified. Nevertheless, the results show that gene upregulation signatures in complex Tn-seq data can be used to identify both antibiotic targets and important mechanisms of resistance.

Several antibacterial compounds tested, including vancomycin, ramoplanin, lysobactin, amphomycin, and bacitracin, do not interact directly with protein targets, so gene upregulation signatures that report on targets would not be expected. Additionally, some antibiotics, e.g., β-lactams, may have several cellular targets, and upregulation of only one of them may not confer a sufficient fitness advantage to be detected. This may explain why we did not detect upregulation signatures for penicillin-binding proteins for the β-lactams we tested. Other antibiotics, e.g., ribosome binders, bind to multicomponent complexes, and upregulation of a single component may not provide a fitness advantage unless it is limiting for complex formation or can titrate the antibiotic away from the complex without creating toxicity. These inherent limitations in single-gene overexpression for target identification prompted us to explore alternative computational approaches to predict mechanism from inactivation mutant fitness data, and we found that a supervised machine learning approach worked well. As Tn-seq data provides information about the fitness of every non-essential gene in the system under a given perturbation, one might expect similar perturbations to give similar fitness profiles. The resulting data is likely more complicated, however, because some antibiotics have more than one target and some mutants may be responsive to a compound’s physical properties rather than its mechanism. For example, changes in cell envelope properties may increase fitness by reducing cell penetration, regardless of the specific molecular target of a compound[46]. Moreover, selection pressures can vary considerably with concentration, and technical challenges in standardizing selection pressures can complicate comparisons of antibiotics having the same mechanism. Having multiple concentrations for each antibiotic is ideal, but for this study we opted for a compromise between the number of antibiotics tested and the number of concentrations tested for each. Our validation rate of 72% compares favorably with a similar algorithm used in yeast[47], but we found that the prediction accuracy varied greatly depending on the target pathway, with predictions for fatty acid biosynthesis inhibitors and for RNA polymerase inhibitors being highly accurate and predictions for DNA synthesis inhibitors being largely incorrect. Evidently, changes in mutant fitness profiles are more diagnostic for some pathways than others and it may be possible to extract this diagnostic information from the data.

In a remarkable example of the fitness profiling approach, we predicted that the lysocins, a family of bacteriolytic lipodepsipeptides, perturb Lipid II pools, and follow-up experiments showed that these compounds are Lipid II binders. Our training set contained three cyclic peptides that bind Lipid II, and was generally overweighted in compounds that perturb Lipid II pools, but we do not think bias in the training data explained the lysocin predictions. The prediction accuracy for fatty acid and RNA polymerase inhibitors was unaffected by the large numbers of samples that affect Lipid II pools. Moreover, we had seven datasets in total for daptomycin and polymyxin, which are also cyclic peptides, and only one of the seven samples was predicted to perturb Lipid II pools. Therefore, the high prediction accuracy for the lysocins again suggests that the data contain information diagnostic for this mechanism.

Our studies indicate that selectivity of the lysocins for bacterial membranes derives from their ability to target Lipid II rather than menaquinone. Lipid II levels are very low in the menaquinone and heme pathway nulls, which may partly explain their lysocin resistance. However, we cannot exclude the possibility that the membrane potential plays another role in the mechanism of action of these compounds. There are many known Lipid II binders, but few are rapidly lytic except at concentrations well in excess of their MIC[48]. The sequelae to Lipid II binding vary depending on structural features of the antibiotics. Vancomycin, for example, acts as an enzyme inhibitor that prevents Lipid II utilization and glycan strand crosslinking, whereas nisin forms 2:1 complexes with Lipid II that assemble into higher order structures in the membrane, resulting in depolarization[49]. For nisin and some other Lipid II binders, especially those having a net positive charge, it is possible that there is an interaction with the membrane potential that facilitates their mechanistic effects[48,49]. Further study will be required to understand what happens after Lipid II binding that explains the lytic behavior of the lysocins.

To summarize, we have demonstrated that Tn-seq libraries containing both upregulation and inactivation mutants provide a rich source of information about compound mechanism of action. While we have treated upregulation and inactivation mutants separately here, we suggest that prediction accuracy can be improved by developing a computational approach that combines upregulation and inactivation mutant fitness data to generate a prediction that factors in both. We also note that the ability to detect resistance mechanisms in addition to targets, either through upregulation signatures or increased fitness of inactivation mutants, may be useful in guiding antibiotic development.

ONLINE METHODS

Materials

Native Lipid II was prepared from Staphylococcus aureus as previously described[41,42]. S. aureus SgtB Y181D and S. aureus PBP4 were purified as reported previously[40,50]. Biotinylated D-Lysine (BDL) was prepared as described[39]. Non-stick conical vials and pipet tips for enzymatic reactions and Fmoc-D-Lys(biotinyl)-OH were purchased from VWR. All other chemicals were purchased from Sigma Aldrich unless otherwise noted.

Treatment and sequencing of the transposon library

We created a transposon library in S. aureus HG003 by phage-based transposition of six different transposon constructs as previously described[11]. Library aliquots stored at −80°C were thawed, diluted to OD600 = 0.2, and incubated for 1 hour with shaking at 30°C. This culture was diluted to 2 × 105 cfu/mL in cation-supplemented TSB (TSB with 25mg/L Ca2+ and 12.5mg/L Mg2+) with the desired antibiotic as described previously[13]. Samples at several antibiotic concentrations above and below the MIC were incubated with shaking at 37°C to OD600 of 1.

For each antibiotic, we processed at least two samples, one that displayed a moderate (three to five hours) and one that displayed an extreme (~20 hours) delay in growth to stationary phase. Samples were prepared for NGS as previously described[11] and sequenced at the Harvard Biopolymers Facility or at the Tufts University Core Facility for Genomics on a High-Seq 2500. The sequencing data was split by transposon construct barcode, filtered for quality, trimmed down to the genomic sequence, and mapped to the genome using the Galaxy software suite[51,52]. SAM files were converted using Tufts Galaxy Tn-Seq software (http://galaxy.med.tufts.edu/) into tab-delimited hopcount files, which were converted into IGV-formatted files using in-house python scripts as previously described[11,13].

Identification of upregulated genes

To identify upregulated genes, the IGV data sets were analyzed using custom python scripts (Supplementary Fig. 1). We separated the sequencing reads by the genomic strand they mapped to and identified the 100–200 TA sites that had the greatest increase in reads in one strand in the treated sample compared to the control. We defined an upregulation signature as one 135 bp genomic window or two consecutive windows with at least three TA sites enriched in the same strand with no intervening sites enriched in the opposing strand. An upregulation signature was considered a hit if there was a proximal gene in the same orientation as the promoter within the transposon

Using machine learning to predict mechanism of action

To create a tool to predict the mechanism of action of unknown antibiotics, we first selected training datasets. We excluded compounds without specific targets or with uncertain mechanisms (daptomycin, polymyxin, CCCP, and gramicidin A). We also removed data files that had too little or too much antibiotic selective pressure, defined as those with fewer than 10 genes or more than 1000 genes at least 4-fold depleted in reads compared to the control file. For these calculations, the control data was an average of untreated datasets from 14 separate experiments to minimize the effects of jackpots. After this pruning, the training dataset had 43 files representing 25 antibiotics. Next, we developed gene fitness scores to compare the relative importance of genes under the different treatments. To avoid highlighting cases in which small read differences yield large fold-changes, we adjusted the read counts of genes with fewer than 0.0001x the total number of reads in the sample up to that value. After this correction, we calculated an initial fitness value (f) by dividing the number of reads in a gene by its length (l) and multiplying by the the total treated: control reads ratio. We then ordered genes by f and ranked them from 0 to 1, with uniform intervals to normalize the fitness values. We have included a supplemental comma-separated file with the fitness values for each gene under each antibiotic treatment, along with files tabulating the gene-by-gene sequencing read counts for the antibiotic-treated and untreated samples (Supplementary Data).

We then set up a K-nearest neighbors algorithm. From each file in the training set, our custom python program selected the 10 genes most depleted of reads compared to the control, excluding genes with fewer than 100 reads in the control,,and the 15 genes most enriched in sequencing reads, provided they were at least 4-fold enriched and had at least 100 reads. Once the gene set had been curated, the program arrayed the normalized fitness values for each gene from each training set file. This “fingerprint” array and the mechanism of action categories listed in Supplementary Table 3 were then fed into the K-nearest neighbors (KNN) function from the Scikit-learn Python library[53]. We chose to use the three nearest neighbors for the prediction, as a smaller k led to more ties between categories and a higher k biased the tool to favor larger categories.

For leave-one-out validation, each dataset in the training set was re-labeled as an unknown, the genes contributed solely by the left-out file were removed from the “fingerprint” array, and the mechanism of action was predicted using the KNN algorithm (Supplementary Table 4). The algorithm reports up to three possible mechanisms for the antibiotic and their respective probabilities. We classified the validation results as correct or incorrect based on the most probable predicted mechanism. However, if the highest probability was less than 0.4, the prediction was labeled as uncertain.

Semi quantitative RT-PCR validation of gene upregulation

Selection for resistant mutants

~107 HG003 cells were plated on TSB agar containing antibiotics: 1 μg/mL DMPI, 1 μg/mL CDFI, 50 μg/mL fosfomycin, or 50 μg/mL D-cycloserine. After incubation at 37°C for 16–20 hours, ~20 colonies were selected from each plate and their resistant phenotypes were reconfirmed. The colonies were screened by PCR to identify transposon insertions in the target genes: murJ (SAOUHSC_01871), fosB (SAOUHSC_02609), and ddL (SAOUHSC_02318). The promoter regions of two mutants from the CDFI/DMPI, fosfomycin, and D-cycloserine plates were sequenced to identify the transposon insertions.

Semi quantitative RT-PCR to measure expression of target genes

Wildtype HG003 and one confirmed mutant from each selection were grown in TSB overnight at 37°C, diluted 100-fold in fresh TSB, and then incubated to OD600 ~ 0.5–0.7. Then, total RNA from each culture was isolated using the RNeasy Kit (Qiagen) and quantified by NanoDrop. Using specific primers for each target gene and gyrA (as a control), semi-quantitative qRT-PCR was performed with One-step real-time RT-PCR (Thermo Fisher Scientific). The PCR settings were the following: 1) reverse transcription (48°C/30min), 2) PCR amplification (94°C/1min, (94°C/15s, 50–65°C/30s, 68°C/1min per 1kb)x40, 68°C/5min). Primer sequences (5′→3′) are as follows: murJ-for: GACAGCAGGTGTACCATTAG, murJ-rev: CTTGGAATATCCCTCTCCATGTC, fosB-for: GCAGGCCTATGGATTGCTTTA, fosB-rev: GTTCTCAAGTGTGCCAGTATGT, ddL-for: CCAGCTGACTTAGACGAAGATG, ddL-rev: CTGTTTATCCTGGTGACGTTCT, gyrA-for: CGGTGTCATACCTTGTTC, gyrA-rev: GTGTTATCGTTGCTCGTG. Reaction products were visualized on a 2% agarose gel with ethidium bromide staining and analyzed using ImageJ 1.48v.

Bacterial growth curves

For validation of mutants found to confer resistance to daptomycin in the TnSeq data, transposon mutants from the Nebraska transposon library[54] were transduced into HG003 using Φ11. For the upregulation constructs, SAOUHSC_00969, SAOUHSC_02149, and SAOUHSC_02164 were amplified using PCR and cloned into the pLOW plasmid using standard methods[55], followed by transduction into HG003. Overnight cultures of all HG003 strains were diluted to an OD600 = 0.001. 150 μL of diluted culture was mixed with 1.5 μL of antibiotic at 100X final concentration in a well of a 96 well plate. The plate was incubated at 37°C with shaking for at least 16 hours, and bacterial growth was monitored by measuring the OD600 in a plate reader (Spectra Max plus384, Molecular Devices).

Production of Lysocins

Lysobacter enzymogenes was isolated and cultivated using previously reported methods[56]. The strain was sequenced by the Illumina paired-end technology on a MiSeq instrument at the Helmholtz Centre for Infection Research and assembled with Abyss-pe 1.3.6 assembler to produce 223 contigs with a total length of 6.1 Mbp. The strain was subsequently assigned to the L. enzymogenes species based on 100% identity across 100% query coverage to the 16S gene of the complete genome sequence of L. enzymogenes M497-1 (accession ID: AP014940). The preculture was used to inoculate 300 mL shake flasks containing 50 mL CY medium (0.3% casiton, 0.1% yeast extract, 0.1% CaCl2, 50 mM HEPES, adjust pH to 7.2 with 10 N KOH) that was incubated at 30°C for 18 hours. Six 5 mL aliquots were used to inoculate six times 2 L of CY medium at 30°C for 3 days at 160 rpm. 5% of XAD adsorber resin (Amberlite-XAD-16, Sigma) was added, and the flasks were shaken for a further 18 hours. The crude extract prepared from a small-scale cultivation (8 L) in Cy-medium displayed strong antibacterial activity against several Gram-positive bacteria.

Isolation of lysocin compounds

The culture broth containing cells and XAD resin was centrifuged at 25°C, 7000 rpm for 30 minutes. The freeze-dried resin-cells mixture was extracted with methanol several times. The methanol solution was then defatted with hexane and the methanol layer was evaporated to dryness. The crude residue was re-suspended in DMSO:MeOH (1:1) and separated by preparative HPLC. A Phenomenex Kinetix 5 μm biphenyl 100 Å 250 × 21.2 mm was used and a gradient from 5 to 95% B in 40 minutes with (A) H2O + 0.1% FA and (B) ACN + 0.1% FA at a flow rate of 25 mL/minute at room temperature. Fractions were collected by MS fractionation. The fractions containing compounds of interest were repurified by semipreparative RP-HPLC (Dionex Ultimate 3000) equipped with a Phenomenex Synergi Polar RP 80 Å 250×10 mm, 4 μm and using a gradient of 5 – 95% B over 55 minutes with (A) H2O + 0.1% FA and (B) ACN + 0.1% FA at a flow rate of 5 mL/minute at 40°C. Elution was monitored at 220 nm.

General derivatization protocol for Advanced Marfey’s analysis

Advanced Marfey’s analysis[35] was used to determine the The absolute configuration (D/L) of the amino acids.

General method for identification of lysocins

The measurements to detect all lysocin derivatives were performed on a Dionex Ultimate 3000 RSLC system using a BEH C18, 100 × 2.1 mm, 1.7 μm dp column (Waters, Germany). Separation of 1 μL sample was achieved by a linear gradient from (A) H2O + 0.1% FA to (B) ACN + 0.1% FA at a flow rate of 600 μL/minute and 45°C. The gradient was initiated by a 0.5 minute isocratic step at 5% B, followed by an increase to 95% B in 18 minutes to end up with a 2 minute step at 95% B before reequilibration under the initial conditions. UV spectra were recorded by a DAD in the range from 200 to 600 nm. The LC flow was split to 75 μL/minute before entering the maXis 4G hr-ToF mass spectrometer (Bruker Daltonics, Germany) using the Apollo ESI source. Mass spectra were acquired in centroid mode ranging from 150 – 2500 m/z at a 2 Hz scan rate. Settings for MS/MS measurements were: minimum precursor intensity is set to 10000. Full scan spectra are acquired at 2 Hz followed by MS/MS spectra acquisition at variable scan speed ranging from 1 to 3 Hz, as a function of precursor intensity. CID energy varies linearly from 30, 35, 45, to 55 eV with respect to the precursor m/z from 300, 600, 1000, to 2000 m/z. The collision cell is set to ramp collision energy (80–120% of the set value with equal weights of both values), collision RF (700 to 1000 Vpp with equal weights of both values) and transfer time (90–110 μs) for every MS/MS scan. The number of precursor was set to 2 and precursors were moved to an exclusion list for 0.2 minute after two spectra were measured (typical chromatographic peak width was 0.10–0.15 minute). Precursors were reconsidered if their intensity changed fivefold.

HPLC-MS analysis of lysocin

Measurements were performed using a maXis system (Dionex Ultimate 3000 RSLC system and maXis 4G hr-ToF mass spectrometer (Bruker Daltonics, Germany) using the Apollo ESI source) or an Orbitrap system (Dionex UltiMate 3000 RSLC system coupled with an Advion Triversa Nanomate nano-ESI system attached to a Thermo Fisher Orbitrap). LC flow was split to 500 nL/minute before entering the ion source. Mass spectra were acquired in centroid mode ranging from 150–1000 m/z, resolution R = 30000. A Waters BEH C18, 100 × 2.1 mm, 1.7 μm dp column was used, injection volume = 1 μL. A gradient of A) H2O + 0.1% FA and B) MeCN + 0.1% FA at a flow rate of 0.55 mL/minute was used to achieve separation. Gradient conditions: start at 5% B increasing to 10% B in 1 minute, increase to 35% B from minute 1→15, increase to 50% B from minute 15→22, increase to 80% B from minute 22→25. After 1 minute hold at 80% B, the system was reequilibrated for 5 minute with the initial conditions. UV data was acquired at 340 ± 8 nm, MS detection was performed simultaneously.

Antimicrobial screening

Minimum inhibitory concentration (MIC) and MIC50 were determined as described previously[57]. Three laboratory strains of S. aureus (Newman, USA300, COL) and 12 clinical S. aureus isolates (kindly provided by Prof. Dr. Markus Bischoff; Saarland University Hospital, Homburg, Germany) with various resistance phenotypes were used for determining MIC50 values (defined as lowest concentration inhibiting the growth of 50% of all tested organisms).

Bactericidal assessment of lysocin compounds

To assess bactericidal activity of the lysocin compounds, HG003 cultures grown in Mueller Hinton II broth at 37°C to mid-log phase (OD600=0.5) were treated with 2 μg/mL vancomycin (2xMIC), 2 μg/mL lysocin compounds (2xMIC), or no antibiotic. Bacterial viability was monitored by CFU counting hourly for 24 hours. To compare bactericidal activity with other known antibiotics, time-kill experiments were performed as described[58] with minor modifications. Briefly, overnight cultures of S. aureus or M. smegmatis were diluted in their respective growth medium to achieve an initial inoculum of 109 or 107 CFU/mL, respectively. Lysocins were added at the assigned concentration and cell viability was assessed at time points as noted. For S. aureus, daptomycin and vancomycin were chosen as controls. Clofazimine served as control for M. smegmatis.

Spot dilution assay to determine antibiotic susceptibility

Transposon mutants (qoxA, qoxB, qoxC, ndh, hemB, hemX, and menH) in USA300 from the Nebraska Transposon Mutant Library and ΔmenA and ΔmenB Newman strains[59] were used in antibiotic susceptibility assays. Cultures of mutant and wildtype strains grown at 37°C to mid-log phase (OD600=0.5) were diluted in ten-fold increments as indicated. 5 μL from each dilution was spotted on TSB agar containing the following antibiotics: 1 μg/mL of gentamicin A; 28 μg/mL of nisin for USA300 strains, 38 μg/mL for Newman strains; 2 μg/mL of lysocin E. The plates were incubated at 37°C overnight before imaging. To test susceptibility of lysocin E in liquid media, overnight cultures of the mutants were inoculated in 150 μL TSB containing lysocin E at various concentrations as indicated (0 to 40 μg/mL) in a 96-well plate, the plate was incubated at 37°C overnight, and cell growth (OD600) was measured.

Lipid II accumulation assay

The Lipid II assay was performed as previously described[39,43]. 2 mL of S. aureus cultures grown at 37 °C to mid-log phase (OD600= 0.4~0.5) were treated with the following antibiotics for 10 min at 37 °C: 0.3 μg/mL moenomycin, 4 μg/mL vancomycin, 0.2 μg/mL lysocin E. Samples were normalized based on OD600 and immediately harvested by centrifugation, and total lipid was extracted as described and resuspended in 20 μL DMSO. S. aureus Lipid II in the lipid extract was labeled with biotin-D-lys(BDL) by S. aureus PBP4 as described[39,43]. After 1 hr at room temperature, reactions were quenched by adding an equal amount of 2x SDS loading buffer, and then resolved on a 4–20% gradient SDS polyacrylamide gel followed by western blot analysis as described.

Determining the inhibitory activity of lysocin E on Lipid II polymerization

The S. aureus SgtBY181D mutant was purified and used in the Lipid II polymerization assay as previously described[40,50]. This mutant makes shorter polymers than wildtype SgtB, which allows for better quantification of products. For Lipid II polymerization, native Lipid II (at concentrations shown elsewhere) and 2 μM SgtB were added to PGT buffer (12.5 mM HEPES pH 7.5, 20 mM MnCl2, 2.5 mM Tween-80) containing 3 mM biotin-D-Lys (BDL). Various concentrations of lysocin E (0, 4, or 8 μM) were added, and the reactions were incubated for 10 minutes at room temperature. Then, 4 μM S. aureus PBP4 was added and incubated for 60 minutes to label the polymers with BDL. Reactions were quenched by adding 10 μL 2xSDS loading buffer, and then resolved on a 4–20% gradient SDS polyacrylamide gel followed by western blot analysis as described above. The densitometric quantitation for each lane was performed using Image J (http:rsbweb.nih.gov/ij). To calculate percent conversion of Lipid II in each reaction, the intensity from each lane was divided by the intensity from the reaction with 16 μM Lipid II.

Lipid II bead binding assay

Lipid extraction

An overnight culture of S. aureus RN4220 was diluted 100-fold in 500 mL TSB broth and grown at 37°C to exponential phase. At exponential phase (OD600=0.6), 5 μg/mL (final concentration) moenomycin was added and the culture was further incubated at 37°C for 30 minutes. Bacterial cells were harvested by centrifugation at 4000 rpm for 10 minutes and resuspended in 15 mL PBS buffer pH 7.4. Total lipid was extracted by adding 75 mL of 1:2 chloroform:methanol to the cell suspension. The mixture was stirred for one hour at room temperature. Cell debris was separated by centrifugation at 4000 rpm for 10 minutes, and the supernatant was transferred to 25 mL of 1:1 chloroform:PBS buffer pH 7.4. The mixture was stirred for one hour at room temperature and centrifuged at 4000 rpm for 10 minutes to separate the phases. The organic phase was transferred into a new glass tube, dried using rotary vacuum, and resuspended in DMSO[39,43].

Lipid II labeling with BDL

In a total reaction volume of 1 mL, 100 μL of the total lipid extract was incubated with 10 μM PBP4 and 3 mM BDL in the PGT buffer for 2 hours at room temperature. To remove free BDL, 1 mL chloroform was added to the reaction and vortexed for 10 minutes at room temperature. The organic phase containing BDL-labeled Lipid II was separated by centrifugation at 3000 rpm at room temperature and transferred into a glass vial, then dried and dissolved in 100 μL of DMSO[39,43].

Lysocin E capture and release from Lipid II bead

100 μL of streptavidin magnetic beads (Dynabeads M-280 Streptavidin, Fisher scientific) was washed with 1 mL of PBS buffer pH 7.4 twice and resuspended in 100 μL PBS buffer pH 7.4. 10 μL of biotinylated-Lipid II was added to the bead mixture and incubated at room temperature for 1 hour to allow capture of BDL-Lipid II on the streptavidin beads. The beads were then washed with 1 mL PBS buffer pH 7.4 twice and resuspended in 100 μL PBS buffer pH 7.4 supplemented with 0.01% bovine serum albumin (BSA). 10 μL of 1 mg/mL lysocin E was added to this bead mixture and incubated at room temperature for 1 hour, followed by washing with 500 μL of PBS buffer pH 7.4 three times. To recover lysocin E, beads were treated with 50 μL of 6 M guanidinium thiocyanate, S. aureus Lipid II (40 μM or 80 μM), or menaquinone-4 (80 μM) as indicated and incubated at room temperature for 10 minutes, after which the supernatant was subjected to LC/MS analysis. LC/MS analysis was performed with ESI-MS operating in positive ion mode using an Agilent Technologies 6120 Quadrupole LC-MSD instrument. The collected fractions were separated using a 16-minute HPLC method with an SB-C18 column (1.8 μm, 2.1 × 50 mm, Agilent) at a constant flow rate of 0.2 mL/minute. The gradient method was set as follows: 0–10 minutes, 95% solution A (0.1% formic acid in H2O)/5% solution B (0.1% formic acid in acetonitrile); 10–11 minutes, 50% solution A/50% solution B; 11–14 minutes, 10% solution A/90% solution B; 14–16 minutes, 95% solution A/5% solution B. Ions corresponding to expected compound were identified in the chromatograms.

Red blood cell (RBC) lysis

To assess hemolytic activity of lysocin compounds, human red blood cells were treated with compounds, and lysis of the red blood cells was quantified as described[60]. Briefly, isolated red blood cells (RBCs) from defibrinated whole human blood were washed three times with a buffer containing 10 mM Tris-HCl (pH 7.4) and 0.9% NaCl, and then diluted to 5% in the same buffer. Lysocin compounds and ramoplanin were added to 300 μL of the buffer with 1% RBC as previously described, and these mixtures were incubated for 30 minutes at room temperature. The mixtures were centrifuged at 1300xg for 5 minutes and 250 μL of the supernatant was transferred to a 96 well plate. The released hemoglobin from lysed RBCs was measured at OD540 using a plate reader (SpectraMax, Molecular Dyanmics).

Data availability

All data generated or analyzed during this study are included in this published article, its supplementary information files, and the associated SRA entry or are available from the corresponding author on reasonable request.