Adaptation to the cervical environment is associated with increased antibiotic susceptibility in Neisseria gonorrhoeae

Neisseria gonorrhoeae is an urgent public health threat due to rapidly increasing incidence and antibiotic resistance. In contrast with the trend of increasing resistance, clinical isolates that have reverted to susceptibility regularly appear, prompting questions about which pressures compete with antibiotics to shape gonococcal evolution. Here, we used genome-wide association to identify loss-of-function (LOF) mutations in the efflux pump mtrCDE operon as a mechanism of increased antibiotic susceptibility and demonstrate that these mutations are overrepresented in cervical relative to urethral isolates. This enrichment holds true for LOF mutations in another efflux pump, farAB, and in urogenitally-adapted versus typical N. meningitidis, providing evidence for a model in which expression of these pumps in the female urogenital tract incurs a fitness cost for pathogenic Neisseria. Overall, our findings highlight the impact of integrating microbial population genomics with host metadata and demonstrate how host environmental pressures can lead to increased antibiotic susceptibility.

N eisseria gonorrhoeae is the causative agent of the sexually transmitted disease gonorrhea. Antibiotics have played a key role in shaping gonococcal evolution [1][2][3] , with N. gonorrhoeae gaining resistance to each of the first-line antibiotics used to treat it [4][5][6] . As N. gonorrhoeae is an obligate human pathogen, the mucosal niches it infects-most commonly including the urethra, cervix, pharynx, and rectum-must also influence its evolution 7 . The gonococcal phylogeny suggests the interaction of these factors, with an ancestral split between a drug-susceptible lineage circulating primarily in heterosexuals and a drug-resistant lineage circulating primarily in men who have sex with men 3 .
Despite the deeply concerning increase in antibiotic resistance reported in gonococcal populations globally 8 , some clinical isolates of N. gonorrhoeae have become more susceptible to antibiotics 9,10 . This unexpected phenomenon prompts questions about which environmental pressures could be drivers of increased susceptibility and the mechanisms by which suppression or reversion of resistance may occur. To address these questions, we analyzed the genomes of a global collection of clinical isolates together with patient demographic and clinical data to identify mutations associated with increased susceptibility and define the environments in which they appear. We find that loss-of-function mutations in the efflux pump component mtrC are significantly associated with both antibiotic susceptibility and cervical infections, demonstrating how antibiotic and mucosal niche selective pressures intersect. In support of a model in which efflux pump expression incurs a cost in this niche, we also observe enrichment of loss-of-function mutations in cervical isolates in another efflux pump in N. gonorrhoeae and in urogenitallyadapted N. meningitidis. Our findings demonstrate how shifts in environmental pressures experienced by pathogenic Neisseria can lead to loss of efflux pump function and suppression of antibiotic resistance.

Results
Unknown genetic loci influence antibiotic susceptibility. We first assessed how well variation in antibiotic resistance phenotype was captured by the presence and absence of known resistance markers. To do so, we assembled and examined a global dataset comprising the genomes and minimum inhibitory concentrations (MICs) of 4852 isolates collected across 65 countries and 38 years ( Fig. 1, Supplementary Table 1). We modeled log-transformed MICs using multiple regression on a panel of experimentally characterized resistance markers for the three most clinically relevant antibiotics 5,11,12 (Supplementary Data 1). This enabled us to make quantitative predictions of MIC based on known genotypic markers and to assess how well these markers predicted true MIC values. For the macrolide azithromycin, we observed that 434/4505 (9.63%) isolates had predicted MICs that deviated by two dilutions or more from their reported phenotypic values. The majority (59.4%) of these isolates had MICs that were lower than expected, indicative of increased susceptibility unexplained by the genetic determinants in our model. Overall MIC variance explained by known resistance mutations was relatively low (adjusted R 2 = 0.667), in agreement with prior studies that employed whole-genome supervised learning algorithms to predict azithromycin resistance 13 . MIC variance explained by known resistance mutations was also low for ceftriaxone (adjusted R 2 = 0.674) but higher for ciprofloxacin (adjusted R 2 = 0.937), with 2.02% and 2.90% of strains, respectively, exhibiting two dilutions or lower reported MICs compared to predictions, similarly indicating unexplained susceptibility. The predictive modeling results, therefore, suggested unknown modifiers that promote susceptibility for multiple classes of antibiotics in N. gonorrhoeae.
GWAS identifies a susceptibility-associated variant in mtrC. To identify novel antibiotic susceptibility loci in an unbiased manner, we conducted a bacterial genome-wide association study (GWAS). We used a linear mixed model framework to control for population structure, and we used unitigs constructed from genome assemblies to capture SNPs, indels, and accessory genome elements (see "Methods") [14][15][16] . Unitigs are a flexible representation of the genetic variation across a dataset that are constructed using compacted de Bruijn graphs and have been previously applied as markers for microbial GWAS 16 . We performed a GWAS on the sequences of 4505 isolates with associated azithromycin MICs using a Bonferroni-corrected significance threshold of 3.38 × 10 −7 . The linear mixed model adequately controlled for population structure ( Supplementary Fig. 1), and the proportion of phenotypic MIC variance attributable to genotype (i.e., narrow-sense heritability) estimated by the linear mixed model was high (h 2 = 0.97). In line with this, we observed highly significant unitigs with positive effect sizes corresponding to the known resistance substitutions C2611T and A2059G (E. coli numbering) in the 23S ribosomal RNA gene (Fig. 2) 17 . The next most significant variant was a unitig associated with increased susceptibility that mapped to mtrC (β, or effect size on the log 2 -transformed MIC scale = −2.82, 95% CI [−3.06, −2.57]; p-value = 2.81 × 10 −108 ). Overexpression of the mtrCDE efflux pump operon has been shown to decrease gonococcal susceptibility to a range of hydrophobic acids and antimicrobial agents 4,18 , and conversely, knockout of the pump results in multidrug hypersusceptibility 19 . To assess whether this mtrC variant was associated with increased susceptibility to other antibiotics, we performed GWAS for ceftriaxone (for which MICs were available from 4497 isolates) and for ciprofloxacin (4135 isolates). We recovered known ceftriaxone resistance mutations including recombination in the penA gene and ciprofloxacin resistance substitutions in DNA gyrase (gyrA). In agreement with the known pleiotropic effect of the MtrCDE efflux pump 19 , we observed the same mtrC unitig at genome-wide significance associated with increased susceptibility to both ceftriaxone (β = −1.18, 95% CI [−1.34, −1.02]; p-value = 2.00 × 10 −44 ) and ciprofloxacin (β = −1.29, 95% CI [−1.54, −1.04]; p-value = 1.87 × 10 −23 ) (Fig. 2). Across all three drugs, heritability estimates for this mtrC variant were comparable to that of prevalent major resistance determinants (azithromycin h 2 : 0.323; ceftriaxone h 2 : 0.208; ciprofloxacin h 2 : 0.155), indicating that unexplained susceptibility in our model could be partially addressed by the inclusion of this mutation.
Annotation of the mtrC unitig revealed that it represented a two base pair deletion in a 'GC' dinucleotide hexarepeat, leading to early termination of mtrC translation and probable loss of MtrCDE activity 20 (Fig. 2 inset). We also checked whether the two base pair deletion would affect recognition by any of the gonococcal methylases 21 , but no methylase target motif sites mapped to the hexarepeat or its direct surrounding sequences. A laboratory-generated gonococcal mutant with a four base pair deletion in this same mtrC dinucleotide hexarepeat exhibited multi-drug susceptibility 20 , and clinical gonococcal isolates hypersensitive to erythromycin were shown to have mutations mapping to this locus 22 . To directly test the hypothesis that the two base pair deletion also contributed to increased susceptibility for the panel of antibiotics we examined, we complemented the mutation in a clinical isolate belonging to the multidrug-resistant lineage ST-1901 23 and observed significant increases in MICs for all three antibiotics, as predicted by the GWAS (Supplementary  Table 2).
We searched for additional mtrC loss-of-function (LOF) mutations and found six clinical isolates with genomes encoding indels outside of the dinucleotide hexarepeat that also were associated with increased susceptibility (Supplementary Fig. 2A). Ten isolates that had acquired the two base pair deletion also have a two base pair insertion elsewhere in mtrC that restores the original coding frame, suggesting that loss of MtrC function may be reverted by further mutation or recombination (Supplementary Fig. 2A). In line with this, mtrC LOF mutations have emerged numerous times throughout the phylogeny ( Supplementary Fig. 3), indicative of possible repeated losses of a dinucleotide in the hexarepeat region due to DNA polymerase slippage, which may occur at a higher rate than single nucleotide nonsense mutations 24 . In total, including all strains with mtrC frameshift mutations and excluding revertants, we identified 185 isolates (3.82%) that encoded a LOF allele of mtrC (Supplementary Table 3). Presence of the mtrC LOF mutation in isolates with known resistance markers was correlated with significantly reduced MICs ( Supplementary Fig. 4 25 . We observed an additional 10 strains with the 'GC' two base pair deletion ranging in isolation date from 1951 to 2000, indicating that mtrC LOF strains have either repeatedly arisen or persistently circulated for decades. Our results demonstrate that a relatively common mechanism of gonococcal acquired antibiotic susceptibility is a two base pair deletion in mtrC and that such mutations are globally and temporally widespread.  Table 1) was annotated with binarized resistance (ciprofloxacin) or decreased susceptibility (azithromycin, ceftriaxone) values. Annotation rings are listed in order of ciprofloxacin, ceftriaxone, and azithromycin from innermost to outermost. For ciprofloxacin, MIC < 1 μg/ml is light purple, and MIC ≥ 1 μg/ml is dark purple. For ceftriaxone, MIC < 0.125 μg/ml is light blue, and MIC ≥ 0.125 μg/ml is dark blue. For azithromycin, MIC ≤ 1 μg/ml is light pink, and MIC > 1 μg/ml is dark pink. Branch length represents total number of substitutions after removal of predicted recombination.
the female murine genital tract model of gonococcal infection, and overexpression of mtrCDE results in substantial fitness benefits for dealing with both antimicrobial and environmental pressures [26][27][28][29] . The relative frequency of the mtrC LOF mutations we observe (occurring in~1 in every 25 isolates) thus seems unusual for a mutation predicted to be deleterious for human infection. mtrC LOF strains do not grow more or less quickly in vitro than mtrC wild-type strains, indicating that this mutation does not confer a simple fitness benefit due to reduced energetic cost 22,26,30 . Instead, we hypothesized that there are unique environments that select for non-functional efflux pump. We aggregated patient-level metadata across included studies on sex partner preferences and anatomical site of infection. Sexual behavior and mtrC genotypic information was available for 1975 isolates from individual patients. There was a significant association between mtrC LOF and sexual behavior (p-value = 0.04021; two-sided Fisher's exact test) (Fig. 3a 20-11.78], p-value = 0.0173; two-sided Fisher's exact test). These results indicate that environmental pressures unique to female urogenital infection may select for loss of the primary gonococcal efflux pump resulting in broadly increased susceptibility to antibiotics and host-derived antimicrobial peptides. mtrA LOF offers an additional level of adaptive regulation. The association of mtrC LOF mutations with cervical specimens suggests that other mutations that downregulate expression of the mtrCDE operon should also promote adaptation to the cervical niche. The MtrCDE efflux pump regulon comprises the MtrR repressor and the MtrA activator ( Fig. 2 inset), the latter of which exists in two allelic forms: a wild-type functional gene capable of inducing mtrCDE expression and a variant rendered nonfunctional by an 11-bp deletion near the 5' end of the gene 32 ( Supplementary Fig. 2B). Knocking out mtrA has a detrimental effect on fitness in the gonococcal mouse model, and epistatic mtrR mutations resulting in overexpression of mtrCDE compensate for this fitness defect by masking the effect of the mtrA knockout 27 . Prior work assessing the genomic diversity of mtrA in a set of 922 primarily male urethral specimens found only four isolates with the 11-bp deletion (0.43%) 33 , seemingly in agreement with the in vivo importance of mtrA. However, in our global meta-analysis dataset, 362/4842 isolates (7.48%) were predicted to be mtrA LOF, of which the majority (357/362, 98.6%) were due to the 11-bp deletion. Of the 4842 isolates, 268 (5.53%) had mtrA LOF mutations in non-mtrCDE overexpression backgrounds (as defined by the absence of known mtrR promoter or coding sequence mutations or mtrCDE mosaic alleles) and therefore not epistatically masked (Supplementary Table 8  dataset, 332/4838 (6.86%) of isolates were predicted to have farA LOF mutations, of which the majority (316/332; 95.2%) were due to indels in a homopolymeric stretch of eight 'T' nucleotides near the 5' end of the gene (Supplementary Fig. 2C). farA LOF mutations were associated with patient sexual behavior (p-value = 5.06 × 10 −10 ; two-sided Fisher's exact test) and site of infection (p-value = 1.78 × 10 −12 ; two-sided Fisher's exact test) and overrepresented in cervical isolates: 33 Fig. 5, Supplementary Tables 16 and 17). The phylogenetic distribution of farA LOF indicated sustained transmission ( Supplementary  Fig. 3) and overlapped with that of mtrA LOF mutations (48.9% of isolates with mtrA LOF mutations also had farA LOF mutations), potentially indicating additive contributions to cervical adaptation. Furthermore, MtrR activates farAB expression by repressing the farR repressor 37 . This cross-talk between the two efflux pump operons indicates that in mtrCDE overexpression strains where MtrR activity is impaired, the effect of farA LOF-like mtrA LOF-may be masked. We did not observe frequent LOF mutations in the sodium gradientdependent MATE family efflux pump NorM 38 or in the ATPdependent ABC family pump MacAB 39 (Supplementary Table 3). The prevalence and cervical enrichment of farA LOF mutations and the relative rarity of LOF mutations in other non-proton motive force-driven pumps suggest that cytoplasmic acidification may be a mechanism by which the female urogenital tract selects for efflux pump loss.
Meningococcal mtrC LOF is driven by urogenital adaptation. N. meningitidis, a species closely related to N. gonorrhoeae, colonizes the oropharyngeal tract and can cause invasive disease, including meningitis and septicemia 40 . We characterized mtrC diversity in a collection of 14,798 N. meningitidis genomes, reasoning that the cervical environmental pressures that select for efflux pump LOF in the gonococcus will be rarely encountered by the meningococcus. In agreement with this, the 'GC' hexarepeat associated with most gonococcal mtrC LOF mutations was less conserved in N. meningitidis; only 9684/14798 (65.4%) isolates contained an intact hexarepeat compared to 4644/4847 (95.8%) of N. gonorrhoeae isolates (p-value < 2.2 × 10 −16 ; two-sided Fisher's exact test). In this same collection, we observed mtrC LOF due to deletions in the hexarepeat region in only 82 meningococcal isolates (0.55%), with a similar frequency of 25/4059 (0.62%) in a curated dataset comprising all invasive meningococcal disease isolates collected in the UK from 2009 to 2013 41 . The observed interruption of 'GC' dinucleotide repeats, predicted to result in a lower mutation rate 42 , and the relative rarity of mtrC LOF mutations suggests that efflux pump loss is not generally adaptive in N. meningitidis. However, a urogenitally-adapted meningococcal lineage has recently emerged in the US associated with outbreaks of non-gonococcal urethritis in heterosexual patients [43][44][45] . In isolates from this lineage, the prevalence of mtrC LOF mutations was 18/207 (8.70%), substantially higher than typical N. meningitidis and comparable to the prevalence of gonococcal mtrC LOF mutations in MSW in our global dataset (4.47%). We compared the frequency of mtrC LOF mutations in the urogenital lineage to geographically and genetically matched isolates (i.e., all publicly available n = 456 PubMLST ST-11 North American isolates) and observed a significant difference in prevalence (18 out of 207 or 8.70% versus 2 out of 249 or 0.80%; OR = 11.71, 95% CI [2.75-105.37], p-value = 3.31 × 10 −5 ; two-sided Fisher's exact test). Most mtrC LOF mutations occurred due to the same hexarepeat two base pair deletion that we previously observed for N. gonorrhoeae, and in line with this, mtrC LOF in this urogenital lineage arose multiple times independently similarly to gonococcal mtrC LOF mutations (Fig. 4, Supplementary  Fig. 3). farA LOF mutations were not observed in this meningococcal lineage. We conclude that MtrCDE efflux pump LOF is rare in typical meningococcal strains that inhabit the oropharynx but elevated in frequency in a unique urogenitally-adapted lineage circulating in heterosexuals, indicative of potential ongoing adaptation to the cervical niche. Our results suggest that efflux pump loss is broadly adaptive for cervical colonization across pathogenic Neisseria.

Discussion
In an era in which widespread antimicrobial pressure has led to the emergence of extensively drug-resistant N. gonorrhoeae 46 , isolates that appear to have reverted to susceptibility still arise 9,10 , demonstrating that antibiotic and host environmental pressures interact to shape the evolution of N. gonorrhoeae. Here, we showed that frameshift-mediated truncations in the mtrC component of the MtrCDE efflux pump are the primary mechanism for epistatic increases in antibiotic susceptibility across a global collection of clinical gonococcal isolates, as suggested by prior work 20,22 . mtrC LOF mutations are enriched in cervical isolates and a frameshifted form of the pump activator MtrA exhibits similar trends, supporting a model in which reduced or eliminated mtrCDE efflux pump expression contributes to adaptation to the female genital tract. We hypothesized that the mechanism by which this occurs is through increased cytoplasmic acidification in anaerobic conditions 34 and demonstrated that LOF mutations in farA, encoding a subunit of the other proton motive force-driven pump FarAB, were likewise enriched in cervical isolates. The LOF mutations we observed in mtrC and farA primarily occurred in short homopolymeric sequences (though with low numbers of repeated units) and thus may occur at higher rates than insertions or deletions in non-repetitive regions or nonsense mutations, similar to other resistance suppressor mutations 47 , though this will need to be confirmed in future experiments. In total, 42.6% of cervical isolates in the global dataset and 32.6% in the validation dataset contained a LOF mutation in either mtrC, farA, or mtrA, indicating that efflux pump downregulation via multiple genetic mechanisms is prevalent in cervical infection. These results complement prior studies suggesting that mtrR LOF resulting in increased resistance to fecal lipids plays a critical role in gonococcal adaptation to the rectal environment 48,49 and taken together suggest a model in which the fitness benefit of efflux pump expression is highly context dependent.
Other selective forces could also have contributed to the observed enrichment of LOF mutations in cervical isolates. For instance, iron levels modulate mtrCDE expression through Fur (the ferric uptake regulator) and MpeR 50 . Iron limitation results in increased expression of mtrCDE, and conversely, iron enrichment result in decreased expression, suggesting a fitness cost for mtrCDE expression during high iron conditions. Variation in environmental iron levels, such as in the menstrual cycle, may provide another selective pressure for LOF mutations to arise particularly when MtrR function is impaired through active site or promoter mutations. Differing rates of antibiotic use for gonorrhea in men and women due to increased asymptomatic infection in women might also select for mtrC LOF mutations, but this would not explain the associations we observed for the non-antibiotic substrate efflux pump farAB or the increased frequency of mtrC LOF mutations in urogenitally-adapted meningococci. RNA sequencing from men and women infected with gonorrhea demonstrated a 4-fold lower expression of mtrCDE in women, re-affirming the idea that efflux pump expression in the female genital tract incurs a fitness cost 51 .
Despite significant associations, only a proportion of cervical isolates exhibited these LOF genotypes, suggesting variation in cervix-associated pressures or indicating that cervical culture specimens were obtained before niche pressures could select for pump downregulation. This variation could also lead to mixed populations of efflux pump wild-type and LOF strains; however, because only one clonal isolate per site per patient is typically sequenced in clinical surveillance studies, we would be unable to detect this intra-host patient diversity. Targeted amplicon sequencing of LOF loci directly from patients in future studies would help to assess whether this intra-host diversity plays a role in infection and transmission. In particular, this intra-host pathogen diversity could facilitate transmission from the female genital tract to other sites of infection, where efflux pump activity incurs less of a fitness cost. In those new sites, isolates with wildtype efflux pump loci in the mixed population could selectively expand relative to LOF efflux pump strains and also serve as possible recombination donors of wild-type alleles. This standing genetic variation would, therefore, facilitate gonococcal adaptation across different mucosal niches. Additionally, while the cervix is the primary site of infection and source for culture in women, the selective pressures at play may include other sites more broadly in the female genital tract and may be influenced by the presence of other microbial species both pathogenic and commensal.
Our model extended to the other pathogenic Neisseria species, N. meningitidis, in that a urogenital clade transmitting in primarily heterosexual populations appeared to be undergoing further urogenital adaptation via the same mtrC frameshift mutation that was most commonly observed for N. gonorrhoeae.
In the absence of data on cases of cervicitis, we hypothesized that for this meningococcal lineage, efflux pump LOF emerged in the female urogenital tract and was transmitted to heterosexual men resulting in the enrichment we observed. Efflux pumps are common across Gram-negative bacteria 52 , and their loss may be a general adaptive strategy for species that face similar pressures as N. gonorrhoeae and urogenitally-adapted N. meningitidis. In support of this, clinical isolates of Pseudomonas aeruginosa with truncations in genes homologous to mtrC 53,54 and exhibiting antibiotic hypersensitivity have been obtained from cystic fibrosis patients, in whom the thick mucus in airway environments can likewise exhibit increased acidity and decreased oxygen availability 55,56 .
Our results also suggest potential therapeutic avenues for addressing the emergence of multidrug-resistant gonococcal strains. Selective knockdown of MtrCDE homologs in other bacteria via antisense RNA 57 and bacteriophages 58 has successfully re-sensitized resistant strains and enhanced antibiotic efficacy, and ectopic expression in N. gonorrhoeae of the mtrR repressor in a cephalosporin-resistant strain enhances gonococcal killing by β-lactam antibiotics in the mouse model 59 . Our population-wide estimated effect sizes for mtrC LOF mutations provide a prediction for the re-sensitization effect of MtrCDE knockdown across multiple genetic backgrounds and suggest particularly strong effects for the macrolide azithromycin (Supplementary Fig. 4). Because the correlation between MIC differences and clinical efficacy is still not well understood 60,61 , follow up studies to assess treatment efficacy differences in patients with and without mtrC LOF strains can help to quantify the expected effect of MtrCDE knockdown in the clinical context.
In summary, by analysis of population genomics and patient clinical data, we have shown that pathogenic Neisseria can use multiple avenues of efflux pump perturbation as an adaptive strategy to respond to host environmental pressures and illustrate

Methods
Genomics pipeline. Reads for isolates with either associated azithromycin, ciprofloxacin, or ceftriaxone MIC metadata were downloaded from datasets listed in Supplementary Table 1. Reads were inspected using FastQC version 0.11.7 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and removed if GC content diverged notably from expected values (~52-54%) or if base quality was systematically poor. We mapped read data to the NCCP11945 reference genome (RefSeq accession: NC_011035.1) using BWA-MEM (version 0.7.17-r1188) 62,63 and deduplicated reads using Picard (version 2.8.0) (https://github.com/ broadinstitute/picard). BamQC in Qualimap (version 2.2.1) 64 was run to identify samples with less than 70% of reads aligned or samples with less than 40X coverage, which were discarded. We used Pilon (version 1.16) 65 to call variants with mindepth set to 10 and minmq set to 20 and generated pseudogenomes from Pilon VCFs by including all PASS sites and alternate alleles with AF > 0.9; all other sites were assigned as 'N'. Samples with greater than 15% of sites across the genome missing were also excluded. We created de novo assemblies using SPAdes (version 3.12.0 run using 8 threads, paired end reads where available, and the --careful flag set) 66 and quality filtered contigs to ensure coverage greater than 10X, length greater than 500 base pairs, and total genome size approximately equal to the FA1090 genome size (2.0-2.3 Mbp). We annotated assemblies with Prokka (version 1.13) 67 , and clustered core genes using Roary (version 3.12) 68 (flags -z -e -n -v -s -i 92) and core intergenic regions using piggy (version 1.2) 69  Resistance allele calling. Known resistance determinants in single-copy genes were called by identifying expected SNPs in the pseudogenomes. For categorizing mosaic alleles of mtr, we ran BLASTn (version 2.6.0) 73 on the de novo assemblies using a query sequence from FA1090 (Genbank accession: NC_002946.2) comprising the mtr intergenic promoter region and mtrCDE. BLAST results were aligned using MAFFT (version 7.450) 74 and clustered into distinct allelic families using FastBAPS (version 1.0.0) 75 . We confirmed that horizontally-transferred mtr alleles associated with resistance from prior studies 5 corresponded to distinct clusters in FastBAPS. A similar approach was used to cluster penA alleles after running BLASTn with a penA reference sequence from FA1090. Variant calling in the multi-copy 23S rRNA locus was done by mapping to a modified NCCP11945 reference genome containing only one copy of the 23S rRNA and analyzing variant allele frequencies 76 . We identified truncated MtrR proteins using Prokka annotations, and mutations in the mtr promoter region associated with upregulation of mtrCDE (A deletion and TT insertion in inverted repeat, mtr 120) using an alignment of the mtr promoter from piggy output.
Phenotype processing and linear models. We doubled GISP azithromycin MICs before 2005 to account for the GISP MIC protocol testing change 77 . Samples with binary resistance phenotypes (i.e., "SUS" and "RES") were discarded. For samples with MICs listed as above or below a threshold (indicated by greater than or less than symbols), the MIC was set to equal the provided threshold. MICs were log 2transformed for use as continuous outcome variables in linear modeling and GWAS. We modeled transformed MICs using a panel of known resistance markers 11,12 and included the recently characterized mosaic mtrCDE alleles 5 and rplD G70D substitution 78  To visualize the continuous MICs using thresholds as on Fig. 1, we binarized MICs using the CLSI resistance breakpoint for ciprofloxacin, the CLSI nonsusceptibility breakpoint for azithromycin, and the CDC GISP surveillance breakpoint for ceftriaxone.
GWAS and unitig annotation. We used a regression-based GWAS approach to identify novel susceptibility mutations. In particular, we employed a linear mixed model with a random effect to control for the confounding influence of population structure and a fixed effect to control for isolate country of origin. Though the outcome variable (log 2 -transformed MICs) is the same, in contrast to the linear modeling approach described above, which models the linear, additive effect of multiple, known resistance mutations, regression in a GWAS is usually run independently and univariately on each variant for all identified variants in the genome, providing a systematic way to identify novel contributors to the outcome variable. Linear mixed model GWAS was run using Pyseer (version 1.2.0 with default allele frequency filters) on the 480,902 unitigs generated from GATB (version 1.3.0); the recombination-corrected phylogeny from Gubbins was used to parameterize the Pyseer population structure random effects term and isolate country of origin was included as a fixed effect covariate. To create the Manhattan plot, we mapped all unitigs from the GWAS using BWA-MEM (modified parameters: -B 2 and -O 3) to the pan-susceptible WHO F strain reference genome (Genbank accession: GCA_900087635.2) edited to contain only one locus of the 23S rRNA. Significant unitigs were annotated using Pyseer's annotation pipeline. Unitigs mapping to multiple sites in the genome and in or near the highly variable pilE (encoding pilus subunit) or piiC (encoding opacity protein family) genes were excluded, as were unitigs less than twenty base pairs in length. Due to redundancy and linkage, variants will be spanned by multiple overlapping unitigs with similar frequencies and p-values. For ease of interpretation, we grouped unitigs within 50 base pairs of each other and represented each cluster by the most significant unitig. Unitigs with allele frequency greater than 50% were also excluded as they represented the majority allele. Unitig clusters were then annotated by gene or adjacent genes for unitigs mapping to intergenic regions and further analyzed for predicted functional effect relative to the WHO F reference genome in Geneious Prime (version 2019.2.1, https://www.geneious.com).
Identifying LOF and upregulation alleles. To identify predicted LOF alleles of efflux pump proteins, we ran BLASTn on the de novo assemblies using a query sequence from FA1090 (reference genome FA19 was used for mtrA). Sequences that were full-length or approximately full-length (+ /−5 nucleotides) beginning with expected start codons were translated using Python (version 3.6.5) and Biopython (version 1.69) 79 . Peptides shorter than 90% of the expected full-length size of the protein were further analyzed using Geneious Prime (version 2019.2.1, https://www.geneious.com) to identify the nucleotide mutations resulting in predicted LOF by alignment of the nucleotide sequences. We called mtrCDE overexpression status by identifying the presence of any of the known mtrR promoter mutations, MtrR coding sequence mutations, and mosaic mtrCDE alleles.
Experimental validation. N. gonorrhoeae culture was conducted on GCB agar (Difco) plates supplemented with 1% Kellogg's supplements 80 at 37°C in a 5% CO 2 atmosphere. Antimicrobial susceptibility testing was conducted on GCB agar supplemented with 1% IsoVitaleX (Becton Dickinson) using Etests (bioMérieux) at 37°C in a 5% CO 2 atmosphere. We selected a clinical isolate (NY0195 81 ) from the multidrug-resistant lineage ST-1901 23 that contained an mtrC LOF mutation mediated by a two base pair hexarepeat deletion and confirmed via Etests that its MIC matched, within one dilution, its reported MIC. Isolate NY0195 contained mosaic penA allele XXXIV conferring cephalosporin reduced susceptibility and the gyrA S91F substitution conferring ciprofloxacin resistance 9 . We complemented the mtrC LOF mutation in this strain by transforming it via electroporation 80 with a 2 kb PCR product containing a Neisserial DNA uptake sequence and an in-frame mtrC allele, obtained by colony PCR from a neighboring isolate (GCGS0759). After obtaining transformants by selecting on an azithromycin 0.05 μg/mL GCB plate supplemented with Kellogg's supplement, we confirmed successful transformation by Sanger sequencing of the mtrC gene. No spontaneous mutants on azithromycin 0.05 μg/mL plates were observed after conducting control transformations in the absence of GCGS0759 mtrC PCR product. We conducted antimicrobial susceptibility testing in triplicate using Etests, assessing statistical significance between parental and transformant MICs by a two-sample t-test.
Metadata analysis. Patient metadata were collected from the following publications from Supplementary Table 1 that had information on site of infection: Demczuk et al. 82 , Demczuk et al. 83 , Ezewudo et al. 84 , Grad et al. 85 , Grad et al. 9 , Kwong et al. 86 , Lee et al. 15 , and Mortimer et al. 81 . Sites of infection were standardized across datasets using a common ontology (i.e., specified as urethra, rectum, pharynx, cervix, or other). Two-sided Fisher's exact test in R (version 3.5.1) was used to infer whether there was nonrandom association between mtrC LOF presence and either anatomical site of infection or sexual behavior. For sexual behavior analysis, isolates cultured from multiple sites on the same patient were counted as only one data point.