Evidence of Horizontal Gene Transfer of 50S Ribosomal Genes rplB, rplD, and rplY in Neisseria gonorrhoeae

Horizontal gene transfer (HGT) in the penA and multidrug efflux pump genes has been shown to play a key role in the genesis of antimicrobial resistance in Neisseria gonorrhoeae. In this study, we evaluated if there was evidence of HGT in the genes coding for the ribosomal proteins in the Neisseria genus. We did this in a collection of 11,659 isolates of Neisseria, including N. gonorrhoeae and commensal Neisseria species (N. cinerea, N. elongata, N. flavescens, N. mucosa, N. polysaccharea, and N. subflava). Comparative genomic analyses identified HGT events in three genes: rplB, rplD, and rplY coding for ribosomal proteins L2, L4 and L25, respectively. Recombination events were predicted in N. gonorrhoeae and N. cinerea, N. subflava, and N. lactamica were identified as likely progenitors. In total, 2,337, 2,355, and 1,127 isolates possessed L2, L4, and L25 HGT events. Strong associations were found between HGT in L2/L4 and the C2597T 23S rRNA mutation that confers reduced susceptibility to macrolides. Whilst previous studies have found evidence of HGT of entire genes coding for ribosomal proteins in other bacterial species, this is the first study to find evidence of HGT-mediated chimerization of ribosomal proteins.


INTRODUCTION
Neisseria gonorrhoeae is a sexually transmitted pathogen that causes the disease gonorrhea (Tapsall et al., 2009;Quillin and Seifert, 2018). Treatment guidelines typically recommend dual antimicrobial therapy, including ceftriaxone and azithromycin . The rapid emergence of antimicrobial resistance to azithromycin in numerous populations threatens this treatment approach and has recently led to the dropping of azithromycin from treatment guidelines in countries such as the United States Wan et al., 2018;St Cyr et al., 2020). Azithromycin [9-deoxo-9a-aza-9a-methyl-9a-homoerythromycin], a 15-membered macrolide, binds to the bacterial large (50S) ribosomal subunit, which consists of 23S rRNA, 5S rRNA, and ribosomal proteins. It interferes with protein synthesis via various mechanisms, including inhibiting the transpeptidation/translocation step, blocking the 50S peptide exit tunnel, and causing ribosomes to release incomplete peptides (Douthwaite and Champney, 2001;Champney and Miller, 2002;Unemo and Shafer, 2014).
The molecular mechanisms that cause azithromycin resistance in N. gonorrhoeae include an alteration of the drug target either by methylation of a single adenine in 23S rRNA (rrl) or chromosomal mutations (Roberts et al., 1999;Chisholm et al., 2010). Three such mutations have been found in domain V of the 23S rRNA: C2597T, A2045G, and A2046G (Chisholm et al., 2010;Pham et al., 2020). In addition, overexpressed efflux pumps [multidrug efflux pump (MtrCDE), MacAB, and mef -encoded pumps] may reduce azithromycin susceptibility . The MtrCDE efflux pump plays an essential role in azithromycin resistance. The mtr locus contains three genes (mtrC -mtrD -mtrE) within an operon that forms a MtrCDE and contains genes for a transcriptional repressor (mtrR) upstream of mtrCDE (Golparian et al., 2014). Mutations upstream of mtrC, including those within the MtrR binding region and mutations in the promoter region and mosaic-like sequences within the mtrD gene, contribute to reduced susceptibility to azithromycin (Hagman and Shafer, 1995;Zarantonelli et al., 1999;Shafer, 2018).
The above well-characterized resistance-associated mutations (RAMs) are, however, unable to explain all the variation in gonococcal susceptibility to macrolides (Ma et al., 2020). A point mutation in the 50S ribosomal protein, L4 (G70D), has been implicated in macrolide resistance in N. gonorrhoeae and other bacteria (Gregory and Dahlberg, 1999;Tait-Kamradt et al., 2000;Ma et al., 2020). Recently, a genome-wide association study conditioned on the known resistance mechanisms, with experimental validation, confirmed that the G70D mutation resulted in reduced susceptibility to azithromycin. In addition, the study identified other L4 mutations at amino acid positions 68 (G68D, G68C), 69 (T69I), and 70 (G70S, G70A, G70R, and G70duplication) associated with reduced azithromycin susceptibility (Ma et al., 2020). These mutations are all situated at the end of the L4 loop close to the macrolide binding site (Gregory and Dahlberg, 1999;Tait-Kamradt et al., 2000;Ma et al., 2020). In a recent in vitro study investigating the molecular pathways to high-level azithromycin resistance in N. gonorrhoeae, we found mutations in 50S ribosomal proteins L4, L22, and L34 that play a prominent role in the genesis of azithromycin resistance (Laumen et al., 2020). Similar mutations in L4 and L22 have been shown to cause macrolide resistance in several other bacterial species (Gregory and Dahlberg, 1999;O'Connor et al., 2004;Halling and Jensen, 2006;Diner and Hayes, 2009).
Horizontal gene transfer (HGT) occurs via conjugation, transformation and transduction. In Neisseria, natural transformation is the primary mode of HGT which occurs with remarkable frequency and efficiency and is dependent on type IV pilus complex (Biswas et al., 1977;Koomey, 1998;Chen and Dubnau, 2004). DNA uptake sequence (DUS), 5'-GCCGTCTGAA-3' is required for efficient transformation and is both strain and strand specific in N. gonorrhoeae (Goodman and Scocca, 1988;Smith et al., 1999;Ambur et al., 2007;Duffin and Seifert, 2010).
Both chromosomal mutations and HGT have been shown to play an important role in the genesis of aminocyclitol (spectinomycin), macrolide (azithromycin), and cephalosporins (ceftriaxone and cefixime) resistance (Ilina et al., 2013;Grad et al., 2016). HGT has been particularly important in the genesis of resistance at the rpsE, mtrCDE, and penA loci (Bowler et al., 1994;Ilina et al., 2013;Wadsworth et al., 2018;Yahara et al., 2018). We hypothesized that HGT in the 50S ribosomal genes might play a similar role in generating macrolide resistance in N. gonorrhoeae. Several studies suggest that ribosomal proteins that are both essential to function (e.g., L12 and L2) and not essential to function are amenable to HGT in other bacterial species (Weglöhner et al., 1997;Ühlein et al., 1998;Garcia-Vallvé et al., 2002;Miyazaki and Tomariguchi, 2019). As such, we assessed if there was any evidence of HGT in the gonococcal 50S ribosomal genes in N. gonorrhoeae. We considered all N. gonorrhoeae isolates (11,659 genomes available from Pathogenwatch at the time of data collection during April 2020). We set-out by analyzing these species that may have undergone HGT in the 50S ribosomal genes with the Neisseria commensal species (N. cinerea, N. elongata, N. flavescens, N. mucosa, N. polysaccharea, and N. subflava).

Data Extraction
As of April 2020, all 27 collections of N. gonorrhoeae whole genome sequences (WGS) that were available from the global platform for genomic surveillance of microbial pathogens, Pathogenwatch 1 , were downloaded. These comprised a total of 11,993 genomes. Plausibility checks were carried out to remove duplicates (n = 349), and in total, 11,659 WGS [N. gonorrhoeae (n = 11,644) from pathogenwatch and commensal Neisseria spp. (n = 15) from NCBI (GenBank)] were used in the study (Tables 1, 2 and Supplementary Table 1). The following outcomes were extracted from the metadata: azithromycin MIC, azithromycin SIR (susceptible/intermediate/resistant), collection label, country, genogroup, MLST, NG-STAR type, SNPs (23S and mtrR), and year. When the azithromycin MIC was only reported as S (susceptible) or R (resistant), we assigned an azithromycin value of 0.25 mg/L for S, and 2 mg/L for R isolates. Also, we defined isolates as azithromycin susceptible, intermediate, and resistant if their MICs were <0.5 mg/L, ≥0.5 and <1 mg/L, and ≥1 mg/L, respectively. The isolates were divided into three distinct eras: "pre-antibiotic" (pre-1950s), "golden" (1950-1970s), and "post-modern" (1980-twenty-first century) following the topology of Golparian et al., 2020.

Gene by Gene Analysis
A gene-by-gene approach on whole-genome sequences (WGS; n = 11659) followed by recombination analysis was carried out (Figure 1). A study-specific Neisseria scheme was created from 35 complete Neisseria genomes [N. gonorrhoeae reference genomes (n = 20) and commensal Neisseria sps. (n = 15), Table 2] using the Blast Score Ratio Based Allele Calling Algorithm (chewBBACA; Silva et al., 2018). The complete genome of N. gonorrhoeae FA1090 was used to create a training file using Prodigal, which was used in subsequent steps (Hyatt et al.,  2010). Firstly, coding sequences (CDS) for each genome were defined and compared in a pairwise manner to generate a single FASTA file containing the selected CDS. Secondly, AlleleCall implemented in chewBBACCA was carried out and followed by filtering out paralogous alleles, thus creating the wgMLST profiles used to define the cgMLST profiles. The SchemaEvaluator option allows for multiple sequence alignments of the alleles of each locus using MAFFT (Katoh et al., 2002) and the construction of a neighbor-joining tree using ClustalW2 (Larkin et al., 2007). The UniProtFinder 2 was used to retrieve the functional information of the CDS. 50S ribosomal and mtrCDE genes were selected based on the UniProt identifier, and the corresponding multiple sequence alignments of the alleles and neighbor-joining trees were extracted from the schema evaluator. The trees and the corresponding metadata were visualized using microreact (Argimón et al., 2016).

Allelic Profiling
Alleles were identified for the ribosomal genes. The FASTA files of the ribosomal genes that were generated from cgMLST were used as input sequences in MEGAX (Kumar et al., 2018). The CDS were translated, and all the variable amino acid sites that were singleton or parsimony-informative were exported. If a non-synonymous substitution (variant) site was present only in a commensal, it was excluded from variant calling and further analysis. The presence or absence of a variant for each site was denoted as "0" and "1, " respectively. The geometric mean (GM) AZM MIC was calculated for each allele.

Recombination Analysis
To identify the donors of mutation-harboring ribosomal genes (n = 23) in N. gonorrhoeae isolates, complete genomes of commensal Neisseria spp., including N. cinerea, N. elongata, N. flavescens, N. lactamica, N. mucosa, N. polysaccharea, and N. subflava were analyzed. The nucleotide alignment of each ribosomal gene derived from cgMLST were screened for the presence of recombination events using the Recombination Detection Program (RDP4) program . Recombinant events supported by at least two of the seven algorithms implemented in RDP software: RDP, GENECONV, Bootscan, Maxchi, Chimera, SiSscan, and 3Seq were used with default settings except for window size that was increased to 60 nt in RDP, to 120 in MaxChi and Chimera, and to 500 in BootScan and SiScan (Lemey et al., 2009;Gomez-Valero et al., 2011). Neighbor-joining phylogenetic trees were constructed using 1,000 bootstrap replicates (Salminen et al., 1995). We defined the minor parent as the one contributing the smaller fraction of the recombinant, and the major parent as the one contributing the larger fraction of the recombinant . Additionally, sequences in the order of 5' to 3': 100 bp upstream of rplD (621 bp), rplW (321 bp), rplB (834 bp) and 100 bp downstream of rplB (length −1,971 nt) from Neisseria commensals (n = 15) and subset of N. gonorrhoeae sequences, ECDC isolates (n = 1,054) were aligned and followed by recombination analysis using RDP4 program.

Statistical Analysis
Statistical analyses were performed using JMP Pro V14.0.0 (SAS Institute, NC, United States) and XLSTAT (Statistical and data analysis solution, New York, United States). MIC variables were log 2 -transformed for use as continuous outcome variables (Ma et al., 2020). Differences between variant and wildtype were analyzed by the Mann-Whitney test. Only mutations present in >10 isolates were used in these analyses. To evaluate the correlation between multiple pairs of variables (SNPs) and to assess the correlations between phenotypic and genotypic patterns of resistance to azithromycin, statistical correlation tests including Kendall's tau-b and Spearman's rank correlation were used (Elhadidy et al., 2020). Differences in categorical variables were evaluated by Chi-square/Fischer exact tests. A p-value of <0.001 was considered statistically significant.

Prevalence of Known RAMs
The MICs of azithromycin and the molecular determinants of known RAMs of the N. gonorrhoeae isolates analyzed in this study are summarized in Figure 2. Of the 11,644 N. gonorrhoeae isolates, AZM MIC value was not available for 2,520 (21.64%) isolates; 6,124 (52.59%) isolates were susceptible, 3,000 (25.76%) isolates were non-susceptible to azithromycin, including 1,593 (13.68%) resistant, and 1,407 (12.08%) intermediate isolates. The MIC's of azithromycin ranged from 0.008 to 512 mg/L (Figure 2). Mutations in domain V of 23S rRNA were detected in 540 (4.63%) isolates. Azithromycin-resistant N. gonorrhoeae isolates contained A2045G (n = 124) and C2597T (n = 416) substitutions with MICs of azithromycin ranging from 1 to 512 mg/L. Of note, two isolates with C2597T mutation were susceptible to azithromycin with an MIC of 0.125 mg/L (ERR3325414) and 0.25 mg/L (SRR8560311).
Additionally, recombination analysis were carried out for the closely linked rplD and rplB genes. We found two unique recombination events. The first event was supported by 6 out of 7 detection methods. The major and minor parent were N. gonorrhoeae ERR1528156 (region -1 to 277 and 959 to 1,971 nt) and N. cinerea NZ_LS483369 (region -278 to 958 nt), respectively, with N. gonorrhoeae ERR1560943 as the recombinant. The second event was supported by all the 7 methods. N. gonorrhoeae ERR155164 was the recombinant, and the sequences with the same recombinant event were present in 180 N. gonorrhoeae ECDC isolates and had corresponding minor parent N. gonorrhoeae ERR1528156 (region -1 to 4 and 1,423 to 1,971 nt) and major parent N. lactamica NZ_LS483369 (region -5 to 1,422 nt) and minor parent N. gonorrhoeae ERR1528156 (Figure 6). The non-synonymous SNPs in HGT regions of rplB, rplD, and rplY were further investigated to assess the population-wide prevalence and association with azithromycin MIC. Mutations were identified in 16 amino acid positions in L2, 5 in L4 and 1 in L25. The following mutations at amino acids positions R13C, N45H/S, S75P, A112V, V114A, A160V, and D187G/N for L2; A118T, V125A, T130I, A147G, and R157Q for L4; and L100P for L25 were analyzed. All the above mutations except L2 (R13C) and L25 (L100P) were associated with significantly higher azithromycin MICs (Figures 4D-F) across varied genetic backgrounds (Figures 4A-C). Some of these mutations were associated with established macrolide RAMs, meaning the association between the ribosomal protein mutations and elevated MICs may not be causal. Statistically significant associations were found between the presence of C2597T in 23S rRNA and mutations in L2 (N45H/S, S75P, A112V, and V114A) and L4 (A118T, V125A, A147G, and R157Q; all p < 0.001; Supplementary Tables 4, 5). Likewise, significant correlations were found between mutations in MtrR (-35A del, A39T,

Distribution of Horizontally Acquired Ribosomal Genes
A total of 66 (0.56%), 59 (0.5%), and 11,432 (98.17%) isolates belonged to pre-antibiotic (pre-1950s), golden (1950-1970s), and post-modern era (1980s to 21st century), respectively. Eightyseven (0.74%) isolates did not have any information on the isolation year. Of these, AZM MICs were available for 49 isolates of the pre-antibiotic and golden era and 8988 isolates from the post-modern era (Supplementary Table 9). The proportion of isolates with HGT in rplB increased from 7.5% in the preantibiotic era to 8.4% in golden era and 20.2% in post-modern era. A similar trend was seen for rplD with 1.5%/1.6% in the pre-antibiotic/golden eras to 20.3% in the post-modern era. The opposite trend was observed for isolates with HGT in rplY with 51%/50.8% in the pre-antibiotic/golden eras and 9.08% in postmodern era.

DISCUSSION
We found evidence of HGT in three gonococcal 50S ribosomal genes rplB, rplD, and rplY. Our analyses suggest that the horizontally transferred DNA was acquired from three commensal Neisseria species: N. cinerea, N. subflava, and N. lactamica. In Neisseria, the genetic organization of 50S ribosomal protein in the order from 5' to 3' is rplD (621 bp), rplW (321 bp) and rplB (834 bp). In order to discern the linkage between rplD and rplB about 1970 nt including rplD, rplW, rplB and 100 bp downstream and upstream of the genes were analyzed. Recombination with N. cinerea generated the rplD recombinant, followed by another recombination event with N. lactamica generating the rplB recombinant. The lengths of the recombined sequence were approximately 678 bases and 544 bases, respectively, which is within the recent estimate (2.5 kb) of the mean of the geometrically distributed DNA tract lengths transferred between donors and recipients in N. gonorrhoeae Yahara et al., 2021). It should be noted that the AT-DUS was present in N. cinerea NZ_LS483369 at 815 bp upstream of rplD (Figure 6). The transferred DNA contained mutations associated with reduced susceptibility to azithromycin: 7, 5 and 1 mutations in rplB, rplD, and rplY, respectively. These mutations were, however, also associated with other mutations that are well established determinants of macrolide resistance. Confounding could thus explain the association we observed between these novel mutations in rplB, rplD, rplY, and azithromycin MIC. It is also possible that these associations could be explained by the ribosomal mutations compensating for the fitness costs of the RAMs such as those in the 23S rRNA (Melnyk et al., 2015).
The strong association between HGT in L2/L4 and the 23S C2597T mutation could be explained by such a fitness restoring effect. Of note, studies of clinical isolates of Mycoplasma genitalium have found associations between the presence of mutations in L4 and L22 and macrolide resistance conferring mutations in 23S rRNA (Cao et al., 2010;Jensen et al., 2014). Likewise, mutations in the S5 and S12 ribosomal protein can confer resistance to spectinomycin and streptomycin, respectively, via structural changes in its rRNA target (Allen and Noller, 1989;Ilina et al., 2013). Mutations in L2 in Bacillus subtilis have also been found to confer resistance to bactobolin via altering its interaction with the 23S rRNA (Chandler et al., 2012). Whilst these hypotheses will require experimental validation, recent studies have illustrated the complexity of the effects of ribosomal mutations on bacterial phenotypes. Multiple studies have, for example, found various mutations in L4 and L22 in a range of bacteria, including N. gonorrhoeae, leading to macrolide resistance (Gregory and Dahlberg, 1999;Tait-Kamradt et al., 2000;Ma et al., 2020). Many of these studies have assumed that these mutations reduce macrolide susceptibility by narrowing the peptide exit tunnel so that macrolides are not able to bind to  the 50S ribosomal subunit (Figures 3A,B) (Moore and Sauer, 2008). However, further experimental evidence in Escherichia coli showed that the mechanism was not simple blockage of the exit-tunnel but rather due to reduced intracellular erythromycin concentration via increased expression of AcrAB-TolC efflux pumps. This effect was, in turn, thought to be mediated by increased translation of the relevant mRNAs by reducing programmed ribosome stalling. Similarly, in Mycobacterium smegmatis, antibiotic-induced mutations in four ribosomal proteins resulted in large alterations in the transcriptome and proteome, which facilitated the acquisition of mutations in other genes that in turn conferred higher-level resistance (Gomez et al., 2017). These findings suggest that mutations in the ribosomal proteins and their interactions with mutations in other genes (such as the 23S rRNA) may be more important and complex than previously appreciated. Previous studies in other organisms have found evidence of HGT in the ribosomal proteins, but in all cases, this involved uptake of the entire ribosomal gene (Makarova et al., 2001;Garcia-Vallvé et al., 2002). HGT-mediated chimerization has been previously reported for 16S rRNAs but not for ribosomal proteins (Miyazaki and Tomariguchi, 2019). To the best of our knowledge, this is the first study showing chimerization of ribosomal proteins. This finding provides further support for Miyazaki and Tomariguchi's "cradle model" of ribosomal evolution which posits that the functionally essential interactions between rRNA and ribosomal proteins (the cradle) were lockedin early on in ribosomal mutation and therefore amenable to HGT events (Miyazaki and Tomariguchi, 2019). Gonococci have a fundamentally non-clonal population structure due to high rates of intraspecies HGT (O'Rourke and Stevens, 1993). In addition, HGT from a number of commensal species has led to mosaicism of a number of gonococcal and meningococcal genes such as penA and mtrCDE (Bowler et al., 1994;Wadsworth et al., 2018;Yahara et al., 2018). HGT in rplB, rplD, and rplY genes was first detected between 1930 and 1932, which means the original events cannot plausibly be linked with antimicrobial pressure. In fact, the prevalence of L25 HGT alleles declined over time. The increase in the prevalence of isolates with L2/L4 HGT over time could, however, be explained by antimicrobial pressure.
An important limitation of our study is that it only included 15 commensal Neisseria isolates and no N. meningitidis isolates. The inclusion of a greater number of commensals would be expected to increase the probability of detecting HGT events.
Moreover, HGT events that resulted in no or few nucleotide changes would not have been detected by our methodology (Nielsen et al., 2013;Overballe-Petersen et al., 2013). Our methodology also did not allow us to confirm what phenotypic effects the variations in ribosomal proteins had. Nevertheless, we were able to characterize the extensive allelic variations present in Neisserial ribosomal proteins. For the first time, we also established that HGT can contribute to this variation via chimerization of ribosomal proteins between different Neisseria species.

AUTHOR CONTRIBUTIONS
SSM-B and CK conceptualized the study, interpreted the data, performed the statistical analysis, and wrote the first draft. SSM-B was responsible for data collection, bioinformatic analysis, and Figure 3A. All the authors read and approved the final draft.

FUNDING
The study was funded by SOFI 2021 grant -"PReventing the Emergence of untreatable STIs via radical Prevention" (PRESTIP).