Genomic characterization and phylogenetic analysis of Salmonella enterica serovar Javiana

Sequencing, preprocessing, and genome assembly of TN isolates

BioSample numbers and metadata for Salmonella ser. Javiana (n = 111) isolates from patients in TN from January 2017 through October 2018 were obtained from the Tennessee Department of Health (TDH) (Data S1). Tennessee population data (2018) used for calculating incidence rates (IR) was obtained from the U.S. Census Bureau (U.S. Census Bureau, 2020) and IR per county were mapped using Tableau Desktop Public Edition (v2019.2.1) (Tableau Software, 2019). PFGE and whole-genome sequencing were performed by the TDH Division of Laboratory Services according to PulseNet protocols (Centers for Disease Control and Prevention, 2016; Centers for Disease Control and Prevention, 2017b). For PFGE, XbaI was used as the primary restriction enzyme. Genomic DNA was extracted using Qiagen DNeasy Blood & Tissue kits, libraries were prepared using Nextera XT kits, and sequencing was performed on an Illumina MiSeq platform using Illumina MiSeq v2 chemistry (500 cycle) to produce 250bp paired-end reads. Raw reads were downloaded from the NCBI SRA database, trimmed using Trimmomatic v0.35 (Bolger, Lohse & Usadel, 2014) (with the following parameters: ILLUMINACLIP: NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36), and quality checked using FastQC v0.11.7 (Andrews, 2010) and MultiQC v1.5 (Ewels et al., 2016) to combine the results. The trimmed reads were assembled into contigs using SPAdes v3.12.0 (Bankevich et al., 2012) with the careful option. Assembly statistics were generated by BBMap v38.88 (Bushnell, 2018), SAMtools v0.1.8 (Li et al., 2009), and QUAST v4.6.3 (Gurevich et al., 2013). SeqSero (Zhang et al., 2015) was used to confirm serovar designations.

SNP detection and phylogenetic analyses

A reference-free SNP detection analysis was initially performed with the TN isolates to determine the overall population structure free of reference choice bias. The assemblies were analyzed using KSNP3.1 (Gardner, Slezak & Hall, 2015) and the resulting core SNP matrix fasta file was then used to construct a phylogenetic tree in Mega7 (Kumar, Stecher & Tamura, 2016) with 100 bootstrap replicates (Felsenstein, 1985). The evolutionary distances were computed using the number of differences method (Nei & Kumar, 2000) and the evolutionary history was inferred using the Neighbor-Joining method (Saitou & Nei, 1987). The final tree was visualized and annotated using iTOL (Letunic & Bork, 2016). Isolates that weren’t serovar Javiana (based on SeqSero results) and were very divergent based on the KSNP analysis were removed from the analysis. TN clades (defined as groups of three or more isolates that were all within 500 SNPs of each other) were identified. Next, reference-based hqSNP analyses were performed for each TN clade independently to determine high-resolution SNP differences between isolates. For the hqSNP analyses, an appropriate internal reference genome assembly (with adequate assembly quality and expected assembly size and G+C content) for each clade was identified (SRS2420927 for TN clade I, SRS2822480 for TN clade II, and SRS3010019 for TN clade III). Additionally, the Salmonella enterica subsp. enterica serovar Javiana str. CFSAN001992 assembly (GCF_000341425.1) was downloaded from the NCBI RefSeq database for use as an external and closed reference genome. The hqSNP analyses were performed, both with the internal and external references and for the 111 isolates together and for each TN clade individually. For each analysis, high quality single nucleotide polymorphisms (hqSNPs) were identified using the CFSAN SNP Pipeline v1.0.1 (Davis et al., 2015). The resulting hqSNP matrix fasta files were then used to construct phylogenetic trees as described above. The matrices were sorted and clustered using the hclust function (gtools package) in R studio. For the individual clade analyses using internal references, clusters of two or more related isolates were identified at hqSNP distance threshold levels of 5, 10 and 25; isolation date and other epidemiological information were not considered.

Genome annotation and pan-GWAS

TN isolate assemblies were annotated using Prokka v1.14-dev (Seemann, 2014) and RASTtk (Brettin et al., 2015). A pan-genome-wide association study (pan-GWAS) was performed to compare gene content among the isolates using Roary v3.12.0 (with Prokka annotation output files, previously described, used as input files) (Page et al., 2015) and statistical analysis was done using Scoary v1.6.16 (with the following arguments: -c I B BH PW EPW P -p 0.05 -e 100) (Brynildsrud et al., 2016) to identify genes or markers associated with inclusion in each TN clade. Genes predominantly present or absent among isolates in each clade were identified. Genes or clusters of genes (loci) were considered significantly associated with a clade if they had a Bonferroni-corrected P-value <0.05. From the pan-GWAS, the positively and negatively associated genes were classified by having an Odds Ratio of >1 and <1, respectively. To further analyze the results in a genomic context, loci that were adjacent or located in close proximity were combined into a single region.

In silico identification of genomic features and genes

Phage regions were identified in a diverse representative subset of the TN isolates (n = 31) using Phaster (Arndt et al., 2016; Zhou et al., 2011). Potential plasmids were predicted and classified using PlasmidFinder v2.0.2 (database version 2019-05-16) (Carattoli et al., 2014), Unicycler v0.4.8-beta (Wick et al., 2017), and plasmidSPAdes (Antipov et al., 2016) and visualized using Bandage (Wick et al., 2015). They were further confirmed by examining the associated assembly contigs for plasmid-associated genes, comparing to known plasmids using PLSDB (Galata et al., 2018) and BLASTn, and comparing coverage and G+C content to the whole assembly. Antibiotic resistance (ABR) determinants in genomes were predicted using ResFinder (90% threshold for identity and 60% for minimum length) (Zankari et al., 2012) to identify acquired ABR genes and PointFinder (90% threshold for identity and 60% for minimum length) (Zankari et al., 2017) to identify point mutations conferring ABR. The representative subset of the isolates (n = 31) was also examined for virulence factors using VirulenceFinder (Joensen et al., 2014) and VFDB VFanalyzer (Liu et al., 2018) and Salmonella Pathogenicity Islands (SPIs) using SPIFinder (v1.0; 95% threshold for identity and 60% for minimum length) (Roer et al., 2016).

Global phylogenetic analysis

The Salmonella database (Alikhan et al., 2018) on EnteroBase (Zhou et al., 2019) was queried for isolates with “human” listed as source niche in the strain metadata and “Javiana” listed as serovar in the strain metadata or experimental data (SISTR1 (Yoshida et al., 2016) or SeqSero2 (Zhang et al., 2019)). Only strains with country (and state for strains from the United States) included in the metadata were retained. A cgMLST + HierCC minimal spanning tree (RapidNJ algorithm) was created with GrapeTree (Zhou et al., 2018) on EnteroBase using all of the resulting strains. Strains that were likely not Javiana (conflicting serovar designations and distant on the tree) were removed from the dataset, leaving 466 strains (Data S2). The dataset was further refined to select representative strain(s) for each HC100 level cluster (≤100 cgMLST allelic differences). For each HC100 cluster, a single strain from each country/state was retained. If there were more than one strain from a country/state, the representative strain was chosen based on assembly quality (N50; and if N50 values were identical or similar, coverage and number of contigs were also considered). At least one TN strain representing each HC100 cluster (if available) was chosen to be included in the final dataset of genomes representing global diversity of of Salmonella ser. Javiana clinical isolates. The final dataset consisted of 162 strains: 29 TN isolates (11 from TN clade I, 3 from TN clade II, 10 from TN clade III, and the five isolates that didn’t fall into the main clades), 43 strains from other states in the US, and 90 isolates from other countries (Data S2). Collection date years ranged from 2002 to 2020. Assemblies for the non-TN strains were downloaded from Enterobase. All assemblies were analyzed using KSNP3.1 (Gardner, Slezak & Hall, 2015) and the resulting core SNP matrix fasta file was then used to construct a phylogenetic tree in Mega7 (Kumar, Stecher & Tamura, 2016) with 100 bootstrap replicates (Felsenstein, 1985). The evolutionary distances were computed using the number of differences method (Nei & Kumar, 2000) and the evolutionary history was inferred using the Neighbor-Joining method (Saitou & Nei, 1987). The final tree was visualized and annotated using iTOL (Letunic & Bork, 2016).

Comparison to polyphyletic serovars

Eight strain datasets were created: one for serovar Javiana and one each for other serovars that have been reported as polyphyletic (Derby, Kentucky, Mississippi, Montevideo, Newport, Saintpaul, and Senftenberg). The Salmonella database (Alikhan et al., 2018) on EnteroBase (Zhou et al., 2019) was queried for isolates with the specified serovar listed as serovar in the experimental data (SISTR1 (Yoshida et al., 2016) or SeqSero2 (Zhang et al., 2019)). For each, a cgMLST + HierCC minimal spanning tree (RapidNJ algorithm) was created. Strains that were likely not the serovar of interest (conflicting serovar designations and distant on the tree) were removed from the datasets. The final datasets (Data S3) were used to create cgMLST + HierCC minimal spanning trees (improved minimal spanning tree algorithm, MSTree V2) using GrapeTree (Zhou et al., 2018) on EnteroBase. The branch lengths between the cgMLST eBurstGroups (ceBGs) of the other polyphyletic serovars were used for comparison to the branch length between the two Javiana ceBGs. ceBG designations associated with each serovar were retrieved from the EnteroBase documentation (EnteroBase Team, 2018).

Tennessee Salmonella ser. Javiana population structure

This analysis included a diverse set of 111 Salmonella ser. Javiana clinical isolates from TN (Data S1). On average, the assembled genomes from this study had 74.6x coverage, contained 71.7 contigs (34.8 contigs ≥ 1 kb), were 46.45 kb in length, and had 52.11% GC content (Data S1). Based on the KSNP analysis, the Salmonella ser. Javiana isolates from TN displayed a population structure with three main clades (Fig. 1). TN clade I contained 54 isolates, TN clade II contained four, TN clade III contained 48 isolates, and five isolates didn’t fall into the main clades (Fig. 1 and Table 1). Isolates in TN clades I, II, and III had average hqSNP distances of 119.4 (range 0 to 631), 210.3 (range 3 to 396), and 66 (range 0 to 361), respectively (Table 1 and Data S4).

The 111 TN Salmonella ser. Javiana clinical isolates represented 47 different PFGE patterns (Fig. 2; Data S1). The most common PFGE patterns were JGGX01.0012 (n = 18), JGGX01.0065 (n = 17), and JGGX01.0072 (n = 10). TN clade I isolates represented 28 different PFGE patterns, with the most common being JGGX01.0012 (n = 18). TN clade II isolates represented three different PFGE patterns. TN clade III isolates represented 11 different PFGE patterns, with the most common being JGGX01.0065 (n = 17) and JGGX01.0072 (n = 10). Each of the five isolates that did not belong to a clade had a distinct PFGE pattern. All PFGE patterns were unique to only one TN clade.

High-quality single nucleotide polymorphism (hqSNP) analysis for cluster detection

Using a 5, 10, and 25 hqSNP distance threshold for cluster identification, 9, 12, and 10 potential clusters were identified, respectively (Data S4). The number of clusters decreases from thresholds of 10 to 25, as with the larger threshold, some of the clusters contain multiple subclusters identified at the lower threshold. Within TN clade I, five clusters were identified at each hqSNP distance thresholds of 5 and 10, and four clusters at 25 (Data S4). Only one cluster was identified at each of the distance thresholds for TN clade II (Data S4). Within TN clade III, three clusters were identified at a distance threshold of 5 hqSNPs, six clusters at 10 hqSNPs, and five clusters at 25 hqSNPs (Data S4).

To evaluate the effects of reference choice and isolate diversity, we ran our hqSNP analyses on the entire TN dataset (111 isolates) and on the three clades individually and with both internal draft assemblies and an external closed assembly as reference genomes. When all isolates were analyzed together using different reference genomes, the average hqSNP distance was lowest when using the internal reference from TN clade I and highest when using the external closed genome and the average percentage of reads mapped differed by up to 1.57% (Table 1 and Data S4). It should be noted that the closed external reference assembly we used (GCF_000341425.1) was most closely related to TN clade II in our original KSNP analysis (Fig. 1) and not representative of the overall population (at the time that this analysis was performed, there was only one closed Salmonella ser. Javiana genome available on NCBI RefSeq). For all three clades, regardless of the reference genome used, hqSNP distances were higher when they were analyzed independently, with the differences being only slightly higher for clades I and III. The average hqSNP distances for clades I and II were lower when using the internal references, but were slightly higher for TN clade III. For all three clades, the average percent reads mapped was higher when using the internal reference genome than the external reference genome, which is to be expected.

Epidemiological trends

The TN Salmonella ser. Javiana isolates were sourced from patients with an average age of 40.0 (range of 1 month to 90 years; standard deviation of 29.1). The highest incidence was in patients ≤4 (6.64 per 100,000) and ≥85 (4.16) years-old. Previous studies have reported that Salmonella ser. Javiana infections are more prevalent in infants and young children than for other serovars (Jones et al., 2008; Shaw et al., 2016; Srikantiah et al., 2004). Overall, 51.4% of isolates were from male patients (incidence of 1.73) and 46.8% from female (incidence of 1.50; Table 2). TN clades I and II contained more isolates from male patients than female, 57.4% and 75%, respectively. Conversely, 52.1% of TN clade III isolates were from female patients.

Most of the TN Salmonella ser. Javiana isolates were taken from stool samples (84.7%, n = 94), followed by urine (6.3%, n = 7) and blood (5.4%, n = 6) (Table 2). The portion of isolates taken from blood samples exceeds the 2.8% invasive disease outcome (defined by isolation from blood, cerebrospinal fluid, bone or joint fluid, or another sterile site; does not include urine, wound, abscess cultures) reported for this serovar in FoodNet states (Jones et al., 2008). Of the urine isolates, most (n = 6) belonged to TN clade III and all were from females, with an average patient age of 55.1 (standard deviation of 26.2). All of the isolates recovered from blood samples belonged to TN clades I (n = 3) and III (n = 3) and most were from males (n = 4), with an average patient age of 71.8 (standard deviation of 15.2). The majority of isolates recovered from extraintestinal sites were collected from elderly patients, which may indicate a correlation between invasive infection and age.

Geographical and temporal distribution

A geographical distribution can be seen for the TN isolates, with 65.8% (n = 73) isolated in counties in the western region of TN (Table 2; Data S1). In contrast, only 17.1% (n = 19) and 16.2% (n = 18) were isolated in counties in east and middle TN, respectively. Incorporating county population data, the west region had an IR of 4.69 clinical isolates per 100,000 population, while the east had an IR of 0.79 and the middle had an IR of 0.64, with an overall IR of 1.62 per 100,000 for the state (Fig. S2). Three counties had noticeably higher IR: Madison with 31.0, Crockett with 27.9, and Carroll with 25.0. As was the trend for all isolates, the majority of TN clade I and III isolates originated in the west region (59.3% and 83.3%, respectively; Table 2). However, a sizable amount of TN clade I isolates also originated in the east region (25.9%). Additionally, a temporal distribution was also clear, with 68.5% (n = 76) of isolates collected in July through September (Table 2; Data S1). This trend was also seen within the three TN clades.

Identification of mobile genetic elements (MGEs)

All of the TN isolates were examined for plasmids and 32 putative plasmids (19 unique) were identified in 30 isolates (Table 3). They ranged in size from 23 to 108 kb and included replicon types (a plasmid typing scheme based on replication control regions (Carattoli et al., 2014)) IncFIB, IncFII, IncI1, IncN3, and IncX4. In the representative subset of TN isolates (n = 31), Phaster predicted an average of 8.90 [range of 5-13] prophage regions per isolate (2.29 intact, 4.94 incomplete, and 1.68 questionable) (Data S5). The TN clade I isolates had the highest number of predicted prophage regions (average of 9.64 and range of 8–11), followed by TN clade II (average of 8.25 and range of 8–9) and TN clade III (average of 8.09 and range of 5–10).

Identification of virulence factors and pathogenicity islands

All of the representative subset of TN isolates (n = 31) contained pathogenicity islands C63PI, SPI-13, and SPI-14 (Data S6). Most of the TN isolates examined contained SPI-2 (except for SRS2442409 [TN clade II] and SRS2998834 [TN clade III]) and SPI-4 (except SRS3453943 [TN clade I], SRS3643364 [TN clade III], and SRS2998834 [TN clade III]).

Most identified virulence genes were present in all of the TN isolates analyzed (Data S6), including the three genes (cdtB, pltA, and pltB) encoding the subunits of the cytolethal distending toxin (CDT) (Miller et al., 2018). Only the TN clade II isolates contained the pefC and pefD genes, which were on plasmid-associated contigs, that are part of the pef (plasmid-encoded fimbriae) operon and associated with fimbrial adherence (Bäumler et al., 1996). Some of the isolates (two from TN clade I and three that were not part of a clade) contained genes from the saf (Salmonella atypical fimbria) operon on putative plasmid-associated contigs (Folkesson et al., 1999). One isolate (SRS2442415, no clade) contained 11 genes associated with the yersiniabactin iron uptake system on a plasmid-associated contig (Carniel, 2001).

Identification of antibiotic resistance genes

All 111 TN Salmonella ser. Javiana isolates analyzed in the present study contained the aac(6′)-Iaa gene, which has been previously reported to confer aminoglycoside resistance (Shaw et al., 1993). One isolate (SRS2783476; TN clade I) contained the aph(3′)-Ia and sul3 genes on a contig (contig 32) associated with a putative plasmid. The aph(3′)-Ia gene has been found to confer resistance to aminoglycosides (Shaw et al., 1993) and sul3 gene has been shown to confer resistance to sulfonamides/sulfones through antibiotic target replacement and has been shown to be associated with resistance to sulfamethoxazole (Perreten & Boerlin, 2003).

Additionally, one isolate (SRS2628542; no clade) contained the qnrB19 gene, which has been shown to confer resistance to fluoroquinolones through physical protection of the antibiotic target (Correia et al., 2017). The qnr gene is plasmid-associated and has been linked with reduced susceptibility to ciprofloxacin (Casas et al., 2016; Crump et al., 2015; Redgrave et al., 2014). It is unclear if this gene in SRS2628542 is on a plasmid, as it is located on a small (703 bp) contig and no plasmids were predicted in this isolate. However, the gene showed 100% identity over only 72.6% of the alignment with the reference gene (accession EU432277), so it may not be functionally capable of conferring the quinolone resistance phenotype.

Clade-enriched genes

A pan-genome analysis of the TN isolates revealed that the core genome (genes contained in ≥99% of isolates) consisted of 4,022 genes and the accessory genome consisted of 3,920 genes (Table 1). The difference in gene content between the identified clades were mostly found in mobile genetic elements (91.5%).

TN clade I isolates had a core genome of 4,106 genes and an accessory genome of 2,513 genes (Table 1). This clade has a much larger accessory genome than the other two clades identified in this study. This is likely due to the large variety of mobile genetic elements (e.g., plasmids, prophages) present in isolates from this clade, which is also reflected in the much larger number of PFGE patterns (n = 28) displayed by these isolates as compared to other clades (n = 3 for TN clade II and n = 11 for TN clade III). The pan-GWAS revealed 153 loci (genes or groups of genes) to be significantly associated with inclusion in this clade (54 positively associated and 99 negatively associated) (Table 1 and Data S7). These loci consisted of 338 total genes (94 positively associated and 243 negatively associated) (Table 1). The positively associated loci are found in 29 distinct genomic regions and the majority of the genes over-represented in TN clade I (73 genes) were located in eight prophage regions (Table 4).

Twelve of the overrepresented genes in TN clade I correspond to pathogenicity-related protein families (as identified by PathogenFinder): DNA damage-inducible protein I, phenolic acid decarboxylase subunit D, small toxic polypeptide LdrD, PTS system fructose-specific EIIB’BC component, PTS system mannose/fructose/sorbose family IID component, prepilin peptidase dependent protein A precursor, phage DNA binding protein, and other hypothetical proteins. Overexpression of ldrD, which is part of a chromosomal toxin–antitoxin gene system (Alix & Blanc-Potard, 2009; Kawano et al., 2002), is toxic to the cell and leads to growth inhibition and rapid cell killing (Fozo et al., 2008; Kawano et al., 2002). ldrD homologs have not been found in plasmids, but may be involved in cellular response to environmental stress (Kawano et al., 2002). Prepilin peptidase dependent protein A precursor is known to be plasmid-associated (Raspoet et al., 2019; Zhang, Lory & Donnenberg, 1994) and is involved in processing of the major pilus subunit (Filloux, Michel & Bally, 1998; Zhang, Lory & Donnenberg, 1994). Other genes of interest enriched in isolates from TN clade I include those encoding autotransporter adhesin SadA, which is associated with pathogenesis, and virulence protein MsgA, which is involved with survival within macrophage (Skyberg, Logue & Nolan, 2006). Mezal, Stefanova & Khan (2013) identified msgA virulence gene in 7 (out of 50) Salmonella ser. Javiana isolates, all of which were clinical.

TN clade II isolates had a core genome of 4,290 genes and an accessory genome of 322 genes (Table 1). The pan-GWAS revealed 22 loci to be significantly associated with inclusion in this clade (16 positively associated and 6 negatively associated) (Table 1 and Data S7). These loci consisted of 221 total genes (207 positively associated and 14 negatively associated) (Table 1). The positively associated loci are found in 17 distinct genomic regions and many of the genes over-represented in TN clade II (93 genes) were located in the 87.5 kb IncFII type plasmid identified in SRS2922480 (Table 4 and Fig. S3). Among the over-represented genes contained on this plasmid are ccdAB, which are part of a toxin/antitoxin system. This system contributes to stability of the plasmid through post-segregational killing (killing new cells that do not inherit a plasmid copy during cell division) (Van Melderen, 2001). Additionally, 72 over-represented genes were located in four predicted prophage regions and 33 in three other putative MGE regions (indicated by gene annotations, clustering of genes, and/or close proximity to predicted prophage regions; Table 4). One of the overrepresented genes, alpha-xylosidase, corresponds to a pathogenicity-related protein family (as identified by PathogenFinder). VFDB identified pefC and pefD, fimbrial adherence determinants, on the plasmid (contig 13); sinH, a nonfrimbrial adherence determinant; and pipB, a TTSS-2 translocated effector.

TN clade III isolates had a core genome of 4,115 genes and an accessory genome of 889 genes (Table 1). The pan-GWAS revealed 155 loci to be significantly associated with inclusion in this clade (101 positively associated and 54 negatively associated) (Table 1 and Data S7). These loci contained 332 total genes (238 positively associated and 94 negatively associated) (Table 1). The positively associated loci are found in 29 distinct genomic regions and the majority of the genes over-represented in TN clade III were located in four predicted prophage regions (88 genes) or nine other putative MGE regions (134 genes) (Table 4). Four of the overrepresented genes correspond to pathogenicity-related protein families (as identified by PathogenFinder: arginine/lysine/ornithine decarboxylase) and other hypothetical proteins.

Global population structure

The KSNP analysis of the diverse set of global clinical Salmonella ser. Javiana strains revealed three major clades (Fig. 3; Fig. S1). Major clade I contained 107 strains, including TN isolates (from TN clades I, II, and III and the five isolates that didn’t fall into the main clades). Major clade II contained 23 strains and major clade III contained 31 strains. Strains from major clades I and II belong to the 590 cgMLST eBurstGroup (ceBG) and strains from major clade III belong to the 204 ceBG; both of these are associated with this serovar (EnteroBase Team, 2018).

To further explore if Salmonella ser. Javiana is polyphyletic, we constructed minimal spanning trees based on cgMLST allele distances of all available Salmonella ser. Javiana strains and of other Salmonella serovars previously described as polyphyletic (Derby, Kentucky, Mississippi, Montevideo, Newport, Saintpaul, and Senftenberg (Achtman et al., 2012; Alikhan et al., 2018; Banerji et al., 2020; Cao et al., 2013; Sangal et al., 2010; Sévellec et al., 2018; Tang et al., 2019; Timme et al., 2013; Vosik et al., 2018; Zhang, Payne & Lan, 2019)) (Fig. 4; Fig. S2). The branch length between the two ceBG clusters on the ser. Javiana tree was 1,423 allelic differences. The branch lengths between ceBG clusters on the other trees ranged from 2,280-2,769 allelic differences for ser. Derby, 2,452-2,756 for ser. Kentucky, 929-2,850 for ser. Mississippi, 1,188-2,844 for ser. Montevideo, 957-2,636 for ser. Newport, 1,513-2,686 for ser. Saintpaul, and 1,844-2,743 for ser. Senftenberg.