A genomic overview of the population structure of Salmonella

1S. enterica has long been subdivided by differential antibody reactions into what are now called serovars. The use of specific antibodies that could identify distinct cell-surface antigens within lipopolysaccharide and flagella have led to the distinction of >2,500 serovars that differ in their antigenic formulas. Some of these serovars cause invasive diseases of specific hosts, including the human-specific typhoidal serovars such as Typhi, Paratyphi A, and Paratyphi C, whereas others cause gastroenteritis of varying severity in multiple hosts (see Fig 1 in). Salmonella is also subdivided taxonomically into S. enterica, which contains multiple subspecies, and a separate species, S. bongori.

2In 2012, Achtman et al. described multilocus sequence typing (MLST) (Box 1) of S. enterica subspecies enterica on the basis of sequences from seven housekeeping gene fragments (legacy MLST). That publication also defined 138 genetic populations called eBurst groups (eBG1 to eBG138), which consisted of single-linkage clusters of sequence types (STs) that differed sequentially by only one of the seven sequences (single-locus variants). Most isolates within each eBG shared the same serological type, i.e., they were of the same serovar. However, many serovars were associated with multiple genetically unrelated eBGs, which indicates they are polyphyletic (derived from multiple independent ancestors). That publication also posed the question, ‘Will these assignments remain stable as 10,000’s of additional isolates are investigated by MLST and genomic analyses?’ The data summarised here demonstrate that these assignments are indeed extremely stable after comparing >100,000 genomes.

3MLST (multilocus sequence typing): MLST was designed to mimic multilocus enzyme electrophoresis for its ability to reveal population structures and to reduce the impact of homologous recombination on genetic distances, which resulted because each recombination event tends to introduce multiple clustered SNPs.

4SNP (single nucleotide polymorphism): a measure of genetic diversity and a common result of mutations.

5ST (sequence type): a combination of numbers, each denoting a unique sequence variant at a particular gene locus.

6Legacy MLST: MLST based on seven housekeeping genes. First described for Salmonella on the basis of serovar Typhi by Kidgell et al. in 2002. Extended to all of Salmonella and illustrated for S. enterica subspecies I in 2012 by Achtman et al..

7eBG (eBurst group): a group of related STs based on single-linkage clustering, in which the distance between nodes is only one allele. EnteroBase automatically adds new related STs to the most similar eBG unless they are one or two alleles distant from a second eBG.

8rMLST: MLST based on 53 genes encoding ribosomal proteins. For Salmonella, only 51 genes are used because 2/53 have paralogous sequences elsewhere in the genome.

9rST: ST based on rMLST.

10reBG: eBG in which single-linkage clustering is used to identify chains or clusters of rSTs that differ by ≤2 alleles. EnteroBase automatically adds new rSTs to the most similar reBG unless they are ≤5 alleles away from a second reBG.

11cgMLST (core genome MLST): MLST based on the core genome. In the case of Salmonella, a soft core of 3,002 genes was selected as a subset from a whole genome MLST set.

12wgMLST (whole genome MLST): MLST based on a nonredundant set of genes across a species, similar to a ‘pan-genome’.

13Due to the polyphyletic nature of serovars and an error rate of serotyping of at least several percent, Achtman et al. also recommended ‘phasing out the routine use of serovars, and replacing it with a classification that is based on population genetic groupings such as eBGs and STs’. This recommendation has largely been ignored. Instead, reference laboratories continue to use serovar designations, although several now simply predict the most likely serovar from a genomic sequence on the basis of MLST. Alternatively, serovars are predicted from genomes with the help of software such as Sistr, which uses MLST-based predictions as a tiebreaker when genomic sequences are not sufficiently discriminatory.

14Finally, Achtman et al. stated, ‘the evaluation of SNPs and genomic sequences from WGS takes much more time than the evaluation of paired traces from seven gene fragments… The S. enterica MLST database will probably contain data for >10,000 isolates in the near future, as do three other MLST databases today, whereas it would currently be difficult to extract information with comparable certainty from that many genomes’. These statements were correct at the time, albeit wordy. However, those interpretations and predictions were nevertheless pessimistic, backward-looking, and wrong. In particular, the number of Salmonella genomic sequences has been increasing exponentially since 2012 (Fig 1A), due in large part to the adoption of high throughput genomic sequencing of Salmonella and other pathogens by the Wellcome Trust Sanger Centre, Public Health England, the FDA, and PulseNet International. Academics are also beginning to sequence larger numbers of bacterial genomes. However, most genomic sequences have been deposited as sets of short reads, and >100,000 short-read sets from Salmonella are currently available. More and more are being deposited in the public domain as genomic assemblies, but their numbers are still far fewer than the numbers in the short-read archives.

15Talented bioinformaticians have analysed population structures or epidemiological transmission chains within a single eBG with the help of phylogenetic trees based on single nucleotide polymorphisms (SNPs) from hundreds or even thousands of Salmonella genomes. However, this approach is time consuming and resource intensive and does not scale adequately to deal with the entire diversity of Salmonella that is represented by >100,000 genomes. Other current approaches are also problematic. PulseNet International has recommended the use of genotypes based on whole genome MLST (wgMLST), but much of the genomic diversity among closely related strains represents variability in selfish mobile genetic elements, which can distort phylogenetic relationships (see below). An alternative approach is to use a core genome MLST (cgMLST) scheme. However, the strict core genome of 361 Salmonella genomes, in which all genes are present and intact in all genomes, consisted of only 360 genes.

16We have pursued still another alternative, a soft core genome: in Salmonella, this consists of 3,002 genes (Table 1) that were found to be present in ≥98%, intact in ≥94%, and of unexceptional diversity in 3,144 representative Salmonella genomes. Publicly accessible websites based on the software framework BigsDB present cgMLST schemes and their alleles for a variety of bacterial species, but not Salmonella, Escherichia, or Yersinia. We therefore developed a website, EnteroBase, which scours all short-read archives for sets of Illumina short reads from Salmonella, Escherichia, Yersinia, Moraxella, and Clostridiodes and supports uploading short reads by registered users. It assembles and polishes genomic contigs from the short reads within 2 hours and calculates MLST assignments from those genomes at the levels of legacy MLST, ribosomal gene MLST (rMLST), soft cgMLST, and wgMLST (Table 1; Box 1). EnteroBase also supports calling SNPs from up to 1,000 genomes against a reference genome as well as the graphic evaluation of genetic relationships between entries. Here, we re-examine the population structure of Salmonella based on the genomic contents of more than 110,000 Salmonella genomes in EnteroBase (Fig 1B).

17We have recently terminated our MLST website, which supported legacy MLST based on Applied Biosystems (ABI) sequencing. This decision was triggered because rMLST and cgMLST from genomic assemblies were much more informative than legacy MLST, and only a few new entries based on ABI sequencing were being submitted in recent years (Fig 1A). All of the historical MLST data are being maintained within EnteroBase in order to allow comparisons between data based on ABI sequences and data based on genomic assemblies.

18As of 1 November 2017, 651 common Salmonella legacy STs had been identified by both genomic and ABI sequencing (Fig 1B). A further 2,261 less common STs were only identified by genomic sequences, and 1,017 STs are currently only defined by ABI sequences (Fig 1B). The lack of genomic support for the 1,017 STs partially reflects a strong selection bias on the strains submitted to genomic sequencing, due to a focus on human disease in Europe or North America or driven by concerns about food safety in those areas. The STs defined by ABI sequences were less biased because they represented a global effort by academics plus reference laboratories and originated from a more diverse set of sources and environments. Indeed, over the last two years genomic sequencing has identified several hundred STs that were previously only known from ABI sequences, and we anticipate that over the next few years, additional historical STs will be confirmed by genomic data. However, other STs may never be confirmed because some ABI sequence-based STs represent ABI sequencing mistakes rather than rare STs. In recent years, the MLST website curators were rejecting about 50% of supposed new alleles and new STs based on ABI sequencing that were submitted by users. Similarly, we have identified multiple false allele calls produced by our own robotic and bioinformatic pipelines. Due to this experience, EnteroBase only accepts draft genomes that were generated by its own pipeline and that meet minimal criteria for contig size (N50 ≥ 20 KB), genomic size (≥4 MB), and unambiguous nucleotide calls (≥97%). Of such draft genomes, 99.3% yield sequences for all seven genes in legacy MLST. We therefore distrust the validity of the 1,017 legacy STs that are currently based exclusively on ABI sequencing (Fig 1B) until they are independently confirmed by genomic sequences.

19eBGs were defined by Achtman et al. as two or more legacy STs in a linear chain or cluster of single-locus variants, or at least 10 isolates within a singleton ST without close relatives. EnteroBase automatically adds new STs to existing eBGs when they match these criteria, unless they would result in the merging of eBGs. Additional eBGs are added at sporadic intervals with the aid of semiautomated curation, for a current total of 360 legacy eBGs. Almost all eBGs are relatively uniform for serovar. However, several percent of the publicly available metadata from the short-read archives include serovar assignments that are implausible according to their eBG assignments and in some cases even differ from published data. EnteroBase therefore provides independent serovar predictions based on the predominant serovar within each eBG and also provides the serovar predicted by Sistr. The two predictions are largely congruent. However, this approach enshrines the current composition of eBGs and ignores the existence of rare exceptions. We therefore tested the stability of eBG composition by comparing them with genetic clusters defined by a totally independent set of genes, namely rMLST, based on 51 genes that encode ribosomal proteins (Table 1).

20It might have been anticipated that 51 genes would have resulted in much higher resolution than seven genes, but the ribosomal genes are more conserved than housekeeping genes. As a result, legacy MLST of the genomes in EnteroBase defined 2,912 STs (Fig 1) and rMLST defined less than twice as many, namely 5,454 rSTs (Table 1). EnteroBase also defines reBGs as clusters of rSTs that differ in stepwise fashion by either one or two alleles. Almost all of the reBG assignments were consistent between legacy MLST and rMLST (Fig 2), with an Adjusted Rand similarity index of 0.992 for the eBG and reBG assignments of STs and rSTs among 101,254 genomes (Table 2). In general, reBGs were also consistent with serovar assignments. For example, most eBG14 isolates are serovar Saintpaul, but several are serovar Haifa, which differs at the H1 flagellar antigen. eBG14 splits into reBG14.1 and reBG14.2, both of which are pure Saintpaul, as well as reBG14.3, which contains all the Haifa variants plus some Saintpaul strains. This consistent assignment to eBGs and reBGs by independent genes provide very strong evidence that almost all eBGs/reBGs represent natural genetic populations within all of Salmonella and can be used to predict serovar with fewer mistakes than is currently possible with serotyping.

21We note that five exceptional sets of comparisons between eBG and reBG composition yielded conflicting clustering (Table 3). A closer examination of these five sets of conflicts with cgMLST identified complicated clustering patterns that are suggestive of the partial disruption of individual genetic populations by homologous recombination.

22Phylogenetic analyses of >100,000 genomes are currently very difficult, and it would be simpler to perform phylogenetic analyses on a representative subset that is more tractable, for example, representatives of all the diversity revealed by rMLST. Fig 3 shows a maximum-likelihood tree of core SNPs from one representative genome from each of the 297 reBGs in subspecies enterica plus one representative from each of 593 other rSTs in S. enterica and 36 in S. bongori. The results clearly delineate S. bongori as well as all seven known subspecies within S. enterica. They also indicate that at least three more undefined subspecies exist within Salmonella, which we have temporarily named novel subspecies A, B, and C, pending confirmation by standard taxonomical techniques.

23Elsewhere, in a second analysis dedicated to subsp. enterica itself, we calculated a maximum-likelihood tree from 2,964 genomes, one from each of the rSTs found within the first 50,000 assembled genomes. That tree revealed the existence of 48 ‘super-lineages’ within subspecies enterica, a higher level of population structure than that of eBGs. One of the super-lineages, the Para C Lineage, included reBGs 6 (Choleraesuis; a swine pathogen), 20.1 (Paratyphi C; human specific), and 20.2 (Typhisuis; pigs) as well as the single isolate in rST27817 (Lomita; a very rare serovar). That analysis showed that Paratyphi C was implicated in historical invasive human disease in Europe 800 years ago and subsequently in Mexico in 1545, and dated the time to the most recent common ancestor of reBGs 6, 20.1, and 20.2 to about 4,000 years ago. This analysis warrants further investigations of the properties of the other super-lineages.

24The literature consists of multiple conflicting claims regarding the optimal methods for pursuing fine-scale epidemiology and transmission chains. PulseNet International now recommends wgMLST as the preferred tool for fine-scale epidemiology, replacing pulsed field gel electrophoresis (PFGE), which was previously its primary tool for this purpose. EnteroBase implements a wgMLST scheme for Salmonella based on 21,065 orthologs from 537 complete genomes. This scheme was invaluable for defining the soft core cgMLST genome of 3,002 genes that EnteroBase uses. However, we do not recommend wgMLST for fine-scale epidemiology for several reasons. Firstly, Salmonella has an open pan-genome, as do many other recombining bacterial genera, and the number of orthologs in an open pan-genome grows with the number of sequenced genomes. A growing pan-genome results in a continuously growing wgMLST scheme, which in turn requires continuous calling of all genotypes as the scheme grows. Continuous reannotation of databases that already include >100,000 genomes would impose very high CPU requirements simply for maintenance. Alternatively, it might be sufficient to call genotypes on the basis of a frozen scheme, as is the case with EnteroBase. However, a frozen wgMLST scheme would lack any newly imported genes encoding antimicrobial resistance and other examples of horizontal gene transfer that were not included in the starter set of genomes from which the scheme was developed.

25The accessory genome is extremely volatile due to the repeated acquisition and loss of mobile, selfish DNA within individual strains. Thus, the primary phylogenetic signals in a wgMLST scheme are related to presence or absence of genes, and their evolutionary dynamics are currently much less well understood than are those of nucleotide diversification. To illustrate this phenomenon, we revisit serovar Agona, which caused a fast food-associated outbreak of human gastroenteritis in Ireland and the United Kingdom in 2008. The outbreak was caused by inadequately sterilised meat products due to a defective autoclave in a meat processing factory. PFGE analyses falsely indicated close relationships between the Agona from the food outbreak and other Agona from domesticated animals in Ireland in 2005. Furthermore, PFGE falsely indicated that Agona from the wastewater in the food plant were less closely related to the outbreak strain and that isolates from humans in Ireland in 2005 were only distantly related to the strains from domesticated animals. However, the PFGE patterns reflected the independent sequential acquisition by the two groups of strains of two bacterial prophages rather than changes in core restriction sites. According to core genome SNPs, the outbreak strains and the strains from the factory’s wastewater (group D) were very closely related but distant from the 2005 isolates from domesticated animals and from humans (group A) (Fig 4A). Fig 4 also shows the results of analysis of all 1,082 Agona genomes in EnteroBase by cgMLST and wgMLST. At first glance, both phylograms look similar and are dominated by the radiation of large numbers of independent branches from a common root. However, groups A and D are totally distinct according to cgMLST (Fig 4B), similar to the core SNP phylogeny (Fig 4A), whereas they represent two related subbranches according to wgMLST (Fig 4C). We therefore recommend investigating the utility of cgMLST instead of wgMLST for epidemiological purposes because it can identify very closely related Salmonella from a clonal outbreak as well as slightly more distantly related isolates that might represent longer term chains of transmission.

26In earlier publications, SNP typing was recommended for recognising outbreaks of salmonellosis by Public Health England and GenomeTrakr. SNP typing can be extremely powerful within an eBG, but the current existence of >20,000 genomes of each of Typhimurium or Enteritidis in EnteroBase render such an approach computationally difficult. Our experience with cgMLST until now is that it approaches core SNP typing in sensitivity. Hence, we also recommend investigating the utility of cgMLST rather than SNP typing for tracking outbreaks.

27We note that the graphic tool GrapeTree (Fig 2) was designed to allow the depiction of relationships based on character data such as MLST from up to 100,000 genomes and visualise the cgMLST genetic distances from all of the current Enteritidis or Typhimurium isolates in EnteroBase within minutes. GrapeTree has been implemented as a visualisation tool for EnteroBase or BigsDB (personal communication, K. Jolley) and can also be used in stand-alone mode. However, although GrapeTree can also represent SNPs (Figs 3 and 4), it is too slow to visualise all core SNPs from 20,000 genomes.

28Although we recommend investigating the utility of cgMLST for epidemiological purposes, we note that cgMLST is more demanding for analyses than are legacy MLST or rMLST, because cgMLST STs often contain low levels of missing data. Draft genomes based on assemblies from short reads contain multiple contigs, i.e., multiple genes are incomplete and cannot be used for calling MLST alleles. And the missing alleles can differ between independent draft genomes of the same organism. As a result, STs of natural isolates may have identical information content but differing missing data and will be assigned differing ST designations. This problem is not generally appreciated. However, EnteroBase includes approaches that minimise the missing data problem. Entries that only differ through missing data are collapsed into a single node by GrapeTree, and missing data are ignored for the calculations of genetic distances. We will soon release a form of hierarchical clustering within EnteroBase that ignores missing data and allows the recognition of natural groupings of isolates at any desired level of resolution within the cgMLST scheme, ranging from the very few loci/SNPs that distinguish isolates during a chain of transmission through to the large number that distinguish full species.

29In summary, eBGs represent natural populations, and we continue to recommend that they should replace serotyping. The data also indicate that multiple additional subspecies of Salmonella may exist, and we recommend that epidemiologists test the abilities of cgMLST for tracing chains of disease transmission.

30Until now, it has been necessary to use a plethora of techniques to work with the genetic diversity of Salmonella. EnteroBase now provides a one-stop solution to analyses at all scales, ranging from the genus level (Fig 3) down to epidemiological tracing (Fig 4). EnteroBase also provides comparable opportunities for the other genera it supports, which are currently Escherichia, Yersinia, and Clostridiodes. Within these genera, EnteroBase provides users with the opportunity to quickly identify the genetic relatives of the bacteria whose genomes they have recently sequenced within the framework of the currently sequenced global diversity at previously unknown scales (>100,000 Salmonella genomes; >50,000 Escherichia genomes). Finally, as demonstrated here, EnteroBase supports teamwork through the facile sharing of work spaces and trees with other selected users and with the general public.

31EnteroBase: http://enterobase.warwick.ac.uk.