Independent evolution of the core and accessory gene sets in the genus Neisseria: insights gained from the genome of Neisseria lactamica isolate 020-06

Genealogical relationships among species

Relationships among N. meningitidis, N. lactamica and N. gonorrhoeae were investigated by an analysis of nucleotide sequences from 19 housekeeping gene loci. The loci included were those used for Neisseria MLST (abcZ, adk, aroE, fumC, gdh, pdhC, and pgm) [29, 33] supplemented with 12 additional loci (aspA, carB, dhpS, glnA, gpm, pilA, pip, ppk, pykA, rpiA, serC, talA) with alleles generated as described previously [34]. The analysis included data from 52 N. lactamica isolates, comprising 46 unique seven locus sequence types (STs) sampled from a population of 271 isolates used to analyse diversity in N. lactamica[30], 20 isolates representing 20 unique STs (one ST from each clonal complex and two STs not belonging to clonal complexes) from the 107 N. meningitidis isolates originally used to define MLST [29], and subsequently analysed at additional loci [34], and seven complete Neisseria genomes downloaded from publicly available databases. These genomes include: meningococcal isolates Z2491 (serogroup A, ST-4 complex) [35], MC58 (serogroup B, ST-32 complex) [36], FAM18 (serogroup C, ST-11 complex) [37], 053442 (serogroup C, ST-4821 complex) [38], α14 [13], and gonococcal isolates FA1090 (The University of Oklahoma, U.S.A), and NCCP11945 [39]. All isolates used (Additional file 1) have been deposited in the Neisseria PubMLST database: http://pubmlst.org/neisseria/, which is searchable by publication record.

Genealogies were constructed with CLONAL FRAME version 1.1 [34, 40] using 50000 burn-in iterations and 50000 Monte-Carlo Markov Chain iterations. The results from eight runs were combined to achieve maximum robustness and the 75% consensus tree was viewed using MEGA version 4.0 [41]. To add support to the species-specific clustering, neighbour-joining trees were drawn from the same data using MEGA version 4.0 and Split Decomposition was carried out using SplitsTree version 4.10 [42]. MEGA version 4.0 was used to determine p distances.

Determination of the complete genome sequence of N. lactamica 020-06

The genome was sequenced to an approximately 10-fold shotgun sequence, totaling 42232 end sequences, from pUC19 (with insert sizes of 1.4-2 kb; 2-2.8 kb and 3-3.3 kb respectively. 2.8-3.3 kb), and pMAQ1Sac_BstXI (with insert sizes of 5.5 kb-6; 9-10 and 10-12 kb respectively) genomic shotgun libraries using big-dye terminator chemistry on ABI3730 automated sequencers. End sequences from large insert BAC libraries in pBACe3.6_BamHI with an insert size 15-18 kb were used as a scaffold. This generated an approximate 1-fold coverage from 4670 end sequences). All repeat regions and gaps were bridged by read-pairs or end-sequenced polymerase chain reaction (PCR) products again sequenced with big dye terminator chemistry on ABI3730 capillary sequencers. The sequence was manipulated to the 'Finished' standard [43].

Annotation and genome comparison

An automatic gene prediction program (Glimmer3) [44] was used to identify coding regions. Putative orthologues were identified by reciprocal-best-match FASTA searches between N. lactamica 020-06 and meningococcal strain FAM18 amino acid sequences with cut-offs of 80% sequence length and 30% identity. The complete genome was annotated manually using the genome viewer Artemis [45]. The N. lactamica genome sequence was compared to genomes from N. meningitidis (Z2491, MC58, FAM18, 053442 and α14) and N. gonorrhoeae (FA1090, NCCP11945) using the Artemis comparison tool (ACT) [46]. The Dot plot figure was generated using MUMmer version 3.22 [47] and indicates matching sequences, with codirectional and reversed regions of synteny shown in red and blue, respectively. The figure is a plot for multiple query sequences (other Neisseria) and one reference sequence (N. lactamica) where each reference/query comparison gives its dotplot. The genome sequences were aligned to start/finish at the origin of replication. The comparative analysis of the core gene dataset was carried out using an in-house pipeline of all-against-all best reciprocal FASTA searches that outputs the respective percentage (pair-wise) identity. The analysis and annotation of repeat families in N. lactamica, was carried out implementing a hidden Markov model based methodology, as previously described [37].

Genealogical relationships of N. lactamica, N. meningitidis, and N. gonorrheoae

The nucleotide sequences of 19 housekeeping genes from 52 N. lactamica, 25 meningococcal and two gonococcal isolates were used to generate a genealogy with the CLONAL FRAME algorithm [40]. ClonalFrame uses a statistical algorithm that infers clonal relationships while taking into account the effect of homologous recombination. This algorithm clustered the isolates into three groups, each corresponding to one of the microbiological species (Figure 1). These groups were also evident from other clustering algorithms, including neighbour joining trees and split decomposition (Additional file 2). Neighbour joining trees for each of the individual loci mostly replicated the same groups with some inconsistencies (Additional files 3, 4 and 5); for example, the genealogy for the dhps locus did not show species-specific clustering for N. meningitidis and N. lactamica, which might suggest horizontal genetic exchange between species. As dhps encodes dihydropteroate synthase, and sulphonamide resistance is mediated by altered forms of this enzyme [48], positive selection may promote recombination at this locus. However, no dhps alleles were shared among the three species.

The nucleotide sequence diversity for the 19 loci among the species groups as assessed by pairwise p-distances was similar, ranging from a p-distance of 0.053 between the meningococci and the gonococci to a p-distance of 0.073 between N. lactamica and N. meningitidis. The diversity within each of the groups, measured by the same means, was also similar among the N. lactamica isolates (average p-distance of 0.026) and the meningococcal isolates (average p-distance of 0.032) and appreciably lower between the two gonococcal isolates (p-distance of 0.06). These data support the conclusions reached from seven-locus multilocus sequence typing (MLST) and multilocus enzyme electrophoresis studies that the three microbiological species represent valid functional and taxonomic groups [3, 28] notwithstanding their close genetic similarity [2], which is consistent with recent divergence.

For the isolates included in this analysis, the number of fixed nucleotide differences greatly exceeded the number of shared polymorphisms in comparisons between the gonococci and meningococci (147 fixed differences, 33 shared polymorphisms) and the gonococci and N. lactamica isolates (316 fixed differences, 11 shared polymorphisms). This contrasted with comparisons of the meningococcal and N. lactamica isolates, where the number of fixed differences (118) was much smaller than the number of shared polymorphisms (435). Taken together these data indicate that all three species groups share an ancestral population, with the gonococcus likely to be descended from a single member of that population. N. gonorrhoeae does not exchange genetic material with other Neisseria frequently but, unlike other pathogens that have emerged by a similar process [49, 50], genetic material is regularly exchanged among gonococci such that their populations do not have a clonal structure [51, 52]. Present day meningococcal and N. lactamica populations, on the other hand, have diverged from the ancestral population without undergoing a severe bottleneck and now represent distinct species groups, with the nucleotide sequence diversity in their housekeeping genes exhibiting extensive shared ancestry.

In the current analysis there were only two cases of alleles shared among the different species groups, both at the pykA locus. The pykA-7 allele was present in 5 meningococcal isolates and one N. lactamica isolate and the pyKA-24 allele was present in one meningococcal isolate and one N. lactamica isolate. In the N. lactamica 020-06 genome this locus is immediately adjacent to an orthologue of gene NMB0088 in the meningococcal MC 58 genome, which encodes a putative immunogenic outer membrane protein orthologous to the FadL potential vaccine component of Escherichia coli[53]. This observation complements previously published evidence for the hitch-hiking of housekeeping alleles with antigen alleles promoting genetic exchange among meningococci and N. lactamica[54], although such hybrids as do occasionally arise appear to be less fit and are purged by purifying selection [55]. In conclusion, while genetic exchange between meningococci and N. lactamica has been reported for a number of genes under positive selection [55, 56], the two groups remain distinct as a result of limited genetic exchange among housekeeping genes [3], which is consistent with the species status of these organisms. Frequent horizontal genetic exchange among groups would prevent speciation, or if it occurs after the speciation event, would lead to a 'despeciation' process, as has been proposed recently for the gastric pathogens Campylobacter jejuni and Campylobacter coli[57]. There is no convincing evidence of such a process occurring among the Neisseria.

As with meningococci, N. lactamica isolates are currently assigned to clonal complexes on the basis of their seven locus sequence types [30]http://neisseria.org/nm/typing/mlst/. These clonal complexes comprise a central genotype ST and STs that share identical alleles for at least four of the seven MLST loci in the central genotype. At the time of writing, six N. lactamica clonal complexes had been identified, each named after the central genotype of the complex (the ST-595, ST-613, ST-624, ST-640, ST-1494 and ST-1540 clonal complexes), five of which were represented in the genealogical analysis reported here. Members of the same clonal complexes were consistently clustered together in the 19 locus genealogy, indicating the general robustness of the assignment of clonal complexes on the basis of seven loci. The genealogy did, however, suggest that some of the genealogical groups were larger than the currently defined clonal complexes and also indicated the existence of several other clonal complexes, currently represented by one or a few sequence types (Figure 1). The star phylogeny of N. lactamica is redolent of that of the meningococcus [34] and consistent with both species exhibiting a clonal complex population structure, comprising groups of related isolates that cannot be genealogically linked to each other as a consequence of high rates of intra-species horizontal genetic exchange [33, 34, 58].

On the basis of these analyses there was no a priori reason to choose a member of one particular clonal complex over any other for genome sequencing, although it was considered desirable to choose a recent bacterial isolate that had not been maintained in the laboratory for long periods of time and one that represented a common and widespread genotype. For these reasons, isolate 020-06, obtained from a six-week-old child from the UK in 1997 was chosen as the most appropriate candidate for genome sequencing. This isolate has ST-640, the central genotype of a clonal complex from which members have been isolated in the UK [30] and four other countries (unpublished data). The MLST profiles for these isolates can be viewed at: http://pubmlst.org/neisseria/[59].

General features and gene order

Whole genome comparisons of gene order showed that, with the exception of an inversion of 404215 bp around the terminus of genome replication, the genome sequence of N. lactamica 020-06 was essentially co-linear with that of meningococcal isolate Z2491 (Figure 2A). The larger inversion was also evident in comparisons of the N. lactamica genome with the genome sequences of the remaining meningococcal genomes examined, those of isolates MC58, FAM18, α14, and 053442, one of the two gonococcal isolates, NCCP11945, but not the genome of isolate FA1090. Other than these changes in gene order, all of which maintained features such as GC skew (Additional file 6), there were no consistent changes in co-linearity among any of the Neisseria genomes examined (Figure 2B). Although the total number of completely assembled genomes available remains small, these comparisons suggest that there is a consensus genome order, described previously for N. meningitidis[37], for all three Neisseria species examined here. This is essentially the gene order seen in meningococcal isolates Z2491, α14, 053442, and FAM18. Most of the large-scale departures from this consensus are present in only one meningococcal genome examined so far and may be particular to that isolate alone or that isolate and its close relatives. Such rearrangements, frequently observed in bacteria [62] may well be biologically neutral or artifacts of growth in the disease state [63] or pure culture. One exception to this was the large rearrangement around the terminus that is present in the N. lactamica isolate and one gonococcal isolate (FA1090). The very close relationship of all gonococci (Figure 1), suggests that this inversion may have occurred more than once in the history of the genus, but more genome data are required to definitively answer this question. There is consequently no evidence to suggest that gross genome rearrangements have played a role in the emergence of pathogenicity in the Neisseria as there is no consistent pattern of genomic rearrangements among non-pathogenic and pathogenic Neisseria, as have been observed in some other groups containing bacterial pathogens [64], although a role of such rearrangements in the emergence of clonal complexes cannot be excluded on the basis of information currently available.

Gene content and core genome

Of the 2018 N. lactamica CDSs identified, 1447 (72%) were assigned to functional categories. Of the remaining CDSs, 421 (21%) were not assignable to a known function and 150 (7%) were categorized as conserved hypothetical proteins. Comparisons with the genomes of N. meningitidis Z2491, MC58, FAM18, 053442, α14, and N. gonorrhoeae genomes FA1090 and NCCP11945 showed a similar distribution of functions (Figure 3). Changes in gene content among the different genome sequences were confined to CDSs classified as involved in transport and metabolism. These were relatively overrepresented among those CDSs present in other Neisseria, but absent in N. lactamica. They were also underrepresented among those only present in the N. lactamica genome. CDSs of unknown function were overrepresented in the N. lactamica genome and underrepresented in the core genome. A total of 1190 (59%) of the N. lactamica CDSs were present in all the seven completed genomes available at the time of writing (Table 1) and these can be regarded as an estimate of the core genome for these three species (Additional file 7). This may actually be an underestimate of the core genome as this dataset was automatically produced and not manually curated. Therefore mis-annotations in any of the genomes would mean that the genes involved might automatically be excluded. The core genome of N. meningitidis alone has been estimated to consist of 1706 genes when the genomes of Z2491, MC58 and FAM18 were examined [8] and approximately 1300 genes when seven meningococcal genomes were examined [65]. The latter estimate is very similar to the core genome proposed for the three species here and, given that the number of core genes in purely meningococcal comparisons is likely to reduce as more genomes become available, suggests that the core genome, amounting to some 60% of the total genome, is common to all three species.

The draft genomes sequences of eight human Neisseria commensal species have been compared to eleven previously sequenced Neisseria genomes [66] and a Neisseria core genome of 896 genes was defined, consisting mainly of housekeeping genes. This analysis, however included Neisseria elongata, which is more distantly related to the meningococcus than the other Neisseria commensals analysed, and is an unusual member of the genus in that it is rod-shaped, in contrast to other Neisseria spp. which are diplococcic. This species also shares fewer core genes with the pathogenic species than the other commensals. Further studies, involving more commensal genomes will help define further the core Neisseria genome and may also alter the species definitions of organisms currently defined as Neisseria.

For the core genome defined here CDSs ranged from 36.4% to 100% amino acid identity. However, 1180 of these CDSs had greater than 70% amino acid identity among all the genomes examined. The ten CDSs that had less than 70% amino acid identity to the corresponding CDSs in isolate Z2491 included seven that were classified as transport/membrane proteins; one unknown protein; a 3-oxoacyl-(acyl carrier protein) synthase II protein; and putative ribonuclease BN. However, the lower identity was not consistent across the genomes. For example, the transferrin-binding protein B (TbpB) in Z2491 was greater than 70% identical to the corresponding proteins in the other genomes except for the genomes of the two serogroup C meningococci (FAM18 and 053442). Whereas the TbpB in isolate 053442 was 67.3% identical to the Z2491 TbpB, the FAM18 TbpB was only 41.7% identical. Previous work has determined that there are two tbpB isotypes, with isotype 1 solely identified among N. meningitidis isolates belonging to the ST-11 clonal complex [67], of which isolate FAM18 is a member. The CDS encoding the porin protein PorB also showed differing degrees of identity among the genomes, with the FAM18 porin having the least identity to Z2491 (61.8%). The PorB proteins are known to be diverse [68] with the FAM18 porin belonging to class PorB2 and the Z2491 porin belonging to class PorB3.

Pair-wise comparisons of the 1190 core CDS sequences from meningococcal isolate α14 against each of the remaining genomes generated a profile of relatedness that was concordant with the microbiological species groups (Figure 4). All five meningococcal isolates were identical at 72 loci, with many of these loci encoding ribosomal proteins (Additional file 7). By contrast, a total of 30 N. lactamica 020-06 and 29 gonococcal FA1090 CDSs were identical with their orthologues in Z2491. As the threshold of identity was lowered, these relationships were retained, with all four meningococcal genomes showing close relationships to the genome of isolate Z2491 and the N. lactamica and gonococcal genomes showing lower numbers of CDSs at each level of sequence identity, with all genomes converging at >70% identity for 1180 CDSs. These relationships were consistent with the CLONAL FRAME genealogy obtained from the 19 locus analysis (Figure 1) and provided no support for the idea that the acapsulate ST-53 meningococcal isolate α14 was more closely related to N. lactamica and N. gonorrhoeae than the encapsulated meningococci [13, 65]. Therefore, although 60% of the CDSs present in a given Neisseria genome appear to be common to all three species, and there is a high level of shared ancestry among meningococci and N. lactamica, the core genomes of each species are characterized by particular sequence poylmorphisms. Again, this is consistent with the results of the 19 locus analysis and provided more evidence for lack of frequent inter-species genetic exchange among the core genomes of these organisms.

The accessory genome

There were several genes not present in meningococci and gonococci which may have been lost by these organisms on speciation, including those encoding: a putative toxin-antitoxin system (RelEB) [70] and seven CDSs encoding hypothetical proteins not found in the pathogenic Neisseria upstream of these loci; a putative TonB-dependant receptor that is also present in a number of other Neisseria species (Neisseria cinerea ATCC14685, Neisseria sicca ATCC29256, Neisseria subflava NJ9703, Neisseria flavescens NRL30031, Neisseria mucosa ATCC25996 sequenced at the Genome Sequencing Center at the University of Washington and available at http://www.ncbi.nlm.nih.gov/). The N. lactamica genome appears to have two different copies of genes encoding putative L-lactate permeases. One (NLA16970), encodes an L-lactate permease similar to those found in N. gonorrhoeae and N. meningitidis and the other (NLA5610) is present only in N. lactamica 020-06 and other genomes of commensal Neisseria (N. lactamica ATCC 23970, N. cinerea ATCC14685, N. sicca ATCC29256 , N. mucosa ATCC25996, sequenced at The Genome Sequencing Center, University of Washington). Downstream of NLA5610 are genes that encode two putative TonB dependant receptors: NLA5590 which is present in pathogenic and commensal Neisseria, and NLA5600, which has only been found in N. lactamica 020-06 and commensals sequenced at The Genome Sequencing Center, University of Washington (N. lactamica N. cinerea, N. sicca, N. mucosa, N. subflava, and N. flavescens). These two genes may have been deleted from the meningococcal and gonococcal genomes, as the GC content of these sequences in N. lactamica 020-06 did not suggest horizontal genetic exchange.

Although N. lactamica is distinct from N. gonorrhoeae, there are some similarities between the two species. The most obvious is the lack of capsule in both N. lactamica and N. gonorrhoeae. In addition, the putative toxin-antitoxin system FitAB is present in N. lactamica and N. gonorrhoeae but not N. meningitidis. A gonococcal mutant that lacks fitAB grows normally extracellularly, but has an accelerated rate of intracellular replication with a concomitant increase in the rate at which this mutant traverses a monolayer of polarized epithelial cells [71]. As N. meningitidis lacks both FitA and FitB, intracellular replication and intracellular trafficking may be greater in meningococci. NLA13470 and NLA13480 (putatively azlC and azlD respectively) are also present in N. gonorrhoeae. In N. meningitidis, azlC is truncated and azlD is not present, suggesting this region has been disrupted in N. meningitidis. A list of genes present in meningococcal genomes Z2491, MC58, FAM18, and the gonococcal genome FA1090, but absent from N. lactamica 020-06 is included in Additional file 8.

Out of a total of 134 candidate meningococcal virulence genes proposed previously, 115 (86%) have been shown to be present in the genomes of three meningococcal isolates from asymptomatic carriage (α14, capsule null; α153, serogroup 29E; and α275, serogroup W-135) [13]. A total of 78 (58%) of these genes were present in the N. lactamica 020-06 genome (Additional file 9). The remaining 56 'N. meningitidis specific' genes were not found in all meningococcal genomes and some were present as putative pseudogenes [13]. Functional genes found in all the meningococcal genomes analysed to date but not found in N. lactamica 020-06 included: lgtA, opcB, porA, fetB2, frpA, frpC, iga1, sodC, gna1870, natC and nlp; however, gna1870 appeared to be present in other N. lactamica isolates (our unpublished data) and many of the other genes may yet be found in N. lactamica isolates that have not been sequenced. There is an iga2 gene present in N. lactamica (NLA19200); however, this gene has a ~708 bp insertion relative to the iga2 gene in the pathogenic Neisseria and bears no resemblance to the gene that encodes the IgA protease, which is implicated in pathogenesis.

These analyses confirm that the majority of the accessory genome is shared among all three Neisseria species, probably by horizontal genetic exchange, and that a 'pathogenome' of accessory genes, which explains differences in the pathogenic phenotypes among these organisms cannot readily be defined. Some of the shared CDSs have previously been described as 'putative virulence genes' but these are better thought of as genes that confer fitness benefits to transmission only in certain circumstances or genetic backgrounds and which affect the probability of causing invasive disease only incidentally [13]; indeed, many of these genes have now been identified in the complete genome of N. lactamica 020-06 and other N. lactamica isolates [7]. As meningococcal-specific sequences are thought to represent only about 5% of the genome of N. meningitidis[72], differences in virulence among isolates and species is likely to be due to subtle changes in their genetic organisation and in their alleles, making the determination of what constitutes a pathogenic phenotype complex for these organisms. This observation is consistent with the conclusions of a previous micro-array-based study, which suggested that a complex integrated network of genes, regulation and diversity of function in common genes could be responsible for the behavioural differences among N. lactamica, N. meningitidis and N. gonorrhoeae. Consequently, for these organisms it is not possible to attribute the emergence of an invasive phenotype to a single event, a single gene set, or loss or gain of a single function [7].