Whole genome sequencing analysis of multiple Salmonella serovars provides insights into phylogenetic relatedness, antimicrobial resistance, and virulence markers across humans, food animals and agriculture environmental sources

Infection with antimicrobial resistant Salmonella in humans and animals is a global threat that has caught the public attention worldwide [1–3]. Human foodborne salmonellosis causes an estimated 100,000 domestic cases and 40 deaths annually in the United States [1]. The U.S. Department of Health and Human Services reported an increase in Salmonella infections from 13.6 to 16.4 cases per 100,000 population, which represented a 17.1% increase from 1997 to 2011 [4]. In the European Union, Salmonella-infected gastroenteritis was the second most frequently reported foodborne illness with 91,408 clinical cases reported by thirty EU/EEA countries, and a confirmed case rate of 25.4 cases per 100,000 population in 2014 compared to 21.4 cases per 100,000 population in 2013, which represented a 19% increase in the notification rate [3].

Inappropriate use of antimicrobials in livestock production and the association to resistant Salmonella infection in humans are a growing concern to public health agencies, and have led to the rise of new multidrug resistant (MDR) bacteria and transferable genetic loci, such as colistin resistance mediated by the MCR-1 gene [5, 6]. Given the ever-growing requirement to maintain the efficacy of antimicrobials as well as decrease the emergence of antimicrobial resistance in human infections, the antimicrobial use in veterinary and agricultural practices is being extensively re-evaluated [7–9]. Humans and animals are linked to each other through the environmental reservoirs which have long been implicated as a source of Salmonella and antimicrobial resistance found in human and animals [8–10]. The selection pressure on Salmonella is created by antimicrobial use in human health and food animal production leading to development and potential spread of antimicrobial resistance [8–11]. Our previous studies reported the persistence and dissemination of multiple resistant Salmonella serovars along with their determinants in the environment of commercial swine operation due to the manure application on land [12, 13].

Multiple Salmonella serovars, including Agona, Anatum, Derby, Heidelberg, Infantis, Kentucky, Muenchen, Newport, Schwarzengrund, and Typhimurium are commonly detected in food animals, food products, and agricultural environments, and are associated with resistant Salmonella infections in humans [14–17]. The Centers for Disease Control and Prevention (CDC) reported that the incidence of human Salmonella infections caused by monophasic 4,[5],12:i:-, which is in the top 4 of the most frequently reported Salmonella serovars, continue to rise while the incidence of the other serovars is decreasing [16]. The increase in the incidence of this serovar in human cases is paralleled by a similar increase in swine and environmental detection of this serovar variant [12, 13, 18]. However, there are gaps that still exist in our understanding of the temporal and spatial connection of resistant Salmonella transmission within humans, animals, and the environment sources.

A number of studies have used the classical molecular typing methods such as pulsed-field gel electrophoresis (PFGE), multilocus sequence-based typing (MLST), and multilocus variable-number tandem repeat analysis (MLVA) to assess the relatedness and the subsequent transmission of antimicrobial resistant (AMR) Salmonella in human, animals, and environment [19–21]. However, the limitation of these methods lies in insufficient discriminatory power to separate closely related Salmonella isolates in outbreak investigations and to differentiate between the intra-serovar isolates from different hosts [20–22]. The use of whole genome sequencing (WGS) has had a major impact on the study of the molecular epidemiology of AMR bacterial pathogens associated and transmitted between human, animal and environmental sources. A WGS study in Denmark reported that SNP, pan-genome, k-mer and nucleotide difference trees were superior to the classical typing method and evaluated the association of the isolates to specific outbreaks of S. Typhimurium [23]. Additionally, WGS has been used to identify known AMR determinants among strains of Escherichia coli and Salmonella [24, 25]. The objectives of this study were to use WGS to analyze multiple Salmonella serovars isolated from human, food-animals and environments in the two states of the US and to clarify the epidemiological transmission of AMR Salmonella within these studied populations. In addition, the capability of WGS to predict antimicrobial resistance and virulence genes in antimicrobial resistant Salmonella retrieved from different sources was evaluated.

The objective of this study was to characterize Salmonella serovar, AMR determinants and virulence genes using whole genome sequencing. The 200 Salmonella enterica genomes were isolated from different sources of origin including human, swine, poultry, and environment, and were analyzed using the core genome SNP-based analysis and the alignment-free analysis method FFP. The phylogenetic trees obtained from parSNP and FFPry showed that the clusters observed matched Salmonella serovars (Fig. 1). The branch length in FFP-based trees is more representative of differences over the whole genome, which may be due to differential plasmid, prophage content, or other accessory genome [28], while SNP-based trees use the core genome derived from whole-genome alignment and read mapping for phylogeny construction [27, 28]. The major difference of the SNP- and FFP-based analyses was in the order of the serovar clusters within the tree, however, the overall approach selected had relatively little effect on the topology of the phylogenetic trees (Fig. 1). A number of studies have reported the use of SNP-based analysis as a potential molecular subtyping tool for outbreak investigation in multiple Salmonella serovars including Dublin [29], Enteritidis [22, 30], Heidelberg [31], Manhattan [21], Montevideo [32, 33], and Typhimurium [34–36].

The phylogenetic analysis based on WGS-derived SNPs has been shown to provide greater cluster resolution than the gold standard subtyping method, pulsed-field gel electrophoresis (PFGE), resulting in discrimination of outbreak-related human clinical isolates and food or environmental origins [21–23]. ParSNP was run for all the 200 Salmonella genomes and separately for each of the serovars. Although the separate comparisons indeed increased the core genome component as would be predicted, this was not a major increase, and we do not expect that the removal of the outliers will have a significant effect on the trees. In concordance with those prior literatures, the current study found that the core genome SNP-based trees of individual Salmonella serovar including Typhimurium and 4,[5],12:i:- (Fig. 2), Schwarzengrund (Fig. 4), and Rissen (Fig. 5) were mostly clustered based on source of origin. However, there were some exceptions in each individual tree. As shown in Fig. 2, some clinical S. Typhimurium isolates (HS71549, HS51537, and HS51628) were closely related to the environmental, swine, and chicken isolates, respectively. The Schwarzengrund cluster in Fig. 4 showed that the Salmonella isolates from chicken feces clustered with the isolates obtained from environmental isolates which were derived from the same farm. There was no isolation source-dependent clustering in S. Derby (Fig. 3), with genomes from human, swine and environmental isolates clustering together. These findings can point towards the potential transmission of Salmonella among humans, animals and the environment and support the idea of zoonotic transmission, while independent human sub-clustering in each serovar might be referred to human-to-human transmission. However, the human Salmonella isolates included in our analysis were only from the North Carolina State Public Health Laboratory which might not be represent all the human clinical cases. According to the same timeline as human Salmonella outbreaks belonged to the independent human sub-clusters, the animal/environmental sources which may have a chance to group with those human sub-clusters were not scheduled for sampling. In addition, S. 4,[5],12:i:- has been defined as a monophasic variant of serovar Typhimurium because of their antigenic and genetic similarities, and the characterization of S. 4,[5],12:i:- using the typical molecular approaches revealed that S. Typhimurium is the direct ancestor of S. 4,[5],12:i:- [37]. Even the two serovars were clustered together (Fig. 2), parSNP-based subtyping could be a suitable analysis applied to differentiate these serovars. In contrast to S. Derby (Fig. 3), the sources of origin cannot be differentiated using parSNP analysis. This serovar has a highly homogeneous genetic composition and can be carried by different hosts [38]. Moreover, the MLST database (http://mlst.ucc.ie/mlst/mlst/dbs/Senterica/) indicates that Derby is a polyphyletic serovar, having originated from more than one common ancestor, and possesses several distantly related sequence types (ST) [39]. Thus parSNP-based analysis might not be an appropriate method for this serovar. However, a recent study in China reported that the clustered regularly interspaced short palindromic repeats (CRISPRs) could be a useful subtyping tool for S. Derby in molecular epidemiological investigations [40]. Though the SNP typing is the reliable tool for genomic and epidemiologic studies, it is not without limitations. SNP-based analysis requires alignment of whole genome sequences and only utilizes the core genome, which may be less sensitive as a result. In addition, this method is still limited to the intragenus analysis of closely related species and strains [27, 41].

The FFP phylogenetic clustering is an effective tool that relies on an alignment-free approach for genomic evolution study. Unlike parSNP-based method which focuses on core genome, the phylogenetic trees acquired from FFP-based method are affected by recombinant genes and/or horizontal gene transfer including plasmid, prophage, and other accessory gene contents [28]. However, the main clustering structure is not significantly different (Fig. 1). The advantages of FFP-based analysis are that it is independent of a reference genome, and has lower hardware requirements. Additionally, FFP analysis can be performed with whole genome shotgun samples as it is not affected by contig orientation, and contig order [28]. FFP-based analyses has been reported in a number of bacterial genomic studies, including Helicobacter pylori [28], Bacillus spp. [42], Escherichia coli [43, 44], and Shigella [44]. These studies have revealed that the FFP method can contribute to the phylogenetic clusters based on geographic relation and outbreak detection, and could provide a complementary analysis approach. Our study is the first to utilize the alignment-free FFP analysis in Salmonella and compare it to the core genome SNP-based analysis. We found that the phylogenetic clusters from these methods were similar in term of serovar characterization, but the branching varied due to differences in analysis approaches (Fig. 1). While SNP- and MLST-based methods are likely to continue to be the default choice for subtyping and comparative genomics in Salmonella, the FFP method can serve as a useful alternative method requiring relatively low-powered computing resources [28].

Antimicrobials are reported extensively used in food animal production to treat clinical disease, to prevent and control common diseases, and to enhance animal growth [45]. Tetracycline and tylosin are frequently mixed in animal feed for disease prevention and growth promotion purposes in commercial swine and poultry systems [45, 46]. In accordance to the high percentages of phenotypic tetracycline and sulfisoxazole resistance were reported in our result. The WGS revealed a number of tetracycline and sulfisoxazole resistance genes such as tetA, tetB, tetC, tetR, sul1, and sul2 (Additional file 2: Table S2). Of interest, gene mechanisms of tetracycline resistance including the efflux genes, the ribosomal protection and enzymatic genes were suggesting a possible ecological role for specific wide spread of tetracycline resistance [47]. However, AMR genes especially tetracycline and sulfonamide were also detected in livestock production surrounding even when antimicrobial drugs were not administered to animals [47, 48].

In this study, we have shown that WGS is an excellent tool for accurately predicting antimicrobial resistant phenotype in human, animal, and environment associated multiple Salmonella serovars, as WGS predictions and phenotypic resistance matched well with high sensitivity and specificity in our study. Overall, the resulting resistance genotypes correlated with 87.61% sensitivity and 97.13% specificity to the resistance phenotype (Table 2). Among the discordant results in our study, the lowest specificity of AMR prediction was observed for streptomycin which accounted for the presence of streptomycin-resistance genes but lacked phenotypic resistance. This finding was in concordance with the previous studies in Salmonella [24, 49] and E. coli [25, 50]. The streptomycin discrepancies have been commonly detected in other studies too because streptomycin is not used to treat enteric infections, and as such, results in the absence of precise clinical breakpoint for streptomycin susceptibility in Salmonella and E. coli [24]. Although the strA/strB and aadA genes were detected, the strA/strB genes conferred higher resistance than aadA genes [25, 51]. Thus, the presence of aadA genes by in silico method may not result in streptomycin resistance phenotypically. In addition, the mechanism of streptomycin resistance is frequently due to lacking of the gene expression as well as mutations in the 16S rRNA gene leading to difficulty of phenotypic prediction [50, 52]. Our results suggest that the refinement of WGS-based AMR prediction could be beneficial and can definitely enhance the monitoring of AMR strains and determinants detected in humans, foods, animals, and environment.

The Salmonella serovars significantly correlated with the presence/absence of AMR genes, plasmid replicons, and virulence genes. We observed specific AMR genes in each Salmonella serovar (Table 3). This result along with the phylogenetic relatedness revealed that the type of serovar in discussion had the greatest impact on AMR characterization. Previous studies reported the presence of AMR genes has been shown to be primarily associated with a particular host and is not frequently transmitted among different species which in accordance to our finding (Additional file 2: Table S2) [20, 49]. Multiple plasmid replicons were detected in this study using WGS method (Additional file 2: Table S2). Plasmids were observed specific to Salmonella serovar that was very similar to the AMR genes (Table 3). This is in accordance to our previous study that the plasmid profiles were correlated to Salmonella serovar and incompatibility (Inc) groups [13]. The IncF (both FI and FII) family found across the different Salmonella serovars in our study is known to be a well-adapted and commonly distributed plasmid among members of the Enterobacterceae family [53, 54]. Although our data cannot fully explain the transmission of AMR determinants among various species, they are in line with previous studies that reported on the role of animals and environment as important sources of multiple AMR genes as well as plasmids, and that transmission can occur by horizontal gene transfer [13, 55, 56].

Multiple virulence genes were identified among the several Salmonella serovars across different sources by WGS (Additional file 2: Table S2). These genes have been described to be involved in several processes important for Salmonella transmission and infection, including adhesion, type III secretion system (T3SS), host recognition/invasion, filamentous formation, magnesium uptake, iron acquisition, and regulation of stress factors. Our data showed that Salmonella isolates recovered from animal or/and environmental sources contained the same virulence genes as carried by human clinical isolates. Along with the phylogenetic analysis, these findings support our view that the high frequency of virulence genes detected in food animal and environment may be transmitted and cause infections in humans, a suggestion that has been previously made in prior studies [57–59]. Figueira et al. (2013) reported that the lack of sseJ, a particular virulence gene makes S. Typhimurium strain became more heterogenous [60]. In our study, this gene was only present in S. Derby (Table 3) which may relate to the non-source-dependent clustering found in this homogenous serovar as mentioned previously (Fig. 3). One of the typhoid-associated virulence factors, the cytolethal distending toxin cdtB, was detected in all isolates of Schwarzengrund, Johannesburg, and Muenster (Additional file 2: Table S2), which was similar to a previous study that detected this gene in S. Schwarzengrund [61]. The cdtB encodes the typhoid toxins of S. Typhi and is not reported from a wider distribution among non-typhoidal Salmonella serovars (NTS) [57, 61]. However, there were a few reports of the prevalence of this virulence gene in several NTS, including Javiana [57], Montevideo, Schwarzengrund, and Bredeney [61]. This data suggested that the cdtB toxin may contribute to the pathogenicity in human and animal.

WGS is a helpful tool to assess the phylogenetic relations among multiple serotypes, AMR and virulenuce gene evaluation and assist in the molecular epidemiological studies of foodborne pathogens. The SNP-based and FFP-based analysis provided the higher resolution Salmonella phylogenetic trees that could differentiate the isolates recovered from human, animal, and environment. In addition, WGS is a useful tool for AMR prediction, plasmid replicon, and virulence gene detections. Our study shows the close relationship between Salmonella isolates associated with different hosts, which is supportive of possible zoonotic transmission. This is seen among multiple serotypes, and the prevalence of AMR genes, plasmid replicons and virulence genes that were identical in different species and could potentially highlight exchange of serovars across different hosts.