Azithromycin Resistance through Interspecific Acquisition of an Epistasis-Dependent Efflux Pump Component and Transcriptional Regulator in Neisseria gonorrhoeae

1The causal agent of gonorrhea, Neisseria gonorrhoeae, is a Gram-negative diplococcus and an exclusively human pathogen. The prevalence of N. gonorrhoeae with resistance to azithromycin has dramatically increased in recent years from just 0.6% in 2013 to 3.6% in 2016 in the United States, 0.8% in 2013 to 4.7% in 2016 in England and Wales, and 5.4% in 2013 to 7.1% in 2015 across Europe. Additionally, reports in 2015 from both China and Japan have documented resistance in as high as 30% of the gonococcal population in some regions. This spike in resistance is alarming, as azithromycin is one of the two first-line drugs, in conjunction with ceftriaxone, recommended as dual-antimicrobial therapy for uncomplicated cases of gonococcal infection by the Centers for Disease Control and Prevention (CDC). Azithromycin is a macrolide antibiotic that inhibits protein synthesis by binding to the 23S rRNA component of the 50S ribosome, and while the majority of resistance can be explained by mutations in the 23S rRNA azithromycin binding sites (C2611T and A2059G), the genetic basis of a large fraction of resistance in the U.S. population is still unexplained, thus limiting the potential for development of molecular biology-based resistance diagnostics and restricting our understanding of the evolutionary paths to resistance.

2Gonococci are adept at acquiring chromosomally encoded antimicrobial resistance determinants as a result of their natural competence for transformation, allowing for the spread of resistance and other adaptively advantageous alleles between lineages and even across species boundaries. Extensive intragenus gene exchange has led to the concept of Neisseria as a consortium of species interconnected by allele sharing, with “fuzzy” borders permitting rapid access to new adaptive solutions. This produces genetic mosaicism within particular lineages, whereby some loci are the products of horizontal gene transfer and homologous recombination from other species. In gonococci, intragenus recombination is an important source of novel genetic variation with many observations of mosaic loci gained from other Neisseria species. However, aside from horizontal gene transfer facilitating the evolution of resistance to third-generation cephalosporins through acquisition of mosaic penA, allelic exchange has not yet been demonstrated to be the basis for resistance to any other antibiotic class in this species.

3Recent epidemiological studies from the United States, Canada, and Australia have reported an association between mosaic multiple transferable resistance (mtr) efflux pump alleles and increased resistance to azithromycin. Mosaics at mtr appear to have originated through horizontal gene exchange from other Neisseria and have been identified by high sequence homology of the repressor of the pump (MtrR) to Neisseria meningitidis and Neisseria lactamica. Mosaic mtr alleles have previously been associated with an outbreak of azithromycin resistance in Kansas City, Missouri, from 1999 to 2000, and with the majority of azithromycin resistance reported in New South Wales, Australia. While correlation between mosaic mtr and azithromycin resistance suggests causality, there is little experimental evidence to confirm the association.

4The Mtr efflux pump is comprised of the MtrC-MtrD-MtrE cell envelope proteins, which together export diverse hydrophobic antimicrobial agents such as antibiotics, nonionic detergents, antibacterial peptides, bile salts, and gonadal steroidal hormones from the cell. Overexpression of mtrCDE has been linked to Mtr-mediated resistance to diverse antimicrobial agents through increased drug export and also higher in vivo fitness in the murine genital tract. Known mutations that alter expression of the pump include the mtrC120 substitution, an adenine-to-guanine transition located 120 bp upstream of the mtrC start codon which acts as an alternative promoter for mtrCDE; deletion of a single adenine (A deletion) in the mtrCDE promoter that has been shown to repress the transcription of mtrR while simultaneously enhancing transcription of mtrCDE; and mutations that abrogate the function of MtrR by inducing premature stop codons or radical amino acid substitutions in the DNA-binding motif. Finally, the presence or absence of a Correia element with an integration host factor (IHF) binding site inserted between mtrCDE and the mtrCDE promoter has also been shown to negatively impact mtrCDE transcription in N. meningitidis, perhaps through the formation of DNA secondary structure, which could act as a weak transcriptional terminator. However, it is unclear if resistance in mtr mosaics is derived from any of these mechanisms and alteration of mtrCDE transcription.

5Here, we used a combination of population genomic and experimental approaches to dissect the mechanism of resistance in mosaics. We first assessed patterns of allelic diversity within the gonococcal population to define the boundaries of horizontal gene transfer at mtrRCDE and found that the entire mtr region is a hot spot of interspecies recombination which has introduced multiple rare and divergent mosaic alleles, likely from N. meningitidis and N. lactamica, into the gonococcal population. Strong linkage disequilibrium at mtrD and the mtr promoter region suggested the maintenance of epistatic allelic combinations; therefore, we tested for interaction effects within and between mtr loci via transformation. We discovered epistatic interactions across almost the entire mtrD gene, and also between mosaic mtrD and mosaic mtr promoter regions, that synergistically increased resistance to azithromycin and other substrates of the Mtr efflux pump. Furthermore, patterns of nucleotide diversity in this region coupled with experimental in vitro evidence suggest antibiotic-mediated selection may be acting to maintain these epistatic interactions at the population level in nature. Finally, we tested for regulatory evolution of pump components, as previous mechanisms of azithromycin resistance through the Mtr efflux pump have been demonstrated to be driven by expression. Our results support that inheritance of mosaic promoter regions increases the expression of mtrCDE while gaining mosaic mtrD alone does not. Thus, the likely mechanism of resistance in mosaics is a structural change to MtrD, which enhances the capacity of the protein to recognize or transport diverse antimicrobial agents, coupled with increased efflux through the amplified production of pump components.

6To gain insight into the evolutionary history of the mtrR transcriptional repressor and the mtrCDE pump, we analyzed patterns of diversity using the 1,102 Gonococcal Isolate Surveillance Project (GISP) isolates described by Grad et al.. A significant increase in allelic diversity was observed across mtrRCDE compared to the entire genome, with the highest diversity at mtrD (Fig. 1A and B; see Table S1 in the supplemental material). We also detected a significant enrichment of rare alleles in the population across mtrRCDE (Fig. 1A and C; Table S1). Linkage disequilibrium was strongest at mtrD and the mtr promoter region in a comparison of all pairs of single-nucleotide polymorphisms (SNPs) within mtrRCDE, with the highest linkage observed at pairs of variant sites within each of these loci (Fig. 1E and F; Table S1).

7To define interspecific admixture events within Neisseria, we characterized the genealogical sorting index (gsi) to explore gene tree topology measures of species-specific phylogenetic exclusivity. gsi ranges from 0 (no exclusivity) to 1 (monophyletic) and serves as a metric to assess allele sharing that may arise through interspecific recombination or incomplete divergence from a recent split between species. We calculated gsi values for genes in a 25-kb window encompassing mtrRCDE with conserved microsynteny between N. gonorrhoeae and other Neisseria species (see Fig. S1 in the supplemental material). This region included 29 core genes excluding mtrRCDE. To define the region-specific gsi background, we calculated gsi values across 100 bootstrap replicates for each gene in the mtrRCDE flanking region by species. The mean gonococcal gsi values for flanking genes were 0.95 for N. meningitidis, 0.99 for N. lactamica, and 0.92 for Neisseria polysaccharea (Fig. 1D; Table S2). Significant reductions in gsi values were detected across mtrRCDE compared to the 29 loci within the surrounding 25-kb region (Fig. 1D). Significant allele sharing between N. gonorrhoeae and N. meningitidis was exclusively observed at mtrR (Fig. 1D; Table S2), while significant allele sharing between N. gonorrhoeae and both N. lactamica and N. polysaccharea occurred across all mtr loci (Fig. 1D; Table S2).

8There were multiple recombined mosaic haplotypes present in the gonococcal population carrying the full-length mtrD (n = 80), mtrRCD (n = 9), and mtrRCDE (n = 20) genes and some isolates with partial mosaic mtrD with the majority of the gene homologous to native gonococcal sequence (n = 13) (Fig. 2). Of the 109 isolates with full-length mosaic mtrD, 4 isolates had mosaic mtrD genes that were 99% identical to the N. meningitidis gene, 5 isolates had mosaic mtrD alleles with 94 to 96% identity to the N. lactamica gene, and the remainder had alleles that aligned equally well to N. meningitidis and N. lactamica with identities ranging from 91 to 92%. Of the 29 isolates with mosaic promoter regions identified by Grad et al., 24 isolates had mosaic promoter regions that were 96 to 98% identical to those in N. lactamica, 4 had regions that were 99% identical to the N. meningitidis region with the presence of a 153-bp Correia element insertion, and 1 isolate had a region that was 92% similar to the N. meningitidis region but lacked the Correia element that was present in the other 4 N. meningitidis-like isolates. While we find the highest sequence similarity of mosaic sequences to N. meningitidis and N. lactamica in the NCBI database, the paucity of sequences representing commensal Neisseria presents an issue with conclusive species annotation. However, it is clear that mosaic alleles are highly divergent from the native gonococcus-like alleles, which are relatively monomorphic. All isolates with full-length mosaic mtrD had azithromycin MICs of ≥0.25 µg/ml, while all isolates with full-length mosaic mtrD and a mosaic mtr promoter had azithromycin MICs of ≥1 µg/ml (Fig. 2).

9Of the 29 mtr mosaics described by Grad et al. and identified by a mosaic mtr promoter sequence, none had the mtrC120 substitution, A deletion, 23S rRNA mutation A2059G, mutations in rplD, rplV tandem duplications, or variants of the rRNA methylase genes ermC and ermB that have been associated with or experimentally confirmed to be involved in azithromycin resistance. However, four mosaic isolates had the premature stop codon mutations in mtrR, and five had the C2611T 23S rRNA mutation. While most of the isolates with full-length mosaic mtrD and a mosaic mtr promoter had azithromycin MICs ranging from 1 to 4 µg/ml, those with the mtrR premature stop codon or the C2611T 23S rRNA mutation had higher MICs ranging from 8 to 16 µg/ml, most likely due to the additive or epistatic effects of inheriting multiple mechanisms of resistance.

10We exploited the natural competence of Neisseria to explore the potential for mosaic mtr alleles to produce an increase in resistance to azithromycin by transforming susceptible strains with either genomic DNA (gDNA) or PCR-amplified products from mosaic donors. Susceptible recipient strains for transformations included 28Bl, GCGS0353, and GCGS0465 (MIC ≤ 0.125 µg/ml; Table 1). Three strains with reported mosaic mtr alleles and azithromycin MICs of ≥1 µg/ml were selected as donors for DNA transfer (Table 1). These isolates included GCGS0276, GCGS0834, and GCGS0402. Of these mosaics, GCGS0276 had an mtrR sequence with the highest percent sequence identity to N. meningitidis isolates in the NCBI database, while GCGS0834 and GCGS0402 had mtrR genes most similar to N. lactamica isolates. Thus, we refer to these isolates as N. meningitidis-like and N. lactamica-like mtr mosaics, respectively. None of the selected donor strains had premature stop codons in mtrR or the C2611T mutation.

11Genomic DNA from Neisseria GCGS0276, GCGS0402, and GCGS0834 donors was able to transform multiple susceptible recipient isolates to resistance on selective media containing azithromycin (Table 2, genomic DNA transformant strains). To identify the locus responsible, we sequenced the genomes of 28Bl cell lines transformed with gDNA from mosaic donors and characterized SNPs that had been inherited from donor strains that were not present in the 28Bl recipient. The only common region that had been inherited across all transformants contained the mtrRCDE locus (Fig. 3). Likewise, the susceptible GCGS0353 and GCGS0465 cell lines transformed with gDNA from mosaic donors and subsequently selected for resistance to azithromycin were also confirmed to have inherited mosaic mtrD and mtr promoter alleles (Table 2, genomic DNA transformant strains). Thus, we further characterized this region using both population genomic and experimental approaches to carefully detail the genetic basis of resistance at mtr in mosaic isolates.

12Primers were designed to amplify each locus in mtr for both a N. meningitidis-like mosaic (GCGS0276) and a N. lactamica-like mosaic (GCGS0402), and the resultant amplicons were then used to transform strain 28Bl (Fig. 4A; Table S3). For strain GCGS0276, the only locus within the mtrRCDE operon that was found to increase resistance to azithromycin independently was mtrD (Fig. 4A; Table 2, PCR product transformant strains). GCGS0276 mtrD in the 28Bl background raised the azithromycin MIC threefold from 0.125 to 0.5 µg/ml, yet no single region of mtrD was able to produce the 0.5 µg/ml phenotype independently (Fig. 4B). Population genomics results indicated the presence of high linkage disequilibrium across mtrD (Fig. 1E and F). Thus, to test for possible within-gene epistatic interaction effects, we designed a series of amplicons that incorporated increasingly larger fragments of the mtrD locus (Fig. 4C and D). Amplicons that contained both the 5′ end (18 to 356 bp) and 3′ end (2356 to 2724 bp) of GCGS0276 mtrD were together sufficient to increase the 28Bl azithromycin MIC to 0.5 µg/ml (Fig. 4E). There were four changes at the amino acid level between GCGS0276 and 28Bl in these regions, two in the PN1 domain of MtrD (I48T and G59D) and two in the PC2 domain (K823D and F854L) (Fig. 5), and in total, 20 amino acid changes between the GCGS0276 and 28Bl proteins (Fig. 5). Unlike the N. meningitidis-like mtrD allele from GCGS0276, the N. lactamica-like GCGS0402 mtrD allele was not able to produce resistance to azithromycin in the 28Bl background above the 0.38-µg/ml-concentration selection plates.

13Though the mosaic GCGS0276 mtrD allele was clearly instrumental in enhancing resistance to azithromycin in strain 28Bl, it was not sufficient on its own to reproduce the level of resistance observed in the GCGS0276 donor strain. However, transformants that inherited the entire mtrRCDE operons of strains GCGS0276 and GCGS0402 had azithromycin MICs of 1 and 2 µg/ml, respectively, perfectly mirroring the donor strain resistance phenotypes (Table 2; Fig. 6A). Thus, in order to determine whether between-gene epistatic interactions resulted in donor strain resistance, amplicons were designed for each of these strains to amplify the mtrD locus in combination with other regions of the operon. Using this strategy, we found that for both strains GCGS0276 and GCGS0402, the donor strain azithromycin MICs of 1 and 2 µg/ml could be reproduced in strain 28Bl by transforming both mtrD and the mtr promoter region together (Fig. 6B; Table 2, PCR product transformant strains). Multiple mutations were present between isolates with mosaic and wild-type gonococcal mtr promoters; however, the only conserved mutation in all 29 mosaics was an A-to-C transversion in the mtrR promoter inverted repeat (Fig. 7).

14The Mtr efflux pump exports diverse antimicrobial agents from the cell. To determine whether mosaic mtr mutations were specific to enhancing azithromycin recognition or transport, or were capable of enhancing resistance to multiple substrates of the pump, we tested the susceptibilities of all strains and transformant lines to the dye crystal violet and the detergents Triton X-100 and polymyxin B.

15All azithromycin-resistant mosaic mtr isolates had higher MICs for the additional antimicrobials tested than the azithromycin-susceptible GCGS0353, GCGS0465, and 28Bl isolates (Table 1). Genomic DNA transformants screened for the presence of both a mosaic mtrD and a mosaic mtr promoter region from either N. meningitidis-like or N. lactamica-like donors and transformants with mosaic mtrRCDE or mtrpCD from GCGS0276 or GCGS0402 had elevated MICs for these drugs compared to susceptible recipient strains (Table 2). GCGS0276 mtrD in the 28Bl genomic background increased MICs for crystal violet and polymyxin B by twofold (from 4 µg/ml to 8 µg/ml) and fourfold (from 60 µg/ml to 250 µg/ml), respectively, while addition of the mosaic mtr promoter increased MICs 1 dilution higher in both cases (Table 2).

16We tested for the contribution of transcript regulatory variation to the mechanism of resistance by profiling gene expression via transcriptome sequencing (RNA-seq) of strains 28Bl, 28BlΔGCGS0276-mtrD, and 28BlΔGCGS0276-mtrRpCDE. As expression of the mtr efflux pump is inducible by exposure to antimicrobial agents, we evaluated expression before azithromycin exposure and 120 min after the addition of a sub-MIC dose of azithromycin (0.125 µg/ml) to the culture media. Across 24 libraries, a total of 106 million 50-bp paired-end reads mapped to the FA1090 reference genome. Each library had on average 4.44 ± 3.49 million mappable reads.

17We assessed the impact of mosaic mtrD on mtrRCDE mRNA expression by comparing 28BlΔGCGS0276-mtrD transformants to 28Bl and found no significant differential regulation of transcripts encoding mtr efflux pump components (Fig. 8; Table S4). To determine the effect of the mosaic mtr promoter on pump expression, we compared 28BlΔGCGS0276-mtrD and 28BlΔGCGS0276-mtrRpCDE transformants (Fig. 8; Table S4). Here, the presence of a mosaic mtr promoter resulted in the significant upregulation of mtrC, mtrD, and mtrE across conditions (false-discovery rate [FDR] of <0.0001) and upregulation of mtrR in the absence of azithromycin (FDR of 0.003).

18Genomic and surveillance efforts from across the globe have shown an association between mosaic mtr alleles and azithromycin resistance, but the causal role of these alleles in generating resistance has been unclear. Mtr-mediated resistance to macrolides has been attributed exclusively to increased mtrCDE expression. However, a majority of mtr mosaic isolates in the population genomic data set have none of the variants that alter regulation of mtrCDE, such as the mtrC120 substitution, deletion of a single adenine (A deletion) in the mtrR promoter inverted repeat, radical mtrR amino acid changes, or premature stop codon mutations in mtrR. Thus, using a combination of experimental and population genomic approaches, we set out to better characterize the mechanism of resistance in isolates with mosaic mtr sequences.

19Population genomics and phylogenetic reconstruction results support mtrRCDE as a hot spot of interspecies allele sharing. Expectations for genomic signatures of recent horizontal gene transfer events include the following: (i) local elevation of allelic diversity and an excess of rare mutations prior to the fixation, rise to intermediate frequencies, or loss of new alleles from the population; (ii) a departure from species-specific phylogenetic tree topology, indicating the admixture of alleles for particular loci; and (iii) high linkage within the recombination tract as the combined effects of recombination and mutation do not have sufficient evolutionary time to break down local associations between sites. Our results demonstrate that each of these parameters departed from the genomic background at mtrRCDE, suggesting recent horizontal gene transfer as opposed to historical admixture in introducing novel haplotypes at this locus (Fig. 1). While historical admixture can also produce patterns of polyphyly between species, the rarity of mosaic alleles and their high linkage suggests a more recent origin in the gonococcal population.

20At least 12 independent acquisitions of mosaic mtr alleles are supported by a whole-genome-based phylogeny of the U.S. gonococcal population samples (Fig. 2A), which have introduced multiple rare and divergent mtr haplotypes from N. meningitidis and N. lactamica into the gonococcal population. Multiple mosaic haplotypes were present in the population, some spanning exclusively over mtrD, while others extended over both mtrD and the mtr promoter region (Fig. 2C and E). Isolates with mosaic sequences were associated with increased azithromycin MIC, where isolates with full-length mosaic mtrD sequences had MICs of ≥0.25 µg/ml (Fig. 2D) and isolates with both mosaic mtrD and mtr promoters had MICs of ≥1 µg/ml (Fig. 2F). These results mirror reports of isolates with mosaic mtrR or mtr promoter regions having azithromycin MICs of ≥1 µg/ml; however, prior studies did not screen for mosaic mtrD. Interestingly, while admixture between N. gonorrhoeae and N. polysaccharea is also supported by depressed gsi values across mtrRCDE, we find no evidence of N. polysaccharea-like mtr alleles in the gonococcal population sample. This may suggest unidirectional transfer of N. gonorrhoeae mtr alleles into the N. polysaccharea population via horizontal gene transfer, or one ancestrally shared mtr allele in both species has risen to high frequency in the N. gonorrhoeae population.

21Though mosaic mtr alleles have previously been associated with resistance, our results conclusively demonstrate that they are responsible for resistance to diverse antimicrobial agents and show that the mechanism of resistance involves multiple synergistically acting epistatic interactions, including epistasis between multiple sites within mtrD and epistasis between MtrD and the mtr promoter region. Surprisingly, despite the sequence divergence (8%) between N. meningitidis- and N. lactamica-like mtr mosaics, mtr sequences from both generated azithromycin resistance through the same epistatic mechanism. For the N. meningitidis-like GCGS0276 mosaic, mtrD alone was capable of raising the azithromycin MIC of strain 28Bl to 0.5 µg/ml, independent of transcriptional changes to the pump’s regulation, with the donor azithromycin resistance phenotype of 1 µg/ml reproduced only by adding the N. meningitidis mtr promoter region, which increased expression of mtrCDE (Table 2; Fig. 4, 6, and 8). Similarly, although acquisition of the N. lactamica-like GCGS0402 mtrD was not sufficient on its own to enhance azithromycin resistance, transformation of both the mtrD and the mtr promoter region yielded the donor’s azithromycin MIC of 2 µg/ml (Table 2, PCR product transformant strains). All transformants that inherited mosaic mtr alleles in addition to becoming resistant to azithromycin also became resistant to diverse antimicrobials (Table 2). Thus, the mechanism of resistance in mosaics is likely derived from both structural changes to MtrD, which nonspecifically export antimicrobials more efficiently, coupled with promoter mutations that alter regulation of mtrCDE, enhancing drug export as a result of the increased availability of Mtr pump complexes.

22Addition of the GCGS0276 mosaic mtr promoter region in the 28Bl background increased expression of mtrCDE (Fig. 8); however, this isolate did not have the mtrC120 substitution, A deletion, or mutations in mtrR that have been previously demonstrated to increase mtrCDE transcription. While we were unable to pinpoint the precise mutation(s) responsible for increased mtrCDE expression in this mosaic, two polymorphisms present in strain GCGS0276 have previously been implicated in altered mtr regulation. The GCGS0276 mosaic has a 153-bp Correia element inserted between mtrC and the mtrCDE promoter with an IHF binding site, which has been demonstrated in meningococci to decrease mtrCDE transcription, possibly by acting as a weak transcriptional terminator. As we show that the addition of the mosaic promoter sequence increases mtrCDE expression in N. gonorrhoeae, the Correia element is either (i) not having the same phenotypic effect in N. gonorrhoeae and N. meningitidis or (ii) decreases expression of mtrCDE as in N. meningitidis while another mutation increases the transcription of mtrCDE to an extent that the level of transcriptional termination caused by the Correia element has no significant impact. The second polymorphism, a conserved A-to-C transversion in the mtrR promoter inverted repeat, was present in all 29 isolates with mosaic mtr promoters (Fig. 7). Interestingly, this mutation is located in the same inverted repeat as an A deletion that has been shown to both increase transcription of mtrCDE and increase resistance to diverse antimicrobials. Additionally, a previously described isolate with an azithromycin MIC of 2 µg/ml was noted to have an A-to-C transversion of another site within this repeat.

23Nearly the full-length mtrD was required for azithromycin resistance in GCGS0276 transformants, indicating the role of within-gene epistasis, rather than a single point mutation within mtrD, in generating resistance. Transformants that inherited two regions at the 5′ and 3′ ends of mtrD from strain GCGS0276 increased the azithromycin MIC of strain 28Bl to 0.5 µg/ml (Fig. 4). These two regions are part of the central pore (PN1) that stabilizes the trimeric organization of the protein and the outer periplasmic region (PC2) of the protein that may interact with MtrC to form a functional pump complex. Though we show that the 5′ and 3′ regions were both required to increase resistance in the 28Bl background, we are currently unable to exclude the possibility that additional SNPs between these two regions are involved in the interaction. Twenty additional amino acid changes were present between the GCGS0276 and 28Bl MtrD proteins (Fig. 5). Of note, none of the mutations observed in strains 28Bl and GCGS0276 have been shown to contribute to macrolide resistance in the orthologous proteins AcrB and MexB in other species, nor are they located in the direct contact site (residue 616) for macrolide recognition, the multidrug binding site of the pump (residues 608 to 619), or the sites involved in forming the proton relay network of the pump. Additional studies will be needed to determine the combination of mutations in mtrD responsible for enhanced resistance, though it is clear that whatever the causal SNPs are, they are acting to enhance resistance nonspecifically to multiple substrates of the pump (Table 2).

24The selective maintenance of epistasis between domains of MtrD in natural populations of gonococci is supported by our experimental and genomics results. Across mtrD, we observed local increases in linkage disequilibrium coupled with an increased frequency of rare mutations (Fig. 1). These signatures could be explained by the recent acquisition of neutral diversity from closely related species, with too little time for the effects of recombination and mutation to break down linkage of sites across imported DNA tracts or the spread of these mutations to higher frequencies. However, our experimental results demonstrate strong purifying selection on plates containing azithromycin after inheritance of partial mtrD mosaic haplotypes (Fig. 4), suggesting that some of the linkage within mtrD observed in natural gonococcal populations may be driven by selection maintaining allelic combinations within the gene that increase resistance to antimicrobials. We speculate that this is due to the effects of selection against partial gene recombination events between highly divergent mtrD alleles, which could produce defective lethal or inefficient protein products when under antimicrobial selection.

25Overall, our results defining the role of mosaic mtr in azithromycin resistance affirm the importance of other Neisseria species as an antibiotic resistance reservoir for N. gonorrhoeae. Moreover, whereas mosaic penA genes arise from interspecies recombinations within a single gene and confer cephalosporin resistance through novel structural forms, our findings of horizontally acquired, epistatically interacting structural and regulatory variants in mtr point to the potential complexity by which antibiotic resistance can arise through the interactions of multiple loci. Interspecies mosaicism will be an important consideration for future development of sequence-based molecular resistance diagnostics, as markers designed to amplify gonococcus-specific sequence will overlook or incorrectly diagnose resistance phenotype. Thus, as the number of commensal neisserial genome sequences increase, analyses that map the patterns and extent of interspecies recombination may be a valuable guide in understanding pathways to resistance and in designing appropriate diagnostic tools.

26Sequencing libraries were prepared using a modification of Illumina’s Nextera XP protocol. Samples were dual-indexed and pooled (n = 15 per pool). Paired-end 150-bp sequencing was conducted on an Illumina MiSeq (Illumina Corp., San Diego, CA) platform located at the Harvard T. H. Chan School of Public Health to an average depth of 40×. Previously sequenced read libraries were obtained from the NCBI’s Short Read Archive (project PRJEB2090) and the European Nucleotide Archive (projects PRJEB2999 and PRJEB7904).

27To determine the impact of interspecific recombination at mtrRCDE, we assessed patterns of allelic diversity across the mtrR transcriptional repressor and the mtrCDE pump compared to the rest of the genome for the 1,102 GISP gonococcal isolates. Reads were aligned to the FA1090 reference using Bowtie2 v.2.2.4, and variants were called using pilon v.1.16. Vcftools v.0.1.12 was used to merge resultant vcf files and calculate genome-wide values of π and Tajima’s D over 100-bp sliding windows, and r2 linkage by site. Gubbins v.2.2.0 was used to predict regions of elevated single-nucleotide polymorphism (SNP) densities. For each isolate, BLASTn was used to identify the top hit and highest percent sequence identity for mtrD and the mtr promoter to all Neisseria within the NCBI database (E value of <10−40).

28The extent of exclusive ancestry between Neisseria gonorrhoeae, N. meningitidis, N. lactamica, and N. polysaccharea was assessed using the genealogical sorting index (gsi) for each gene across a 25-kb window surrounding mtrRCDE. In brief, we downloaded de novo assemblies and raw sequencing reads from NCBI for N. meningitidis (n = 431), N. lactamica (n = 326), N. polysaccharea (n = 37), and N. gonorrhoeae (n = 1102). Raw reads were assembled with SPAdes v.3.7.0, and assemblies were aligned to the N. gonorrhoeae FA1090 reference genome (AE004969.1) using progressiveMauve (snapshot 2015-02-13 for linux-x64), since a multigenome alignment for all genomes was not computationally tractable. The sequences that were aligned to each gene within 25 kb of mtrRCDE in strain FA1090 were then extracted with custom Perl scripts (available on request) and realigned with MAFFT v.7.309. This method of pairwise alignments of de novo assemblies to the FA1090 reference identifies orthologs between N. gonorrhoeae and the other Neisseria species using both sequence identity and microsynteny, which is conserved in the genomic region surrounding mtrRCDE (see Fig. S1 in the supplemental material). We then used RAxML v.8.1.4 to reconstruct the phylogeny for each gene, using 50 bootstrap replicates and the GTRCAT substitution model. With these multispecies phylogenies, we calculated gsi with the genealogicalSorting R package.

29N. gonorrhoeae isolates were provided by the CDC (Table 1). Isolates were cultured on GCB agar medium supplemented with 1% IsoVitaleX (Becton Dickinson Co., Franklin Lakes, NJ). After inoculation, the plates were incubated at 37°C in a 5% CO2 atmosphere incubator for 16 to 18 h. Antimicrobial susceptibility testing was conducted using the agar dilution method at a range of concentrations of azithromycin from 0 µg/ml to 16 µg/ml, crystal violet from 0.25 µg/ml to 32 µg/ml, Triton X-100 from 250 µg/ml to 16,000 µg/ml, and polymyxin B from 30 µg/ml to 1,000 µg/ml. MICs were recorded after 24 h of growth. All isolate stocks were stored at −80°C in Trypticase soy broth containing 20% glycerol.

30Genomic DNA was extracted from isolates by lysing growth from overnight plates in TE buffer (10 mM Tris [pH 8.0], 10 mM EDTA) with 0.5 mg/ml lysozyme and 3 mg/ml proteinase K (Sigma-Aldrich Corp., St. Louis, MO). DNA was purified using the PureLink Genomic DNA minikit, treated with RNase A (Thermo Fisher Corp., Waltham, MA), and stored in water. Primers were designed to amplify regions of mtrRCDE. For primer pairs that did not amplify over a region containing a DNA uptake sequence, to enhance transformation efficiency, the 12-bp AT-DNA uptake sequence (AT-DUS) was added to the forward primer (5′-ATGCCGTCTGAA-3′). PCR was conducted in 50-µl volumes using Phusion High-Fidelity DNA polymerase (New England Biolabs, Ipswich, MA) using the conditions listed in Table S3 in the supplemental material. Amplified products were run on a 0.8% agarose gel, excised, and purified with the QIAEX II gel extraction kit (Qiagen Inc., Valencia, CA) to remove genomic DNA (gDNA) contamination.

31Transformations were conducted in liquid culture as described previously at 37°C in a 5% CO2 atmosphere incubator. In brief, we inoculated GCP broth (7.5 g protease peptone 3, 0.5 g soluble starch, 2 g dibasic K2HPO4, 0.5 g monobasic KH2PO4, 2.5 g NaCl, and double-distilled water [ddH2O] to 500 ml; Becton Dickinson Co., Franklin Lakes, NJ) supplemented with 1% IsoVitaleX and 10 µM MgSO4 (Sigma-Aldrich Corp., St. Louis, MO) with naturally competent piliated cells to an optical density (OD) of ~0.5. Suspensions were then incubated with approximately 1 µg of gDNA or purified PCR product for 10 min to allow DNA uptake and homologous recombination and subsequently spread on GCB with 1% IsoVitaleX and incubated for 4 h to allow for expression of novel alleles. Transformants resistant to azithromycin were selected by plating cells on GCB supplemented with 1% IsoVitaleX and 0.38 to 1 µg/ml azithromycin and picking single colonies after 18 h.

32All transformants were screened for successful transformation by testing the levels of resistance to known substrates of the Mtr efflux pump system as well as sequencing to couple the acquisition of novel mosaic alleles with gain of resistance. Transformants generated with targeted PCR products from mosaic donors were screened for mosaic alleles by sequencing mtrD and the mtr promoter region via Sanger sequencing using the GeneWiz sequencing service (GeneWiz Inc., Cambridge, MA) (Table 2). Undirected gDNA-transformed cell lines were sequenced via Illumina whole-genome sequencing as described above. To detect the alleles gained from mosaic donors, we searched for SNP homology between the transformants and mosaic donors that was not present in the 28Bl recipient strain (Fig. 3). Genomes were de novo assembled using SPAdes v.3.7.0 and annotated with prokka v.1.11. Alignment of transformant libraries to the 28Bl and mosaic references was conducted in Bowtie2 v.2.2.4 using the “end-to-end” and “very-sensitive” options, which resulted in a higher percentage of uniquely mapped reads than any other combination of preset options. SNPs and indels between read libraries and the reference genomes were called with pilon v.1.16, removing variant calls with qualities of <15.

33Cells harvested from overnight plates were suspended in GCP supplemented with 1% IsoVitaleX and 0.042% sodium bicarbonate. Cultures were incubated at 37°C for 2 h to mid-log phase and then exposed to a sublethal dose of azithromycin (AZI) (0.125 µg/ml). RNA was extracted at 0 min (pre-AZI) and 120 min (post-AZI) exposure using the Direct-Zol kit (Zymo Research, Irvine, CA). Transcriptome libraries were prepared at the Broad Institute at the Microbial ‘Omics Core using a modified version of the RNAtag-seq protocol. Five hundred nanograms of total RNA was fragmented, depleted of genomic DNA, dephosphorylated, and ligated to DNA adapters carrying 5′-AN8-3′ barcodes of known sequence with a 5′ phosphate and a 3′ blocking group. Barcoded RNAs were pooled and depleted of rRNA using the RiboZero rRNA depletion kit (Epicentre, Madison, WI). Pools of barcoded RNAs were converted to Illumina cDNA libraries in two main steps: (i) reverse transcription of the RNA using a primer designed to the constant region of the barcoded adapter with addition of an adapter to the 3′ end of the cDNA by template switching using SMARTScribe (Clontech, Mountain View, CA) as described previously; (ii) PCR amplification using primers whose 5′ ends target the constant regions of the 3′ or 5′ adapter and whose 3′ ends contain the full Illumina P5 or P7 sequences. cDNA libraries were sequenced on the Illumina Nextseq 500 platform to generate 50-bp paired-end reads.

34Barcode sequences were removed, and reads were aligned to the FA1090 reference genome. Read counts were assigned to genes and other genomic features using custom scripts. For the FA1090 reference genome, we mapped reads to either the sense or antisense strand for coding domain sequences (CDSs) (n = 1,894), tRNAs (n = 55), and rRNAs (n = 12). For intergenic regions (IGRs) (n = 1,722), we mapped to each antiparallel strand. Differential expression analysis was conducted in DESeq2 v.1.10.1.

35All genomic and transcriptomic read libraries generated in this study are available from the NCBI SRA database (accession numbers PRJNA475134 and PRJNA475139); all Sanger sequencing files are available in the NCBI database (accession numbers MH576918 to MH576947).