Isolates were chosen for the study based on similarity or dissimilarity by PFGE to S. Heidelberg isolated from a large ground turkey-associated outbreak in 2011 (supplementary fig. S1, Supplementary Material online) (CDC 2011; Folster et al. 2012). Forty-four isolates of S. Heidelberg were included in the study (table 1). Isolates were obtained from animal (n = 9), retail meat (n = 27), and human clinical (n = 7) sources in the United States, and one of the isolates was obtained from an unknown sample source in Brazil. Four isolates were collected from 1982 to 1987; 38 of the isolates were collected from 2002 to 2011, and collection dates for two of the isolates were unknown. Five of the isolates were from the Salmonella Reference Collection A (SARA) (Beltran et al. 1991), and nine of the isolates were collected in the course of the investigation of a ground turkey-associated outbreak in 2011 (Folster et al. 2012). Table 1

Metadata of S. Heidelberg Isolates Used in This Study

Salmonella isolates were cultured on trypticase soy agar (TSA; Becton, Dickinson, NJ) and in trypticase soy broth (TSA; Becton) overnight at 37 °C. All isolates used for WGS were serotyped by conventional methods and tested for antimicrobial susceptibility according to the NARMS standard protocol as previously described (Zhao et al. 2008). Antimicrobial susceptibility testing, using a panel consisting of 15 antimicrobial agents (amikacin, ampicillin, amoxicillin-clavulanic acid, cefoxitin, ceftiofur, ceftriaxone, chloramphenicol, ciprofloxacin, gentamicin, kanamycin, nalidixic acid, streptomycin, sulfisoxazole, tetracycline, and trimethoprim-sulfamethoxazole) was performed according to the NARMS methodology with the Sensititre automated antimicrobial susceptibility system (Trek Diagnostic Systems, Westlake, OH). Isolate antimicrobial resistance was determined by comparison of MICs to values established MICs by the Clinical and Laboratory Standards Institute (CLSI). PFGE was performed according to the CDC PulseNet protocol (http://www.cdc.gov/pulsenet/pathogens/index.html, last accessed April 23, 2014). Genomic DNA of each strain was isolated from overnight cultures using DNeasy Blood and Tissue Kit (Qiagen, CA). All S. Heidelberg isolates were stored in TSB containing 15% glycerol at −80 °C.

We performed shotgun sequencing of the 44 S. Heidelberg using the Genome Sequencer FLX 454 (Roche, Branford, CT) and the GS FLX Titanium Sequencing Kit XLR70 according to the manufacturer’s protocol to generate an average genome coverage of 23×. De novo assemblies were performed using Roche’s Newbler software (v.2.6) with the resulting contigs being annotated using the NCBIs Prokaryotic Genomes Automatic Annotation Pipeline (Klimke et al. 2009).

Using DNA from the 2011 clinical outbreak isolate 41578, we also sought to close the genome using the Pacific Biosciences (PacBio) sequencing platform. Specifically, we prepared a single 10-kb library that was sequenced using the C₂ chemistry on eight single-molecule real-time cells with a 90-min collection protocol on the PacBio RS. The 10-kb continuous-long-read data were de novo assembled using the PacBio hierarchical genome assembly process/Quiver software package, followed by Minimus 2, and they were polished with Quiver.

For comparative analysis, in addition to the 44 isolates sequenced in this study, we included sequence data available for S. Heidelberg SL476 complete genome (NC_011083), two plasmids (NC_011081 and NC_011082), and the whole-genome shotgun sequence for S. Heidelberg SL486 (ABEL00000000) at the NCBI genome database.

Phylogenetic informative SNP sites (i.e., SNPs shared by two or more strains in the alignment) were identified by two different methods. The first involved mapping the 454 reads to the complete reference genome of S. Heidelberg strain SL476 using Roche Newbler software GS Reference Mapper (v.2.8.). SNP sites were then found using the SNP calling function of GS Reference Mapper. SNP positions were defined, where one or more isolates differed from strain SL476 with coverage ≥10× and with ≥95% of the reads containing the SNP after having filtered out the SNPs from homopolymer artifacts, followed by a custom pipeline to construct a SNP matrix (Allard et al. 2013). S. Heidelberg SL486 could not be included in this analysis as its raw sequence data were not available. The second SNP detection method was a reference free k-mer-based approach implemented in the program kSNP (Gardner and Slezak 2010). This is a collection of Perl scripts that aggregate the results from Jellyfish v1.1.3 (used for k-mer counting) (Marcais and Kingsford 2011) and MUMmer v3.22 (used to align k-mers and detect variable positions) (Kurtz et al. 2004). Analyses with kSNP were based on a k-mer length of 25.

To construct evolutionary relationships among the isolates, the maximum-likelihood (ML) method was implemented in the Genetic Algorithm for Rapid Likelihood Inference (GARLI) software (Zwickl 2006b). GARLI analyses were performed using a web service (Bazinet and Cummings 2011) that uses a special programming library and associated tools (Bazinet et al. 2007). All ML trees were constructed with GTR+I+G nucleotide substitution model. The 100 replicate runs of the nonbootstrapped data set were conducted to identify the most probable phylogeny based on our observed SNP matrix. We determined branch support by performing 1,000 bootstrap replicates.

We elucidated the relationships among isolates using the Bayesian clustering method implemented in STRUCTURE v2.3.2 (Pritchard et al. 2000; Falush et al. 2003). The analyses were based on the kSNP matrix, and we ran ten replicate analyses at K = 2–9 under the admixture model with correlated allele frequencies being performed and visualized with DISTRUCT v1.1. (Rosenberg 2004). Each independent run consisted of 50,000 generations serving as burnin followed by 100,000 generations. The ΔK statistic (Evanno et al. 2005) was used to identify the k value that best fits the data. Pairwise genetic distances, calculated as the number of nucleotide differences, were generated in MEGA5 (Tamura et al. 2007). Core genes were identified using the UCLUST algorithm (Edgar 2010) using 95% sequence identity as the cutoff. Alignments of core genes were accomplished using MUSCLE with default settings; these alignments were then used to identify variable core genes.

Prophages were identified using PHAST (Zhou et al. 2011). Sequences for intact phages were extracted from the original contig, in which they were found and mapped against the entire assembled data to determine whether the phage was present in an isolate at multiple sites on different contigs. Contigs that could not be aligned to the reference genome were evaluated using BLAST to identify plasmids. The plasmids were analyzed by comparative analysis using MAUVE algorithm (Darling et al. 2010) and with the comparative analysis tools of RAST (Aziz et al. 2008). For plasmid closure, a concatenated sequence was generated from a 500-bp sequence from each end of the contig. This artificial sequence from the single contig of the plasmid in study was used as a reference to map all the 454 raw reads, using runMapping from Newbler. If there were mapped reads that covered the conjunction point in the reference, the contig was a closed plasmid.

To determine the incompatibility (Inc) groups for plasmids, we used BLAST to find sequences described by Johnson and Nolan (2009) for specific Inc groups that would produce theoretical PCR amplicons for known Inc group sequences. Resistance and virulence genes were identified by mapping sequence data available at an in-house database, consisting of 1,379 resistance genes, and 107 previously characterized virulence genes (84 pathogenicity and 23 fimbrial markers; Huehn et al. 2009). Using a presence/absence matrix of genes conferring resistance to antibiotics and disinfectant agents, we constructed a similarity tree based on binary distances under neighbor-joining algorithm for tree construction; topological support was assessed based on 100 bootstrap replicates. The analyses were performed using the R package ape (Paradis et al. 2004; Liu et al. 2011).

Pathogenicity was evaluated using the Caenorhabditis elegans survival assay. Caenorhabditis elegans strain SS104 glp-4 (bn2) was acquired from the Caenorhabditis Genetics Center. Worm cultures were maintained at 16 °C, which is the permissive temperature for this temperature-sensitive sterile mutant of C. elegans. Strains were cultured in C. elegans habitation media (CeHM) in tissue culture flasks on a platform shaker (Sprando et al. 2009). Nematodes were bleached (0.5 M NaOH, 1% hypochlorite) to collect eggs, which were incubated in M9 media for 24 h to bring them to synchronized L1 stage and then transferred to CeHM. To produce synchronized L4 stage worms, L1 worms were grown for an additional 72 h in CeHM.

Pathogen lawns for survival assays with test strains, as well as the food bacteria OP50, were prepared by inoculating Nematode Growth Media (in 6-cm Petri plates) with 50 μl of an overnight bacterial culture. Plates were incubated overnight at room temperature before worms were added. We used 60–80 synchronized L4 worms for each treatment group. Worms were scored every 24 h for survival (Aballay et al. 2000).

Animal survival was plotted using Kaplan–Meier survival curves and analyzed by log rank test using GraphPad Prism (GraphPad Software, Inc., La Jolla, CA). Survival assays were repeated at least two times. Survival curves resulting in P values < 0.05 relative to control were considered significantly different.

Two representative isolates from the 2011 ground turkey-associated outbreak were previously reported to be resistant to ampicillin, gentamicin, streptomycin, and tetracycline (Folster et al. 2012). Antimicrobial susceptibility testing of seven additional representative isolates showed resistance to ampicillin, gentamicin, and tetracycline among six; and streptomycin, tetracycline, and kanamycin resistance in one of the isolates. The susceptibility and PFGE (Xbal and BlnI) results for all nine of the outbreak representative isolates included in this study and the reference isolates (table 1) are shown in supplementary figure S1, Supplementary Material online. The nine outbreak associated isolates exhibited indistinguishable, or nearly indistinguishable (one band difference), Xbal PFGE patterns JF6X01.0032 or JF6X01.0058, and BlnI pattern JF6A26.0076 and one food isolate had Xbal PFGE pattern JF6X01.0058 and BlnI pattern JF6A26.0017. Of the 35 isolates not associated with the outbreak, 15 had composite Xbal/BlnI PFGE patterns indistinguishable from ground turkey outbreak isolates whereas 20 isolates had composite XbaI/BlnI PFGE patterns ≤98% similar to isolates from the 2011 outbreak (supplementary fig. S1, Supplementary Material online). Nine reference isolates showed the same PFGE pattern, as well as the same antibiotic resistance profile, noted for most of the outbreak isolates (supplementary fig. S1, Supplementary Material online).

Shotgun sequencing produced 44 draft genomes with an average coverage of 23× and a minimum of 4,700 genes. The genomic size of the 44 S. Heidelberg isolates (including extrachromosomal DNA) ranged from 4.671 to 5.130 Mb (fig. 1). Differences in genomic size were due primarily to the presence or absence of mobile genetic elements, including phages and plasmids, especially the carriage of IncHI2 or IncA/C plasmids (fig. 1). Next to the genome size variation, figure 1 shows the estimated N50 sizes within S. Heidelberg draft genome sequences. The estimated N50 value is a rough estimate of the quality and coverage of the draft genomes, and it represents the average contig size after assembly with the Newbler software. Fig. 1.—

The number of assembled bases (Mb) and N50 contig size (kb) for each sequenced S. Heidelberg isolate. Samples are colored according to the presence of antimicrobial resistance plasmids. No antimicrobial resistance plasmid, filled triangles; antimicrobial resistance Inc-I plasmid, filled diamonds; antimicrobial resistance Inc-H1/2 plasmid, filled circles; antimicrobial resistance plasmid Inc-I and H1/2, filled circles with borders; antimicrobial resistance plasmid Inc-AC, filled squares; antimicrobial resistance plasmid IncI and Inc-AC, filled squares with borders.

Phylogenetically informative SNP sites were identified using two independent methods. Figure 2 shows a ML tree based on SNP analysis of 44 S. Heidelberg isolates mapped to strain SL476. Among the collective 40,716 variable SNPs identified, 12,187 were determined to be “informative” (i.e., SNPs shared by at least two isolates). The ML tree shown in figure 2 clearly demonstrates that S. Heidelberg formed a monophyletic group distinct from the outgroup comprised Salmonella enterica supsp. enterica serovar Newport (S. Newport) and S. Typhimurium. The SNP divergence between S. Heidelberg and the two outgroups, S. Newport and S. Typhimurium, is 24,724 and 24,305 SNPs, respectively. The resultant tree shows two more interesting points. First, isolates having a similar Xbal PFGE pattern (supplementary fig. S1, Supplementary Material online) clustered together. Second, there is a tendency for isolates to cluster based on date of collection, such as noted for isolates from 2003, 2008, and 2011. In addition, the 2011 clinical and food isolates that were originally thought to be associated with the ground turkey outbreak isolated in different states clustered together with 80% bootstrap support (fig. 2), which could indicate that they are all members of a clonal community likely derived from the same source. Consistent with the trend for temporally isolated strains to cluster closely together, one isolate, N29341, that was isolated in Minnesota in 2011 from ground turkey by NARMS, clustered tightly together with the outbreak isolates received from CDC, whereas all other reference isolates having the same PFGE pattern not associated with the ground turkey outbreak are distinct from the isolates of the 2011 outbreak (fig. 2). These results suggest that N29341 might be epidemiologically related to the 2011 outbreaks strains, as Minnesota was one of the states involved in that outbreak. Fig. 2.—

ML tree based on SNP analysis of 44 S. Heidelberg isolates and previously reported S. Heidelberg genome sequences SL476 (12). A total of 40,716 variable SNPs with 12,187 being informative were found using GS Reference Mapper followed by a custom pipeline. The ML tree was generated in GARLI v.2.0 (Zwickl 2006a) under the GTR+Γ model of nucleotide evolution and visualized using Figtree v1.3.1. Parameter space was searched for the best tree with simultaneous estimation for model parameters using a ML search. The best tree was identified from 100 runs on the nonbootstrapped data set. Measures of clade confidence are reported below each node in the form of bootstrap values (1,000 iterations). Bootstrap values <70% were not shown. The tree was rooted using S. Newport 637564 and S. Typhimurium AZ057. The taxa of source for each isolate, geographic location, and date were mapped onto the tree. Prophage observations are further depicted on the tree using colored bars, shown on the right.

To further investigate the relationships among outbreak isolates, we constructed the Distruct plot for group 3. Within that group, there is support for four subgroups; these four groups also are found in the phylogenetic analysis, each with a ≥75% bootstrap support. Pairwise SNP variation between the four sublineages also support the phylogenetic partitions described above (table 2 and fig. 3). Mean divergence among the four subclades ranged from 36 to 153 nt differences, whereas the mean intragroup nucleotide differences ranged from 15 to 30 (table 2). The 2008 strains seem to have quite a unique SNP profile when compared with the remaining isolates of group 3. The ML tree, the Distruct plot, and the pairwise SNP variation show the evolutionary changes occurring among isolates from 2003 to 2009 and 2011, revealing the apparent emergence of the genetically unique S. Heidelberg strain that was responsible for the 2011 outbreak.

Table 3 lists variable genes having unique nucleotide substitutions which defined specific sublineages that could be used to separate them from other highly clonal S. Heidelberg isolates in group 3. Table 3 also lists their SNP change and, if applicable, the gene carrying the SNPs and the resulting amino acid change. Twenty-three unique SNPs, with 14 being nonsynonymous (i.e., causing an amino acid change) SNPs, were found to be unique to group 3 and clustering together isolates with XbaI pattern (JF6X01.0032, JF6X01.0034, JF6X01.0058). These nucleotide substitutions were found on 20 different genes that could be useful for subtyping of strains belonging to group 3 (table 3 and fig. 3). Furthermore, 24 unique SNPs (11 being nonsynonymous) were found at the node that characterizes the 2008 ground turkey isolates, which seemed to have their own unique SNP profile distinct from that of the remaining group 3 isolates. Twelve of the 24 SNPs are located on the insertion sequence (IS) 2 element with five SNPs being within the IS2 TnpA transposase and seven SNPs on the IS2 TnpB transposase genes (table 3 and fig. 3). Table 3

Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study

Importantly, two nonsynonymous SNPs were found to be only present in all isolates associated with the 2011 outbreak. These two SNPs were located within two different variable core genes: a putative hydrolase with amino acid change P653Q, and the gluconate transporter GntP gene with amino acid change L187V (table 3 and fig. 3). Five SNPs, with four being nonsynonymous, were unique to isolates 41563, 41566, 41567, 41573, 42578, 41579, and N29341 (associated to the outbreak) and not observed in 41565 (table 3 and fig. 3). From these five SNPs, three were in plasmid sequences, whereas two were found on variable core genes: the flagellin CDS resulting in an amino acid change from threonine to serine at position 282 (ΔT282S), and within a putative hypothetical protein (possibly functioning as H+ gluconate symporter related to permeases) producing an amino acid change from glycine to serine at position 358 (ΔG358S) (table 3).

At a minimum, 273 genes involved in a variety of different functions, such as DNA replication and repair, cell division, transcription, metabolism, and virulence, vary among the S. Heidelberg isolates observed in clinical, animals, and/or food samples, including the isolates associated with the 2011 outbreak. The variable core genes, including the number of SNPs and haplotypes are listed in supplementary table S1, Supplementary Material online. Figure 4 shows a histogram comparing strains based on: 1) the number of SNPs according to the number of core genes and 2) the haplotype diversity according to the number of core genes. Only a few mutation “hotspot” genes, such as tailspike (35 SNPs), scaffolding protein (22 SNPs), DNA polymerase V subunit umuD (20 SNPs), phage lysozyme (18 SNPs), gifsy-2 prophage protein (15 SNPs), terminase small subunit (14 SNPs), outer membrane lipoprotein Blc (7 SNPs), and alpha amylase (6 SNPs), were found to have a disproportionate number of SNPs that could be useful for rapid subtyping the serovar Heidelberg and/or as targets of resequencing. Interestingly, although these genes have several SNPs, the haplotype diversity is fairly low, with only a few haplotypes (2–4) among 46 isolates being identified. This fact suggests that, even though there are genes with up to 35 SNPs, there are only two to four unique sequences among S. Heidelberg that are congruent with our tree topology for three different groups of S. Heidelberg. Fig. 4.—

(A) Histogram showing the number of SNPs per core genes. (B) Histogram showing haplotype diversity for all variable, core genes.

In the course of performing this study, we demonstrated that it was possible to rapidly determine the complete sequence for the chromosome and three plasmids from a clinical isolate associated with the 2011 outbreak using the Pacific Biosciences (PacBio) system. A single contig of 4,793,479 bp (GC content 52.2%) representing the complete chromosome and three contigs of 4,773 bp, 35,297 bp, and 117,929 bp representing three plasmids were generated. The chromosome contains 12 different Salmonella pathogenicity islands and four prophages (prophage P22, lambda-like prophage Gifsy-2, prophage P4, and prophage P2-like; fig. 5). The largest plasmid pSEEH1578_01 is a resistance plasmid with incompatibility group IncI1 encoding the resistance to gentamicin (aacC), streptomycin (aadA1), tetracycline (tetA, tetR(A)), and ampicillin (bla_tem1) (fig. 5). The plasmid pSEEH1578_02 is a VirB/D4 plasmid that carries the type IV secretion system (T4SS; fig. 5), and plasmid pSEEH1578_03 is a mobilization plasmid that carries genes capable of promoting plasmid mobilization such as mbeA, mbeB, mbeC, and mbeD. Fig. 5.—

Chromosome and plasmids features of a clinical S. Heidelberg isolate 41578. The circular map was drawn using dnaplotter. Different features are shown in different colored bars. (A) Chromosome: The coding sequence are shown in dark blue, rRNA is shown in green, tRNA is shown in yellow, prophages are shown in blue, and Salmonella pathogenicity islands are shown in red. Track 7 represents the GC content while Track 8 shows the GC skew [(G − C)/G + C]. (B) IncI1 antimicrobial resistance plasmid: The coding sequences are shown in dark blue and resistance genes are shown in red. Track 5 shows the GC skew [(G − C)/G + C]. Regions of GC content above average of the plasmid are drawn outside the ring in yellow, whereas regions below average are inside the ring in purple. (C) VirB/D4 virulence plasmid: The coding sequences are shown in dark blue, genes that carry the T4SS are shown in red, genes responsible for plasmid stability, and replication is shown in green. Track 6 shows the GC skew [(G − C)/G + C]. Regions of GC content above average of the plasmid are drawn outside the ring in yellow, whereas regions below average are inside the ring in purple.

PHAST analysis identified six different phages among the S. Heidelberg isolates. Only those phages that PHAST determined as being intact were further analyzed. All 44 isolates contained prophage P22 and lambda-like prophage Gifsy-2. Group 1, which was composed of three historical isolates (SARA 39, SARA 37, and 82-2052), also contained prophage Fels-2, not observed in any of the group 2 or group 3 isolates. Additionally, SARA 37 and 82-2053 also contained prophage ST64B (table 1 and fig. 2).

Heidelberg isolates belonging to group 3, which includes the 2011 outbreak isolates, contained prophage P4 and P2-like, neither of which was seen in any of the isolates from group 1 and 2, excepting isolate N19992, a group 2 isolate that was confirmed to contain the P2-like prophage (table 1 and fig. 2). The phage presence and absence matrix shows correlation with the tree phylogeny shown in figure 2.

From the 107 known and recognized Salmonella-associated virulence genes, we identified 62 pathogenicity and 13 fimbrial markers that were highly conserved among the 44 S. Heidelberg isolates (supplementary table S2, Supplementary Material online). One outer membrane fimbrial usher gene, safC, was only found in the three historical isolates (SARA 39, SARA 37, and 82-2052), all belonging to group 1. All other detected virulence genes, observed to be located on the Salmonella pathogenicity islands (SPI-1, SPI-2, SPI-3, SPI-4, and SPI-5), the lambda-like prophage Gifsy-2, and/or islets, were present in all S. Heidelberg isolates, suggesting that all isolates are pathogenic (supplementary table S2, Supplementary Material online). For example, all isolates carried a chromosomal msgA gene, which is essential to Salmonella mouse virulence (Gunn et al. 1995), and is common in S. Typhimurium (Huehn et al. 2010). Other genes present among all isolates encoded proteins involved in type III secretion system (T3SS) and adhesins. Remarkably, the sopE1 virulence gene, involved in the translocation of effector proteins into host cells, was found in all of our strains. Moreover, although sopE1 is known to be carried by P2-like phages, in this study we also identified this gene on lambda-like phage Gifsy-2, which is present in all S. Heidelberg isolates. Consequently, all S. Heidelberg isolates have the sopE1 virulence gene. Group 3 isolates, which also contain the P2-like phage, have the sopE1 gene twice in their chromosome; one gene is carried by lambda-like phage Gifsy-2 and the other one is located on the P2-like phage. Analysis confirmed that these two genes shared a high degree of similarity having no more than eight mutations between the two different sopE1 alleles. In contrast, virulence determinants known to be associated with SPI-7, prophage Gifsy-1, Gifsy-3, Fels-1, and Salmonella plasmid virulence (spv) were not identified in any of the 44 isolates examined in this study, confirming sequencing results that none of the strains carried any of these mobile genetic elements.

Several plasmids were found among the 44 S. Heidelberg isolates sequenced in this study (table 1). Electrophoresis of samples obtained using a basic alkaline plasmid purification procedure demonstrated the presence of endogenous plasmids that had not integrated into the host genome (data not shown). The plasmid sizes varied from 3.7 to over 200 kb. Twenty-one S. Heidelberg isolates, including all five human isolates and two ground turkey isolates associated with the 2011 outbreak, contained a plasmid ranging in size from 33 to 40 kb, which carried genes (virB1, virB2, virB3-4, virB5, virB6, virB7, virB8, virB9, virB10, virB11, virD2, and virD4) associated with the VirB/D4 T4SS (table 1). Three specific roles for the bacterial T4SS have been identified: 1) to facilitate translocation of DNA via a conjugative mechanism to recipient cells; 2) to facilitate translocation of effector molecules, such as proteins, to eukaryotic target cells; and 3) to function in DNA uptake or dissemination from or into the environment milieu (Christie et al. 2005; Alvarez-Martinez and Christie 2009).

To further study the diversity among all 22 identified VB/D4 plasmids, we identified SNPs that might be useful in differentiating these plasmids. The SNP matrix from the 22 plasmids and the reference plasmid built using the k-mer strategy of kSNP consisted of 338 informative SNP positions. Figure 6 shows a ML tree based on SNP analysis. The phylogeny placed the VB/D4 plasmids into three well-supported clades (95%, 100%, and 100% bootstrap support). Clade I contained two plasmids isolated from the historical isolates, SARA 37 and 82-2052, belonging to group 1 on the ML tree shown in figure 3. Clade III is composed of ten plasmids, including the reference plasmid pSARA30 and 9 isolates belonging to group 2 in figure 4. Clade II contains 11 plasmids, including 7 plasmids from 7 isolates associated with the 2011 outbreak. Clade II can be further separated into two subclades, one of which contains nine plasmids derived from isolates belonging to group 3 in figure 4, while subclade 2 contains one clinical isolate and isolate N653, the unique isolate found to carry two VB/D4 plasmids (the other plasmid partitioned into clade III). Moreover, using the primers designed by Johnson et al. (2012) for the relaxase gene (taxC) (VirD2 component) to screen for the four different IncX subgroups, we determined that the VB/D4 plasmids from clade I, subclade 1 from clade II, and clade III belonged to the IncX1 plasmid subgroup, whereas the two plasmids from clade II subclade 2 were deemed to be members of the IncX4 subgroup. Fig. 6.—

ML tree based for the 22 identified VB/D4 plasmids, including the reference S. Heidelberg plasmid pSARA30. A total of 338 SNPs with all being informative, were found based on k-mer analysis using kSNP. The numbers of unambiguous substitutions that mapped to the tree only once are given above each node. ML trees were generated as described in figure 3. Bootstrap values (1,000 iterations) are reported below each node.

To answer the question whether the isolates that contain the VirB/D4 plasmid are more virulent, we performed a survival assay using a C. elegans system to compare the virulence potential of our strains. This system has been used successfully as an invertebrate host model to assess the virulence determinants of human pathogens, such as Vibrio cholerae (Sahu et al. 2012) and S. Typhimurium (Aballay et al. 2000). We tested six different S. Heidelberg isolates (four ground turkey isolates, one swine isolate, and one clinical outbreak isolate) having different plasmid profiles. Strain 29169 singly carries a small mobilization plasmid; strain SARA 31 carries only the VirB/D4 T4SS plasmid; strain 418 carries a small mobilization plasmid and an IncA/C plasmid; strain 41565 carries a small mobilization plasmid and an IncHI2 plasmid; strain 41578 has a small mobilization plasmid, the VirB/D4 T4SS plasmid, and an IncI1 plasmid; and N30678 has four plasmids: the small mobilization plasmid, the VirB/D4 T4SS plasmid, the IncI1, and the IncHI1. Except SARA 31, in which phage P4 and P2-like were not present, all other isolates had the same phage profile.

The results, which showed that S. Heidelberg isolates are significantly (P < 0.0001) more pathogenic than the E. coli OP 50 control strain, failed to demonstrate any appreciable difference in pathogenicity between strains carrying the VirB/D4 T4SS plasmid and those strains that did not (fig. 7). However, we also identified some VirB/D4 T4SS components on the resistance plasmids IncA/C, IncHI2, and IncI1. Interestingly, isolate 29169, which did not carry any large plasmid and any of the T4SS components, was found to be significantly (P < 0.0001) less pathogenic using this assay than any of the other isolates evaluated. Isolate 41565, which carried the small mobility plasmid as well as additional IncHI2 plasmid, demonstrated the second weakest pathogenic potential. Isolate SARA 31, which carried the VB/D4 T4SS plasmid but no other resistance plasmid, showed a survival curve profile similar to that of isolate N418, which did not have the VB/D4 T4SS plasmid but carried the IncA/C plasmid. Fig. 7.—

Caenorhabditis elegans survival data from six S. Heidelberg isolates. The figure shows that the six S. Heidelberg isolates (29169, SARA 31, 418, 41565, 41578, and N30678) are significantly (P < 0.0001) more pathogenic than the Escherichia coli OP 50 control strain. Further the figures show that isolates not carrying any T4SS components tend to be significant less pathogenic than those isolates that do carry them.

Analysis of WGS data showed 27 plasmid-associated antibiotic resistance genes among the 44 isolates and included genes expected to produce resistance to aminoglycosides, beta-lactams, phenicols, folate pathway inhibitors, tetracyclines, and cephems. These genes were not detected in the eight pan susceptible isolates that also did not carry any of the antibiotic resistance plasmids. Resistance phenotypes correlated with genotypes in all strains examined (table 4). However, some strains with the same genotype exhibited different antimicrobial resistance patterns. For example, both strains N18440 and N18413 contain the aadA1 gene but only isolate N18440 showed resistance to streptomycin. In this study, streptomycin resistance was defined using the MIC value >32 mg/l determined by CLSI. However, studies have shown that strains harboring an aadA gene cassette can have MICs of 16 mg/l, which would classify these strains as streptomycin-susceptible (Sunde and Norström 2005). Table 4

Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study

Among the 36 S. Heidelberg isolates resistant to one or more antibiotics, 24 contained sequences indicative of presence of an IncI1 plasmid, 5 had an IncHI2 plasmid, and 1 carried an IncA/C plasmid. Three isolates contained two plasmids, one with IncI1 + IncHI2, and two isolates each having IncI1 + IncA/C replicons (table 4). Besides isolate 41565, which carries an IncHI2 plasmid, all Heidelberg isolates associated with 2011 outbreak carry an IncI1 plasmid, which was detected the most overall among our strains. The IncI1 plasmid carries several antibiotic resistance genes and components of the T4SS. The IncHI2 plasmid was found in eight isolates and carries several antibiotic resistance genes and genes conferring resistance to quaternary ammonium compounds, as well as genes encoding the conjugative T4SS. Furthermore, unique to the IncHI2 plasmid is the presence of genes encoding resistance to heavy metals, including tellurium, copper, cadmium, mercury, and silver.

The three isolates resistant to the most (ten or more) antibiotics carried the IncA/C multidrug resistance plasmid. In addition to carrying the genes encoding resistance to different antibiotics and quaternary ammonium compounds, this plasmid also carries the conjugative T4SS transfer system. All three IncA/C plasmids carried the resistance element bla_CMY-2-blc-sugE2, which was only found in one other isolate, N189, on an IncI1 plasmid.

SARA 35, which was only resistant to ampicillin, did not carry a multidrug resistance plasmid but, instead, carried a smaller plasmid (∼8.0 kb) having only a single determinant that correlated with resistance to β-lactams (tnpA-bla_TEM-1). This plasmid shared a high degree of sequence similarity to S. Typhimurium pAnkS (accession number: DQ916413), which also carried the transposon resistance element (tnpA-bla_TEM-1). The Inc group for this plasmid could not be determined from the sequence.

Figure 8 shows a neighbor-joining tree generated from a presence/absence matrix of resistance genes for the 30 group 3 Heidelberg isolates having a common PFGE profile. As expected, resistance tended to segregate according to the particular Inc group plasmid they carried and not by clonality or evolutionary relationship. This observation indicated that the presence of specific resistance genes correlates with the particular Inc group of the plasmids. For example, the resistance genes (aphA1, bla_CMY2, sul2, cmlA1, floR) were only detected in isolates that carry an IncA/C plasmid, whereas only isolates carrying an IncHI2 plasmid were found to contain the resistance genes aph, sph, tetC, tetD, and ble. Similarly, the aacC gene only appears to be present in those isolates identified as carrying an IncI1 plasmid. The association of different resistance genes with specific Inc groups is well known. Fig. 8.—

Absence and presence tree of resistance genes among S. Heidelberg isolates associated with paraphyletic group 3. It is a similarity tree based on binary distances under neighbor-joining algorithm for tree construction. The resistance genes and incompatibility group of the resistance plasmid are mapped to the tree.

Plasmids below 8 kb did not carry any antimicrobial resistance or virulence-associated genes. Instead, these smaller plasmids tended to carry genes capable of promoting plasmid mobilization, such as the mbeA, mbeB, mbeC, mbeD genes. S. Heidelberg isolates N189, N653, and N4403, which were originally isolated from chicken breast, carry a small, 3,772 bp plasmid that is identical to the S. Heidelberg SL476 plasmid pSL476_3 (GenBank accession number: CP001119) (Fricke et al. 2011). All group 3 isolates (fig. 3), except isolate 24393, have similar XbaI patterns (JF6X01.0032, JF6X01.0033, JF6X01.0034, JF6X01.0058), and carry a 4,473 bp plasmid that shares high sequence similarity with the S. Heidelberg plasmid pSH1148_4.8 (GenBank accession number: JX494965). SARA 35, which was isolated in Brazil, was found to carry a 6,647 bp ColE1 plasmid (pCFSAN000443_01) that is similar to E. coli plasmid ColE1 (GenBank accession number: J01566).