A total of 214 Salmonella complete genomes and genome assemblies of highest quality (Table S1) were retrieved from the NCBI (National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov/). All the Salmonella genomes with annotations of different serotypes were downloaded as of May, 2016, including most of complete genomes as well as draft genomes of type strains or strains of some serovars for which a complete genome was not available. Five subspecies of S. enterica (I), salamae (II), arizonae (IIIa), diarizonae (IIIb), and houtenae (IV), and a total of 87 serovars of enterica (I), were under consideration. MLST (Multilocus Sequence Typing) (Achtman et al., 2012), together with checking the original studies was implemented to identify their serovars. Housekeeping genes, thrA, purE, sucA, hisD, aroC, hemD, and dnaN were targeted for the MLST scheme (Achtman et al., 2012). The nucleotide sequences were then aligned against seven housekeeping genes. Outcome of sequences were subsequently submitted to the MLST database (http://mlst.warwick.ac.uk/mlst/dbs/Senterica) and assigned existing or novel allele type numbers. The composite sequence STs (sequence types) were defined by the database based on the set of allelic profiles derived from each of the seven loci.

OrthoMCL (version 2.0.9) (Li et al., 2003) was used to determine homologous clusters of the pan-genome for subsequent analyses. All extracted protein sequences were adjusted to a prescribed format and were grouped into homologous clusters based on sequence similarity. The BLAST reciprocal best hit algorithm (Moreno-Hagelsieb and Latimer, 2008) was employed with 70% match cutoff, and Markov Cluster Algorithms (MCL) (Enright et al., 2002) were applied with an inflation index of 1.5.

To determine the phylogenetic relationships among Salmonella subspecies and different serovars, a phylogenetic tree was inferred from single nucleotide polymorphisms (SNPs) of core genes. According to the identification of homologous clusters, a total of 127,544 genes presented in single-copy orthologous genes shared by all 214 Salmonella strains were selected for phylogenetic analysis. For each orthologous cluster, nucleotide sequences were aligned and then concatenated. Custom Perl scripts were used to parse the concatenated sequence and extract SNPs. The organismal phylogeny of SNPs and other genes were inferred using the maximum likelihood (ML) algorithm in RAxML version 8.2.9 (Stamatakis, 2014). The tree reconstruction is repeated 1,000 times, and the percentage of times each interior branch which is noted as bootstrap-value is given. Each tree of genes in the two flagellin gene clusters was constructed by MEGA (Kumar et al., 2016), using Neighbor-Joining method with 1,000 bootstraps.

To characterize the evolution pattern of flagellin genes, reference sequences of fliC, fljB, and the flanked genes, fliA, fliB, fliD, fliS, fljA, and hin, were extracted from the S. Typhimurium LT2 genome (GenBank accession number AE006468) (McClelland et al., 2001). The relevant homologous clusters in the pan-genome were determined by BLAST with 70% match cutoff and 50% coverage cutoff. In cases of shared homologous genes related to flagellin genes, phylogenetic trees were inferred by ML algorithms in RAxML version 8.2.9 (Stamatakis, 2014). Together with the phylogenetic trees, the flanked gene contents were taken in account to confirm the gene loci of fliA, fliB, fliC, fliD, fliS, fljA, fljB, and hin.

PAUP version 4.0a152 (Swofford, 1998) was used to conduct maximum-likelihood (ML) analysis of the congruence of gene trees, as described by Feil et al. (2001). The HKY85 model of DNA substitutions was used to estimate α parameter, and the transition/transversion (Ti/Tv) ratio, assuming gamma distribution. The difference between the log likelihood scores of the trees of two genes Δ-lnL was calculated, generating from the comparison of one gene tree against another based on the 99th percentile of the distribution of log likelihood scores for 200 trees from random topology. If the Δ-lnL values for the comparisons is less than that reference gene comparisons within the 99th percentile of this distribution, then these two gene trees should be considered significantly congruent (Octavia and Lan, 2006).

A key element of our study was to determine whether the evolution of flagella gene was associated with the evolutionary rate estimate. To test for this pattern we conducted a least squares linear regression for the genome-wide rate as a function of sampling distance. The MEGA version 7 (Kumar et al., 2016) was used for pairwise distance calculating. A potential shortcoming of this analysis is that the data do not necessarily represent independent samples because some correspond to closely related lineages. This can be addressed by using phylogenetic independent contrasts, or phylogenetic generalized least squares. For the regression of the rate as a function of sampling evolutionary consensus the degree of association was estimated when fitting a regression for a ratio as a function of its denominator.

Potential recombination within fliB, fliC, and fliD genes was screened using seven methods (RDP, GENECONV, MaxChi, Bootscan, Chimera, SiScan, and 3Seq) implemented in the Recombination Detection Program version 4.88 (RDP4) (Martin et al., 2015). The breakpoints were also defined by RDP4. Similarity between the recombinants and their possible major and minor parents was estimated using BootScan, embedded in RDP4 (Stamatakis, 2014). The recombination event evaluated by RDP4 was considered significant if it satisfied at least two criteria when the P-value (p) < 0.05 and the RDP recombination consensus score (RDPRCS) was >0.4. When p < 0.05 and the RDPRCS was between 0.4 and 0.6, the recombination event was possible. An RDPRCS < 0.4 with p < 0.05 indicated the rejection of the recombination event (Martin et al., 2015).

A total of 214 Salmonella genomes were analyzed in this study, consisting of 205 S. enterica subsp. enterica genomes and 9 genomes of other four S. enterica subspecies, which covered 87 serovars. More than one genome (maximum twelve) was included for some famous serovars (e.g. S. Enteritidis and S. Typhi), while only one genome was considered for most other serovars because no more genomes were available. More detailed information on these genomes was summarized in Table S1. After homologous gene clustering, the total 947,793 proteins were grouped into 14,583 homologous clusters, including common genes represented by 596 single-copy groups of 712 core homologous clusters. A total of 57,946 SNPs in the “core” regions were identified that were present and aligned with high confidence in all taxa. These SNPs were used to construct the phylogeny (Figure 1).

The role of the phylogenetic tree based on SNPs in genome variation and evolution has been documented previously (Filliol et al., 2006; Holt et al., 2008). The core genome tree can support the major genetic group phylogeny and reveal much more phylogenetic detail, including the divergence of major phylogenetic groups. Moreover, the core genome tree also revealed a number of strong concordances with serotyping. From Figure 1, core genome tree was able to categorize S. Typhimurium isolates but unable to ultimately separate the S. Stanleyville from the S. Gallinarum with high resolution. In most cases, the genomes within a cluster belonged to the same serovar, as for S. Typhi or Enteritidis. Although some isolates were mixed up with each other, for example, it was hard to differentiate S. Give, S. Schwarzengrund, and S. Bredney, the core genome tree generally separated the serovars from each other.

The prevalence of phase I and phase II related flagelllin gene locus of the 214 genomes was also indicated in Figure 1. In total, fliC gene cluster was found in 207 genomes (96.7%) of 91 serovars (97.8%), including subspecies, while fljB gene cluster was found in 138 genomes (64.5%) of 61 serovars (65.6%). The detailed list (Table S2) was dominated by fliC and fljB genes and included their flanked genes and hin. Some well-studied serovars (Enteritidis, Dublin, Typhi etc.) are monophasic the same as the previous study (Selander et al., 1992; Baker et al., 2007). Nevertheless, the fliC gene was not found present in the genomes of strain S. Poona 2010K-2244, S. Heidelberg B182, S. Abortusovis SS44, S. enterica subsp. salamae 3588/07, S. Bovismorbificans A1608, S. Hadar ABBSB1020-2, and S. London CFSAN001081, of which only a draft genome was obtained for these strains. Duplicated annotations and misprediction of some genes that span contig boundaries may occur in draft genomes. The homologs (Table S2) were able to mitigate but not eliminate these limitations. We further checked the contigs and related annotation of seven genomes of fliC gene cluster absent, and found some residue of the gene cluster at some contigs' ends in each genome. On the other hand, there was no report on the absence of fliC in these serovars, for example, fliC should be present in all S. Heidelberg (Macnab, 1992). Therefore, the absence of fliC gene in these genomes should be a result of misassembly or the related fragment not fully sequenced.

Interestingly, the monophase and biphase being linked to the evolution of S. enterica subsp. enterica showed discrete phylogenetic clusters on the core genome tree (Figure 1). A few serovars with identical phase, clustered together in the core genome tree. In one larger monophasic clade, with 25 genomes, containing S. Dublin, Gallinarum, Pullorum, Stanleyville, Blegdam, Enteritidis, and Berta, all the genomes of each serovar were located together in one distinct cluster in the core genome tree. For the other monophasic serovars, the 45 genomes were identified at different locations, of which another larger clade of 12 genomes, containing S. enterica subsp. enterica serovars Cubana, Rissen, Tennessee, Senftenberg, and Albany were located adjacent to Derby, Havana, and Agona, but were disrupted by Mbandaka and Kentucky of biphase.

Evolutionary analysis of flagellin gene fliC and fljB in Salmonella has been performed previously (Mortimer et al., 2007). However, the co-evolution of these two genes and their flanking genes has never been studied. To evaluate evolutionary congruence between and within the fliC and fljB clusters, ML analysis was carried out, and the results were summarized in Table 1. As is shown, none of the genes in fliC gene cluster were congruent to all the genes in fljB gene cluster, and none of the genes in fljB gene cluster were congruent to all the genes in fliC gene cluster. The gene trees with the largest number of congruencies were those of fliC and hin, which were congruent to four other gene trees. More concretely, fliC gene tree was similar to those of fliB, fliD, fljA, and fljB gene, and hin gene tree was similar to those of fliB, fliD, fljA, and fljB gene. As for fliB (equaling to fljA and fljB gene trees) and fljB (equaling to fliB and fljA trees), they were congruent to two gene trees. And fliA and fliD trees were congruent to fljA and fliB respectively. Overall, only 29.2% (14/48) of the flagellar gene tree comparisons were congruent among these S. enterica strains. As for the fliC and its neighbored genes, fliB and fliD, they were closely related, while for fljA and fljB genes, their relationship was more distant. Moreover, there was generally higher level of congruence between fliC and fljB/fliD genes than between fliA and fliB genes and between fliD and fliS genes. However, correlation calculations were not available between fliC and fliS genes, because approximation limit dynamically readjusted too low.

Besides the relationship among the genes in the fliC and fljB gene clusters, we further studied their congruence to the core genome using linear regression analysis instead of ML analysis, because datasets of the core genome were too large to calculate -lnL score. The mutual information analysis indicated a weak correlation between each flagellin gene and core genome. The low values of both trees for each linear regression indicate a high level of incongruence between individual tree of each gene and the core genome tree, as shown in Figures 2, 3. The topology of fljB, fliC, and fliD gene trees was much distinct from that of the core genome tree, while those of the fliA and fliS genes were more similar to the core genome tree than those of fljB, fliC, and fliD genes. However, the fljA tree for serovars displayed a similar topology to that of the core genome tree, while fljB did not.

To infer evolutionary relationships among different serovars, based upon similarities and differences in their genetic characteristics, we constructed phylogenetic trees of fliA, fliB, fliC, fliD, fliS, and fljA, fljB, hin, and fliC/fljB (Figure 4). Diagrammatic representation of fliB, fliC, and fliD indicated bifurcation to form two distinct lineages, where other Salmonella enteric subspecies joined together with subspecies enteric serovars. The other four subspecies could be easily recognized from fliA and fliS gene tree because of their distant branches, which was more consistent with core genome tree. Subspecies salamae was close to three neighbors, serovar Weltervreden, Sloterdijk and Kentucky, in fljB tree, while these serovars were scattered to different branches in fljA tree. The results of phylogenetic trees were consistent to the congruence analysis (Table 1, Figure 2, 3). Dramatically, the individual gene trees of fliB, fliC (fliC of S. Typhi was an exception) and fliD were consistent with a clustering dividing monophasic and biphasic serovars. This seemed to imply two distinct evolutionary paths of fliC/fljB corresponding to monophasic and biphasic serovars.

Because different genes of the flagellin diverge under different evolutionary pressures, we manually curated the phylogenetic pattern of genes in flagellin gene clusters. For all flagellin and their flanked genes, the individual gene trees (Figure 4) can be used to demonstrate the evolutionary story on individual gene basis. To study the genetic evolution of the 207 fliC and 138 fljB allelic variants, phylogenetic tree based on nucleotide sequences was constructed (Figure 4I). There is a clear separation of monophasic Salmonella groups from the biphasic groups based on the fliC gene, while all the fljB alleles clustered with the biphasic fliC alleles. This may either reflect that monophasic were acting as incipient serotypes with a barrier to recombination with biphasic serotypes, or indicate that there were insufficient incidents in which recombination could have taken place between monophasic and biphasic serotypes. However, S. Typhi was not clustered into other monophasic genomes as one branch together. Similarly, previous study had also identified that monophasic S. Typhimurium existed, of which the serotype was always biphasic (Laorden et al., 2010).

Furthermore, all biphasic fliC and fljB alleles were analyzed based on the amino acid sequences (Figure S1). The genetic distance between fliC gene and fljB gene sequences showed no relationship within/among serovars, which was reflected in the phylogenetic tree. For example, serovar Schwarzengrund, Liverpool, and Livingstone were closely related based on fliC gene while they were distant based on fljB. S. Heidelberg and S. Newport were clustered together according to the fljB gene, while they were clustered separately according to the fliC gene. It suggested that the evolutionary histories were independent between fliC and fljB genes and frequent homologous recombination referring fliC and fljB occurred between biphasic serovars. The phylogenetic tree of 5′ and 3′ conserved segment of fliC and fljB genes was also constructed (data not shown), showing the similar division.

Recombination is a significant source of genetic diversity in Salmonella (Silverman et al., 1979). This phenomenon, which occurs during co-infection with different strains, is usually observed at the fliC and its flanked genes. It is apparent that there were different levels of recombination within flagellar genes of S. enterica. Comparison of recombination of monophasic genes with biphasic genes, which represents the unique ability to change flagellar composition, by both ML (Table 1) and recombination detective analysis (Figure 5), showed that recombination occured far more frequently within monophasic or biphasic gene clusters than between them. For example, the serovar Heidelberg, whose major parent (a sequence closely related to that from which the greater part of the recombinant's sequence may have been derived, Martin et al., 2015) was Virchow, and minor parent (a sequence closely related to that from which sequences in the proposed recombinant region may have been derived, Martin et al., 2015) was Newport, did not recombine to any monophasic genes. It was identical to Hartford, whose major parent was Pomona, and minor parent was Anatum. Even though it was found that Senftenberg (monophase) was minor parent for Newport (biphase), and other similar recombinant events of referring small segments were also found, recombination was believed mainly to have occurred within the same phase. Moreover, there were lots of breakpoints and the breakpoints tended to cluster far more around the edges of genes than they did within the central parts of genes (Figure 5).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.