Phylogenomic Analysis Identifies Gene Gains That Define Salmonella enterica Subspecies I

Comparative methods for analyzing whole genome sequence (WGS) data enable us to assess the genetic information available for reconstructing the evolutionary history of pathogens. We used the comparative approach to determine diagnostic genes for Salmonella enterica subspecies I. S. enterica subsp. I strains are known to infect warm-blooded organisms regularly while its close relatives tend to infect only cold-blooded organisms. We found 71 genes gained by the common ancestor of Salmonella enterica subspecies I and not subsequently lost by any member of this subspecies sequenced to date. These genes included many putative functional phenotypes. Twenty-seven of these genes are found only in Salmonella enterica subspecies I; we designed primers to test these genes for use as diagnostic sequence targets and data mined the NCBI Sequence Read Archive (SRA) database for draft genomes which carried these genes. We found that the sequence specificity and variability of these amplicons can be used to detect and discriminate among 317 different serovars and strains of Salmonella enterica subspecies I.


Introduction
Recently, we applied whole genome phylogenetic analysis to the epidemiological trace-back of an outbreak of Salmonellosis [1]. However, analyses of this type can only give information about past outbreaks, and cannot prevent outbreaks from happening in the first place. In order to prevent outbreaks, we must be able to rapidly identify tainted foods before they come to market.
Some researchers have questioned whether it is possible to reconstruct an accurate evolutionary history of bacteria, given ongoing debates about the influence of horizontal gene transfer [2][3][4][5][6][7][8][9]. However, we believe that a phylogenetic analysis of whole genome sequence (WGS) data can solve these problems and account for HGT. In fact, using a whole genome tree of life hypothesis, we were recently able to corroborate the hypothesis that there is a vertical history of life for bacteria [8]. We expect these techniques will enable us to better understand the genomic evolutionary history of finer scale taxonomic classes of bacteria, including serovars of S. enterica subspecies I. As a step toward this goal, we have applied the comparative method of WGS phylogenetic analysis to discover diagnostic biomarkers [2] capable of identifying and discriminating among forms of Salmonella enterica (S. enterica). We suggest that phylogenetic analysis of WGS data can provide a solution to the problem of effective detection and identification of S. enterica serovars and some strains.
The Salmonella enterica subspecies I Salmonella infection is currently the most common foodborne illness in the United States (US), resulting in thousands of infections per year. These rates have not declined in over a decade, demonstrating the high fitness level of S. enterica. To reduce the human and financial costs of this pathogen, it is imperative to quickly and cheaply detect Salmonella contamination before it enters the food distribution system [10,11].
The current classification of the genus Salmonella divides it into two species: Salmonella bongori and Salmonella enterica [12]. Baumler [13] suggested that a gain of the genetic elements fim, the Salmonella Pathogenicity Island 1 (SPI1), and lpf, are responsible for the ability of this genus to invade intestinal epithelial cells. In the same study, Baumler [13] went on to postulate that the evolutionary transition from the common ancestor of S. bongori and S. enterica to S. enterica occurred, in part, by the acquisition of SPI1, and that the divergence of S. enterica subspecies I from the other subspecies is due to the acquisition of several genes by subspecies I, and loss of the lpf operon by subspecies II, III, IV, and VI. Later, Baumler et al. [14] developed the hypothesis that the complex lymphoid systems of mammals and some bird species drove the evolution of virulence among all of the members of subspecies I. Later research from the same group reported that invA dependent SPI1 is responsible for the ability of non-typhoidal Salmonella to enter gut lymphoid systems [15]. Several approaches have been used to classify the serovars within Salmonella enterica subspecies I and some of the perceived disagreements among researchers may be attributable to differences in methodology. For example, one recent study showed that gene presence-absence data from DNA microarray analyses produced an un-weighted pairwise-distance tree that clusters most serovars together; however, multi-locus-sequence-typing (MLST) analysis showed more variability [16]. One study aimed at classifying serovars within S. enterica subspecies I using WGS information concluded that there is little correspondence of serotype with evolutionary history [17], although this analysis did not address any possible HGT. Another analysis explored gene gains in different subspecies of S. enterica from a functional point of view, noting abundant recombination events between lineages [18]. Another recent analyses with draft and complete genome sequences using Ribosomal 16s and weighted gene presenceabsence matrices came to different conclusions based on the data type and weighting scheme used to correlate serotype and genomic evolutionary history [19]. An MLST and whole genome alignment analysis, using serotypes of both S. bongori and S. enterica that rooted the genus with S. enterica arizonae, found that serovars of S. enterica and S. bongori underwent HGT from other species [20]. Another Salmonella population genetics study, that sequenced 146 regions of 2 to 2.5 kb for 114 strains of Salmonella enterica, found there was significant homologous recombination in the species. Each of these analyses has provided a wealth of information that furthers our understanding of the evolutionary history and function of these important pathogens.
In the current project, we have used different draft genome sequences, in conjunction with complete genome sequences, to further test the evolutionary relationships within Salmonella enterica subspecies I to derive a better-corroborated history of these foodborne pathogens ( Table 1). As draft genome data are only able to describe gene sequences that are present in, but not those absent from, a genome, we focused our analyses on those genes that were present in all samples used in our phylogenetic analysis.

Results/Discussion
The Salmonella enterica subspecies I We used gene presence-absence data and the phylogenetic methods of Lienau et al. [21,22] as heuristic searches to empirically define the Salmonella enterica subspecies I homologous genes. Briefly, these searches define gene similarity thresholds and select the threshold resulting in the most resolved and consistent gene presence-absence phylogeny that also provides the most consistent character statements as measured by the combined corroboration metric (CCM) [21]. Our phylogenetic analysis and homology search showed that a similarity value of 70% yielded the most congruent and resolved gene presence-absence phylogeny.
We only considered open reading frames (ORFs) that are longer than 120 base pairs and matched at least one other ORF across more than 80% of the nucleotide sequence. The heuristic analysis identified 937 genes common to all organisms in our analysis; this is about half the number found by Jacobsen et al. [19], due to the stringency with which we defined our homologous groups. We used these commonly-held genes to create a multi-locus DNA sequence evolution phylogenetic matrix of 937 genes that contained 204,753 total characters. Out of this number, 82,147 characters were constant, 74,345 were variable characters and parsimony-uninformative, and 48,261 characters were parsimonyinformative. Parsimony analysis of this matrix produced 321 most parsimonious trees with a score of 201,886 and CI of 0.653; the strict consensus tree gave a normalized consensus fork component information index of 0.887. We rooted this strict consensus tree with the S. bongori genome and chose to depict this rooting as a polyphyletic relationship. This root placed S. enterica arizonae at the base of the S. enterica subspecies I lineage. We used broken branches, as seen at the top of Figure 1, to denote that S. bongori and S. enterica arizonae had long internal branch lengths compared to the S. enterica subspecies I, a result shared with Fookes, et al. [20].
The strict consensus tree postulates an evolutionary hypothesis classifying most operational taxonomic units (OTUs) that belong to the same O-antigen serovar as monophyletic groups ( Figure 1A).     Figure 1) and by S. Agona (shown in bold toward the bottom of the tree). It is tempting to speculate that the S. 1,4,5,12:i:-phenotype came into these taxa via one or more HGT events conferring genes responsible for generating the somatic antigens of S. enterica. This tree also resolves the S. enterica subspecies I relationships for 12 of 35 strains of S. Montevideo [1] and places S. Javiana and S. Schwarzengrund at the base of the monophyletic S. Montevideo group. S. Enteritidis is monophyletic and clusters with a monophyletic S. Galinarum and a monophyletic S. Dublin. In contrast to the findings of Achtman, et al. [16], S. Kentucky and S. Tennessee appear to be monophyletic. These discrepancies may be due to the different phylogenetic and/or sampling methods and isolates used in these analyses.

Genes gained at the base of the Salmonella enterica subspecies I
We then mapped the gene presence-absence data onto the core gene phylogeny to identify genes gained by S. enterica subspecies I. Previous researchers have identified from 216 to 249 genes gained by the S. enterica subspecies I. [13,[23][24][25] We examined all the genes gained by the common ancestor of all S. enterica subspecies I in our analysis. We identified 377 total genes gained. Of these, 71 genes were gained once and not lost by any members of S. enterica subspecies I included in this analysis ( Table 2, see methods).
We then identified ORFs at the nucleotide level that are both specific to and able to discriminate among the genetically known serovars of Salmonella enterica subspecies I included in this study. We used nucleotide blasts to compare members of these genes against the NCBI non-redundant (nr) nucleotide database. Thirty-one gene sequences had a BLASTn total score of greater than or equal to 200 vs. non-Salmonella organisms (Table 3). We also retrieved 13 significant hits using the default parameters of megaBLAST to S. enterica houtenae, (subspecies IV) of Salmonella enterica in the SRA http://www.ncbi.nlm.nih.gov/Traces/sra).
This left 27 genes completely specific to S. enterica subspecies I and 31 genes with homologs that appear closely related to other non-Salmonella species. The functional classes of these gained genes were varied, indicating that the evolutionary advantage of gene gains to S. enterica subspecies I is not limited to a specific biological pathway (Tables 2 and 3). The relatively equal number of gene gains with similarity to non-Salmonella species (n = 31), compared to those with only weak similarity to other organisms (n = 27), may indicate that recent HGT from distantly related species as well as duplication and divergence within Salmonella enterica subspecies I may play a major role in the evolution of S. enterica subspecies I.

Biomarkers to create improved detection tools for Salmonella enterica subspecies I
We also assessed the suitability of these gene sequences for use as biomarkers in a PCR and/or sequence detection system. The system was designed to quickly identify and/or discriminate among different serovars of S. enterica subspecies I. First, we tested the 27 S. enterica subspecies I-specific genes for utility as diagnostic biomarkers for S. enterica subspecies I. We designed PCR primers for all 27 biomarker candidates (Table 3, methods) and blasted all 27 sequences against members of S. enterica in the SRA [26]. These sequences showed significant similarity to 317 strains of S. enterica subspecies I, suggesting that these amplicons could be used as a gene presence-absence diagnostic test for S. enterica subspecies I. The primers were also tested in the laboratory on a limited set of isolates and amplicons were successfully generated as designed both for positive and negative expectations (data not provided).
We then tested the resolution capabilities of a PCR and sequencing-based test to discriminate among the 317 available  strains of S. enterica subspecies I by doing phylogenetic analysis on the predicted PCR fragment sequences. We made 27 alignments for 317 S. enterica subspecies I serotypes and concatenated them into a MLST super-matrix composed of 12055 characters, 862 of which were parsimony-informative. We analyzed the matrix using the parsimony ratchet strategy. This analysis yielded 99 most parsimonious trees with scores 3720; the majority rule bootstrap consensus is shown in Figure 1B and the strict consensus is shown in Figure S1. While this tree is based on less character information and contains nearly 5 times the number of taxa than the tree derived from the 937 gene analysis, it is resolved enough to discriminate among most of the different serotypes and strains of S enterica subspecies I included in our analysis ( Figure 1B). In order to test whether we could increase the resolution of this MLST approach, we also selected gene sequence fragments for the 13 genes determined to be significantly similar to S. enterica houtenae in the SRA megaBLAST test. We designed and tested PCR primers to amplify these genes (Table 3). We included this additional genetic data to the 27 gene alignments for 317 S. enterica subspecies I serotypes and concatenated them into a MLST supermatrix composed of 319 taxa (including 2 S. enterica houtenae strains). This matrix had 16221 characters, 1361 of which are parsimony-informative. Phylogenetic analysis yielded 256 most parsimonious trees with scores 5820, the strict consensus of which is pictured in Figure S2. This tree is rooted with S. enterica houtenae. It is resolved enough to discriminate among all of the different serotypes and many strains of S. enterica subspecies I tested. This gene set is an improvement to the smaller MLST analysis that did not include the genes also present in S. enterica houtenae ( Figure S2). These results indicate this set of PCR amplicons should be excellent candidate biomarkers for use in MLST and single nucleotide polymorphism (SNP) detection and diagnostic tools for S. enterica subspecies I serovars and some strains.
The first set of 27 genes and the PCR amplicons derived from them are sufficient for a gene presence-absence diagnostic test for S. enterica subspecies I. This set of amplicons can also be used to reliably determine most serotypes and some strains of S. enterica subspecies I in a PCR/SNP detection system. However, because some isolates did not resolve all serotypes into monophyletic clades when using the 27 gene set, it would be better to use the 40 gene set which includes the out-group S.enterica houtenae as a SNP-based MLST diagnostic system, since that method fully resolves all serotypes and many strains of S enterica subspecies I tested ( Figure S2).

Conclusion
Reconstructing the evolutionary history among lineages provides an approach for both identifying the qualities that make pathogens dangerous and detecting those organisms in settings where they pose potential threats to human health. We devised a scheme that would generate an evolutionary hypotheses to test which genes were unique among Salmonella to S. enterica subspecies I. Using the BLAST algorithm, we tested the origin of these gene sequences against the non-redundant nucleotide database and found that some genes were very similar to distantly related organisms, and that others were only weakly similar to distantly related species. We used these unique gene sequences to generate diagnostic biomarkers that can detect the presence and determine the serotype of S. enterica subspecies I. This method of identifying diagnostic characters for a clade of organisms provides a future framework to generate and test hypotheses about genetic variations that may be correlated with disease phenotypes.

Genome Sequences
We used Roche 454 sequencing technology to sequence 34 new Salmonella enterica draft genomes from various sources (Table 1). We assembled the shotgun sequenced genomes using Glimmer (http://cbcb.umd.edu/software/glimmer/) and had them annotated using the National Center for Biotechnology Information (NCBI) Prokaryotic Genome Automatic Annotation Pipeline (PGAAP) [27]. We also downloaded 30 publicly-available Salmonella species genome sequences from NCBI's GenBank, including Salmonella bongori (Accession NC_015761). This yielded a total of 71 Salmonella enterica genomes and one Salmonella bongori genome. The sequences reported in this paper have been deposited in the GenBank and SRA databases, with accession numbers listed in Table 1.

Phylogenetic analysis
An empirical homology cluster search was performed per Lienau et al. [8] to determine the similarity value to generated the gene clusters that yielded the most congruent and best-resolved gene presence-absence phylogenetic tree. We tested similarity value thresholds for gene clustering of 60%, 70%, 80% and 90% at length limit of 120 bp and minimum match length of 80% using megaBLAST, as implemented in a computer program called PathGenome, currently being developed by the FDA and the Food and Environmental Research Agency (FERA). We generated gene presence-absence matrices for each of the similarity values tested and performed tree searches using the Phylogenetic Analysis using Parsimony and Other Methods (PAUP*) 4.10 b portable version [28] with a ratchet search of 9 iterations: 3 iterations each at perturbations of 15% 17% and 21%, respectively, using command files generated by Parsimony Analyses using PAUP* (PRAP) [29]. After establishing which genes in our study resulted in the gene presence-absence tree with the highest CCM score, we then aligned each of those gene sequences using MUSCLE and constructed a multi-locus DNA sequence evolution matrix [30]. We searched for the most optimal tree using the parsimony ratchet searches as described above. All characters were equally weighted. We did a bootstrap analysis on all phylogenetic matrices using PAUP*4.10 b portable version [28] at 100 replicates, holding a maximum of 1000 trees per replicate.

Biomarker identification
We used the method of Lienau, 2012 [2] to identify likely candidate genes for use as diagnostic biomarkers. We defined the node of interest as the node that led to all of the Salmonella enterica subspecies I and then used accelerated transformation of parsimony character reconstruction to identify genes gained at that node. We then selected 71 of these gained genes that showed perfect consistency with the phylogenetic hypothesis as measured by the consistency index of Kluge and Farris [31]. We manually checked these genes for presence in all members of the ingroup (S.enterica subspecies I) and absence in the outgroup (S. bongori and S. enterica arizonae) to rule out symplesiomorphies. All 71 genes met these criteria. We then checked other organisms to see whether these genes were present (see next section).
Blast to identify potential false positive markers and primer design We took an example sequence from each of the 71 genes and blasted them against the nr nucleotide database at NCBI. Using a lower bound of 70% identity, we separate 40 genes with positive hits to only Salmonella, leaving aside the 31 genes that had hits to other organisms. We designed and tested primers to the conserved regions of the 40 Salmonella-only genes using MacVector with Assembler 11.0.2 via Primer Design (Primer3). We set the ideal amplicon length to between 500 and 600 base pairs. (Table 3). Using phylogenetic analysis ( Figure S2), we further tested the potential utility of these 40 predicted sequence amplicons to discriminate among 317 serotypes of S. enterica subspecies I. We extracted the sequence for 27 genes on 317 strains of S. enterica subspecies I using phylogenetic analysis (Figure 1 B, Figure S1). All 40 primers also were tested in the laboratory with a limited set of isolates. Expected amplicons were generated (or not for negative controls and E. coli) for these limited experiments (data not provided).
All NCBI Salmonella genomes are linked to NCBI Sequence Read Archive (SRA) files, and accession numbers. Cultures included in this study are also available upon request. Please direct any queries for isolates to our strain curator Dwayne Roberson, at Dwayne.Roberson@fda.hhs.gov. Figure S1 Majority Rule Bootstrap consensus tree of 317 S. enterica subspecies I serotypes MLST supermatrix composed of 12055 characters from 27 alignments derived from predicted PCR products made from the S. enterica subspecies 1 specific biomarker sequences in Table 3.

Supporting Information
(TIFF) Figure S2 Majority Rule Bootstrap consensus tree of a MLST super-matrix composed of 319 taxa (including 2 S. enterica houtenae strains) 16221 characters from 40 alignments derived from the predicted PCR products made from the S. enterica subspecies 1 and S. enterica houtenae specific biomarker sequences in Table 3.