Association of Neisseria gonorrhoeae Plasmids With Distinct Lineages and The Economic Status of Their Country of Origin

Abstract Plasmids are vehicles for horizontal gene transfer between bacteria, and in Neisseria gonorrhoeae plasmids can mediate high-level antimicrobial resistance (AMR). Using genomic and phylogenetic analyses, we show that plasmids are widespread in a collection of 3724 gonococcal isolates from 56 countries, and characterized the conjugative, β-lactamase and cryptic plasmids. We found that variants of the conjugative plasmid (which can mediate tetracycline resistance) and the β-lactamase plasmid expressing TEM-135 are associated with distinct gonococcal lineages. Furthermore, AMR plasmids are significantly more prevalent in gonococci from less wealthy countries, highlighting the need for further studies. More than 94% of gonococci possess the cryptic plasmid, with its absence correlated with the presence of a novel chromosomal type IV secretion system. Our results reveal the extent of plasmid-mediated AMR in the gonococcus, particularly in less wealthy countries, where diagnostic and therapeutic options can be limited, and highlight the risk of their global spread.

Plasmids are important vehicles for horizontal gene transfer (HGT) in bacteria, and frequently harbor genes encoding virulence factors, antimicrobial resistance (AMR), and properties that allow bacteria to survive in diverse niches [1]. Plasmids spread through bacterial populations by transformation and conjugation, resulting in the rapid dissemination of traits. Characterization of plasmids, including understanding their distribution in bacterial populations, is therefore key to understanding bacterial evolution, and in particular the spread of AMR.
Neisseria gonorrhoeae has developed resistance to all classes of antimicrobials [2]. This bacterium, the causative agent of the sexually transmitted infection gonorrhea, is a major global public health concern and a World Health Organization (WHO) priority antibiotic-resistant pathogen [3]. Complications from gonococcal disease include infertility, pelvic inflammatory disease, ectopic pregnancy, and neonatal conjunctivitis [2]. Gonococcal infection is also a cofactor for the acquisition and transmission of HIV [2].
Gonococcal AMR determinants are either chromosomally or plasmid-encoded [2], with recent research largely focusing on mechanisms of chromosomally mediated AMR, such as mosaic penA alleles [4]. However, N. gonorrhoeae can harbor 2 plasmids, the conjugative and β-lactamase, which can confer high-level resistance against tetracycline and β-lactams, respectively. The 2 AMR plasmids were first reported in gonococci in 1970s and 1980s [5], respectively, and had profound effects on the treatment of gonorrhea, resulting in a decreased use of benzyl penicillin and tetracycline. Although treatment failures with penicillin were often due to chromosomal mutations, tetracycline was discontinued as a first-line treatment in late 1980s owing to plasmid-mediated resistance [2]. In addition, the cryptic plasmid is highly prevalent in the gonococcus, even though it has no known function [5].
While chromosomal AMR genes, such as mosaic penA alleles, can lead to resistance against cephalosporins [4], plasmid-mediated AMR continues to pose a significant threat. A particular concern is that the most prevalent β-lactamase allele in N. gonorrhoeae, TEM-1, requires only 2 amino acid substitutions to become an extended-spectrum β-lactamase (ESBL) [12], which would have profound effects on the management of gonorrhea. Moreover, the TEM-135 allele, which carries an M182T substitution, is widespread among gonococci [13] and requires only a single change to become an ESBL [12].
In the current study, we characterized plasmids in a global collection of >3500 gonococcal isolates, and we compared their distribution with the phylogeny of the N. gonorrhoeae core genome. Of note, we found that the presence and type of plasmid are closely associated with the core genome. Our results highlight how cooperation between plasmids can have important implications for the emergence of AMR, with pConj and pbla often found together. In addition, we demonstrate that gonococcal AMR plasmids are particularly prevalent in low-and middle-income countries (LMICs), with this previously unrecognized epidemic of plasmid-mediated AMR highlighting the need for further characterization of gonococcal plasmids worldwide.

Core Genome Annotation and Visualization
The gonococcal core genome consists of 1668 loci [23] and yielded the N. gonorrhoeae core genome multilocus sequence typing (MLST) scheme (Ng_cgMLST), version 1.0 [23,24]. Deposited WGSs were analyzed with the BIGSdb genomics platform [25], which has "autotagger" and "autodefiner" functions that scan WGS against defined loci. Each different allele at every locus is given a unique integer. Distance matrices, based on pairwise allelic mismatches, were generated, and genetic relationships between gonococci visualized with a minimum spanning algorithm using Grapetree software [26].
Because MLST does not discriminate between gonococcal populations owing to frequent recombination [27], clustering was used to identify gonococcal core genome groups [23].
Briefly, each isolate was assigned a core genome sequence type (ST), grouping them based on a threshold of allelic differences (from <5 to <500) with ≥1 member of the same group. The gonococcal population structure was resolved into distinct, reproducible, and stable groups using a threshold of 400 differences to cluster isolates (Ng_cgc_400) [23].

Plasmid Annotation and Analysis
Loci from pConj, pbla, pCryp, and other mobile genetic elements, such as the VirB type IV secretion system (T4SS), are defined in PubMLST [19,28] (Supplementary Table 2). Manual scanning and curation were used to identify loci and alleles, and their absence validated by examining relevant regions using Artemis version 17.0.1 (Wellcome Sanger Institute) [29].
Gene-by-gene comparison of plasmid loci was undertaken using the BIGSdb Genome Comparator tool enabling assembly of allelic profiles based on all loci in each plasmid. Each profile was given a unique number (plasmid ST). Owing to extensive variation in pConj, a new clustering scheme, N. gonorrhoeae conjugative plasmid (Ng_cp), was derived using tools for core genome analysis. Ng_cp groups plasmids using a threshold of <3, <5, <10, or less, locus differences on pConj and were visualized using Grapetree software [26]. The distribution of pbla was assessed by analyzing the blaTEM allele, because this plasmid frequently contains truncated loci. Plasmid alignments were built using Easyfig software and basic local alignment search tool (BLAST) Ring Image Generator [30]. Long-read sequence data were available for isolates 36248, 39124, 39114, 39097, and 39089, [19], and WHO isolates L, G, and N [31] (GenBank accession nos. LT591902.1, LT591899.1, and LT591912.1, respectively).

Statistical Analysis
Mann-Whitney and × 2 tests were performed using GraphPad Prism software, version 7.

Core Genome Phylogenetic Analysis and Plasmid Characterization
Phylogenetic relationships among the 3724 gonococci in our collection were determined using the Ng_cgMLST scheme, which clusters isolates according to the alleles of 1668 core genome loci (Figure 1). Previously defined STs from MLST [32] associate with multiple isolate clusters [23], consistent with extensive recombination in N. gonorrhoeae [27]. Therefore, to identify distinct groups of gonococcal isolates, we used algorithms using different thresholds. Clusters containing isolates with ≤400 allelic differences (Ng_cgc_400) reproducibly resolve the gonococcal population into distinct, stable groups, which persist over time [23].

Distribution of Plasmids
Of note, phylogenetic analyses show that plasmids are associated with distinct core genome groups (Ng_cgc_400 groups) (

A Novel T4SS in N. gonorrhoeae
Owing to low variation in pCryp, analysis of plasmid STs allows sufficient resolution to assess the spread of this plasmid ( Figure 5A). In total, we identified 95 pCryp plasmid STs (ST pCryp ), with their distribution generally following the gonococcal population structure. For example, 65% of isolates belonging to Ng_cgc group 16 [14,15]. Of note, strains 32028 and 32076, which harbor pCryp lack vapDX, a putative TA system ( Figure 5B). In addition, 23 gonococci that lack pCryp contain vapDX on the chromosome, located within a putative VirB T4SS locus ( Figure 5C). This is the third T4SS identified in gonococcus, in addition to the T4SSs on the gonococcal genetic island [36] and on pConj [6].

Relationship between AMR Plasmids and LMICs
We also assessed the prevalence of plasmids in bacteria according to the economic status of country from where they had been isolated. Of note, pConj with tetM is significantly more frequent in isolates from LMICs compared with high-income countries (P = .007) ( Figure 6 and Supplementary  Table 6). Conversely, the frequency of pConj without tetM has no relationship with a country's gross domestic product (GDP) (P = .46) (Figure 6 and Supplementary Table 6). Therefore, the presence of genes conferring AMR is necessary for a plasmid to be associated with GDP.  Figure 3. pConj variants in Neisseria gonorrhoeae. A, Alignment of long-read sequences of pConj from the World Health Organization (WHO) [30] and strains from coastal Kenya [19]. Each ring is color coded. Plasmid pEP5289 [6] was used as the reference sequence. Loci are colored as follows: green indicates mating pair formation; red, genetic load; blue, conjugative transfer; orange, inheritance control; and light blue, replication initiation [6].

DISCUSSION
Gonorrhea is a major public health problem because of its impact on reproductive health, increasing AMR and the potential for untreatable gonorrhea in the near future [2]. Of particular concern is reduced susceptibility to third-generation cephalosporins, which are the mainstay of treatment [2], with recent research focused on chromosomally mediated AMR [2,4]. In the current study, we characterized plasmids in N. gonorrhoeae from across the world and assessed their distribution in relation to gonococcal lineages and country of origin. Plasmids play a vital role in disseminating AMR and often encode complex mechanisms to promote their stable inheritance and spread [1]. We provide the first global genomic analysis of plasmids in gonococci. Although many strains were obtained from the United Kingdom and the United States, our collection includes strains from >50 countries. There is a remarkably high prevalence of AMR plasmids in gonococci isolated in LMICs. This indicates that there is extensive spread and maintenance of plasmid-mediated AMR in gonococci circulating in LMICs and highlights the need for further studies to analyze N. gonorrhoeae, particularly from countries in South America, Africa, and Asia.
We analyzed WGS data deposited in the pubmlst.org/ neisseria [25], which allows annotation of core and accessory genomes across the Neisseria genus. We were able to assess the diversity of plasmid-encoded genes and the association of certain plasmids with the gonococcal population structure. Our results demonstrate that distinct gonococcal core genomes are closely associated with particular plasmids. For example, certain types of pConj (variants 2, 4, and 7), TEM-135 and TEM P14S , and pCryp are found in particular Ng_cgc groups. Associations between plasmids and bacterial lineages have been described in other pathogens [38,39]. Potential reasons for this include the presence of compensatory mutations, which mitigate the fitness costs of possessing a plasmid, chromosomally encoded restriction systems or Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR)-Cas, and ancestral plasmid acquisition followed by clonal expansion; further studies are underway to examine these mechanisms.
We also demonstrate an association between Ng_cgc groups and the absence of plasmids. For example, isolates belonging to Ng_cgc_400 group 15 lack plasmids. In the absence of pCryp, isolates can harbor chromosomally encoded vapD/X in a novel VirB T4SS locus. T4SSs export DNA or proteins and can contribute to virulence in human pathogens such as Helicobacter pylori and Bordetella pertussis [40]. vapD/X is a putative TA system [41] usually present on pCryp. Gonococci can possess 2 vapD genes, 1 on pConj and the other associated with pCryp. The genes are 30% identical, and have homology to the CRISPR/ Cas-associated protein Cas2, which forms a complex with Cas1 during the acquisition of foreign DNA into CRISPR loci that provide immunity against invading mobile genetic elements [42]. In addition, VapD homologues in other bacteria have , and markerless pConj (gray circles) were calculated against the total number of Neisseria gonorrhoeae strains from each country (only countries with ≥5 strains were considered; the size of the circle corresponds to the total number of isolates in each country). Gross domestic product (GDP) data for fiscal year 2020 were obtained from the World Bank [37]. Downloaded from https://academic.oup.com/jid/advance-article-abstract/doi/10.1093/infdis/jiaa003/5804211 by guest on 26 March 2020 endoribonuclease activity [41]. The presence of vapD on either pCryp or the VirB T4SS could inhibit the maintenance of other vapD-containing mobile genetic elements in the bacterium. We identified 7 pConj variants, which differ in the genetic load and mating pair formation regions in genes encoding zeta/ epsilon TA, and trbK and trbL. Zeta toxins are part of TA systems and interfere with peptidoglycan synthesis [43]. Of the 3 forms of zeta/epsilon TA found on pConj, only the zeta_1/ep-silon_1 TA has been studied so far [43]. The presence of different zeta toxins might influence retention of pConj variants, and/or the fitness or virulence of gonococci.
Certain Ng_cgc groups are associated with specific TEM alleles. For example, pbla in Ng_cgc_400 groups 20, 25, and 109 have TEM-135 alleles with a M182T amino acid substitution, which could develop into an ESBL with 1 further mutation [12]. The prospect of gonococcal TEM developing into an ESBL is worrying as gonococci can harbor ESBL bla TEM in a laboratory setting [44], and we found evidence of further substitutions in TEM-135 enzymes (eg, blaTEM allele 11). It is important to note that certain bla TEM alleles in N. gonorrhoeae can also harbor mutations that slow down the hydrolysis of penicillins, as in Canada plasmid pJRD20 [45]. Future work will determine the effect of substitutions in bla TEM identified here on susceptibility to β-lactams. Overall, pbla is much more likely to be found in strains with pConj. Mechanistically, this is likely to be a consequence of the ability of pConj to comobilize pbla [7,35]. Strains with TEM-135 also harbor variant 3 or 5 pConj, although the significance of this particular association is not known.
The high frequency of plasmid-encoded AMR among strains from LMICs is a particular concern, given the speed at which plasmids can spread in bacterial populations [5]; the high proportion of plasmidmediated resistance in some of these countries has been observed before [46][47][48]. The prevalence of plasmids in LMICs could be driven by ≥2 factors. First, the lack of association between GDP and pConj without tetM indicates that only plasmids conferring resistance are under positive selection in LMICs. In these countries, diagnostic facilities can be limited, and syndromic management of sexually transmitted infections is often used [19], with doxycycline being prescribed in line with WHO guidelines for urethritis, proctitis, and cervicitis [49]. Therefore, pConj harboring tetM is likely to be maintained and propagated through selective pressure imposed by use of doxycycline. Second, it is possible that fitness costs resulting from a plasmid are mitigated in strains circulating in LMICs. Regardless of the reason for the high prevalence of pConj, this plasmid can facilitate the cotransfer of pbla through bacterial populations [7], which in some countries carry TEM alleles that are a single change away from encoding an ESBL. Therefore, our findings emphasize the need to prioritize the development of rapid diagnostic measures and combination therapy to prevent the emergence of plasmid-mediated resistance to thirdgeneration cephalosporins [2].
In addition to AMR, plasmids can encode virulence factors [1]. Therefore the occurrence of specific plasmid variants could affect host-pathogen interactions. Indeed, vapD from Haemophilus influenzae can promote intracellular survival [50]. Similarly, vapD homologues on pCryp or pConj might enhance gonococcal intracellular survival [51]. Therefore, the characterization of plasmids, together with an understanding their distribution and phylogenetic relationships, should not only enhance our ability to monitor and predict the emergence and spread of AMR but also aid our understanding of the mechanisms underlying gonococcal virulence.

Supplementary Data
Supplementary materials are available at The Journal of Infectious Diseases online. Consisting of data provided by the authors to benefit the reader, the posted materials are not copyedited and are the sole responsibility of the authors, so questions or comments should be addressed to the corresponding author. Notes