RESEARCH ARTICLE Sánchez-­Osuna et al., Microbial Genomics 2020;6 DOI 10.1099/mgen.0.000440 OPEN OPEN DATA ACCESS Exploration into the origins and mobilization of di-­hydrofolate reductase genes and the emergence of clinical resistance to trimethoprim Miquel Sánchez-­Osuna1, Pilar Cortés1, Montserrat Llagostera1, Jordi Barbé1† and Ivan Erill2,*,† Abstract Trimethoprim is a synthetic antibacterial agent that targets folate biosynthesis by competitively binding to the di-­hydrofolate reductase enzyme (DHFR). Trimethoprim is often administered synergistically with sulfonamide, another chemotherapeutic agent targeting the di-­hydropteroate synthase (DHPS) enzyme in the same pathway. Clinical resistance to both drugs is wide- spread and mediated by enzyme variants capable of performing their biological function without binding to these drugs. These mutant enzymes were assumed to have arisen after the discovery of these synthetic drugs, but recent work has shown that genes conferring resistance to sulfonamide were present in the bacterial pangenome millions of years ago. Here, we apply phylogenetics and comparative genomics methods to study the largest family of mobile trimethoprim-­resistance genes (dfrA). We show that most of the dfrA genes identified to date map to two large clades that likely arose from independent mobilization events. In contrast to sulfonamide resistance (sul) genes, we find evidence of recurrent mobilization in dfrA genes. Phyloge- netic evidence allows us to identify novel dfrA genes in the emerging pathogen Acinetobacter baumannii, and we confirm their resistance phenotype in vitro. We also identify a cluster of dfrA homologues in cryptic plasmid and phage genomes, but we show that these enzymes do not confer resistance to trimethoprim. Our methods also allow us to pinpoint the chromosomal origin of previously reported dfrA genes, and we show that many of these ancient chromosomal genes also confer resistance to trimethoprim. Our work reveals that trimethoprim resistance predated the clinical use of this chemotherapeutic agent, but that novel mutations have likely also arisen and become mobilized following its widespread use within and outside the clinic. Hence, this work confirms that resistance to novel drugs may already be present in the bacterial pangenome, and stresses the importance of rapid mobilization as a fundamental element in the emergence and global spread of resistance determinants. DATA SUMMARY INTRODUCTION • Nucleotide and protein sequences analysed in this study Bacterial resistance to antibacterial agents remains an increas- ingly challenging and global problem in modern health care have been downloaded from publicly available National [2, 3]. Bacterial cells display a diverse array of mechanisms Center for Biotechnology Information databases. to cope with exposure to antibacterial compounds, including • The scripts used for data collection and analysis can be modification or overexpression of the antibacterial target, obtained at the GitHub ErillLab repository (https://​github.​ efflux or reduction of antibacterial uptake and the use of com/​ErillLab/). alternate pathways [4]. Constant exposure to non-­lethal concentrations of antibacterial agents may lead to the selec- • The Bayesian phylogenetic tree can be visualized online on tion of partial resistance to antibiotics over relatively short iTOL (https://​itol.​embl.​de/​tree/​8556​7415​9451​5851​33078) time spans [5], and this evolution may be hastened by simul- [1]. taneous exposure to multiple antibacterials [6]. However, Received 27 April 2020; Accepted 08 September 2020; Published 24 September 2020 Author affiliations: 1Departament de Genètica i de Microbiologia, Universitat Autònoma de Barcelona, Bellaterra, Spain; 2Department of Biological Sciences, University of Maryland, Baltimore County, Baltimore, MD, USA. *Correspondence: Ivan Erill, ​erill@​umbc.​edu Keywords: trimethoprim; antibiotics; chemotherapeutic agent; resistance; sulfonamides; phylogenetics; evolution. Abbreviations: MIC, minimum inhibitory concentration; NCBI, National Center for Biotechnology Information; PABA, para-­aminobenzoic acid. †These authors share senior authorship. Data statement: All supporting data, code and protocols have been provided within the article or through supplementary data files. Sixteen supplementary tables and five supplementary figures are available with the online version of this article. 000440 © 2020 The Authors This is an open-­access article distributed under the terms of the Creative Commons Attribution License. 1 Sánchez-­Osuna et al., Microbial Genomics 2020;6 the rapid proliferation of multidrug-­resistant nosocomial pathogens in the last 50 years has not been driven by the Impact Statement independent evolution of resistance traits, but through the Antibiotic resistance is a pressing and global phenom- extensive dissemination of mobile genetic elements carrying enon. It is well established that resistance to conven- resistance genes [4, 7]. It is widely accepted that most genes conferring resistance to antibiotics present in pathogenic tional antibiotics emerged millions of years ago in bacteria were obtained by successive lateral gene transfer either antibiotic-­producing bacteria or their competitors. of homologues that originally evolved in the microbes that Resistance to synthetic chemotherapeutic agents cannot produce the antibiotic or in their natural competitors [7, 8]. be explained by this paradigm, since these drugs are The high plasticity of bacterial genomes, enabled by a large not naturally produced. Hence, resistance is assumed to repertoire of mobile genetic elements, and the availability of have evolved rapidly following the clinical introduction of a large pool of ancient antibiotic-­resistance determinants, these drugs. Recently, we showed that resistance to one hence, set the stage for the rapid proliferation of antibiotic such drug, sulfonamide, evolved not recently, but millions resistance, giving rise to multi-­resistant clinical strains just a of years ago, suggesting that the diversity of bacterial few years after the commercial introduction of antibiotics [7]. genomes may well contain genes conferring resistance to drugs yet to be developed. Here, we analyse the origin Synthetic chemotherapeutic agents predate antibiotics in the clinical setting, and continue to be used synergistically of resistance to trimethoprim, another chemothera- with antibiotics to treat microbial infections [9]. Following peutic agent developed in the 1960s. Using phylogenetic the initial discovery and clinical use of arsphenamine in methods, we identify new variants of the trimethoprim-­ 1907 [10], interest in chemotherapeutic agents quickly took resistance genes that had not previously been reported, off after the development of sulfa drugs in the 1930s [11]. and we trace the chromosomal origins for a number of The discovery of trimethoprim (a di-­aminopyrimidine) was already known resistance variants. Our results show that received with interest because, like sulfonamides, trimetho- resistance to trimethoprim is very diverse, and has origi- prim targets the bacterial synthesis of tetrahydrofolic acid, nated both from recent mutations and from pre-­existing which is a necessary cofactor in the synthesis of thymine ancient variants. These results stress the importance of and purines [12]. Sulfonamides are structural analogues of gene mobilization mechanisms as the main drivers of the para-­aminobenzoic acid (PABA) and inhibit the synthesis current antibiotic-­resistance phenomenon. of di-­hydropteroate by competing with PABA for binding to the di-­hydropteroate synthase (DHPS) enzyme, resulting in sulfonamide-­bound di-­hydropterin [13]. Trimethoprim is a structural analogue of di-­hydrofolic acid, derived from mutations, but to the acquisition of resistance determinants di-­hydropteroate. It acts by competitively binding to the on mobile genetic elements [13]. Parallel to their systematic di-­hydrofolate reductase (DHFR) enzyme; hence, inhibiting combined use in both clinical and agricultural settings, genes the production of tetrahydrofolic acid [13, 14]. The synergistic conferring resistance to sulfonamides and trimethoprim are use of trimethoprim and sulfonamides was expected to have frequently found together on mobile elements, such as class a potent bactericidal effect by producing a serial blockade on 1 integrons [19] or conjugative plasmids [13, 20]. The mobile the tetrahydrofolic acid pathway [12, 15]. genes conferring resistance to sulfonamide are homologues of the chromosomally encoded folP gene and are collectively Unlike antibiotics, chemotherapeutic agents are not produced known as sul genes (for sulfonamide resistance). Mobile genes by natural organisms, yet resistance to these novel drugs arose conferring resistance to trimethoprim are either homologues quickly after their mass-­production and it is today pervasive or functional analogues of the chromosomally encoded among clinical isolates [7]. In the case of sulfonamides and folA gene and are collectively known as dfr genes (for di-­ trimethoprim, which are usually administered in tandem, hydrofolate reductase) [17]. resistance via chromosomal mutations to both chemothera- peutics was reported soon after their clinical introduction In spite of their frequent co-­occurrence on mobile genetic [13]. Chromosomal resistance to sulfonamides can occur elements, there are significant differences between the mobile through mutations yielding increased production of PABA genes conferring resistance to sulfonamides (sul genes) and [16] or, more commonly, via mutations to the chromosomal trimethoprim (dfr genes). To date, only three sul gene classes folP gene (encoding DHPS), which decrease the affinity of have been described in clinical isolates [21], whereas more DHPS for sulfonamide without detriment to PABA binding than 30 different dfr genes have been reported in clinically [13, 17]. Such mutations have been reported in multiple relevant strains [22]. Trimethoprim-­resistance (dfr) genes bacterial groups and target different conserved regions of have been further classified into two families (dfrA and dfrB). DHPS [13]. Similarly, chromosomal resistance to trimetho- These two families encode evolutionarily unrelated proteins prim may arise via mutations that increase transcription of of markedly different sizes. Sequence similarity indicates that the folA gene (encoding DHFR) [18], or through mutations dfrA genes are homologues of the chromosomally encoded that decrease the affinity of DHFR for trimethoprim [13]. The folA genes, whereas dfrB genes are functional analogues of vast majority of resistant clinical isolates to both sulfonamides unknown origin [23, 24]. Most dfrA genes follow a standard and trimethoprim, however, are not due to chromosomal naming convention consisting of dfrA followed by a numerical 2 Sánchez-­Osuna et al., Microbial Genomics 2020;6 value indicating their discovery rank order. However, several Center for Biotechnology Information (NCBI) GenBank/ dfrA genes first identified in Gram-­positive bacteria, and RefSeq sequences corresponding to mobile genetic elements. thought at the time to be unrelated to the Gram-­negative These were defined as sequences containing the keywords dfrA genes, were originally named following an alphabetical ‘plasmid’, ‘integron’ or ‘transposon’ in their title, belonging convention (dfrC–K). The disparity in genetic diversity to complete genome records [28, 29]. Protein records corre- among sulfonamide and trimethoprim mobile resistance sponding to blastp hits matching stringent value (<1e−20) and determinants is suggestive of different evolutionary processes query coverage (>75 %) thresholds were added to the panel leading to the onset and spread of resistance to these two if non-­redundant (amino acid sequence identity <90 % with chemotherapeutic agents [13]. It was suggested that resistance respect to existing panel members), and their classification as to sulfonamide had arisen in a few isolated organisms and mobile elements was validated by assessing that the nucleotide rapidly spread upon the introduction of sulfa drugs, whereas record encoding them contained at least one gene encoding trimethoprim resistance had independently evolved, and had an integrase, transposase or plasmid replication protein, as been subsequently mobilized multiple times [13]. determined by hmmer (hmmscan, E value <1e−05) with refer- Recently, we examined the origins of sul genes through ence pfam models (Table S3) [30–34]. To detect putative chro- comparative genomics, phylogenetic analysis and in vitro mosomally encoded folA genes associated with mobile dfrA assays [25]. We identified a well-­defined mutational signa- genes, we used the sequences in the extended non-­redundant ture in sul-­encoded proteins with respect to chromosomally panel of DfrA homologues as queries for tblastn searches encoded folP genes, and we used this conserved motif to map against NCBI GenBank complete genomes with stringent E the origins of sul genes in bacterial chromosomes. Our work value (<1e−40) and query coverage (>75 %) settings. Hits with revealed that the three groups of sul genes identified in clinical nearby genes annotated as resistance determinants, trans- isolates originated in the Leptospiraceae and were transferred posases or integrases were considered to encode chromosom- to the Rhodobiaceae more than 500 million years ago. These ally integrated mobile DfrA homologues and not considered two ancient resistant determinants were later independently for further analysis. For each mobile DfrA homologue in the mobilized, and rapidly disseminated following the commer- panel, the first, if any, tblastn hit satisfying these thresholds cial introduction of sulfa drugs. By tracing the phylogenetic was considered, for the purposes of this study, to be a proxy lineage of sul genes and demonstrating that these two bacte- for the closest putative chromosomally encoded FolA protein. rial families were resistant to sulfonamides long before their The choice of representative DfrA sequences did not alter the discovery and clinical use, our work indicated that resistance tblastn results. To complete the panel of protein sequences to novel drugs could very well pre-­exist, and be ready for used to reconstruct the evolutionary history of DfrA/FolA mobilization, within the vast bacterial pangenome. Here, we sequences, we used the non-­redundant panel of mobile DfrA apply similar methods to examine the phylogenetic relation- sequences to identify via blastp (E value <1e−20, coverage ships among dfrA and chromosomally encoded folA genes. >75 %) FolA proteins encoded by the chromosomes of NCBI Our aim is to shed light on the evolutionary processes giving RefSeq representative species for all bacterial orders, and for rise to mobile trimethoprim-­resistance genes. Our work illus- each bacterial family in the Proteobacteria. In addition, the trates significant similarities and differences in the processes closest archaeal homologues of bacterial FolA sequences were leading to the emergence and spread of trimethoprim- and determined by searching with blastp the NCBI protein data- sulfonamide-­ resistance determinants, reveals previously base, restricted to Archaea (taxid:2157), with the Escherichia unreported clusters of dfrA genes, and suggests that system- coli FolA protein. A member of each family from the order atic analyses of the bacterial pangenome may be of use in the (Halobacteriales) of the identified best archaeal hit of E. coli design of novel antibacterials. FolA was sampled to populate the outgroup. Phylogenetic analysis METHODS For phylogenetic inference, we performed a t-­coffee Sequence data collection multiple sequence alignment of protein sequences for the To identify homologues of DfrA proteins, we first compiled a complete panel of DfrA and FolA homologues, combining panel of Dfr proteins reported in the literature. Dfr proteins three clustalw (version 2.1) profile amino acid sequence belong to two distinct families of unrelated sequences (DfrA alignments with different (5, 10, 25) gap opening penalties and DfrB; Fig. S1, available with the online version of this and leveraging the E. coli FolA crystal structure (P0ABQ5) to article). We mapped these sequences to pfam models of DfrA adjust gap penalties [35]. The resulting amino acid sequence (PF00186) and DfrB (PF06442) (Table S1) using hmmer alignment was processed with Gblocks version 0.91b (allowed version 3.1b2 (hmmscan) [26], and we discarded sequences gap positions, with half; minimum number of sequences for mapping to the DfrB family, retaining only DfrA proteins a conserved position, 86; minimum number of sequences for analysis (Table S2). We further excluded redundant for a flanking position, 95; maximum number of contiguous DfrA sequences (amino acid sequence identity >90 %) using nonconserved positions, 5; minimum length of a block, t-­coffee version 11.00.8cbe486 seq_reformat command 4) [36]. Bayesian phylogenetic inference on the trimmed [27], and used the resulting non-­redundant panel to identify multiple amino acid sequence alignment was carried out with DfrA homologues in protein records associated with National MrBayes version 3.2.6 [37]. Four Metropolis-­coupled Markov 3 Sánchez-­Osuna et al., Microbial Genomics 2020;6 chain Monte Carlo simulations with four independent chains Sequence analysis were run for 20 000 000 generations, using a mixed four-­ To assess whether the identified chromosomal gene associated category gamma distributed rate plus proportion of invariable with a mobile dfrA gene is the canonical folA gene for the sites model [invgamma] and a JTT (Jones–Taylor–Thornton) genus, and not the product of a subsequent recombination amino acid substitution model [38]. Chains were sampled of the mobile dfrA gene into the chromosome, we computed every 100 iterations and stationarity was analysed with Tracer the pairwise amino acid identity among the products of all version 1.7.1 [39] by monitoring the estimated sample size chromosomal folA homologues and then compared this (ESS). To determine burn-­in, chain results were summarized distribution with the pairwise amino acid identity of the with MrBayes imposing the restriction that ESS be above 200 putative origin versus the chromosomal folA homologues. We and that the potential scale reduction factor be within 0.005 of used a one-­sided Mann–Whitney U test to determine whether 1. Based on summarization results, the burn-­in was set at 20 % the two distributions were significantly different. To analyse of iterations. A consensus tree was generated with the half-­ the mol% G+C content relationship between sul/dfrA genes compat option and visualized using the ggtree version 1.14.6 and their host chromosomes, we used pre-­compiled panels of R library [40]. Clades of reported DfrA proteins were deter- sequences for non-­redundant Sul [25] and DfrA homologues mined as branches with posterior probability values higher to search protein records associated with NCBI GenBank/ than 0.8 containing at least five reported DfrA sequences. RefSeq sequences of mobile genetic elements. The nucleo- Ancestral state reconstruction of a single binary trait (mobile/ tide sequences of the genes encoding these proteins was then chromosomal) was performed with BayesTraits version 3.0.2 retrieved. The mol% G+C content of the corresponding sul [41]. The mobile/chromosomal state of each sequence was and dfrA genes, as well as the overall mol% G+C content in determined through the data collection pipeline outlined both the mobile genetic element and the chromosome of the above. Known states at tree tips were labelled, and ancestral species harbouring it, were computed with custom Python states were reconstructed using the multistate and maximum-­ scripts. To analyse whether mobile dfrA genes with mol% likelihood settings. G+C content close to their hosts’ genomes are more similar to the hosts’ folA genes than expected if dfrA–host associations were arbitrary, we performed a permutation test comparing DNA techniques and in vitro trimethoprim- the mean pairwise amino acid sequence alignment distance susceptibility assay between DfrA proteins and host-­encoded FolA proteins. We With the exception of the Ralstonia solanacearum GMI1000 randomly permuted DfrA-­host assignments 1000 times and (Marc Valls, Center for Research in Agricultural Genomics, computed the corresponding P value as the rank of the non-­ Barcelona, Spain) and E. coli K-12 (CGSC5073) folA genes, permuted mean pairwise alignment distance. The input files which were amplified from genomic DNA, dfrA and folA and scripts used for data collection and analysis are available homologues were adapted to E. coli codon usage, synthesized in public repositories: (i) input files (json, txt and fasta) (ATG:biosynthetics) and then subcloned into a dephospho- and blast database for Python scripts used in data collection rylated pUA1108 vector [42] using an NdeI and BamHI (New and analysis – https://​doi.​org/​10.​6084/​m9.​figshare.​12156891.​ England Biolabs) double digest when possible. Genes with v1, and (ii) GitHub repository containing the Python scripts internal restriction sites for any of these two enzymes were used for data collection and analysis – https://​doi.​org/​10.​ subcloned into the same vector using the HiFi DNA assembly 5281/​zenodo.​3760352. kit (New England Biolabs), following the manufacturer’s protocol. Oligonucleotides used in this work are listed in Table S4. All constructs were validated by sequencing (Macrogen) RESULTS AND DISCUSSION prior to their use in transforming E. coli K-12 (CGSC5073). A large fraction of reported dfrA genes share a The minimum inhibitory concentration (MIC) for trimetho- common evolutionary origin prim (Sigma-­Aldrich) for strains of E. coli K-12 (CGSC5073) To explore the phylogenetic relationship of trimethoprim-­ carrying different versions of pUA1108 encoding folA or dfrA resistance determinants with their chromosomally encoded homologues was determined following the Clinical and Labo- folA counterparts, we used a non-­redundant panel of protein ratory Standards Institute guidelines using microdilution tests sequences encoded by reported dfr genes (Table S2) to detect in Mueller–Hinton (MH) broth (Merck) [43]. All MIC assays putative DHFR homologues in sequenced mobile elements. were performed in triplicate. Colonies were grown on Luria– We discarded sequences associated with the dfrB gene family, Bertani agar for 18 h and then suspended in sterile 0.9 % NaCl and retained for analysis non-­redundant sequences mapping solution to a McFarland 0.5 turbidity level. Suspensions were to the clades defined by dfrA genes reported in the Proteobac- then diluted at 10−2 in MH broth, and 50 µl (5×104) cells were teria and by dfrDEFGK genes associated with Gram-­positive inoculated into microtitre plates that contained 50 µl MH bacteria. For convenience, and in accordance with recent broth supplemented with 1024–0.250 mg trimethoprim l−1. reports on dfr nomenclature [46], we hereinafter designate To determine growth, optical density at 550 nm was measured these two groups (dfrA and dfrDEFGK) with the umbrella after 24 h incubation at 37 °C. The dfrA1 gene was used as a term dfrA. These reference mobile DfrA homologues were positive control [44] and the E. coli folA gene as a negative then combined with FolA homologues sampled from repre- control [45]. sentative genomes of all bacterial orders with complete 4 Sánchez-­Osuna et al., Microbial Genomics 2020;6 Fig. 1. Consensus tree of DHFR protein sequences. Branch support values are provided as Bayesian posterior probabilities estimated after four independent runs of 20 000 000 generations. Support values are only shown for branches with posterior probability values higher than 0.8. For chromosomal DHFR, the species name is displayed. Mobile DHFRs are denoted by their established dfrA name or by their NCBI GenBank accession number. Reported dfrA genes deemed redundant (>90 % identity) are listed next to the corresponding non-­redundant taxon included in the analysis. Next to each tip label, coloured boxes designate trimethoprim-­resistant (orange) and -sensitive (purple) DHFR. Numbers between brackets indicate the mol% G+C content of the sequence for the gene encoding the DHFR. Tip label colouring denotes previously reported (green) and novel (blue) DHFRs. Bold label text indicates that resistance has been experimentally assessed in this work. DHFR variants marked with an asterisk are encoded in megaplasmids (>400 kbp). The internal ring shows the mol% G+C of the gene encoding the DHFR in a yellow−red colour scale, while the external ring displays the ratio between the mol% G+C content of the genome harbouring the DHFR gene and the mol% G+C content of the gene. Dotted lines from the inner ring to tip labels denote genes discussed in the text. Reconstructed mobile/chromosomal states are displayed on ancestral nodes as pink/ black pie charts. genome assemblies, and of each bacterial family within the to two well-­supported (>0.8 posterior probability), distinct Proteobacteria (Table S5). The resulting multiple amino acid clades that likely arose from two different ancestors. The first sequence alignment was used to perform Bayesian phyloge- clade (clade 1), typified by the DfrA1 and DfrA12 proteins netic analysis of FolA/DfrA sequences. [47], includes 22 sequences encoded by previously reported The phylogenetic tree shown in Fig. 1 showcases the genetic dfrA genes with a mean amino acid identity of 51.19 % ±17.63 diversity of DfrA/FolA proteins, which encompass sequences sd, divided into two subgroups (containing 17 and 5 known with pairwise amino acid sequence identities ranging from dfrA genes, respectively) and associated with Gammapro- 99 to 20 % (Table S6). The resulting phylogeny also reveals teobacteria pathogens. This clade also includes a basal set of that the vast majority (70.7%) of known DfrA sequences map taxa composed of the Clostridioides difficile dfrF gene and 5 Sánchez-­Osuna et al., Microbial Genomics 2020;6 which also shows lower mol% G+C content (43.88 mol%±4.91 sd) than expected for the Enterobacteriaceae. To ascertain whether this pattern of mol% G+C heterogeneity extended to other previously reported and putative dfrA genes, we examined the mol% G+C content of dfrA (935 genes) and sul (408 genes) homologues identified in this analysis with respect to the genome mol% G+C content of the host species harbouring these mobile resistance genes. The results shown in Fig. 2(a) and Table S8 reveal that dfrA genes tend to align with host genome mol% G+C content (Pearson ⍴=0.56), whereas sul genes display a two-­tiered distri- bution of mol% G+C content that is essentially independent of host genome mol% G+C (Pearson ⍴=0.14). Available dfrA and sul sequences are dominated by variants of a known dfrA and sul genes that have been isolated predominantly in a select group of bacterial hosts (Fig. 2b). To correct for this skew, we filtered dfrA sequences based on the amino acid identity (<90 %) of their encoded proteins. This filtering resulted in a significantly smaller number of non-­redundant representa- tive dfrA (63 genes) and sul (4 genes) sequences (Table S9). The four representative sul genes correspond to one exemplar of the sul1 and sul2 families, and two exemplars of the sul3 Fig. 2. Correlation between the mol% G+C content of mobile dfrA (red family. Among representative dfrA genes, 14 map to the first circles) and sul (green squares) genes and that of their host genome. clade (clade 1) of Fig. 1 and 4 to the second clade (clade 2). Large open circles/squares denote representatives of clusters of redundant sequences (identity >90 %), and dfrA genes from clade 1 The correlation of dfrA genes with host genome mol% G+C and clade 2 in Fig. 1 are marked with an additional corona. A 0.75 % increases significantly (Pearson ⍴=0.78) when considering jitter to both x- and y-­axis values has been applied for visualization only non-­redundant representative dfrA sequences. The fact purposes. The red line shows the linear regression for representative that the mol% G+C content of representative dfrA sequences dfrA gene values. The Pearson R2 coefficient is superimposed. Vertical aligns well with their host genome mol% G+C could suggest background bars in (a) designate DfrA sequences harboured by that mol% G+C content in dfrA genes has been ameliorated mobile genetic elements (MGEs) identified in E. coli and K. pneumoniae to match the host’s. Alternatively, it could indicate that the isolates, which are heavily overrepresented in the dataset. Sequences from clusters with more than 100 sequences (represented by dfrA12, mobile dfrA gene originated via mobilization of a chromo- dfrA5 and dfrA1) are shown with specific markers, and highlighted somal folA gene from a bacterium in the same clade as the by horizontal background bars. The number of MGEs identified as current host. The later scenario posits that, besides mol% harbouring dfrA genes, before and after filtering DfrA sequence identity G+C content similarity, representative dfrA genes should (>90 %), is shown in (b). also encode proteins with significant sequence similarity to their hosts’ FolA protein. We performed a permutation test to analyse whether representative dfrA gene products show two newly identified mobile dfrA homologues also from significant similarity with their hosts’ FolA protein (Table Firmicutes isolates. The second clade (clade 2), exemplified S10). Our results indicate that this is the case (P<0.001), by DfrA18 [48], comprises a group of six highly diverged suggesting that most mobile dfrA genes are still associated (34.37%±10.15 sd mean amino acid identity) DfrA sequences with species from the same clade they originated in. from Gammaproteobacteria isolates. The filtering of dfrA sequences based on amino acid identity Analysis of the dfrA gene sequences in these two clades reveals brings forward three large clusters (represented by dfrA12, an unexpected degree of heterogeneity in mol% G+C content. dfrA5 and dfrA1, and belonging to clade 1 from Fig. 1) In the first clade, several dfrA homologues, including the containing more than 100 genes with amino acid identity C. difficile dfrF gene, show relatively low mol% G+C content larger than 90 %. The dfrA genes in these clusters show a distri- (Fig. 1, inner ring), matching the Firmicutes species they bution of mol %G+C content that is essentially independent were reported on (Table S7, Fig. 1, outer ring). Similarly, dfrA of the host genome mol% G+C, as in the case of sul genes genes in the dfrA12 group show a mol% G+C content (53.28 (Fig. 2a), and their products show no significant sequence mol%±1.80 sd) that is well in line with that of the Enterobac- similarity with the hosts’ FolA (permutation test P>0.1). teriaceae isolates harbouring them. Conversely, the largest This indicates that the dfrA genes in these large clusters have group in this clade, encompassing dfrA1, dfrA7 and dfrA14, spread across distantly related bacterial clades, primarily shows a mean mol% G+C content of 41.03 mol%±3.99 sd, through their association with sul-­containing integron-­based which is substantially lower than the average mol% G+C transposable elements that are widely disseminated among content of the Enterobacteriaceae harbouring these mobile clinically relevant bacteria [49, 50]. This analysis also brings to elements. The same holds true for the second clade (dfrA18), the fore the presence of multiple dfrA cluster representatives, 6 Sánchez-­Osuna et al., Microbial Genomics 2020;6 with widely divergent mol% G+C, on narrow bands of host homologue was a bona fide FolA protein and not a mobile genome mol% G+C content. These bands correspond to DHFR homologue, and we excluded from further analysis E. coli (50.7 mol% G+C) and Klebsiella pneumoniae (57.4 all other DHFR homologues identified in megaplasmids mol% G+C) isolates, which are heavily oversampled in the (>400 kbp). dataset (Fig. 2b). The marked divergence in mol% G+C Two remaining clades of novel mobile DHFR homologues content (E. coli, 9.55 mol %±6.15 sd; K. pneumoniae, 6.93 from clinically relevant bacteria associated with chromo- mol %±4.63 sd) and amino acid sequence identity (E. coli, somal FolA proteins were shown to confer resistance to 38.57 % ±15.24 sd; K. pneumoniae, 38.39 % ±14.02 sd; Table trimethoprim on E. coli (Table 1). To investigate whether S11, Fig. S2) among these representative dfrA genes suggests the sequence determinants conferring resistance had that they originated via mobilization from a diverse set of originated in the associated chromosomal background, we chromosomal backgrounds. cloned the most closely related chromosomal folA gene as The fact that the mol% G+C of dfrA genes aligns with their well as a gene encoding an additional DHFR homologue host genome’s mol% G+C, and that DfrA proteins display from the same genus, and performed broth microdilution higher amino acid sequence identity when aligned to their assays to determine the MIC of trimethoprim. We also host genome FolA proteins than to other FolA proteins, performed ancestral state reconstruction of the molecule strongly supports the notion that dfrA genes have been encoding the DHFR homologues (chromosomal/mobile mobilized multiple times within different bacterial clades trait), as determined during the data collection process [13]. In a few instances, typified by the large dfrA clusters (Tables S5 and S12). illustrated in Fig. 2, dfrA genes have been captured by highly The combined results of Table 1 and Fig. 1 reveal different efficient mobile elements and dispersed widely across unre- patterns of trimethoprim-­ resistance acquisition. The lated groups of bacteria [49]. These mobile elements often protein KMV08986 is a DHFR homologue harboured by harbour sul genes, which also display a host-­independent a conjugative plasmid from an Acinetobacter baumannii mol% G+C distribution. Many of the dfrA genes identified clinical isolate. Its most closely related chromosomally here are associated with clinical isolates. The divergent mol% encoded DHFR homologue is the FolA protein of Flavo- G+C content and amino acid identity of these dfrA genes indi- bacterium branchiophilum, which confers resistance to cates that pathogenic bacteria have obtained dfrA genes on trimethoprim (Table 1). To ascertain whether this chro- multiple occasions and from different sources, highlighting mosomally encoded DHFR homologue was encoded by a the ability of mobilized resistance determinants to reach clini- bona fide folA gene, instead of a mobile dfrA gene that inte- cally relevant pathogens [17, 51]. grated into the chromosome, we compared the genus-­wide distribution of pairwise amino acid sequence alignment Novel trimethoprim-resistance determinants of distances between FolA proteins to the pairwise distance Acinetobacter clinical isolates identified through of the identified homologue versus all other FolA proteins phylogenetic methods in the genus. The Flavobacterium branchiophilum FolA The phylogenetic tree in Fig. 1 includes reported DfrA sequence is significantly different from other Flavobac- proteins and their putative homologues, as well as FolA terium FolA sequences (Mann–Whitney U P<0.05; Table proteins identified via tblastn as putative DfrA homo- S13), raising the possibility that this chromosomal gene logues or sampled uniformly across bacterial clades. The could be in fact a recombined mobile dfrA gene. However, inferred phylogeny also reveals several groups of previ- phylogenetic analysis with a broader representation of ously unreported mobile DHFR homologues that form Flavobacterium sequences (Fig. S3) confirms the well-­ well-­supported clades in association with chromosomal supported branching of Flavobacterium branchiophilum FolA proteins. Hence, these FolA proteins could constitute FolA with other Flavobacterium species FolA proteins, the chromosomal origins of the associated mobile DHFR and comparative genomics analysis reveals that the genetic homologues, and provide insights into the emergence neighbourhood of the chromosomal folA gene is conserved and dissemination of trimethoprim-­resistance genes. To in the genus Flavobacterium (Fig. S4). Furthermore, the determine whether these mobile DHFR homologues did FolA protein of a prototypical genus member, Flavobacte- confer resistance to trimethoprim, we cloned a subset rium faecale, also confers resistance to trimethoprim on E. of dfrA/folA genes and performed broth microdilution coli (Table 1). These results indicate that the FolA protein assays to determine the MIC of trimethoprim. Consid- was likely resistant to trimethoprim in the ancestor of ering that the clinical breakpoint for trimethoprim in extant Flavobacterium species, which diverged more than E. coli is 4 mg l−1 [52], the results, shown in Table 1, reveal 50 million years ago [53]. The branching of the Acineto- that most of the mobile DHFR homologues identified bacter baumannii protein KMV08986 in the reconstructed here do confer significant resistance to trimethoprim. The phylogeny and the associated ancestral state reconstruction sole exception is the protein AQW32254. Close inspec- indicates that this mobile DHFR homologue likely origi- tion revealed that this DHFR homologue is encoded by nated via mobilization of a chromosomal folA gene within a megaplasmid (1.2 Mb) from a Ralstonia isolate, and the phylum Bacteroidetes. The encoded FolA protein that this is the only DHFR homologue present in its was likely resistant to trimethoprim, but the exact donor complete genome. Hence, we determined that this DHFR species remains to be elucidated. 7 Sánchez-­Osuna et al., Microbial Genomics 2020;6 Table 1. MICs of trimethoprim for wild-­type E. coli K-12 (CGSC5073) and derivatives carrying different versions of dfr/folA or the control empty vector Values are representative of four independent replicates. Strain Mobile / Nucleotide accession Cloned protein ID Trimethoprim (mg l−1) chromosomal no. E. coli CGSC5073 – – – 0.25 E. coli pUA1108 – – – 0.25 E. coli pUA1108::folA E. coli C NC_000913 WP_000624375 4 E. coli pUA1108::dfrA1 M NC_002525 WP_000777554 >512 E. coli pUA1108::folA Flavobacterium branchiophilum C NC_016001 WP_014083133 256 E. coli pUA1108::folA Flavobacterium faecale C NZ_CP020918 WP_108740183 >512 E. coli pUA1108::dfrA38 Acinetobacter baumannii M CP021344 KMV08986 256 E. coli pUA1108::folA Acinetobacter schindleri C NZ_CP025618 WP_004813248 0.25 E. coli pUA1108::dfrA39 Acinetobacter baumannii M NZ_CP021785 WP_031380727 512 E. coli pUA1108::dfrA40 Acinetobacter baumannii M NZ_JEVW01000010 WP_034702334 128 E. coli pUA1108::dfrA41 Acinetobacter defluvii M NZ_CP029396 WP_004729503 >512 E. coli pUA1108::folA Fluviicola taffensis C NC_015321 WP_013685591 >512 E. coli pUA1108::folA 'Candidatus Fluviicola riflensis' C CP022585 ASS49886 >512 E. coli pUA1108::folA Alcanivorax pacificus C NZ_CP004387 WP_008736147 32 E. coli pUA1108::folA Alcanivorax borkumensis C AM286690 CAL17791 16 E. coli pUA1108::folA Bacillus mobilis C NZ_CP031443 WP_000637217 >512 E. coli pUA1108::folA Ralstonia solanacearum C NC_003295 WP_011000898 0.5 E. coli pUA1108::folA blood disease bacterium A2-­HR MARDI M CP019912 AQW32254 1 E. coli pUA1108::folA E. coli O104:H4 M CP003298 AFS59762 2 In contrast to Flavobacterium proteins, the Acinetobacter multiple antibiotic-­resistance determinants (Fig. S5). The schindleri FolA protein does not confer resistance to gene encoding this DHFR homologue is preceded by an trimethoprim, in agreement with previous reports of insertion sequence transposase (Fig. S5), in an arrangement Acinetobacter schindleri susceptibility to trimethoprim that has been reported to drive up expression of the DHFR [54], and with the well-­established susceptibility of Acine- homologue through promoter enhancement [57]. However, tobacter baumannii FolA to trimethoprim [55, 56]. The the MIC determined here corresponds to that of the isolated Acinetobacter schindleri FolA protein is closely related to DHFR ORF, indicating that it confers heightened resistance three mobile DHFR homologues conferring resistance to irrespective of the promoter driving its expression. This trimethoprim and harboured by Acinetobacter baumannii suggests that these DHFR homologues have acquired muta- (protein ID: WP_031380727, WP_034702334) and Acine- tions conferring heightened resistance to trimethoprim in tobacter defluvii (protein ID: WP_004729503) clinical and parallel to their broader dissemination on multi-­resistant environmental isolates. These mobile DHFR homologues mobile elements. Based on their validated trimethoprim-­ branch within a well-­supported clade of chromosomal resistance phenotype and their level of amino acid sequence Acinetobacter FolA proteins, as supported by ancestral identity versus previously reported DfrA proteins (<95 %; state reconstruction (Fig. 1, Table S12). The trimethoprim Table S14) [13], we propose to designate these Acinetobacter susceptibility of Acinetobacter chromosomal folA genes and DHFR homologues as DfrA38 (protein ID: KMV08986), the phylogenetic placement of these DHFR homologues, DfrA39 (protein ID: WP_031380727), DfrA40 (protein ID: hence, indicates that the observed resistance to trimetho- WP_034702334) and DfrA41 (protein ID: WP_004729503). prim was acquired immediately prior to or after mobiliza- tion from an Acinetobacter chromosomal background. This Here, we report the identification of trimethoprim-­ is supported by the observation that these mobile DHFR susceptible chromosomal folA genes that are closely related homologues confer different levels of resistance to trimeth- to mobile dfrA genes, as well as the discovery of chromosom- oprim (Table 1), and that the largest MIC correlates with the ally encoded folA genes conferring resistance to trimetho- location of the DHFR homologue on a plasmid harbouring prim. This indicates that, in contrast to sulfonamides [25], 8 Sánchez-­Osuna et al., Microbial Genomics 2020;6 trimethoprim-­resistance mutations with small or negligible FolA proteins confer resistance to trimethoprim (Table 1). fitness cost must occur frequently enough in natural envi- The most closely associated chromosomal folA genes are ronments. These folA variants can then be selected for not significantly different from other folA genes in their and mobilized upon exposure to trimethoprim. It is well-­ respective genera (Mann–Whitney U P>0.05; Table S13), documented that resistance to trimethoprim, mediated by as reflected also by substantial conservation of the folA mutations in the chromosomal folA gene, develops very genomic neighbourhood (Fig. S4). Together, these data rapidly and in a fairly structured way [58–60], whereas indicate that resistance to trimethoprim was present in resistance to sulfonamides takes much longer to evolve the ancestor of these genera. The dfrA26 gene was iden- in a laboratory setting. Moreover, sulfonamide-­resistant tified in a K. pneumoniae clinical isolate and its most mutants typically show significantly reduced affinity to closely associated chromosomal folA gene is a member of PABA. This results in a net fitness cost in the absence of the genus Alcanivorax. The branching pattern of dfrA26 sulfonamide that is only palliated by the emergence of within this clade and ancestral state reconstruction results subsequent compensatory mutations [61, 62]. Beyond (Fig. 1, Table S12) suggest that it arose via mobilization of structural constraints on the respective binding pockets, a a chromosomal folA gene from the genus Alcalinivorax. The crucial difference between both chemotherapeutic agents dfrDGK genes have been reported in Enterococcus faecalis, lies in their respective targets. While trimethoprim directly Enterococcus faecium and Staphylococcus aureus, and ances- inhibits DHFR, sulfonamides compete with PABA for tral state reconstruction results indicate that these mobile access to DHPS, yielding a non-­productive sulfonamide-­ dfrA genes originated through mobilization of a member of bound di-­hydropterin. For sulfonamides, therefore, it is closely related genus Bacillus, members of which have been the PABA-­to-­sulfonamide ratio that limits the produc- reported to be naturally resistant to trimethoprim [66]. In tion of di-­hydropteroate from a limited pool of pteridine both cases, therefore, the phylogenetic evidence and the di-­phosphate, and this cannot be altered via overexpres- similarity in mol% G+C content among chromosomal and sion of DHPS [63]. Conversely, DHFR overexpression can mobile genes (Fig. 1, Table S15) point towards a mobili- provide partial resistance to trimethoprim, and mutations zation event that has to date remained circumscribed to enhancing DHFR expression have been reported to be the related genera. Conversely, the dfrA20 gene was identified first to appear in directed evolution experiments [60]. The in a Pasteurella multocida isolate, yet the chromosomal folA ability to obtain partial resistance through overexpression gene most closely associated to it is encoded by Fluviicola may provide a stepping stone for the gradual accumula- taffensis, a Bacteroidetes; hence, suggesting a much more tion and refinement of mutations conferring substantial distant mobilization event (Fig. 1, Table S15). In all three resistance with little fitness cost and, hence, facilitate the cases, however, we find evidence that pre-­existing resistant development of trimethoprim resistance [59, 60]. folA genes can be readily mobilized from both close (e.g. dfrDGK) or distant (e.g. dfrA20) species. Trimethoprim resistance in chromosomally The resistance to trimethoprim reported here for the chro- encoded folA genes mosomal folA genes of two different genera of Bacteroi- Besides uncovering novel dfrA genes, the phylogenetic detes, two distinct Alcanivorax species and a Bacillus strain analysis in Fig. 1 also identifies several chromosomal folA underscores the deep ancestry of chromosomal mutations genes associated with previously reported dfrA genes. Two yielding resistance to trimethoprim. The folA genes of of these chromosomal folA genes have already been reported Flavobacterium and Fluviicola were shown here to confer in the literature as putative origins of dfrA genes, and their resistance to trimethoprim. These two genera are thought to identification here provides some degree of validation for have diverged more than 500 million years ago and define the phylogenetic approach implemented in this work. The major lineages within the Flavobacteriales, suggesting putative chromosomal origin for Staphylococcus aureus that resistance to trimethoprim emerged in an ancestor of Tn4003 S1-­DHFR has been identified as the chromosom- this bacterial order. It is worth noting that several of the ally encoded dfrC gene (Staphylococcus epidermidis) and is chromosomal folA genes shown here to be associated with reported to be susceptible to trimethoprim [64]. The Ente- mobile DHFR homologues (Alcanivorax, Flavobacterium rococcus faecalis dfrE gene, identical to the chromosomally and Fluviicola) appear to be resistant at the genus level and encoded folA gene of Enterococcus faecalis, was reported correspond to genera of aquatic bacteria. This parallels our to confer moderate resistance to trimethoprim in E. coli, recent identification of soil and subterranean water bacteria but only when cloned in a multicopy plasmid, which could as the likely originators of clinical sulfonamide-­resistance easily result in overexpression-­mediated resistance [63, 65]. genes [25], and suggests that the intensive use of trimetho- To ascertain whether the chromosomal folA genes found prim/sulfamethoxazole in agriculture, aquaculture and here to be associated with other known dfrA genes (dfrA20, animal husbandry in the last 50 years may have acted as dfrA26 and the dfrDGK cluster) confer resistance to a trigger for the selection and mobilization of pre-­existing trimethoprim, we performed broth microdilution assays folA and folP genes conferring resistance to trimethoprim to determine the MIC of trimethoprim on these chro- and sulfonamides. Conversely, trimethoprim-­susceptible mosomally encoded FolA proteins and on another FolA chromosomal folA genes found here to be associated protein from the same genus. In all cases, both related with dfrA genes belong to clinically relevant genera 9 Sánchez-­Osuna et al., Microbial Genomics 2020;6 Fig. 3. Schematic representation of the two proposed evolutionary processes (based on the results presented in Figs 1 and 2, and Table 1) leading to the dissemination of trimethoprim-­resistance determinants. Left panel: upon the introduction of trimethoprim, mobilization events involving pre-­existing resistant chromosomal folA genes can be favourably selected. Right panel: following the introduction of trimethoprim, mobilization events involving folA genes with novel mutations that confer resistance to this chemotherapeutic agent may be selected for and disseminated among closely related bacteria. (Staphylococcus and Acinetobacter) that may have been of metagenomics [72, 73]. However, recent studies under more direct trimethoprim pressure. This suggests that have shown that many potential resistant determinants among relatively isolated bacterial populations, frequent encoded by phages do not confer resistance against their exposure to high levels of trimethoprim may trigger the putative targets. Furthermore, only a small proportion mobilization of spontaneous folA mutants, whereas longer of complete phage genomes contain putative antibiotic-­ term exposure to sub-­lethal doses of trimethoprim in resistance genes [74]. Enzymes participating in the folate ecological rich habitats might instead rely predominantly biosynthesis pathway, however, are relatively frequent in on the mobilization of naturally resistant folA genes (Fig. 3). phage genomes. These include homologues of the folP gene encoding DHPS, of the thyX gene encoding flavin-­ Phage-encoded folA genes do not confer resistance dependent thymidylate synthase [75–77] and, predomi- to trimethoprim nantly, homologues of the folA gene encoding DHFR, Our phylogenetic analysis also identifies a well-­defined often found in tandem with the thyA gene encoding type clade of Enterobacteriaceae cryptic plasmids derived from 1 thymidylate synthase [71]. Salmonella phage SSU5 and encoding DHFR homologues [67–70]. Genes encoding DHFR homologues occur Early work on Enterobacteria phage T4 showed that the frequently in many bacteriophage families, often in tandem phage-­encoded thyA and folA gene products are functional with thymidylate synthase genes [71], but their functional and also participate in the phage baseplate structure [78], role has not been fully elucidated. We performed broth and thyX has been shown to be functional in a number of microdilution assays to determine the MIC of trimetho- phages [75–77]. It has been proposed that these genes help prim of E. coli O104:H4 DHFR (protein ID: AFS59762). bacteriophages overcome shortages in the deoxynucleotide This phage-­encoded DHFR does not confer resistance to pool during replication, but their potential in conferring trimethoprim (Table 1). The high amino acid sequence resistance to sulfonamides or trimethoprim remains largely identity and neighbourhood conservation among the unexplored. The detection here of DHFR homologues in DHFR enzymes encoded by these Enterobacteriaceae Enterobacteriaceae cryptic plasmids and phages, and the cryptic plasmids and phages (Table S16, Fig. S4) would subsequent assessment of their trimethoprim susceptibility, presumably suggest that these DHFR enzymes are suscep- reinforces the notion that these genes have been function- tible to trimethoprim. ally co-­opted by phages principally for deoxynucleotide Bacteriophages can transfer substantial amounts of synthesis. Nonetheless, these genes may still confer partial genetic material via generalized transduction, and their trimethoprim resistance as a by-­product of folA overexpres- potential as reservoirs of antibiotic-­resistance determi- sion, as recently reported for Stenotrophomonas maltophilia nants has gained increased attention with the advent phage DLP4 [79]. 10 Sánchez-­Osuna et al., Microbial Genomics 2020;6 Conclusions 4. Davies J, Davies D. Origins and evolution of antibiotic resistance. Microbiol Mol Biol Rev 2010;74:417–433. Recent work has shown that resistance to sulfonamide, 5. Baym M, Lieberman TD, Kelsic ED, Chait R, Gross R et al. Spati- a synthetic chemotherapeutic agent, can be present in otemporal microbial evolution on antibiotic landscapes. Science the bacterial pangenome well before the discovery of the 2016;353:1147–1151. agent. Here, we have used a combination of in silico and in 6. Hegreness M, Shoresh N, Damian D, Hartl D, Kishony R. Acceler- vitro techniques to identify novel trimethoprim-­resistance ated evolution of resistance in multidrug environments. Proc Natl Acad Sci USA 2008;105:13977–13981. genes, and to identify chromosomal folA genes that are 7. Aminov RI, Mackie RI. Evolution and ecology of antibiotic resist- strongly associated with novel and previously reported ance genes. FEMS Microbiol Lett 2007;271:147–161. dfrA genes. We find that most of the chromosomal folA 8. Sengupta S, Chattopadhyay MK, Grossart H-­P. The multifaceted genes associated with mobile dfrA genes confer resistance roles of antibiotics and antibiotic resistance in nature. Front Micro- to trimethoprim, but we detect cases of novel mutations biol 2013;4:47. being rapidly mobilized. Hence, our work shows that 9. Kmeid JG, Youssef MM, Kanafani ZA, Kanj SS. Combination therapy the observations from sulfonamide resistance extend to for Gram-­negative bacteria: what is the evidence? Expert Rev Anti Infect Ther 2013;11:1355–1362. trimethoprim, with generalized chromosomal resistance 10. Williams KJ. The introduction of ‘chemotherapy’ using arsphena- determinants predating the origin of several genera and mine – the first magic bullet. J R Soc Med 2009;102:343–348. several clusters of resistance genes disseminated broadly 11. Aminov RI. A brief history of the antibiotic era: lessons learned and among clinical isolates. Moreover, this work also reveals challenges for the future. Front Microbiol 2010;1:134. that, unlike sulfonamides, resistance to trimethoprim is 12. Masters PA, O'Bryan TA, Zurlo J, Miller DQ, Joshi N. Trimethoprim-­ relatively easy to generate and frequently associated with sulfamethoxazole revisited. Arch Intern Med 2003;163:402–410. species from the same clade it originated in. The identifica- 13. Sköld O. Resistance to trimethoprim and sulfonamides. Vet Res tion of ancient resistance determinants for two synthetic 2001;32:261–273. chemotherapeutic agents strongly suggests that resistance 14. Quinlivan EP, McPartlin J, Weir DG, Scott J. Mechanism of the antimi- crobial drug trimethoprim revisited. FASEB J 2000;14:2519–2524. to any novel drugs is likely to be already present in the 15. Hitchings GH. Mechanism of action of trimethoprim-­ bacterial pangenome. Systematic screening of existing sulfamethoxazole-­I. J Infect Dis 1973;128:S433–S436. natural variants could provide, therefore, the means to pre-­ 16. Landy M, Larkum NW, Oswald EJ, Streightoff F. Increased emptively identify derivatives presenting widely distributed synthesis of p-­aminobenzoic acid associated with the develop- natural resistance determinants and, conversely, to engineer ment of sulfonamide resistance in Staphylococcus aureus. Science derivatives that circumvent most, if not all, natural resistant 1943;97:265–267. variants. 17. Huovinen P, Sundström L, Swedberg G, Sköld O. Trimetho- prim and sulfonamide resistance. Antimicrob Agents Chemother 1995;39:279–289. 18. Flensburg J, Sköld O. Massive overproduction of dihydrofolate Funding information reductase in bacteria as a response to the use of trimethoprim. This work was supported by grant BIO2016-77011-­R from the Spanish Eur J Biochem 1987;162:473–476. Ministerio de Economia y Competitividad to J.B. M.S.-O. was the recip- ient of a predoctoral fellowship from the Ministerio de Educación, 19. Shin HW, Lim J, Kim S, Kim J, Kwon GC et al. Characterization of Cultura y Deporte de España. trimethoprim-­sulfamethoxazole resistance genes and their relat- edness to class 1 integron and insertion sequence common region Acknowledgements in Gram-­negative bacilli. J Microbiol Biotechnol 2015;25:137–142. The authors wish to thank Joan Ruiz and Susana Escribano for their 20. Rådström P, Swedberg G, Sköld O. Genetic analyses of sulfona- technical support during some of the experimental procedures, as well mide resistance and its dissemination in gram-­negative bacteria as Ángela Martínez-­Mateos for her continued support. The authors also illustrate new aspects of R plasmid evolution. Antimicrob Agents express their gratitude to Dr Marc Valls for kindly providing the R. sola- Chemother 1991;35:1840–1848. nacearum GMI1000 strain. 21. Perreten V, Boerlin P. A new sulfonamide resistance gene (sul3) in Author contributions Escherichia coli is widespread in the pig population of Switzerland. Conceptualization, J. B., I. E.; data curation, M.S.-O., I. E.; formal analysis, Antimicrob Agents Chemother 2003;47:1169–1172. M.S.-O., I. E.; funding acquisition, M. L., J. B.; investigation, M.S.-O., P. C., I. 22. Tagg KA, Watkins LF, Moore MD, Bennett C, Joung YJ et al. Novel E.; methodology, M. L., J. B., I. E.; project administration, P. C., M. L., J. B., trimethoprim resistance gene dfrA34 identified in Salmonella I. E.; resources, P. C., M. L., J. B.; software, M.S.-O., I. E.; supervision, P. C., Heidelberg in the USA. J Antimicrob Chemother 2019;74:38–41. J. B., I. E.; visualization, M.S.-O., I. E.; manuscript preparation – original 23. White PA, Rawlinson WD. Current status of the aadA and dfr gene draft, M.S.-O., I. E.; manuscript preparation – review and editing, M.S.-O., cassette families. J Antimicrob Chemother 2001;47:495–496. P. C., M. L., J. B., I. E. 24. Toulouse JL, Edens TJ, Alejaldre L, Manges AR, Pelletier JN. Conflicts of interest Integron-­associated DfrB4, a previously uncharacterized member The authors declare that there are no conflicts of interest. of the trimethoprim-­resistant dihydrofolate reductase B family, is a clinically identified emergent source of antibiotic resistance. Anti- References microb Agents Chemother 2017;61:e02665-16. 1. Letunic I, Bork P. Interactive tree of life (iTOL) v4: recent updates 25. Sánchez-­Osuna M, Cortés P, Barbé J, Erill I. Origin of the mobile and new developments. Nucleic Acids Res 2019;47:W256–W259. di-­hydro-­pteroate synthase gene determining sulfonamide resist- 2. Carlet J, Rambaud C, Pulcini C. Save antibiotics: a call for action ance in clinical isolates. Front Microbiol 2018;9:3332. of the World Alliance Against Antibiotic Resistance (WAAAR). BMC 26. Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol Infect Dis 2014;14:436. 2011;7:e1002195. 3. Rossolini GM, Arena F, Pecile P, Pollini S. Update on the antibiotic 27. Taly JF, Magis C, Bussotti G, Chang JM, Di Tommaso P et al. Using resistance crisis. Curr Opin Pharmacol 2014;18:56–60. the T-­Coffee package to build multiple sequence alignments 11 Sánchez-­Osuna et al., Microbial Genomics 2020;6 of protein, RNA, DNA sequences and 3D structures. Nat Protoc 50. Ho PL, Wong RC, Chow KH, Que TL. Distribution of integron-­ 2011;6:1669–1682. associated trimethoprim-­ sulfamethoxazole resistance determi- 28. O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D et al. Refer- nants among Escherichia coli from humans and food-­producing ence sequence (RefSeq) database at NCBI: current status, taxo- animals. Lett Appl Microbiol 2009;49:627–634. nomic expansion, and functional annotation. Nucleic Acids Res 51. Volz C, Ramoni J, Beisken S, Galata V, Keller A et al. Clinical 2016;44:D733–D745. resistome screening of 1,110 Escherichia coli isolates efficiently 29. Benson DA, Cavanaugh M, Clark K, Karsch-­Mizrachi I, Lipman DJ recovers diagnostically relevant antibiotic resistance biomarkers et al. GenBank. Nucleic Acids Res 2017;45:D37–D42. and potential novel resistance mechanisms. Front Microbiol 2019;10:1671. 30. Holm L, Sander C. Removing near-­neighbour redundancy from large protein sequence collections. Bioinformatics 1998;14:423–429. 52. MacGowan AP, Wise R. Establishing MIC breakpoints and the inter- pretation of in vitro susceptibility tests. J Antimicrob Chemother 31. Altschul Set al. Gapped BLAST and PSI-­BLAST: a new genera- 2001;48:17–28. tion of protein database search programs. Nucleic Acids Res 53. Kumar S, Stecher G, Suleski M, Hedges SB. TimeTree: a resource 1997;25:3389–3402. for timelines, timetrees, and divergence times. Mol Biol Evol 32. Kaushik S, Mutt E, Chellappan A, Sankaran S, Srinivasan N et al. 2017;34:1812–1819. Improved detection of remote homologues using cascade PSI-­ 54. Sigala J-­C, Suárez BP, Lara AR, Borgne SL, Bustos P et al. Genomic BLAST: influence of neighbouring protein families on sequence and physiological characterization of a laboratory-­isolated Acine- coverage. PLoS One;8:e56449. tobacter schindleri ACE strain that quickly and efficiently catabo- 33. Brum JR, Ignacio-­Espinoza JC, Roux S, Doulcier G, Acinas SG lizes acetate. Microbiology 2017;163:1052–1064. et al. Patterns and ecological drivers of ocean viral communities. 55. Falagas ME, Vardakas KZ, Roussos NS. Trimethoprim/sulfameth- Science 2015;348:1261498. oxazole for Acinetobacter spp.: a review of current microbiological 34. Chen Q, Zobel J, Verspoor K. Benchmarks for measurement of and clinical evidence. Int J Antimicrob Agents 2015;46:231–241. duplicate detection methods in nucleotide databases. Database 56. Pérez-­Varela M, Corral J, Aranda J, Barbé J. Roles of efflux pumps 2017;2017:baw164. from different superfamilies in the surface-­ associated motility 35. Wallace IM, O'Sullivan O, Higgins DG, Notredame C. M-­Coffee: and virulence of Acinetobacter baumannii ATCC 17978. Antimicrob combining multiple sequence alignment methods with T-­Coffee. Agents Chemother 2019;63:e02190-18. Nucleic Acids Res 2006;34:1692–1699. 57. Leelaporn A, Firth N, Byrne ME, Roper E, Skurray RA. Possible role 36. Castresana J. Selection of conserved blocks from multiple of insertion sequence IS257 in dissemination and expression of alignments for their use in phylogenetic analysis. Mol Biol Evol high- and low-­level trimethoprim resistance in staphylococci. Anti- 2000;17:540–552. microb Agents Chemother 1994;38:2238–2244. 37. Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic 58. Vickers AA, Potter NJ, Fishwick CWG, Chopra I, O'Neill AJ. Analysis inference under mixed models. Bioinformatics 2003;19:1572–1574. of mutational resistance to trimethoprim in Staphylococcus aureus 38. Erill I. Dispersal and regulation of an adaptive mutagenesis by genetic and structural modelling techniques. J Antimicrob cassette in the bacteria domain. Nucleic Acids Res 2006;34:66–77. Chemother 2009;63:1112–1117. 39. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior 59. Watson M, Liu J-­W, Ollis D. Directed evolution of trimethoprim summarization in Bayesian phylogenetics using Tracer 1.7. Syst resistance in Escherichia coli. FEBS J 2007;274:2661–2671. Biol 2018;67:901–904. 60. Toprak E, Veres A, Michel J-­B, Chait R, Hartl DL et al. Evolutionary 40. Yu G, Smith DK, Zhu H, Guan Y, LamTT‐Y. ggtree: an r package for paths to antibiotic resistance under dynamically sustained drug visualization and annotation of phylogenetic trees with their covar- selection. Nat Genet 2012;44:101–105. iates and other associated data. Methods Ecol Evol 2017;8:28–36. 61. Swedberg G, Fermér C, Sköld O. Point mutations in the dihydrop- 41. Pagel M, Meade A, Barker D. Bayesian estimation of ancestral teroate synthase gene causing sulfonamide resistance. In: Ayling character states on phylogenies. Syst Biol 2004;53:673–684. JE, Nair MG, Baugh CM (eds). Chemistry and Biology of Pteridines and Folates. Boston, MA: Springer; 1993. pp. 555–558. 42. Mayola A, Irazoki O, Martínez IA, Petrov D, Menolascina F et al. RecA protein plays a role in the chemotactic response and 62. Griffith EC, Wallace MJ, Wu Y, Kumar G, Gajewski S et al. The chemoreceptor clustering of Salmonella enterica. PLoS One structural and functional basis for recurring sulfa drug resistance 2014;9:e105578. mutations in Staphylococcus aureus dihydropteroate synthase. Front Microbiol 2018;9:1369. 43. Clinical and Laboratory Standards Institute. Methods for Dilution 63. Palmer AC, Kishony R. Opposing effects of target overexpression Antimicrobial Susceptibility Tests for Bacteria that Grow Aerobically, reveal drug mechanisms. Nat Commun 2014;5:4296. Approved Standard, 6th edn. Wayne, PA: Clinical and Laboratory Standards Institute; 2003. 64. Dale GE, Broger C, Hartman PG, Langen H, Page MG et al. Charac- terization of the gene for the chromosomal dihydrofolate reduc- 44. Fling ME, Richards C. The nucleotide sequence of the trimethoprim-­ tase (DHFR) of Staphylococcus epidermidis ATCC 14990: the origin resistant dihydrofolate reductase gene harbored by Tn7. Nucleic of the trimethoprim-­ resistant S1 DHFR from Staphylococcus Acids Res 1983;11:5147–5158. aureus? J Bacteriol 1995;177:2965–2970. 45. Smith DR, Calvo JM. Nucleotide sequence of the E. coli gene coding 65. Coque TM, Singh KV, Weinstock GM, Murray BE. Characterization of for dihydrofolate reductase. Nucleic Acids Res 1980;8:2255–2274. dihydrofolate reductase genes from trimethoprim-­susceptible and 46. Faltyn M, Alcock B, McArthur A. Evolution and nomenclature of the trimethoprim-­resistant strains of Enterococcus faecalis. Antimicrob trimethoprim resistant dihydrofolate (dfr) reductases. Preprints Agents Chemother 1999;43:141–147. 2019:2019050137. 66. Barrow EW, Bourne PC, Barrow WW. Functional cloning of 47. van Hoek AHAM, Mevius D, Guerra B, Mullany P, Roberts AP et al. Bacillus anthracis dihydrofolate reductase and confirmation of Acquired antibiotic resistance genes: an overview. Front Microbiol natural resistance to trimethoprim. Antimicrob Agents Chemother 2011;2:203. 2004;48:4643–4649. 48. Villa L, Visca P, Tosini F, Pezzella C, Carattoli A. Composite inte- 67. Parkhill J, Dougan G, James KD, Thomson NR, Pickard D et al. gron array generated by insertion of an ORF341-­type integron Complete genome sequence of a multiple drug resistant Salmo- within a Tn21-­like element. Microb Drug Resist 2002;8:1–8. nella enterica serovar Typhi CT18. Nature 2001;413:848–852. 49. Grape M, Farra A, Kronvall G, Sundström L. Integrons and gene 68. Ahmed SA, Awosika J, Baldwin C, Bishop-­Lilly KA, Biswas B cassettes in clinical isolates of co-­trimoxazole-­resistant Gram-­ et al. Genomic comparison of Escherichia coli O104:H4 isolates negative bacteria. Clin Microbiol Infect 2005;11:185–192. 12 Sánchez-­Osuna et al., Microbial Genomics 2020;6 from 2009 and 2011 reveals plasmid, and prophage hetero- 74. Enault F, Briet A, Bouteille L, Roux S, Sullivan MB et al. Phages geneity, including shiga toxin encoding phage stx2. PLoS One rarely encode antibiotic resistance genes: a cautionary tale for 2012;7:e48228. virome analyses. ISME J 2017;11:237–247. 69. Kim M, Kim S, Ryu S. Complete genome sequence of bacteriophage 75. Bhattacharya B, Giri N, Mitra M, Gupta SKD. Cloning, charac- SSU5 specific for Salmonella enterica serovar typhimurium rough terization and expression analysis of nucleotide metabolism-­ related genes of mycobacteriophage L5. FEMS Microbiol Lett strains. J Virol 2012;86:10894. 2008;280:64–72. 70. Octavia S, Sara J, Lan R. Characterization of a large novel phage-­ 76. Wittmann J, Gartemann K-­H, Eichenlaub R, Dreiseikel- like plasmid in Salmonella enterica serovar Typhimurium. FEMS mann B. Genomic and molecular analysis of phage CMP1 from Microbiol Lett 2015;362:fnv044. Clavibacter michiganensis subspecies michiganensis. Bacteriophage 71. Asare PT, Jeong T-­Y, Ryu S, Klumpp J, Loessner MJ et al. Putative 2011;1:6–14. type 1 thymidylate synthase and dihydrofolate reductase as signa- 77. Huang S, Zhang S, Jiao N, Chen F. Comparative genomic and ture genes of a novel bastille-­like group of phages in the subfamily phylogenomic analyses reveal a conserved core genome Spounavirinae. BMC Genomics 2015;16:582. shared by estuarine and oceanic Cyanopodoviruses. PLoS One 2015;10:e0142962. 72. Muniesa M, Colomer-­Lluch M, Jofre J. Could bacteriophages 78. Kozloff LM, Lute M, Crosby LK. Bacteriophage T4 virion base- transfer antibiotic resistance genes from environmental bacteria plate thymidylate synthetase and dihydrofolate reductase. J Virol to human-­ body associated bacterial populations? Mob Genet 1977;23:637–644. Elements 2013;3:e25847. 79. Peters DL, McCutcheon JG, Stothard P, Dennis JJ. Novel Steno- 73. Balcazar JL. Bacteriophages as vehicles for antibiotic resistance trophomonas maltophilia temperate phage DLP4 is capable of lyso- genes in the environment. PLoS Pathog;10:e1004219. genic conversion. BMC Genomics 2019;20:300. Five reasons to publish your next article with a Microbiology Society journal 1. The Microbiology Society is a not-for-profit organization. 2. We offer fast and rigorous peer review – average time to first decision is 4–6 weeks. 3. Our journals have a global readership with subscriptions held in research institutions around the world. 4. 80% of our authors rate our submission process as ‘excellent’ or ‘very good’. 5. Your article will be published on an interactive journal platform with advanced metrics. Find out more and submit your article at microbiologyresearch.org. 13