Population structure of Neisseria gonorrhoeae based on whole genome data and its relationship with antibiotic resistance

Neisseria gonorrhoeae isolates

Sixty-one N. gonorrhoeae isolates of diverse origin were obtained. These included isolates from the Gonococcal Isolation Surveillance Program (GISP) site covering Atlanta, Miami, New York city and North Carolina in the United States (n = 21), from Canada (primarily, Saskatchewan) and Chile (n = 1) (Vidovic et al., 2014) (n = 24), and from WHO global collaborations; Sweden (n = 7), Norway (n = 3), Japan (n = 2), Austria (n = 1), Pakistan (n = 1), Philippines (n = 1), and Australia (n = 1). Phenotypic determination of the minimum inhibitory concenterations (MICs) of all isolates was performed using the agar dilution method or the Etest method (bioMerieux), according to the instructions from the manufacturer. The strains sequenced in this study were tested for resistance to primarily three antibiotics, tetracycline, azithromycin and cefixime, with breakpoints for resistance set at 2, 2.0, and 0.25 µg/mL, respectively, based on the CDC MIC (minimum inhibitory concentration) breakpoints for testing in the GISP protocol (http://www.cdc.gov/std/gisp/gisp-protocol07-15-2010.pdf). Antibiotic resistance profiles of the Canadian strains have been previously reported (Vidovic et al., 2014). Details of the different isolates with their NCBI accession numbers are presented in Table 1.

Sequence generation and assembly

The N. gonorrhoeae strains were shotgun (WGS) sequenced using the Illumina HiSeq™ instrument (Illumina, San Diego, California, USA), utilizing libraries prepared from 5 µg of genomic DNA for each sample. The average sequencing coverage was 225. The sequencing reads were filtered using the prinseq-lite algorithm (Schmieder & Edwards, 2011) to ensure only sequence reads with average phred score ≥30 were used. The reads for each project were then assembled de novo, using the velvet assembler program (Zerbino & Birney, 2008). The optimal kmer length for each assembly prior to assembly was determined using the velvet optimizer algorithm (Gladman & Seemann, 2012). Data was deposited in the NCBI Sequence Read Archive public database (Accession # SRA099559) (Table 1). For this study, we included an additional 14 draft genome sequences of N. gonorrhoeae strains, downloaded from the NCBI draft genomes database (NCBI Bioproject numbers: PRJNA55649, PRJNA55651, PRJNA55653, PRJNA55905, PRJNA46993, PRJNA55657, PRJNA55655, PRJNA55659, PRJNA55661, PRJNA55663, PRJNA55665, PRJNA55667, PRJNA55669, PRJNA55671, and the reference genome sequence Ref_FA_1090 (NC_002946.2).

Genome-wide phylogeny and pangenome analysis

The assembled genomes were annotated individually using the NCBI PGAP annotation pipeline to give predicted proteome for each of the strains. The orthologs were determined by OrthoMCL (Li, Stoeckert & Roos, 2003), which uses bi-directional BLASTP scores of all the protein sequences to perform Markov clustering in order to improve sensitivity and specificity of the orthologs. For the OrthoMCL analysis, we used a BLASTP E-value cut-off of 1e-05, and inflation Markov clustering parameter of 1.5. Core genes were defined as the orthologous genes that are shared among all the N. gonorrhoeae strains used in this analysis.

The nucleotide sequences of all the core genes were concatenated together and core-gene nucleotide alignment was conducted using progressive MAUVE (Darling et al., 2004). Similarly, whole genome amino-acid alignment was also generated by concatenating the deduced amino-acid sequences of all the core genes generated using MUSCLE (Edgar, 2004), and to form a super protein alignment. Homoplasious sites were removed from the whole-genome nucleotide alignment using the Noisy software (Dress et al., 2008). The protein alignments were filtered by GBLOCKS (Talavera & Castresana, 2007) using default settings to remove regions that contained gaps or were highly diverged. A maximum likelihood (ML) tree from the same data set was created using the GTR and JTT substitution models for the nucleotide and protein alignment respectively and the GAMMA evolutionary model (Stamatakis, 2014). The majority rule-consensus tree was generated from 200 bootstrap replicates of the model. Linear regression of the root-to-tip distances against the year of isolation was performed using the Path-O-Gen tool (http://tree.bio.ed.ac.uk/software/pathogen/).

Multi-locus sequence typing (MLST) locus analysis

MLST is a genotyping tool for Neisseria based on sequencing of 7 core housekeeping genes (Jolley & Maiden, 2010). There are currently close to 11,000 individual Neisseria sequence profiles in the publicly available MLST database (http://pubmlst.org). We utilized a custom python script mlstBLAST.py (http://sourceforge.net/projects/srst/files/mlstBLAST/) to perform a BLAST search of these genes across all the strains in our data set and identified the sequence type (ST) for each strain. Novel alleles of the locus and STs were submitted to the MLST database. A phylogeny of the concatenated DNA sequences of all the N. gonorrhoeae STs in the MLST public database was created using the neighbor joining distance matrix approach of the PHYLIP (Felstein, 1989). Mean nucleotide distance for the sequence alignments and MLST genes was computed using MEGA software (Tamura et al., 2013).

Estimating population parameters and homologous recombination

ClonalFrame (Didelot & Falush, 2007) utilizes a statistical framework to reconstruct the clonal genealogy as well as identify the regions along the genomes that has been affected both by recombination and mutation. The model uses a Bayesian approach to predict the phylogenetic relationship in the sample set, given the whole genome sequence alignment data. The input genome alignment data was the core genes (n = 1,189) alignment generated from MAUVE. Four independent ClonalFrame runs were performed for 40,000 iterations, with the first 20,000 discarded as burn-in. This allowed the model parameters to converge, and each of the 4 runs were checked for the consistency of the estimated parameters as well as the consistency of the topology of the inferred clonal genealogies.

Population structure analysis

The program BAPS (Bayesian Analysis of Population Structure) version 5.3 (Corander & Marttinen, 2006; Tang et al., 2009) was used to infer the underlying population structure of the 76 N. gonorrhoeae strains in the sample set. SNPs from the core MAUVE alignment, with gaps removed were converted to a BAPS input file, which is a representation of all the polymorphic loci in the multi-sequence alignment. BAPS applied a Bayesian model to predict the likelihood of a population structure, given the input data and non-parametric assumption approach to trace ancestry of the different individuals in the sample set. For the mixture analysis we used the ‘Clustering of individuals’ approach. We ran a preliminary analysis to evaluate the approximate number of genetically differentiated groups using a vector from 2 to 40 K values, where K is the maximum number of groups. Given that 5 groups was the K value with the best log likelihood, we ran a second analysis using from 3 to 7 K values and again the best K value was 5 groups. We used the ‘Admixture based on mixture clustering’ module for the admixture analysis. For the analysis, the minimum population and the admixture coefficient for the individuals was then set to 5. For the reference individuals from each population and the admixture coefficient for reference individuals we used the values as described by Castillo-Ramírez et al. (2012). In addition, population structure analysis of the sample set using the fineSTRUCTURE tool (Lawson et al., 2012) was performed. The fineSTRUCTURE analysis was a two step process-(1) ChromoPainter algorithm was used to generate the co-ancestry matrix from the genome-wide haplotype data using the linkage model. (2) The fineSTRUCTURE algorithm performed a model-based clustering using a Bayesian MCMC approach to predict the likelihood of a population structure given the input data and non-parametric assumption approach to trace ancestry of the different individuals in the sample set. The fineSTRUCTURE approach was used to corroborate the findings from the BAPS population structure analysis.

Mapping the movement of DNA between Neisseria gonorrhoeae clades

We traced the flow of recombination between strains into five different subgroups in the phylogeny determined from the subgroups of the population defined by the BAPS analysis. We created a BLAST database of the whole genome sequence of all 76 strains in the sample set and included 14 whole genome sequences of all other Neisseria species (NC_008767.1, NC_014752.1, NC_017512.1, NC_017516.1, NC_003112.1, NC_010120.1 NC_017501.1, NC_017514.1, NC_017517.1, NC_003116.1, NC_013016.1, NC_017505.1, NC_017515.1, NC_017518.1) that are present in the NCBI database. Next, we performed a BLASTN search for each of the genomic region within the strains identified by ClonalFrame to be under recombination, selecting the best hit within the sequences in the database we created, with an identity of >98%, to be the source of the recombined region. We also removed BLAST matches found in strains from similar subgroups as the source of the recombined region. We used the migest package (http://cran.r-project.org/web/packages/migest/) implemented in the R statistical language to create a circular representation of the matrix of relationship between the subpopulations identified by BAPS based on the purported recombination between strains in the different subgroups. We also supplied migest with the matrix from BAPS admixture analysis and recreated the circular flow of recombination across only the subpopulations as defined by BAPS.

Comparison of nucleotide substitution rates

Amino acid sequences were aligned using MUSCLE sequence aligner (Edgar, 2004). The amino acid sequence alignment was converted to nucleotide alignment based on the corresponding gene sequence using PAL2NAL (Suyama, Torrents & Bork, 2006) and we implemented the YN00 method of the PAML package (Yang, 2007) to calculate the pairwise dN/dS ratios for the strains (Rocha et al., 2006). The contribution of each strain to the overall variation in the dN/dS rates across the sample set was estimated using ANOVA (Analysis of Variance) approaches.

Analysis of positive selection

For the analysis of positive selection within core genes of the strains in the sample set, we first identified and removed core genes that have signals of homologous recombination using three methods of Pairwise Homoplasy Index (PHI), Neighbor Similarity Score (NSS) and the maximum χ² method. The three methods are implemented in the PhiPack package (Sawyer, 1989). A window size of 50 nucleotides was used to run the methods in the package, and genes shown to have significant probability of homologous recombination by a majority of the methods were not used for the positive selection analysis. Next, we identified core genes under positive selection using codeml of PAML tool version 4.7 (Yang, 2007). We applied the branch-sites test for positive selection Model A test 2 of the tool, to identify genes under positive selection population groups. For each of the clades, we performed the Likelihood Ratio Test (LRT) for two hypotheses—the null hypothesis is the existence of neutral selection as implemented in the null model versus the alternative hypothesis implemented in the test model for positive selection. The LRT was performed to a degree of freedom of 1, and we corrected for multiple testing using the False discovery rate approach (FDR) (Benjamini & Hochberg, 1995). We further identified the Gene Ontology (GO) terms and functional characterizations of the genes under positive selection (see Table 2) and performed an enrichment test for functionality of these genes using the blast2go test pipeline (Götz et al., 2011).

Confirming known predictors of antibiotic resistance phenotype

We downloaded from NCBI reference DNA sequences of resistance determinants that have been shown in the literature to underlie the resistance phenotype we have observed in our sample set (see Table 3), and performed a BLASTN search for each of these DNA sequence regions across all the strains in the database of whole genome sequences. For convenience, the contigs for each assembly were ordered into one pseudocontig after tiling to the reference genome FA1090, using the ABACUS tool (http://abacas.sourceforge.net/).

We selected the top hit (with identity match of 98% or more) for each sequence (strain) in the database and parsed the alignment between the query and the subject sequence in the database for the presence or absence of the underlying resistance genetic mutations as suggested in the literature.

Genome-wide homologous recombination in diverseN. gonorrhoeae

We sequenced 61 recent clinical isolates primarily from the US and Canada but also single representatives from other countries, including Japan, Pakistan, Australia, Austria, Philippines, Norway and Sweden using the WGS approach to a high average coverage (average 225-fold read redundancy). De novo assemblies based on these data produced a set of contigs that represented draft, unordered representations of the genomes with high sequence quality. A preliminary phylogeographical analysis of the Canadian isolates (n = 23) was recently published (Vidovic et al., 2014). For the analysis, we included the 14 N. gonorrhoeae draft NCBI genome sequences (12 from the US and 2 from Europe) and the genome sequence of the FA1090 N. gonorrhoeae reference strain. The 76 were assigned into 23 previously described MLST STs and four new STs (10931,10932,10933,10934). The genetic diversity (measured as pairwise nucleotide distances of MLST loci) of the strains in this study was about half that of the N. gonorrhoeae strains as a whole (0.001 substitutions per site in our study compared to 0.002 in the large MLST set), and the STs from our sample set were represented across the different clades of a phylogeny of housekeeping genes of N. gonorrhoeae strains found in the MLST database (see Fig. S5). Alignment of the shotgun assembly to reference genome FA1090 (NC_002946.2), yielded 10,962 SNPs in the core region (conserved in all strains). The average per nucleotide diversity in the core genome regions was 0.003.

Homologous recombination is known to play a role in shaping bacteria populations (Didelot & Maiden, 2010). The ClonalFrame tool (Didelot & Falush, 2007) detected 952 independent recombination events, covering more than 50% of the reference genome. The average size of the recombination regions identified was 360 base pairs. The estimate for the ratio of effects of recombination and mutation (r/m) for our strain set was 2.2, a relatively high value for bacterial species (Didelot & Maiden, 2010) (Fig. S3) and quite similar to the r/m estimate of 1.9 based on the whole genome alignment on a less genetically diverse group of N. gonorrhoeae strains reported by Grad et al. (2014).

We constructed a maximum likelihood phylogeny of the core genome of the 76 strains using the RAxML program (excluding regions identified as potentially recombinant) (Fig. 1). This tree had similar topology to the clonal frame that determined by the eponymous software (Fig. S3). The tree showed multiple clades but the there was no strong signal of genetic isolation by distance at the continental scale. The rate of the molecular clock was estimated to be 8.93 × 10⁻⁶ mutations per year based on the slope of the regression of the root-to-tip divergence with isolation dates (see Fig. S2). This value was similar to those obtained in other bacterial studies, ranging from 8.6 × 10⁻⁹ to 2.5 × 10⁻⁵ (Zhou et al., 2013). However, because the temporal signal was weak in the root-to-tip analysis (Fig. S2), we did not use these data for Bayesian phylogeny analysis using the BEAST phylogeny tool (Drummond & Rambaut, 2007).

Neisseria gonorrhoeae population structure and biogeography

Given that recombination was frequent in these genomes, we sought to evaluate the genetic substructure of the population. We used two complementary methods. BAPS (Tang et al., 2009) predicts the likelihood of a population structure given the input data and uses a non-parametric assumption approach to trace ancestry. fineSTRUCTURE (Lawson et al., 2012), on the other hand, uses similar methods of predicting population substructure, but to a finer detail and does not assume a prior optimum number of subpopulations (K). The BAPS tool identified 5 subgroups within the N. gonorrhoeae population from the strains within the sample set (Fig. 3). As expected, members with the same subgroup ancestry generally were found near each other when mapped on the ML phylogeny constructed using the nonrecombining portion of the genome. On the other hand, fineSTRUCTURE identified 30 genetic subgroups within our sample set (Fig. S4). However, individual members of each of the fineSTRUCTURE subgroups belong to the same BAPS subgroup. Each of the five BAPS subgroups contained strains from multiple continents based on geography or location of isolation (Fig. 1). It was particularly interesting that each BAPS cluster had at least one US strain and one Canadian strain. The BAPS analysis revealed a complex relationship between Group 3 and 5, with the latter group 5 separated into two group 3 clades (Fig. 1). Group 3 strains in clades closely related to group 5 showed significant genetic import from group 5. It is possible that the extent of admixture occurring in group 3 and 5 may have caused misidentification. Also, strains of sequence type (ST) 1901, which is the most abundant ST in our sample set all belong to subgroup 1, hinting at a correlation between BAPS subgrouping and MLST.

We assessed patterns of genetic drift effects in the population by estimating the pairwise substitution rates between all the core gene orthologs for the strains and determining the mean dN/dS ratio for each strain. The mean pairwise dN/dS ratios for each strain are shown in Fig. 4. There was significant variation in the mean dN/dS ratios among the strains (ANOVA p-value = 2.0e-16). The overall mean of the dN/dS estimate was 0.3184, similar to the 0.402 value estimated for the bacterial pathogen Chlamydia trachomatis (Joseph et al., 2012). The mean dN/dS ratio for the five subgroups respectively was (0.32412976, 0.33325164, 0.31273103, 0.30952504 and 0.30990092). There is no significant difference between group mean dN/dS ratios (p-value =0.921, t test of means).

The mean dN/dS ratio for strains from the Canadian region was 0.3279, which was above the overall mean ratio, while that for strains collected in the US was 0.31708, which is below the overall pairwise dN/dS mean ratio for the sample set. This was also a statistically significant difference (p-value 0.0018 for t test of means), suggesting a possible geographical effect within this subset of strains.

Genetic admixture within N. gonorrhoeae and with other Neisseria species

In order to understand the flow of genetic information between the strains from five different subgroups defined by the BAPS analysis (Fig. 3) as well as strains from other Neisseria species, we used two independent approaches. The first was to search each of the 952 recombination regions identified by ClonalFrame for a best BLASTN match from another subgroup or Neisseria species (We created a blast database of the 76 genomes from this study and representative strains from the Neisseria genus.) (Fig. S3). In parallel, we also counted the occurrence of co-ancestry of genetic markers revealed by the BAPS analysis. Both the BAPS and BLAST analyses suggested that group 3 was the most common nexus of homologous recombination between other clades, consistent with its basal phylogenetic status. In the BAPS-based network groups 1 and 2, and to a greater extent, group 5, were primarily DNA donors to group 3 (Fig. 2B). But this pattern was less visible in the BLAST network (Fig. 2A). It is notable that more than 90% of the recombination with strains from other Neisseria species occured in groups 2 and 3. Group 5 stood out as a significant source of genetic exchange into strains in group 3.

The genetic relatedness of the strains in the sample set or the purported sharing of genetic materials across the different subgroups shown by the BAPS figure paralleled the pattern revealed by the BLAST clonal frame analysis (p-value = 0.048, Mantel test for comparing the distance matrices of the five populations between both methods). The exchange of genetic materials from other Neisseria species was not accounted for in the BAPS admixture analysis. Based on the BLAST analysis, the proportion of DNA transferred within N. gonorrhoeae compared to arriving from Neisseria strains the species was 729 out of 849 intra-specific genetic events. This finding is line with the “fuzzy species” concept of Fraser, Hanage & Spratt (2007): while N. gonorrhoeae is not sexually isolated, DNA flow seemed predominantly through intra-specific exchanges.

Genes under positive selection

Of the 1,189 core genes, we identified 352 genes as likely to contain past recombination histories using the PHIPACK tests (Sawyer, 1989). Thirty-one genes within the subset of 837 non-recombining core genes were found to be under positive selection using the tests implemented by the PAML software (‘Materials and Methods’). BAPS subgroup 5 had the highest number of core genes under selection (14) followed by subgroup 3 (7). While we found no significant enrichment of genes under positive selection in any of the functional classes in the Gene Ontology (GO) database, the functions of the best match proteins from genes under positive selection can be broadly classified to genes involved in DNA or RNA synthesis of gene expression, membrane or transport proteins, and, to a lesser extent, genes involved in metabolic pathways in the bacterial cell (Data S2 spreadsheet). Of the 352 genes found to have signals of recombination, we found no significant enrichment of the genes in any of the functional classes in the GO database. The functions of these genes could broadly be classified into 2 groups: genes encoding membrane and transport proteins; and those involved in metabolic pathways in the cell (Data S3).

In regard to antibiotic resistance and selection, the most interesting gene found to be under positive selection was porB (Smith, Maynard Smith & Spratt, 1995), which has been shown to be involved in mechanisms of resistance to penicillins, macrolides, cephalosporins and tetracyclines (Unemo & Shafer, 2014). porB exhibited signals of selection in subgroups 1, 3 and 5—the groups that harbored most of the antibiotic resistant strains in our sample set (Fig. 5). comA, which encodes a membrane protein necessary for competence of N. gonorrhoeae, was also found to be under selection in a handful of strains that make up subgroup 5. This finding is of interest in regards to the potential for DNA uptake in these strains, since they appear to be primarily DNA donors, rather than recipients in genetic exchanges (Fig. 2). Other genes putatively under selection included a stress response gene, a gene encoding a chaperone protein of the HscA family and a number of proteins: ribosyl transferase, RNA polymerase, and an arsR family transcriptional regulator, which were all linked to gene expression. Genes under positive selection in subgroup 3 were also mainly involved in gene expression or DNA metabolism, including DNA helicase and tRNA pseudo uridine synthase. Most of the genes with known functions, identified to be under positive selection in subgroup 1 were either membrane-associated or transport proteins.

Analysis of known genetic predictors for AMR phenotypes

A substantial amount of research effort over the past 10 years has been devoted to understanding the genetic basis of drug resistance in N. gonorrhoeae (Garvin et al., 2008; Unemo & Shafer, 2011; Veal, Nicholas & Shafer, 2002; World Health Organization (WHO), 2012). Since there is a an increasing interest in the direct attribution of resistance phenotypes based on genome sequencing, we attempted to ascertain how knowledge of existing variants could be applied to the N. gonorrhoeae genomes in this study. We searched for variants known to underlie resistance to 3 antibiotics classes within our study (Table 3). In terms of subgroup distribution, tetracycline resistance was found in each of the 5 population subgroups, azithromycin resistance was present in only 2 of the strains tested (SK36809 and MUNG8) and restricted to subgroup 2 (Fig. 5), and cefixime resistance was found in subgroup 1 and subgroup 2. We identified genes responsible for resistance to the drugs tested in this works using literature searches and the CARD antimicrobial resistance database (McArthur et al., 2013).

The tetM resistant determinant, which confers high-level resistance to tetracycline, is borne on plasmids and is transferred either through conjugation or transformation (Knapp et al., 1988; Morse et al., 1986; Turner, Gough & Leeming, 1999). It was found in only 5 of the 10 strains with high-level resistance to tetracyline (MIC equal or greater than 16 µg/ml). Strain SK1902, one of the 5 strains with the tetM determinant, had a significantly higher MIC (>256 µg/ml) than the rest (see attached Data S3). Other strains with reduced susceptibility or chrosomally-mediated resistance to tetracycline, i.e., without the tetM determinant do have other corresponding chromosomal mutations on one or more of the resistance loci: mtrR (including its promoter), penB, rpsJ. Only one strain (ATL0508) within the sample set exhibits resistance to tetracycline in the laboratory, without the presence of any of the known resistance determinants of the tetracycline resistance phenotype.

Different “mosaic” penA alleles are thought to have developed from recombination with portions of DNA transferred horizontally from commensal Neisseria and/or N. meningitidis and underly decreased susceptibility or resistance to cephalosporins by preventing their binding action on the encoded mosaic PBP2 (Ameyama et al., 2002). The mosaic penA XXXIV (Ohnishi et al., 2011; Grad et al., 2014; Unemo & Shafer, 2014) had the best positive predictive value of all the known resistance determinants we searched for within our dataset, being present in 6/7 of the strains resistant to cefixime. This result echoed the observations made by Grad et al. (2014) in their epidemiologic study of N. gonorrhoeae strains. The other loci (i.e., mutations in the mtrR, mtrCDE operon promoter region and penB gene) also proven to enhance the MICs of cephalosporins (Unemo & Shafer, 2011; Warner, Shafer & Jerse, 2008) did not have a similar predictive property within strains in our data set. These variants were seen in 2 out of 7 and 3 out of 7 cefixime resistant strains, respectively. MUNG17 is the only strain in the sample set that has an elevated MIC (0.38 µg/mL) to cefixime that we could not find any of the known resistance determinants within its genome sequence (see attached Data S1).

Resistance to azithromycin can be mediated by mutations in the previously mentioned penB and mtr operon genes as well as mutations found in the 4 different alleles of the 23S rRNA gene that inhibits protein synthesis (Chisholm et al., 2009; Palmer, Young & Winter, 2008; Starnino, Stefanelli & Neisseria gonorrhoeae Italian Study Group, 2009). The 23S rRNA mutation allele was found in one (SK36809) of the two strains with the azithromycin resistance phenotype. The other azithromycin resistant strain, MUNG8, did not have the 23S rRNA resistance determinant or any of the other known mutations in the mtrR or penB loci (see attached Data S1).