Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Comparative Genomic Analysis and Virulence Differences in Closely Related Salmonella enterica Serotype Heidelberg Isolates from Humans, Retail Meats, and Animals

Comparative Genomic Analysis and Virulence Differences in Closely Related Salmonella enterica... Salmonella enterica subsp. enterica serovar Heidelberg (S. Heidelberg) is one of the top serovars causing human salmonellosis. Recently, an antibiotic-resistant strain of this serovar was implicated in a large 2011 multistate outbreak resulting from consumption of contaminated ground turkey that involved 136 confirmed cases, with one death. In this study, we assessed the evolutionary diversity of 44 S. Heidelberg isolates using whole-genome sequencing (WGS) generated by the 454 GS FLX (Roche) platform. The isolates, including 30 with nearly indistinguishable (one band difference) Xbal pulsed-field gel electrophoresis patterns (JF6X01.0032, JF6X01.0058), were collected from various sources between 1982 and 2011 and included nine isolates associated with the 2011 outbreak. Additionally, we determined the complete sequence for the chromosome and three plasmids from a clinical isolate associated with the 2011 outbreak using the Pacific Biosciences (PacBio) system. Using single-nucleotide polymorphism (SNP) analyses, we were able to distinguish highly clonal isolates, including strains isolated at different times in the same year. The isolates from the recent 2011 outbreak clustered together with a mean SNP variation of only 17 SNPs. The S. Heidelberg isolates carried a variety of phages, such as prophage P22, P4, lambda-like prophage Gifsy-2, and the P2-like phage which carries the sopE1 gene, virulence genes including 62 pathogenicity, and 13 fimbrial markers and resistance plasmids of the incompatibility (Inc)I1, IncA/C, and IncHI2 groups. Twenty-one strains contained an IncX plasmid carrying a type IV secretion system. On the basis of the recent and historical isolates used in this study, our results demonstrated that, in addition to providing detailed genetic information for the isolates, WGS can identify SNP targets that can be utilized for differentiating highly clonal S. Heidelberg isolates. Key words: outbreak, antimicrobial resistance, plasmid, SNP analysis, trace-back. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution 2014. This work is written by US Government employees and is in the public domain in the US. 1046 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE to determine whether the outbreak isolates belonged to a Introduction new strain with potentially higher pathogenicity or whether Disease caused by nontyphoidal serotypes of Salmonella these isolates were typical members of a common PFGE type enterica is the leading cause of food-related death in the present in retail poultry. In addition, the data show the role of United States (Scallan et al. 2011). Salmonella enterica sup- transmissible mobile genetic elements in the evolution of vir- species (supsp.) enterica serovar Heidelberg (S. Heidelberg) has ulence and resistance among S. Heidelberg. been among the most frequently isolated serovars in clinical cases of salmonellosis, causing an estimated 84,000 illnesses Materials and Methods in the United States annually (Foley et al. 2011; Han et al. 2011). It is frequently isolated from poultry and poultry Bacterial Strains, Growth Condition, and Characterization meat (Zhao et al. 2008). S. Heidelberg was the seventh Isolates were chosen for the study based on similarity or dis- most common serovar isolated from humans in 2011 (CDC similarity by PFGE to S. Heidelberg isolated from a large 2013) and was found in chicken breast, ground turkey, and ground turkey-associated outbreak in 2011 (supplementary eggs, the main sources of S. Heidelberg infections (Chittick fig. S1, Supplementary Material online) (CDC 2011; Folster et al. 2006; NARMS 2010). S. Heidelberg also tends to be et al. 2012). Forty-four isolates of S. Heidelberg were included associated with invasive diseases such as septicemia and myo- in the study (table 1). Isolates were obtained from animal carditis (Wilmshurst and Sutcliffe 1995). After S. enterica (n¼ 9), retail meat (n¼ 27), and human clinical (n¼ 7) sources supsp. enterica serovar Typhimurium (S. Typhimurium), in the United States, and one of the isolates was obtained S. Heidelberg is the serovar of Salmonella most often associ- from an unknown sample source in Brazil. Four isolates ated with Salmonella-related deaths in the United States were collected from 1982 to 1987; 38 of the isolates were (Kennedy et al. 2004; Patchanee et al. 2008; Crump et al. collected from 2002 to 2011, and collection dates for two of 2011). The National Antimicrobial Resistance Monitoring the isolates were unknown. Five of the isolates were from the System (NARMS), which is responsible for monitoring antimi- Salmonella Reference Collection A (SARA) (Beltran et al. crobial resistance in Salmonella estimated that 65% of the 1991), and nine of the isolates were collected in the course S. Heidelberg isolates from ground turkey in 2010 were resis- of the investigation of a ground turkey-associated outbreak in tant to 3 antimicrobial classes. Presently, the antimicrobial 2011 (Folsteretal. 2012). agents to which this serovar is most commonly resistant are Salmonella isolates were cultured on trypticase soy agar ceftriaxone (a drug of choice for treatment), along with resis- (TSA; Becton, Dickinson, NJ) and in trypticase soy broth tance to streptomycin, tetracycline, sulfamethoxazole, chlor- (TSA; Becton) overnight at 37 C. All isolates used for WGS amphenicol, and trimethoprim-sulfamethoxazole (NARMS were serotyped by conventional methods and tested for anti- 2010). microbial susceptibility according to the NARMS standard pro- The Centers for Disease Control and Prevention (CDC) in- tocol as previously described (Zhao et al. 2008). Antimicrobial vestigated a multistate (34 states) outbreak of antimicrobial- susceptibility testing, using a panel consisting of 15 antimicro- resistant S. Heidelberg infections comprised of 136 confirmed bial agents (amikacin, ampicillin, amoxicillin-clavulanic acid, cases between February 27 and September 13, 2011. Among cefoxitin, ceftiofur, ceftriaxone, chloramphenicol, ciprofloxa- the 94 case patients for which there was available informa- cin, gentamicin, kanamycin, nalidixic acid, streptomycin, sulfi- tion, 37 (39%) had been hospitalized and one patient died. soxazole, tetracycline, and trimethoprim-sulfamethoxazole) Collaborative investigative efforts by state and federal officials was performed according to the NARMS methodology with implicated ground turkey as the source of this outbreak, and the Sensititre automated antimicrobial susceptibility system as a result, 36 million pounds of ground turkey meat were (Trek Diagnostic Systems, Westlake, OH). Isolate antimicrobial recalled (CDC 2011; Folsteretal. 2012). resistance was determined by comparison of MICs to values Recently, our investigative capabilities have been greatly established MICs by the Clinical and Laboratory Standards enhanced with the development and increasing feasibility of Institute (CLSI). PFGE was performed according to the CDC whole-genome sequencing (WGS) as a molecular epidemio- PulseNet protocol (http://www.cdc.gov/pulsenet/pathogens/ logical tool to complement current foodborne outbreak inves- index.html, last accessed April 23, 2014). Genomic DNA of tigation techniques. WGS is of particular interest because it each strain was isolated from overnight cultures using DNeasy provides definitive data for distinguishing outbreak isolates Blood and Tissue Kit (Qiagen, CA). All S. Heidelberg isolates were stored in TSB containing 15% glycerol at 80 C. from nonoutbreak isolates in common and highly clonal pop- ulations (Allard et al. 2012, 2013). In this study, we sought to Genome Sequencing, Assembly, and Annotation determine how effectively WGS and single nucleotide poly- morphisms (SNPs) analysis would differentiate outbreak iso- We performed shotgun sequencing of the 44 S. Heidelberg lates of S. Heidelberg from nonoutbreak isolates that share the using the Genome Sequencer FLX 454 (Roche, Branford, CT) same Xbaland BlnI pulsed-field gel electrophoresis (PFGE) pat- and the GS FLX Titanium Sequencing Kit XLR70 according to terns. Using WGS data and virulence assays, we also wanted the manufacturer’s protocol to generate an average genome Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1047 Hoffmann et al. GBE 1048 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Table 1 Metadata of S. Heidelberg Isolates Used in This Study Salmonella enterica susp. Source Source Source Isolated Accession Closed Plasmid BioProject Plasmid Profile Phage Profile enterica serovar and strain Country State Type Date Heidelberg str. 20752 USA NC Animal 2002 AMNR00000000 78501 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 21381 USA CT Food 2002 AMMZ00000000 80505 mob1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24359 USA MO Animal 2003 AMNP00000000 78497 mob1, IncI1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24388 USA MO Animal 2003 AMNQ00000000 78499 mob1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24390 USA MO Animal 2003 AMMY00000000 78475 mob1, IncI P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24391 USA MO Animal 2003 AMNO00000000 78495 mob1, IncI P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24393 USA MO Animal 2003 AMNM00000000 78491 IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 29169 USA CT Food 2003 APIX00000000 80527 mob1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 32507 USA GA Food 2003 AMNL00000000 80525 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41563 USA OR Human 2011 AJGX00000000 pCFSAN000312_01 78487 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41565 USA WA Food 2011 AJHA00000000 80521 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41566 USA WA Food 2011 AJGZ00000000 78489 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41567 USA WA Food 2011 AMNG00000000 80511 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41573 USA OH Human 2011 AJGY00000000 pCFSAN000316_01 80513 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41576 USA OH Human 2011 AMNH00000000 80515 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41578 USA OH Human 2011 CP004086 pSEEH1578_01,_02,_03 78477 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41579 USA OH Food 2011 AJGW00000000 78479 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41584 USA MI Human 2011 AMNI00000000 80517 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 82-2052 USA ME Human 1982 AMMX00000000 78473 IncI1, VirB/D4 P22, Gifsy-2 like, Fels-2, ST64B Heidelberg str. N1514 USA GA Food 2004 AKYM00000000 69871 VirB/D4 P22, Gifsy-2 like Heidelberg str. N1536 USA GA Food 2004 AKYN00000000 69869 IncI, VirB/D4 P22, Gifsy-2 like Heidelberg str. N15757 USA GA Food 2007 AMND00000000 80509 mob1, IncI1, VirB/D4, IncHI2 P22, Gifsy-2 like Heidelberg str. N18393 USA CT Food 2008 AMNF00000000 78485 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18413 USA GA Food 2008 AMNS00000000 80529 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18440 USA MD Food 2008 AMNT00000000 80531 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18453 USA MD Food 2008 AMNU00000000 80533 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N189 USA CT Food 2004 AMMP00000000 69873 mob2, IncI1, VB/D4 P22, Gifsy-2 like Heidelberg str. N19871 USA NY Food 2008 AMNC00000000 78481 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N19992 USA CA Food 2009 AMMC00000000 66709 ColE1 P22, Gifsy-2 like, P2-like Heidelberg str. N20134 USA GA Food 2009 AMNN00000000 78493 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N26457 USA CO Food 2010 AMNE00000000 78483 mob1, IncI1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. N29341 USA MN Food 2011 AMNJ00000000 80519 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N30678 USA NM Food 2011 AMNK00000000 80523 mob1, IncI1, VirB/D4, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N418 USA NM Food 2004 AMMB00000000 66707 mob1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. N4403 USA CA Food 2005 AMMT00000000 ‘ 69861 mob2, IncI1, VirB/D4 P22, Gifsy-2 like Heidelberg str. N4496 USA CO Food 2005 AMMQ00000000 69867 VirB/D4 P22, Gifsy-2 like Heidelberg str. N4541 USA CT Food 2005 AMMV00000000 69857 VirB/D4 P22, Gifsy-2 like Heidelberg str. N4630 USA GA Food 2005 AMNA00000000 80507 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N653 USA OR Food 2004 AMMW00000000 pCFSAN000412_01, _02 69855 mob2, IncI1, VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA31 USA MD Animal 1987 AMLV00000000 pCFSAN000441_01 66679 VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA32 USA TX Animal 1986 AMLU00000000 pCFSAN000442_01 66677 VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA35 Brazil n/a Unknown n/a AMLT00000000 pCFSAN000443_01 66675 ColE1, 8 kb resistance plasmid P22, Gifsy-2 like Heidelberg str. SARA37 USA CO Animal 1987 AMLR00000000 66673 IncI1, VirB/D4 P22, Gifsy-2 like, Fels-2, ST64B Heidelberg str. SARA 39 USA NC Human n/a AMMJ00000000 69887 IncHI2 P22, Gifsy-2 like, Fels-2 First released (Hoffmann et al. 2012). First released (Timme et al. 2013). Genomic Analysis and Virulence Differences GBE coverage of 23. De novo assemblies were performed using SNP matrix. We determined branch support by performing Roche’s Newbler software (v.2.6) with the resulting contigs 1,000 bootstrap replicates. being annotated using the NCBIs Prokaryotic Genomes We elucidated the relationships among isolates using the Automatic Annotation Pipeline (Klimke et al. 2009). Bayesian clustering method implemented in STRUCTURE Using DNA from the 2011 clinical outbreak isolate 41578, v2.3.2 (Pritchard et al. 2000; Falush et al. 2003). The analyses we also sought to close the genome using the Pacific were based on the kSNP matrix, and we ran ten replicate Biosciences (PacBio) sequencing platform. Specifically, we pre- analyses at K¼ 2–9 under the admixture model with corre- pared a single 10-kb library that was sequenced using the C lated allele frequencies being performed and visualized with chemistry on eight single-molecule real-time cells with a 90- DISTRUCT v1.1. (Rosenberg 2004). Each independent run min collection protocol on the PacBio RS. The 10-kb continu- consisted of 50,000 generations serving as burnin followed ous-long-read data were de novo assembled using the PacBio by 100,000 generations. The K statistic (Evanno et al. 2005) hierarchical genome assembly process/Quiver software pack- was used to identify the k value that best fits the data. Pairwise age, followed by Minimus 2, and they were polished with genetic distances, calculated as the number of nucleotide dif- Quiver. ferences, 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. Comparative and Phylogenetic Genome Analysis Alignments of core genes were accomplished using For comparative analysis, in addition to the 44 isolates se- MUSCLE with default settings; these alignments were then quenced in this study, we included sequence data available used to identify variable core genes. 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 Prediction of Prophages, Plasmids, Resistance, and (ABEL00000000) at the NCBI genome database. Virulence Genes Phylogenetic informative SNP sites (i.e., SNPs shared by two Prophages were identified using PHAST (Zhou et al. 2011). or more strains in the alignment) were identified by two dif- Sequences for intact phages were extracted from the original ferent methods. The first involved mapping the 454 reads to contig, in which they were found and mapped against the the complete reference genome of S. Heidelberg strain SL476 entire assembled data to determine whether the phage was using Roche Newbler software GS Reference Mapper (v.2.8.). present in an isolate at multiple sites on different contigs. SNP sites were then found using the SNP calling function of Contigs that could not be aligned to the reference genome GS Reference Mapper. SNP positions were defined, where were evaluated using BLAST to identify plasmids. The plasmids one or more isolates differed from strain SL476 with coverage were analyzed by comparative analysis using MAUVE algo- 10 and with 95% of the reads containing the SNP after rithm (Darling et al. 2010) and with the comparative analysis having filtered out the SNPs from homopolymer artifacts, fol- tools of RAST (Aziz et al. 2008). For plasmid closure, a con- lowed by a custom pipeline to construct a SNP matrix (Allard catenated sequence was generated from a 500-bp sequence et al. 2013). S. Heidelberg SL486 could not be included in this from each end of the contig. This artificial sequence from the analysis as its raw sequence data were not available. The single contig of the plasmid in study was used as a reference second SNP detection method was a reference free k-mer- to map all the 454 raw reads, using runMapping from based approach implemented in the program kSNP Newbler. If there were mapped reads that covered the con- (Gardner and Slezak 2010). This is a collection of Perl scripts junction point in the reference, the contig was a closed that aggregate the results from Jellyfish v1.1.3 (used for k-mer plasmid. counting) (Marcais and Kingsford 2011) and MUMmer v3.22 To determine the incompatibility (Inc) groups for plasmids, (used to align k-mers and detect variable positions) (Kurtz we used BLAST to find sequences described by Johnson and et al. 2004). Analyses with kSNP were based on a k-mer Nolan (2009) for specific Inc groups that would produce the- length of 25. oretical PCR amplicons for known Inc group sequences. To construct evolutionary relationships among the isolates, Resistance and virulence genes were identified by mapping the maximum-likelihood (ML) method was implemented in sequence data available at an in-house database, consisting the Genetic Algorithm for Rapid Likelihood Inference of 1,379 resistance genes, and 107 previously characterized (GARLI) software (Zwickl 2006b). GARLI analyses were per- virulence genes (84 pathogenicity and 23 fimbrial markers; formed using a web service (Bazinet and Cummings 2011) Huehn et al. 2009). Using a presence/absence matrix of that uses a special programming library and associated tools (Bazinet et al. 2007). All ML trees were constructed with genes conferring resistance to antibiotics and disinfectant GTR+I+G nucleotide substitution model. The 100 replicate agents, we constructed a similarity tree based on binary dis- runs of the nonbootstrapped data set were conducted to tances under neighbor-joining algorithm for tree construction; identify the most probable phylogeny based on our observed topological support was assessed based on 100 bootstrap Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1049 Hoffmann et al. GBE replicates. The analyses were performed using the R package whereas 20 isolates had composite XbaI/BlnI PFGE patterns ape (Paradis et al. 2004; Liu et al. 2011). 98% similar to isolates from the 2011 outbreak (supplemen- tary fig. S1, Supplementary Material online). Nine reference isolates showed the same PFGE pattern, as well as the same Caenorhabditis elegans Assay antibiotic resistance profile, noted for most of the outbreak Pathogenicity was evaluated using the Caenorhabditis elegans isolates (supplementary fig. S1, Supplementary Material survival assay. Caenorhabditis elegans strain SS104 glp-4 (bn2) online). was acquired from the Caenorhabditis Genetics Center. Worm cultures were maintained at 16 C, which is the per- Draft Genome Size missive temperature for this temperature-sensitive sterile mutant of C. elegans. Strains were cultured in C. elegans hab- Shotgun sequencing produced 44 draft genomes with an av- itation media (CeHM) in tissue culture flasks on a platform erage coverage of 23 and a minimum of 4,700 genes. The shaker (Sprando et al. 2009). Nematodes were bleached genomic size of the 44 S. Heidelberg isolates (including (0.5 M NaOH, 1% hypochlorite) to collect eggs, which were extrachromosomal DNA) ranged from 4.671 to 5.130 Mb incubatedinM9media for24 h to bringthem to synchronized (fig. 1). Differences in genomic size were due primarily to L1 stage and then transferred to CeHM. To produce synchro- the presence or absence of mobile genetic elements, including nized L4 stage worms, L1 worms were grown for an addi- phages and plasmids, especially the carriage of IncHI2 or IncA/ tional 72 h in CeHM. C plasmids (fig. 1). Next to thegenome sizevariation, figure 1 Pathogen lawns for survival assays with test strains, as well shows the estimated N50 sizes within S. Heidelberg draft as the food bacteria OP50, were prepared by inoculating genome sequences. The estimated N50 value is a rough esti- Nematode Growth Media (in 6-cm Petri plates) with 50 mlof mate of the quality and coverage of the draft genomes, and it an overnight bacterial culture. Plates were incubated over- represents the average contig size after assembly with the night at room temperature before worms were added. We Newbler software. used 60–80 synchronized L4 worms for each treatment group. Worms were scored every 24 h for survival (Aballay Phylogenetic Relationships among S. Heidelberg Using et al. 2000). SNP Analysis Animal survival was plotted using Kaplan–Meier survival Phylogenetically informative SNP sites were identified using curves and analyzed by log rank test using GraphPad Prism two independent methods. Figure 2 shows a ML tree based (GraphPad Software, Inc., La Jolla, CA). Survival assays were on SNP analysis of 44 S. Heidelberg isolates mapped to repeated at least two times. Survival curves resulting in strain SL476. Among the collective 40,716 variable SNPs P values< 0.05 relative to control were considered signifi- identified, 12,187 were determined to be “informative” cantly different. (i.e., SNPs shared by at least two isolates). The ML tree shown in figure 2 clearly demonstrates that S. Heidelberg Results formed a monophyletic group distinct from the outgroup comprised Salmonella enterica supsp. enterica serovar Antimicrobial Resistance and PFGE Two representative isolates from the 2011 ground turkey-as- sociated outbreak were previously reported to be resistant to ampicillin, gentamicin, streptomycin, and tetracycline (Folster et al. 2012). Antimicrobial susceptibility testing of seven addi- tional 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 FIG.1.—The number of assembled bases (Mb) and N50 contig size associated isolates exhibited indistinguishable, or nearly indis- (kb) for each sequenced S. Heidelberg isolate. Samples are colored accord- tinguishable (one band difference), XbalPFGEpatterns ing to the presence of antimicrobial resistance plasmids. No antimicrobial JF6X01.0032 or JF6X01.0058, and BlnI pattern JF6A26.0076 resistance plasmid, filled triangles; antimicrobial resistance Inc-I plasmid, and one food isolate had Xbal PFGE pattern JF6X01.0058 and filled diamonds; antimicrobial resistance Inc-H1/2 plasmid, filled circles; BlnI pattern JF6A26.0017. Of the 35 isolates not associated antimicrobial resistance plasmid Inc-I and H1/2, filled circles with borders; with the outbreak, 15 had composite Xbal/BlnI PFGE patterns antimicrobial resistance plasmid Inc-AC, filled squares; antimicrobial resis- tance plasmid IncI and Inc-AC, filled squares with borders. indistinguishable from ground turkey outbreak isolates 1050 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE 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. Newport (S. Newport) and S. Typhimurium. The SNP diver- trend for temporally isolated strains to cluster closely together, gence between S. Heidelberg and the two outgroups, one isolate, N29341, that was isolated in Minnesota in 2011 S. Newport and S. Typhimurium, is 24,724 and 24,305 from ground turkey by NARMS, clustered tightly together SNPs, respectively. The resultant tree shows two more inter- with the outbreak isolates received from CDC, whereas all esting points. First, isolates having a similar XbalPFGEpattern other reference isolates having the same PFGE pattern not (supplementary fig. S1, Supplementary Material online) clus- associated with the ground turkey outbreak are distinct tered together. Second, there is a tendency for isolates to from the isolates of the 2011 outbreak (fig. 2). These results cluster based on date of collection, such as noted for isolates suggest that N29341 might be epidemiologically related to from 2003, 2008, and 2011. In addition, the 2011 clinical and the 2011 outbreaks strains, as Minnesota was one of the food isolates that were originally thought to be associated states involved in that outbreak. with the ground turkey outbreak isolated in different states The SNP matrix, including the same S. Newport and clustered together with 80% bootstrap support (fig. 2), which S. Typhimurium outgroups, built using the k-mer strategy of could indicate that they are all members of a clonal commu- kSNP consisted of 41,432 variable SNP positions with 10,873 nity likely derived from the same source. Consistent with the being informative. Using the k-mer strategy, we identified 716 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1051 Hoffmann et al. GBE more SNPs than were noted within the in-house pipeline, present in genes shared by all 46 S. Heidelberg isolates and the where 40,716 variable SNPs were observed. Notably, the two outgroups. Only 718 SNPs, with 284 being informative, SNP diversity within the highly clonal S. Heidelberg isolates were found in genes shared only by S. Heidelberg isolates. The was lower, only 4,053 SNPs were found with 1,394 being SNP matrix of the core genes was used to construct a ML informative. The resulting ML tree (fig. 3) was well resolved phylogenetic tree (supplementary fig. S2, Supplementary and congruent with the tree topology in figure 2.Toshow Material online), which also confirms the tree topology in that the tree topology is not an artifact of mobile elements figures 2 and 3. that are not present in all isolates, we constructed a kSNP In this study, we determined that all 46 S. Heidelberg iso- matrix using SNPs present only on core genes. This produced lates, including S. Heidelberg SL476 and S. Heidelberg SL486, a matrix of 24,447 SNPs, of which 7,630 were informative and partitioned into three paraphyletic groups, each having 100% FIG.3.—ML tree based for the 44 S. Heidelberg isolates and two previously reported S. Heidelberg genome sequences SL476 and SL486 (Fricke et al. 2011). A total of 4,053 SNPs with 1,394 being informative were found based on k-mer analysis using kSNP. ML trees were generated as described in figure 2. Bootstrap values (1,000 iterations) are reported below each node. The numbers of unambiguous substitutions that mapped to the tree only once and are greater than zero are given above each node in blue. The numbers in parenthesis represent the nodes in table 3. To the right of the tree, two Distruct plots were reconstructed with the same SNP matrix—one including all 46 S. Heidelberg isolates and, adjacent to that, another with only those isolates from group 3—to present a fine-scale structure is shown. The Distruct plot was generated using a model-based Bayesian clustering method implemented in Structure v2.3.2 and visualized with DISTRUCT v1.1. 10 replicate analysis at K¼ 2–9 under the admixture model with correlated allele frequencies were performed. Each independent run consisted of 50,000 generations serving as burnin followed by 100,000 generations. Different colors represent the different clusters and each bar represents an individual isolate. The fraction of the bar that is a given color represents the coefficient of membership to that cluster (e.g., multicolored bars indicate membership to multiple groups indicative of admixture). 1052 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE bootstrap support (fig. 3). Group 1 contains three historical 41578 and 41579 are quite distinct from isolates N19871, N18413, N18440, and N18453, which have 44 SNPs that isolates (two clinical and one turkey isolate) with a PFGE pat- allow complete discrimination between the two groups tern, for which there was no assigned Pulse-Net designation (fig. 3). Furthermore, isolate 41565, which has the same (supplementary fig. S1, Supplementary Material online). PFGE and antibiotic resistance profile as the 2003 isolates Group 2 is comprised of 11 isolates, including clinical, food, 32507, 24393, and 20752, is more closely related to the and animal isolates, with XbaI PFGE pattern JF6X01.0022, as other isolates associated with the 2011 outbreak, even with well as the isolates SL476 and N19992 having different XbaI having a different BlnI pattern and antibiotic resistance profile patterns JF6X01.0133 and JF6X01.0326, respectively. Group (fig. 3). 3 is composed of 30 isolates from food, animal, and clinical To provide an additional picture as to the number of distinct sources, including the outbreak isolates, all having a very sim- groups the sampled serovars segregate into, we constructed ilar XbaI pattern (JF6X01.0032, JF6X01.0034, JF6X01.0058). groups using the program STRUCTURE, the results of which The distribution of PFGE patterns on the tree topology were visualized using DISTRUCT (fig. 3). The results are con- based on SNPs differs between the two enzymes XbaIand sistent with the groupings formed through the phylogenetic BlnI. For example, N19992 (group 2) and N18393 (group 3) analyses in that there are three distinct groups (each repre- are quite distinct on the ML tree, even though both have the sented by a different color [fig. 3]). same BlnI pattern JF6A26.0015 and are adjacent to each other on the PFGE dendrogram (supplementary fig. S1, Differences between Outbreak Isolates and Nonoutbreak Supplementary Material online). Pairwise SNP variation be- Isolates in Group 3 tween the three distinct S. Heidelberg lineages seen in figure 3 are presented in table 2 and support the phylogenetic To further investigate the relationships among outbreak iso- partitions revealed above. Mean divergence among the three lates, we constructed the Distruct plot for group 3. Within that clades ranged from 200 to 712 nt differences, whereas the group, there is support for four subgroups; these four groups mean intragroup nucleotide differences ranged from 47 to also are found in the phylogenetic analysis, each with a75% 155 (table 2). bootstrap support. Pairwise SNP variation between the four Notably, the mean SNP variation within the outbreak iso- sublineages also support the phylogenetic partitions described lates is only 17 SNPs. The phylogenetic analysis clustered the above (table 2 and fig. 3). Mean divergence among the four outbreak isolates together even though the PFGE profile (sup- subclades ranged from 36 to 153 nt differences, whereas the plementary fig. S1, Supplementary Material online) is slightly mean intragroup nucleotide differences ranged from 15 to 30 different (XbaI pattern JF6X01.0032, JF6X01.0058; BlnIpat- (table 2). The 2008 strains seem to have quite a unique SNP tern JF6A26.0017 and JF6A26.0076). For example, isolates profile when compared with the remaining isolates of group 41578 and 41579, sharing the same PFGE and antibiotic re- 3. The ML tree, the Distruct plot, and the pairwise SNP varia- sistance profile with isolates N19871, N18413, N18440, and tion show the evolutionary changes occurring among isolates N18453 (isolated in 2008), are more closely related to the from 2003 to 2009 and 2011, revealing the apparent emer- other outbreak isolates, even though they have a different gence of the genetically unique S. Heidelberg strain that was XbaI pattern. We also found that the outbreak isolates responsible for the 2011 outbreak. Table 2 (Mean Pairwise Distance (Number of Nucleotide Differences) between the (A)Three Major S. Heidelberg Groups in figure 3 and the (B) Four S. Heidelberg Subgroups from Group 3 in figure 3 S. Heidelberg Gr. 1 S. Heidelberg Gr. 2 S. Heidelberg Gr. 3 (A) S. Heidelberg Gr. 1 155 S. Heidelberg Gr. 2 535 101 S. Heidelberg Gr. 3 712 200 47 Gr. 3/I Gr. 3/II Gr. 3/III Gr. 3/IV (B) S. Heidelberg Gr. 3/I n/c S. Heidelberg Gr. 3/II 144 30 S. Heidelberg Gr. 3/III 139 42 15 S. Heidelberg Gr. 3/IV 153 46 51 27 Outbreak isolates 151 46 65 36 Intragroup mean pairwise distance (number of nucleotide differences). Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1053 Hoffmann et al. GBE Table 3 lists variable genes having unique nucleotide sub- protein (15 SNPs), terminase small subunit (14 SNPs), outer stitutions which defined specific sublineages that could be membrane lipoprotein Blc (7 SNPs), and alpha amylase used to separate them from other highly clonal (6 SNPs), were found to have a disproportionate number of S. Heidelberg isolates in group 3. Table 3 also lists their SNP SNPs that could be useful for rapid subtyping the serovar change and, if applicable, the gene carrying the SNPs and the Heidelberg and/or as targets of resequencing. Interestingly, resulting amino acid change. Twenty-three unique SNPs, with although these genes have several SNPs, the haplotype diver- 14 being nonsynonymous (i.e., causing an amino acid change) sity is fairly low, with only a few haplotypes (2–4) among 46 SNPs, were found to be unique to group 3 and clustering isolates being identified. This fact suggests that, even though together isolates with XbaI pattern (JF6X01.0032, thereare geneswithup to35SNPs, thereare only twotofour JF6X01.0034, JF6X01.0058). These nucleotide substitutions unique sequences among S. Heidelberg that are congruent were found on 20 different genes that could be useful for with our tree topology for three different groups of subtyping of strains belonging to group 3 (table 3 and fig. 3). S. Heidelberg. Furthermore, 24 unique SNPs (11 being nonsynonymous) were found at the node that characterizes the 2008 ground Complete Genome from Clinical Isolate 41578 turkey isolates, which seemed to have their own unique SNP profile distinct from that of the remaining group 3 isolates. In the course of performing this study, we demonstrated that Twelve of the 24 SNPs are located on the insertion sequence it was possible to rapidly determine the complete sequence for (IS) 2 element with five SNPs being within the IS2 TnpA trans- the chromosome and three plasmids from a clinical isolate posase and seven SNPs on the IS2 TnpB transposase genes associated with the 2011 outbreak using the Pacific (table 3 and fig. 3). Biosciences (PacBio) system. A single contig of 4,793,479 bp Importantly, two nonsynonymous SNPs were found to be (GC content 52.2%) representing the complete chromosome only present in all isolates associated with the 2011 outbreak. and three contigs of 4,773 bp, 35,297 bp, and 117,929 bp These two SNPs were located within two different variable representing three plasmids were generated. The chromo- core genes: a putative hydrolase with amino acid change some contains 12 different Salmonella pathogenicity islands P653Q, and the gluconate transporter GntP gene with and four prophages (prophage P22, lambda-like prophage amino acid change L187V (table 3 and fig. 3). Five SNPs, Gifsy-2, prophage P4, and prophage P2-like; fig. 5). The larg- with four being nonsynonymous, were unique to isolates est plasmid pSEEH1578_01 is a resistance plasmid with incom- 41563, 41566, 41567, 41573, 42578, 41579, and N29341 patibility group IncI1 encoding the resistance to gentamicin (associated to the outbreak) and not observed in 41565 (aacC), streptomycin (aadA1), tetracycline (tetA, tetR(A)), (table 3 and fig. 3). From these five SNPs, three were in plas- and ampicillin (bla )(fig. 5). The plasmid pSEEH1578_02 tem1 mid sequences, whereas two were found on variable core is a VirB/D4 plasmid that carries the type IV secretion system genes: the flagellin CDS resulting in an amino acid change (T4SS; fig. 5), and plasmid pSEEH1578_03 is a mobilization from threonine to serine at position 282 (T282S), and plasmid that carries genes capable of promoting plasmid mo- within a putative hypothetical protein (possibly functioning bilization such as mbeA, mbeB, mbeC,and mbeD. as H+ gluconate symporter related to permeases) producing an amino acid change from glycine to serine at position 358 (G358S) (table 3). Identification of Prophages PHAST analysis identified six different phages among the Genetic Variation among S. Heidelberg S. Heidelberg isolates. Only those phages that PHAST deter- mined as being intact were further analyzed. All 44 isolates At a minimum, 273 genes involved in a variety of different contained prophage P22 and lambda-like prophage Gifsy-2. functions, such as DNA replication and repair, cell division, Group 1, which was composed of three historical isolates transcription, metabolism, and virulence, vary among the (SARA 39, SARA 37, and 82-2052), also contained prophage S. Heidelberg isolates observed in clinical, animals, and/or Fels-2, not observed in any of the group 2 or group 3 isolates. food samples, including the isolates associated with the Additionally, SARA 37 and 82-2053 also contained prophage 2011 outbreak. The variable core genes, including the ST64B (table 1 and fig. 2). number of SNPs and haplotypes are listed in supplementary Heidelberg isolates belonging to group 3, which includes table S1, Supplementary Material online. Figure 4 shows a the 2011 outbreak isolates, contained prophage P4 and P2- histogram comparing strains based on: 1) the number of like, neither of which was seen in any of the isolates from SNPs according to the number of core genes and 2) the hap- lotype diversity according to the number of core genes. Only a group 1 and 2, excepting isolate N19992, a group 2 isolate few mutation “hotspot” genes, such as tailspike (35 SNPs), that was confirmed to contain the P2-like prophage (table 1 scaffolding protein (22 SNPs), DNA polymerase V subunit and fig. 2). The phage presence and absence matrix shows umuD (20 SNPs), phage lysozyme (18 SNPs), gifsy-2 prophage correlation with the tree phylogeny shown in figure 2. 1054 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Table 3 Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study Node on Tree SNP Location NT Change AA Change NT Position Locus Tag 1 D-Mannonate oxidoreductase A!GD!G 839 SEEHRA37_03221 1 Translocation protein TolB T!A P 1152 SEEHRA37_24108 1 Glutathione S-transferase C!AW!Stopp 233 SEEHRA37_20055 1 Nonannotated C!TR!C n/a SEEHRA37_contig25 1 Hypothetical protein/putative transcriptional regulator G!TV!F 376 SEEHRA37_19825 1 Ureidoglycolate dehydrogenase AllD C!TT!I 20 SEEHRA37_20355 1 Starvation-inducible outer membrane lipoprotein C!T S 63 SEEHRA37_07830 1 ATPase recombination and repair protein G!T L 966 SEEHRA37_05779 1 DNA-specific endonuclease I A!CK!T 341 SEEHRA37_19810 1 Scaffolding protein C!A V 135 SEEHRA37_21903 1 Peptidase T G!TV!L 976 SEEHRA37_13506 1 Predicted metal-dependent membrane protease C!A A 108 SEEHRA37_15707 1 Hypothetical protein A!GN!D 418 SEEHRA37_04201 1 Glycine/serine hydroxymethyltransferase GlyA T!GS!A 145 SEEHRA37_00365 1 Preprotein translocase subunit SecF G!A A 534 SEEHRA37_22413 1 Hypothetical protein A!GT!A 1879 SEEHRA37_01714 1 Ribosomal protein L21 RplU A!CI!L 220 SEEHRA37_20845 1 Nonannotated C!GQ!E n/a SEEHRA37_contig19 1 Isoaspartyl peptidase C!TA!V 380 SEEHRA37_16329 1 Lysine-N-methylase C!T G 57 SEEHRA37_07110 1 Hypothetical protein T!C L 162 SEEHRA37_23678 1 Hypothetical protein A!CK!T 227 SEEHRA37_07475 1 Nonannotated T!CI!V n/a SEEHRA37_contig14 2 Hypothetical protein A!GQ!R 851 SEEHRA37_20355 2SensorkinaseDpiB A!TI!F 67 SEEHRA37_10675 2 Notannotated C!T P 63 SEEHRA37_contig6 2 Hypothetical protein G!A F n/a SEEHRA37_contig6 2 Nonannotated G!AG!S n/a SEEHRA37_contig5 2 Proline dehydrogenase PutA C!A G 975 SEEHRA37_04036 2 Electron transport complex protein RnfC T!C I 276 SEEHRA37_02971 2 Protein involved in chromosome partitioning MukB G!CG!A 1514 SEEHRA37_14706 2 N-ethylmaleimide reductase A!T G 351 SEEHRA37_02866 2 Hypothetical protein C!T L 1039 SEEHRA37_16509 2 Multidrug efflux system subunit MdtC A!TQ!H 2067 SEEHRA37_00060 3 Putrescine/spermidine ABC transporter ATPase G!TP!Q 47 SEEHRA37_13511 protein PotA 3 Pilus assembly protein, porin PapC G!TR!S 4 SEEHRA37_18324 3 IS2 repressor TnpA G!A T 93 SEEH8393_18082 3 IS2 transposase TnpB A!G T 195 SEEH8393_18077 3 IS2 transposase TnpB A!GN!D 502 SEEH8393_18077 3 IS2 transposase TnpB G!A A 483 SEEH8393_18077 3 Not annotated (plasmid sequence) T!CA!V n/a SEEH8393_contig56 3 IS2 repressor TnpA C!T Y 309 SEEH8393_18082 3 IS2 repressor TnpA A!G E 195 SEEH8393_18082 3 IS2 transposase TnpB C!A V 525 SEEH8393_18077 3 IS2 transposase TnpB C!A P 192 SEEH8393_18077 3 IS2 transposase TnpB A!G A 789 SEEH8393_18077 3YacA T!G V 123 SEEH8393_01004 3 Hypothetical protein G!AA!T 232 SEEHRA37_09669 3 Multifunctional fatty acid oxidation complex G!TM!I 2094 SEEHRA37_15762 subunit alpha FadJ (continued) Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1055 Hoffmann et al. GBE Table 3 Continued Node on Tree SNP Location NT Change AA Change NT Position Locus Tag 3 tRNA uridine 5-carboxymethylaminomethyl G!AG!S 1199 SEEHRA37_11040 modification enzyme GidA 3 Permease DsdX C!AA!D 1166 SEEHRA37_02334 3 IS2 transposase TnpB T!C R 549 SEEH8393_18077 3 IS2 repressor TnpA C!T A 225 SEEH8393_18082 3 IS2 transposase TnpB C!AP!T 319 SEEH8393_18077 3 IS2 repressor TnpA C!T L 102 SEEH8393_18082 3 Arginyl-tRNA synthetase ArgS A!T G 237 SEEHRA37_07325 3 Nonannotated C!AA!N n/a SEEHRA37_contig18 3 Hypothetical protein G!AD!N 103 SEEHRA_37_05617 4 Not annotated G!AV!L n/a SEERA37_contig2 4 Sodium/glucose cotransporter G!AS!N 1037 SEERA37_14011 5 Phosphoserine phosphatase C!AP!Q 653 SEEHRA37_11265 5 Gluconate transporter GntP C!GL!V 187 SEEHRA37_18800 6 Flagellin and related hook-associated proteins G!CT!S 282 SEEHRA37_18390 6 Electron transfer flavoprotein, beta subunit G!AG!S 358 SEEHRA37_02456 6 Hypothetical protein (plasmid sequence) C!AL!M 562 SEEH8393_17461 6 Not annotated (plasmid sequence) G!AG!R n/a SEEH8393_contig56 6 Not annotated (plasmid sequence) C!TG!S n/a SEEHRA39_contig149 7 Type III secretion protein SopE T!AH!Q 560 SEEHRA37_14371 FIG.4.—(A) Histogram showing the number of SNPs per core genes. (B) Histogram showing haplotype diversity for all variable, core genes. Identification of Virulence Genes among S. Heidelberg S. Heidelberg isolates (supplementary table S2, Supplementary From the 107 known and recognized Salmonella-associated Material online). One outer membrane fimbrial usher gene, virulence genes, we identified 62 pathogenicity and 13 fim- safC, was only found in the three historical isolates (SARA 39, brial markers that were highly conserved among the 44 SARA 37, and 82-2052), all belonging to group 1. All other 1056 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE A B SPI-1 SPI-9 SPI-13 SPI-12 P-4 SPI-3 S. Heidelberg 41578 P-2 CP004086 tetA pSEEH1578_01 tetR(A) SPI-4 SPI-2 aadA1 SPI-11 SPI-5 blaTEM1 SPI-14 P-22 SPI-6 Gifsy-2 aacC SPI-16 virB3-4 virB5 virB2 virB6 virB1 virB8 virB9 virD2 virB10 hicB hicA virB11 virD4 pSEEH1578_02 stbD repA 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 areshownindarkblueand 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. detected virulence genes, observed to be located on the chromosomal msgA gene, which is essential to Salmonella Salmonella pathogenicity islands (SPI-1, SPI-2, SPI-3, SPI-4, mouse virulence (Gunn et al. 1995), and is common in S. and SPI-5), the lambda-like prophage Gifsy-2, and/or islets, Typhimurium (Huehn et al. 2010). Other genes present were present in all S. Heidelberg isolates, suggesting that all among all isolates encoded proteins involved in type III secre- isolates are pathogenic (supplementary table S2, Supplemen- tion system (T3SS) and adhesins. Remarkably, the sopE1 viru- tary Material online). For example, all isolates carried a lence gene, involved in the translocation of effector proteins Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1057 Hoffmann et al. GBE into host cells, was found in all of our strains. Moreover, al- isolates belonging to group 3 in figure 4, while subclade 2 though sopE1 is known to be carried by P2-like phages, in this contains one clinical isolate and isolate N653, the unique iso- study we also identified this gene on lambda-like phage Gifsy- late found to carry two VB/D4 plasmids (the other plasmid partitioned into clade III). Moreover, using the primers de- 2, which is present in all S. Heidelberg isolates. Consequently, signed by Johnson et al. (2012) for the relaxase gene (taxC) all S. Heidelberg isolates have the sopE1 virulence gene. Group (VirD2 component) to screen for the four different IncX sub- 3 isolates, which also contain the P2-like phage, have the groups, we determined that the VB/D4 plasmids from clade I, sopE1 gene twice in their chromosome; one gene is carried subclade 1 from clade II, and clade III belonged to the IncX1 by lambda-like phage Gifsy-2 and the other one is located on plasmid subgroup, whereas the two plasmids from clade II the P2-like phage. Analysis confirmed that these two genes subclade 2 were deemed to be members of the IncX4 shared a high degree of similarity having no more than eight subgroup. mutations between the two different sopE1 alleles. In con- Interestingly, isolate N653 carried two plasmids, one being trast, virulence determinants known to be associated with in the IncX1 and another in the IncX4 incompatibility groups. It SPI-7, prophage Gifsy-1, Gifsy-3, Fels-1, and Salmonella plas- appears that there are significant differences among these mid virulence (spv) were not identified in any of the 44 isolates groups that allow simultaneous maintenance of two plasmids examined in this study, confirming sequencing results that belonging to the same incompatibility group, a heretofore none of the strains carried any of these mobile genetic believed unachievable occurrence. Given the notion that elements. two plasmids in the same Inc group cannot simultaneously coexist within the same bacterium, another plausible explana- Description of VirB/D4 Plasmids tion could be that IncX4 might actually belong to a different Several plasmids were found among the 44 S. Heidelberg incompatibility group other than that of IncX1, which was isolates sequenced in this study (table 1). Electrophoresis of only previously experimentally confirmed to belong to IncX samples obtained using a basic alkaline plasmid purification (Norman et al. 2008). Consistent with this hypothesis was procedure demonstrated the presence of endogenous plas- the fact that a BLAST search with the IncX4 repA gene (pre- mids that had not integrated into the host genome (data not sent only in both putative IncX4 plasmids) sequence identified shown). The plasmid sizes varied from 3.7 to over 200 kb. significant sequence similarity to only the Citrobacter roden- Twenty-one S. Heidelberg isolates, including all five human tium ICC168 plasmid pCROD2 (previously grouped to IncX4 isolates and two ground turkey isolates associated with the [Johnson et al. 2012]) and two plasmids from S. Heidelberg 2011 outbreak, contained a plasmid ranging in size from 33 to (pSH163_34, pSH696_34), suggesting that the IncX4 group is 40 kb, which carried genes (virB1,virB2,virB3-4,virB5,virB6, a new incompatibility group distinct from that of IncX1. If this virB7,virB8,virB9,virB10,virB11,virD2, and virD4) associated is confirmed to the case, coexistence of these two “IncX” withthe VirB/D4T4SS(table 1). Three specific roles for the plasmids would be explained, although additional research bacterial T4SS have been identified: 1) to facilitate transloca- needs to be conducted to confirm this possibility. tion of DNA via a conjugative mechanism to recipient cells; 2) In this study, eight of the plasmids carrying the VirB/D4 to facilitate translocation of effector molecules, such as pro- T4SS were completely sequenced including the two plasmids teins, to eukaryotic target cells; and 3) to function in DNA from isolate, N653 (table 1). The annotated sequence dem- uptake or dissemination from or into the environment milieu onstrated that the eight plasmids carry, in addition to the VirB/ (Christie et al. 2005; Alvarez-Martinez and Christie 2009). D4 T4SS systems, genes encoding other virulence-associated To further study the diversity among all 22 identified VB/D4 determinants and several proteins involved in a variety of met- plasmids, we identified SNPs that might be useful in differen- abolic processes, for example, hemolysin expression-modulat- tiating these plasmids. The SNP matrix from the 22 plasmids ing gene ymoA, DNA-binding protein genes (hns), (stpA), the and the reference plasmid built using the k-mer strategy of toxin/antitoxin stability genes (stbE), (stbD), IncN plasmid Kika kSNP consisted of 338 informative SNP positions. Figure 6 gene (kika), and DNA topoisomerase III (topB). shows a ML tree based on SNP analysis. The phylogeny A progressive Mauve alignment using the data from the placed the VB/D4 plasmids into three well-supported clades eight completely sequenced plasmids and the reference (95%, 100%, and 100% bootstrap support). Clade I con- S. Heidelberg plasmid pSARA30 (supplementary fig. S3, tained two plasmids isolated from the historical isolates, Supplementary Material online) separated the plasmids into SARA 37 and 82-2052, belonging to group 1 on the ML three distinct groups. As expected, the plasmids partitioned in tree shown in figure 3. Clade III is composed of ten plasmids, different clades in the plasmid VirB/D4 ML tree (fig. 6)were including the reference plasmid pSARA30 and 9 isolates be- also found by Mauve analysis to belong to different groups. longing to group 2 in figure 4. Clade II contains 11 plasmids, Group 1 included plasmids isolated from S. Heidelberg SARA including 7 plasmids from 7 isolates associated with the 2011 31, SARA 32, N4403, and N653, all having a size of approx- outbreak. Clade II can be further separated into two sub- imately ~38.0 kb and showing the highest degree of sequence clades, one of which contains nine plasmids derived from similarity to the S. Heidelberg plasmid pSARA30. Unique to 1058 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE 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. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1059 Hoffmann et al. GBE this group are several hypothetical proteins: truncated trans- mobilization plasmid, the VirB/D4 T4SS plasmid, the IncI1, and posase, DEAD/DEAH box helicase domain protein, and puta- the IncHI1. Except SARA 31, in which phage P4 and P2-like tive N-acetyltransferase. The second group is composed of the were not present, all other isolates had the same phage plasmids isolated from S. Heidelberg 41573 (~38 kb) and profile. 41563 (~33 kb)—both involved in the 2011 outbreak. The The results, which showed that S. Heidelberg isolates are best BLAST hit for them is Escherichia coli plasmid pDKX1- significantly (P< 0.0001) more pathogenic than the E. coli OP TEM-52 (GenBank accession number: JQ269336). Attributes 50 control strain, failed to demonstrate any appreciable dif- unique to this group are several hypothetical proteins, a ference in pathogenicity between strains carrying the VirB/D4 replicon initiation protein, and an ATPase central domain- T4SS plasmid and those strains that did not (fig. 7). However, containing protein. The third group contains plasmids isolated we also identified some VirB/D4 T4SS components on the from S. Heidelberg 41578, a 2011 outbreak clinical isolate, resistance plasmids IncA/C, IncHI2, and IncI1. Interestingly, and the chicken breast isolate N653. Both plasmids are ~35 kb isolate 29169, which did not carry any large plasmid and in size and members of the IncX4 subgroup. They share a high any of the T4SS components, was found to be significantly degree of sequence similarity to C. rodentium ICC168 plasmid (P< 0.0001) less pathogenic using this assay than any of the pCROD2 (FN543504). Unique to this group are some meta- other isolates evaluated. Isolate 41565, which carried the bolic genes, such as putative toxin–antitoxin system genes small mobility plasmid as well as additional IncHI2 plasmid, (hicA), (hicB), DnaJ-class molecular chaperone (dnaJ) (heat demonstrated the second weakest pathogenic potential. shock protein), plasmid replication gene (repA), site-specific Isolate SARA 31, which carried the VB/D4 T4SS plasmid but recombinase (xerD), and acetyl-CoA carboxylase beta subunit no other resistance plasmid, showed a survival curve profile (accD). With regards to the VirB/D4 gene clusters found on the similar to that of isolate N418, which did not have the VB/D4 IncX or IncX-like plasmids in this study, a considerable variation T4SS plasmid but carried the IncA/C plasmid. was observed, sometimes having <55% sequence identity The highest degree of pathogenicity was seen in N30678, between homologous genes (supplementary figs. S3 and which carries four plasmids including resistance and VB/D4 S4, Supplementary Material online), a phenomenon that has plasmids, and it is significantly (P< 0.05) different from all been observed elsewhere among the Enterobacteriaceae other isolates except the clinical outbreak isolate 41578, (Johnson et al. 2012). Comparison of the plasmids from the which also consist resistance and VB/D4 plasmids (fig. 7). three groups demonstrated that most of the T4SS-associated The data show that isolates not carrying any T4SS components genes divergence was seen in plasmids derived from group 3 tend to be significant less (P< 0.0001) pathogenic than those isolates. isolates that do carry the VB/D4 T4SS plasmid. Moreover, Noteworthy was the observation that three of these eight those strains that carry both the VB/D4 T4SS plasmid and an plasmids, each isolated from different 2011 clinical outbreak antibiotic resistance plasmid, which also carries the T4SS strains (41563 OR Blood, 41578 OH stool, 41573 OH Urine), genes, demonstrated a pathogenic potential greater than iso- differed in DNA sequence and gene carriage, suggestive of lates carrying only one of these plasmids. extensive horizontal gene transfer (HGT) or a mixed subpop- ulation of S. Heidelberg isolates. Caenorhabditis elegans System 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 po- tential of our strains. This system has been used successfully as an invertebrate host model to assess the virulence determi- nants 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 FIG.7.—Caenorhabditis elegans survival data from six S. Heidelberg small mobilization plasmid; strain SARA 31 carries only the isolates. The figure shows that the six S. Heidelberg isolates (29169, VirB/D4 T4SS plasmid; strain 418 carries a small mobilization SARA 31, 418, 41565, 41578, and N30678) are significantly plasmid and an IncA/C plasmid; strain 41565 carries a small (P< 0.0001) more pathogenic than the Escherichia coli OP 50 control mobilization plasmid and an IncHI2 plasmid; strain 41578 has strain. Further the figures show that isolates not carrying any T4SS com- a small mobilization plasmid, the VirB/D4 T4SS plasmid, and ponents tend to be significant less pathogenic than those isolates that do carry them. an IncI1 plasmid; and N30678 has four plasmids: the small 1060 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Identification of Antimicrobial Resistance Genes and Figure 8 shows a neighbor-joining tree generated from a Plasmids presence/absence matrix of resistance genes for the 30 group 3 Heidelberg isolates having a common PFGE profile. As ex- Analysis of WGS data showed 27 plasmid-associated antibi- pected, resistance tended to segregate according to the par- otic resistance genes among the 44 isolates and included ticular Inc group plasmid they carried and not by clonality or genes expected to produce resistance to aminoglycosides, evolutionary relationship. This observation indicated that the beta-lactams, phenicols, folate pathway inhibitors, tetracy- presence of specific resistance genes correlates with the par- clines, and cephems. These genes were not detected in the ticular Inc group of the plasmids. For example, the resistance eight pan susceptible isolates that also did not carry any of the genes (aphA1, bla , sul2, cmlA1, floR) were only detected CMY2 antibiotic resistance plasmids. Resistance phenotypes corre- in isolates that carry an IncA/C plasmid, whereas only isolates lated with genotypes in all strains examined (table 4). carrying an IncHI2 plasmid were found to contain the resis- However, some strains with the same genotype exhibited dif- tance genes aph, sph, tetC, tetD,and ble. Similarly, the aacC ferent antimicrobial resistance patterns. For example, both gene only appears to be present in those isolates identified as strains N18440 and N18413 contain the aadA1genebut carrying an IncI1 plasmid. The association of different resis- only isolate N18440 showed resistance to streptomycin. In tance genes with specific Inc groups is well known. this study, streptomycin resistance was defined using the MIC value >32 mg/l determined by CLSI. However, studies Mobilization Plasmids have shown that strains harboring an aadA gene cassette Plasmids below 8 kb did not carry any antimicrobial resistance can have MICs of 16 mg/l, which would classify these strains or virulence-associated genes. Instead, these smaller plasmids as streptomycin-susceptible (Sunde and Norstro ¨ m2005). tended to carry genes capable of promoting plasmid mobili- Among the 36 S. Heidelberg isolates resistant to one or zation, such as the mbeA, mbeB, mbeC, mbeD genes. more antibiotics, 24 contained sequences indicative of pres- S. Heidelberg isolates N189, N653, and N4403, which were ence of an IncI1 plasmid, 5 had an IncHI2 plasmid, and 1 originally isolated from chicken breast, carry a small, 3,772 bp carried an IncA/C plasmid. Three isolates contained two plas- plasmid that is identical to the S. Heidelberg SL476 plasmid mids, one with IncI1+ IncHI2, and two isolates each having pSL476_3 (GenBank accession number: CP001119) (Fricke IncI1+ IncA/C replicons (table 4). Besides isolate 41565, which et al. 2011). Allgroup3isolates (fig. 3), except isolate carries an IncHI2 plasmid, all Heidelberg isolates associated 24393, have similar XbaI patterns (JF6X01.0032, with 2011 outbreak carry an IncI1 plasmid, which was de- JF6X01.0033, JF6X01.0034, JF6X01.0058), and carry a tected the most overall among our strains. The IncI1 plasmid 4,473 bp plasmid that shares high sequence similarity with carries several antibiotic resistance genes and components of the S. Heidelberg plasmid pSH1148_4.8 (GenBank accession the T4SS. The IncHI2 plasmid was found in eight isolates and number: JX494965). SARA 35, which was isolated in Brazil, carries several antibiotic resistance genes and genes conferring was found to carry a 6,647 bp ColE1 plasmid resistance to quaternary ammonium compounds, as well as (pCFSAN000443_01) that is similar to E. coli plasmid ColE1 genes encoding the conjugative T4SS. Furthermore, unique to (GenBank accession number: J01566). the IncHI2 plasmid is the presence of genes encoding resis- tance to heavy metals, including tellurium, copper, cadmium, Discussion mercury, and silver. The three isolates resistant to the most (ten or more) anti- Within the last several years, S. Heidelberg has been identified biotics carried the IncA/C multidrug resistance plasmid. In ad- as one of the top serovars responsible for human illness, in- dition to carrying the genes encoding resistance to different cluding a recent multistate outbreak of an antibiotic-resistant antibiotics and quaternary ammonium compounds, this plas- strain associated with consumption of contaminated ground mid also carries the conjugative T4SS transfer system. All three turkey. Using whole-genome sequences, we characterized the IncA/C plasmids carried the resistance element bla -blc- CMY-2 genetic diversity of 46 S. Heidelberg isolated over ~30 years. sugE2, which was only found in one other isolate, N189, on This study illustrates the novel ways in which WGS, combined an IncI1 plasmid. with phylogenetic analysis, can be used to investigate and SARA 35, which was only resistant to ampicillin, did not characterize the diversity of S. Heidelberg. Using this strategy carry a multidrug resistance plasmid but, instead, carried a with outbreak isolates, we were also able to comprehensively smaller plasmid (~8.0 kb) having only a single determinant characterize and differentiate among highly clonal that correlated with resistance to b-lactams (tnpA-bla ). S. Heidelberg isolates. For example, SNP analysis in combina- TEM-1 This plasmid shared a high degree of sequence similarity to tion with phylogenetic analysis is able to determine the phy- S. Typhimurium pAnkS (accession number: DQ916413), logenetic relationships among S. Heidelberg isolates, easily which also carried the transposon resistance element (tnpA- separating them from other S. Heidelberg strains, including bla ). The Inc group for this plasmid could not be deter- isolates having the same XbaIand BlnI PFGE patterns and TEM-1 mined from the sequence. thought to be highly clonal. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1061 Hoffmann et al. GBE 1062 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Table 4 Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study Isolates Antimicrobial Resistance Profile MDR Plasmids Resistance Genes on Plasmids SARA37 STR, TET IncI1 aph, sph, strA, strB, tetC, ble SARA35 AMP ca 8 kb bla TEM1 SARA32 Susceptible None None SARA31 Susceptible None None N418 AMC,AMP,FOX,AXO,CHL,GEN,KAN,STR,SUL,TET,TIO IncA/C aac3-VI, aadA1, aadA2, aadB, aphA1, strA, strB, bla ,bla ,sul1, sul2,cmlA1, CMY2 TEM1 floR, tetA, tetR(A), qacDE, sugE2 N19992 Susceptible None None SARA39 SUL IncHI2 aadA1, sul1, tetG, tetR(G), qacDE N189 AMC,AMP,FOX,AXO,TIO IncI1 bla , sugE2 CMY2 N1514 Susceptible None None N1536 AMP, GEN,STR,SUL IncI1 aac3-VI, aadA1, bla , sul1, qacDE TEM1 N4496 Susceptible None None N4403 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE N4541 Susceptible None None N653 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 82-2052 KAN, STR, TET IncI1 aph, sph, strA, strB, tetC, ble 24390 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 21381 AMC,AMP,FOX,AXO,GEN,KAN,STR,SUL,TET,TIO IncI1; IncA/C aac3-VI, aadA1, aphA1, strA, strB, bla ,sul1, sul2,tetA, tetR(A), qacDE, sugE2 CMY2 N4630 GEN,STR,SUL IncI1 aac3-VI, aadA1,sul1, qacDE N19871 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N15757 AMP,KAN,STR,TET IncI1; IncHI2 aph, sph, strA, strB, bla , tetC, tetD, ble TEM1 N26457 AMC,AMP,FOX,AXO,GEN,KAN,STR,SUL,TET,TIO IncI1; IncA/C aac3-VI, aadA1, aphA1, strA, strB, bla ,sul1, sul2,tetA, tetR(A), qacDE, sugE2 CMY2 N18393 AMP,GEN,TET IncI1 aacC, aadA1, blaTEM1, tetA, tetR(A) N29341 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41578 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 41563 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41567 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41573 AMP,GEN,TET IncI1 aacC, aadA1,bla , tetA, tetR(A) TEM1 41576 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41584 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41579 AMP,GEN,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 41566 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB,bla , tetA, tetR(A) TEM1 41565 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N30678 AMP,GEN,KAN,STR,SUL,TET IncI1; IncHI2 aac3-VI, aadA1, aph, sph, strA, strB, bla , sul1, tetC, tetD, ble, qacDE TEM1 32507 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble 24393 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N20134 AMP,GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 29169 Susceptible None None 24391 KAN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 24359 GEN,KAN,STR,SUL,TET IncI1; IncHI2 aac3-VI, aadA1, aph, sph, strA, strB, sul1, tetC, tetD, ble, qacDE 24388 Susceptible None None 20752 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N18413 AMP,GEN,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N18440 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N18453 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 Genomic Analysis and Virulence Differences GBE 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 plasmidare mapped to the tree. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1063 Hoffmann et al. GBE We also found that the SNP analysis clustered isolates to carry cassette genes (morons) encoding factors that en- together that did not have the same PFGE profiles, such as hance the proliferation and dissemination of the prophage the 2011 outbreak isolates, which were distinct on the PFGE by improving the fitness and /or virulence of the lysogen two-enzyme dendrogram. SNP analysis tightly clustered (Boyd and Brussow 2002; Pelludat et al. 2003). This possibility the clinical and ground turkey outbreak isolates together has led to the assumption that phage-mediated distribution of and distinguished these from other apparently clonal iso- virulence factors and fitness factors is a key driving force in the lates that shared the same PFGE and antibiotic resistance pro- optimization of Salmonella–host interactions and the emer- files. This result confirms the epidemiological data that the gence of new epidemic clones (Brussow et al. 2004). For ex- clinical and ground turkey isolates belong to the same out- ample, lambda-like phage Gifsy-2, which is present among all break, and are separate from other strains grouped with them S. Heidelberg isolates examined in the study, carries the viru- by PFGE. lence factor periplasmic CuZn-superoxide dismutase (SodCI), We found two nonsynonymous SNPs that were present an enzyme that has the potential to enhance the fitness of a only among the outbreak isolates. The first of these SNPs Gifsy-2 lysogen (Brussow et al. 2004). Moreover, the pro- was found to be located on the gluconate transporter gntP phage P4 and prophage P2-like were associated with isolates gene. The biological role of GntP might be to permit the host from group 3 which included the outbreak isolates. The P2- to take up gluconate when only tiny amounts are present, like phage carries the sopE1 gene that encodes part of the thereby conferring an advantage to the host in a competitive T3SS that translocates bacterial effector proteins directly into environment (Klemm et al. 1996). The second of these SNPs the cytosol of host cells. The phage P2-like has also been was located within a gene encoding a putative hydrolase identified in S. Typhimurium strain DT49/DT204 that caused (phosphoserine phosphatase). This gene is involved in glycine a major outbreaks in Europe in 1970s and 1980s. It has been and serine metabolism, and resulted in a proline to glutamine speculated that lysogenic conversion with sopE1 was an im- amino acid change that is known to change the secondary portant step in the emergence of this epidemic strain (Mirold structure of this enzyme. et al. 2001). Interestingly, Mirold et al. reported that sopE1 is SNP-based phylogeny clearly has a distinct advantage in not restricted to a certain bacteriophage as a “vehicle,” iden- subtyping clonal isolates because the method offers nucleo- tifying sopE1 in other serovars present on a Gifsy-like phage tide base resolution. We detected 4,053 SNPs within highly (Mirold et al. 2001). We also saw the same phenomenon clonal S. Heidelberg isolates. Compared with the PFGE analy- among the S. Heidelberg of this study where the sopE1 ses that grouped isolates together having the same Xbaland gene, which was found to be present among all isolates, BlnI patterns, the resolution of the SNP phylogeny conclusively was carried by both the P2-like and the lambda-like Gifsy-2 demonstrated that genome sequencing distinguishes be- phages. tween outbreak and nonoutbreak isolates that shared the The exchanges of genetic cassettes between unrelated same XbaIand BlnI patterns. The results from the ML tree phages further increase the degree of HGT noted among generated from the SNP matrix constructed with kSNP are S. Heidelberg strains, clearly demonstrating the importance concordant with those observed based on the SNP matrix of phages in the emergence and evolution of Salmonella path- developed with our custom pipeline. The study showed that ogens. To date, this is the first time that phages have been SNP analysis is very consistent and reliable, confirming that characterized by name among S. Heidelberg isolates and that even when different methods were used to identify SNPs, a P2-like phage was identified in S. Heidelberg that carries the there are no appreciable changes in the overall phylogeny of sopE1 gene. This phage and, subsequently, the sopE1 deter- the strains included in the analyses. minant, was present in the outbreak isolates, possibly serving Only a few mutation “hotspot” genes were found to have as one of the virulence factors giving rise to these pathogens. a disproportional number of SNPs and that might serve as Furthermore, we identified the presence of 74 virulence de- targets useful for rapidly subtyping serovar Heidelberg. terminants in all S. Heidelberg isolates. Although these genes have several SNPs, the haplotype diver- In the course of performing this study, both newly se- sity is fairly low, indicating there are only few unique quenced plasmids and plasmids sharing a high sequence sim- sequences among S. Heidelberg. This information supports ilarity to previously described plasmids (Han et al. 2012)were the value of SNP analysis as a useful tool for subtyping this identified. Similar to the phage presence and absence tree serovar. topology previously described, the mobilization plasmids ab- WGS data and comparative genomics not only offer a sence/presence seemed to correlate with the tree topology useful subtyping method to differentiate closely related bac- according to the plasmid(s) they carried. For instance, the teria, but they also provide a better understanding of patho- 3,772 bp mobilization plasmid was only identified in isolates genicity and evolution. Our results show substantial loss and belonging to paraphyletic group 2, whereas the 4,473 bp plas- gain of plasmids and phages between these S. Heidelberg mid was found only in isolates which belong to the paraphy- isolates that were consistent with previously generated tree letic group 3 on the ML tree. The mobilization plasmids carry topologies. The phages found among these strains are known the genes necessary to encode the proteins involved with 1064 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE relaxosome formation and processing. Moreover, these plas- stability to assure that the plasmid as well as its associated mids are efficiently transferred by other plasmids, such as IncI1 virulence determinants are inherited by all bacterial progeny and IncX (Francia et al. 2004), two plasmids frequently found (Yamaguchi and Inouye 2011; Unterholzner et al. 2013). As among the S. Heidelberg in this and other studies (Folsteretal. such, these genes serve to stabilize the “pathogenic poten- 2012; Han et al. 2012). The mobilization plasmids tend to be tial” of these isolates. retained and are, continually carried by specific groups of Finding that the T4SS system was present in the majority of S. Heidelberg isolates, suggesting that they must play an es- the pathogenic isolates and carried on plasmids also bearing sential role in this serovar survival, although the precise role is determinants associated with assuring vertical transmission of unclear. Plasmid sequence analysis detected additional plas- the plasmid to progeny bacteria, we sought to assess the role mids that encode several S. Heidelberg virulence determi- of the T4SS in overall pathogenicity. Using the C. elegans vir- nants, including T4SS and antibiotic resistance genes. Unlike ulence assay, we confirmed that S. Heidelberg isolates that the smaller plasmids, the loss and gain of the larger plasmids carried components of the T4SS were more pathogenic in this (ranging from 33 to >200 kb) did not correlate with the tree model than isolates that did not carry the system, suggesting topology. The plasmid analyses performed in this and previous the T4SS plays a role in the pathogenicity of S. Heidelberg. studies (Folsteretal. 2012; Han et al. 2012) demonstrated that Much more inclusive investigations must be performed to fully S. Heidelberg isolates have at least four transmissible plasmids, understand the role of plasmid-encoded genes in pathogen- including an IncA/C, IncI, and IncHI2 and IncX plasmid. The esis, with special emphasis on the role of virulence-associated IncA/C, IncI, and IncHI2 plasmids carry genes that are impor- VirB/D4 T4SS secretory mechanism. tant for virulence, colonization, and persistence, as well as In summary, the data presented in this study support the those encoding resistance determinants to heavy metals, dis- role of transmissible mobile genetic elements in the evolution infection, and antimicrobial agents. We determined that most of virulence and resistance among S. Heidelberg. Understand- of the antimicrobial-resistance phenotypes could be ac- ing the results and the role of these elements is very important counted for by resistance genes carried on plasmids. We did since S. Heidelberg is the second most common serovar, fol- not identify an instance of resistance occurring in the absence lowing S. Typhimurium, responsible for Salmonella-related of a corresponding gene. This shows that genetic prediction of deaths in the United States (Kennedy et al. 2004; Patchanee phenotypic resistance may be possible and may become et al. 2008). Although the mechanisms responsible for necessary as clinical diagnostics moves toward culture- S. Typhimurium pathogenicity and colonization are well stud- independent technologies. More studies are necessary to de- ied, fewer studies are available for S. Heidelberg (Han et al. termine how well specific DNA sequences, and combinations 2011). Presumably, the transmissible plasmids identified in our of sequences, predict resistance or susceptibility to various strain set play a major role in the spread of virulence and antimicrobial agents. resistance genes among S. Heidelberg and, most likely, Besides the IncI1 and the IncHI2 antibiotic resistance plas- other Salmonella serovars. mids, an IncX plasmid that encodes the structural components The results from this and other studies suggest that plas- of the VirB/D4 T4SS apparatus was detected in most of the mids are being exchanged quickly. For example, it was noted outbreak isolates. The VirB/D4 T4SS apparatus is important in that some strains associated with the 2011 outbreak isolated that it assists with “rapid dissemination of antibiotic resistance from the same state had different plasmid profiles. In addition and virulence determinants” (Bhatty et al. 2013). Unlike other to the ecological benefits conferred as a consequence of secretion systems, the T4SS is not only capable of acting in a this rapid plasmid exchange, large plasmids can change conjugative fashion disseminating DNA to other bacteria, PFGE patterns. We noted that, while most of the se- thereby contributing to genome plasticity and the evolution quenced isolates associated with the 2011 outbreak had of infectious pathogens through dissemination of antibiotic BlnI PFGE pattern JF6A26.0076, a single isolate, 41565, had resistance and virulence genes, but also capable of transfer- a different, JF6A26.0017 BlnI pattern, most likely due to ring protein effectors across both the bacterial and eukaryotic the presence of an IncHI2 plasmid rather than IncI1. Other plasma membrane directly into the cytosol of the target cells, closely related isolates not associated with the outbreak car- thereby contributing directly to bacterial pathogenicity (Juhas rying the IncHI2 plasmid also exhibited BlnIPFGE pattern et al. 2008; Gokulan et al. 2013). The T4SS was first described JF6A26.0017. for Agrobacterium tumefaciens and consists of 11 proteins The results determined in this study have convincingly encoded by the virB gene complex. Of particular interest shown that WGS analysis could be a useful tool for identifying was the fact that some plasmids carrying the T4SS, such as the source of contamination and studying the short-term evo- those of the IncX group, also carry systems that confer plasmid lution of these epidemic clones. In the context of the 2011 stability among the strains possessing the plasmid. For in- ground turkey-associated outbreak, SNP analysis provided an stance, the IncX plasmids carry the toxin–antitoxin system excellent tool for subtyping S. Heidelberg and clustering clo- genes (hicA), (hicB), as well as the toxin/antitoxin stability sely related isolates together. Using these data, we deter- genes (stbE), (stbD), both of which contribute to the plasmid mined that the antibiotic-resistant outbreak isolates diverged Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1065 Hoffmann et al. GBE from environmental isolates through the gain of several viru- Disclosure forms provided by the authors will be available lence factors, such as sopE1 gene carried on Phage P2 or the with the full text of this article. This work was supported by T4SS system carried on the IncX plasmid. Of particular interest, the Center for Veterinary Medicine and the Center for Food other isolates in paraphyletic group 3 on the ML tree carried Safety and Applied Nutrition at the US Food and Drug the same virulence factors as the 2011 outbreak strain, sug- Administration and by an appointment of MH by the Joint gesting that this strain has been seen earlier in the environ- Institute for Food Safety and Applied Nutrition (JIFSAN), ment and that the outbreak occurred as a consequence of University of Maryland, College Park, MD. changes in outside influences such as the dose or the handling of the food. Literature Cited CDC projected that for every case of Salmonella reported Aballay A, Yorgey P, Ausubel FM. 2000. Salmonella typhimurium prolifer- that 29.3 cases occur in the community and 28.3 cases goes ates and establishes a persistent infection in the intestine of undiagnosed or unreported (Scallan et al. 2011), suggesting Caenorhabditis elegans. Curr Biol. 10:1539–1542. that the 2011 outbreak actually may have sickened over Allard MW, et al. 2012. High resolution clustering of Salmonella enterica 4,000 people in the United States. The World Health serovar Montevideo strains using a next-generation sequencing ap- Organization has declared that most countries have docu- proach. BMC Genomics 13:32. Allard MW, et al. 2013. On the evolutionary history, population genetics mented a significant increase of microorganism-mediated and diversity among isolates of Salmonella enteritidis PFGE pattern foodborne illnesses (WHO 2013). Therefore, it is paramount JEGX01.0004. PLoS One 8:e55254. to have a robust identification system that can distinguish Alvarez-Martinez CE, Christie PJ. 2009. Biological diversity of prokaryotic closely related strains in an outbreak situation. To date, type IV secretion systems. Microbiol Mol Biol Rev. 73:775–808. DNA “fingerprinting” facilitated by PFGE at PulseNet has Aziz RK, et al. 2008. The RAST Server: rapid annotations using subsystems technology. BMC Genomics 9:75. often been used to establish the source of clinical isolates. Bazinet AL, Cummings MP. 2011. Computing the tree of life—leveraging However, several examples in the past have shown that PFGE the power of desktop and service grids. In: Proceedings of the Fifth may not be as useful for some serotypes that show less di- Workshop on Desktop Grids and Volunteer Computing Systems versitybyPFGE(CDC 2010a, 2010b). Therefore, it is essential PCGrid 2011. Available from: http://lattice.umiacs.umd.edu/lattice- to develop and evaluate other methods for their ability to files/publications/lattice/computing_tree_of_life.pdf. Bazinet AL, Myers DS, Fuetsch J, Cummings MP. 2007. Grid services base effectively and accurately discriminate outbreak isolates library: a high-level, procedural application program interface for writ- from nonoutbreak isolates, like the combination of WGS ing Globus-based Grid services. Future Gener Comput Syst. 23: and computationally efficient methods of analysis. Such a 517–522. combination of approaches also will help to genetically link Beltran P, et al. 1991. Reference collection of strains of the Salmonella the implicated food sources of contamination with farm or typhimurium complex from natural populations. J Gen Microbiol. 137: factories. As more WGS data are gathered and analytical 601–606. Bhatty M, Laverde Gomez JA, Christie PJ. 2013. The expanding bacterial tools developed in combination with epidemiology investiga- type IV secretion lexicon. Res Microbiol. 164(6):620–639. tions, we anticipate significant improvements in timeliness of Boyd EF, Brussow H. 2002. Common themes among bacteriophage- outbreak investigations and attribution of human infections encoded virulence factors and diversity among the bacteriophages to specific sources. involved. Trends Microbiol. 10:521–529. Brussow H, Canchaya C, Hardt WD. 2004. Phages and the evolu- tion of bacterial pathogens: from genomic rearrangements to lyso- Supplementary Material genic conversion. Microbiol Mol Biol Rev. 68:560–602, table of contents. Supplementary tables S1 and S2 and figures S1–S4 are avail- CDC. 2010a. Centers for Disease Control and Prevention. 2010. Multistate able at Genome Biology and Evolution online (http://www. outbreak of human Salmonella Montevideo infections (final ubdate). gbe.oxfordjournals.org/). Atlanta (GA): Centers for Disease Control and Prevention. Available from: http://www.cdc.gov/salmonella/montevideo. CDC. 2010b. Centers for Disease Control and Prevention. 2011. Multistate Acknowledgments outbreak of human Salmonella Enteritidis infections associated with shell eggs (final ubdate). Atlanta (GA): Centers for Disease The authors thank Ruth Timme for submission of the draft Control and Prevention. Available from: http://www.cdc.gov/salmo genomes to NCBI and the NCBI rapid annotation pipeline nella/enteritidis. team, Bill Klimke, Dmitry Dernovoy, Stacy Ciufo, Kathleen CDC. 2011. Centers for Disease Control and Prevention. 2011. Multistate O’Neill, Azat Badretdin, and Tatiana Tatusova, for key outbreak of human Salmonella Heidelberg infections linked to ground turkey. Atlanta (GA): Centers for Disease Control and Prevention. genome annotation service. Further they thank Yuansha Available from: http://www.cdc.gov/salmonella/heidelberg/. Chen for sharing the CVM/DAFM database of antimicrobial CDC. 2013. Centers for Disease Control and Prevention. 2013. resistance genes. No human subjects or animals were used in National enteric disease surveillance: Salmonella annual report, this study. All authors have read the manuscript and agreed to 2011. Atlanta (GA): Centers for Disease Control and Prevention. its contents, subject matter, and author line order. These data Available from: http://www.cdc.gov/ncezid/dfwed/PDFs/salmonella- are novel and have not been previously published elsewhere. annual-report-2011-508c.pdf. 1066 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Chittick P, Sulka A, Tauxe RV, Fry AM. 2006. A summary of national Juhas M, Crook DW, Hood DW. 2008. Type IV secretion systems: tools of reports of foodborne outbreaks of Salmonella Heidelberg infections bacterial horizontal gene transfer and virulence. Cell Microbiol. 10: in the United States: clues for disease prevention. J Food Prot. 69: 2377–2386. 1150–1153. Kennedy M, et al. 2004. Hospitalizations and deaths due to Salmonella Christie PJ, Atmakuri K, Krishnamoorthy V, Jakubowski S, Cascales E. infections, FoodNet, 1996–1999. Clin Infect Dis. 38(Suppl. 3): 2005. Biogenesis, architecture, and function of bacterial type IV secre- S142–S148. tion systems. Annu Rev Microbiol. 59:451–485. Klemm P, Tong S, Nielsen H, Conway T. 1996. The gntP gene of Crump JA, et al. 2011. Antimicrobial resistance among invasive nontyphoi- Escherichia coli involved in gluconate uptake. J Bacteriol. 178:61–67. dal Salmonella enterica isolates in the United States: National Klimke W, et al. 2009. The National Center for Biotechnology Antimicrobial Resistance Monitoring System, 1996 to 2007. Information’s Protein Clusters Database. Nucleic Acids Res. 37: Antimicrob Agents Chemother. 55:1148–1154. D216–D223. Darling AE, Mau B, Perna NT. 2010. progressiveMauve: multiple genome Kurtz S, et al. 2004. Versatile and open software for comparing large alignment with gene gain, loss and rearrangement. PLoS One 5: genomes. Genome Biol. 5:R12. e11147. Liu F, et al. 2011. Novel virulence gene and clustered regularly interspaced Edgar RC. 2010. Search and clustering orders of magnitude faster than short palindromic repeat (CRISPR) multilocus sequence typing scheme BLAST. Bioinformatics 26:2460–2461. forsubtypingofthe majorserovarsof Salmonella enterica subsp. Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters of enterica. Appl Environ Microbiol. 77:1946–1956. individuals using the software STRUCTURE: a simulation study. Mol Marcais G, Kingsford C. 2011. A fast, lock-free approach for efficient Ecol. 14:2611–2620. parallel counting of occurrences of k-mers. Bioinformatics 27: Falush D, Stephens M, Pritchard JK. 2003. Inference of population struc- 764–770. ture using multilocus genotype data: linked loci and correlated allele Mirold S, et al. 2001. Salmonella host cell invasion emerged by acquisition frequencies. Genetics 164:1567–1587. of a mosaic of separate genetic elements, including Salmonella Foley SL, et al. 2011. Population dynamics of Salmonella enterica serotypes pathogenicity island 1 (SPI1), SPI5, and sopE2. J Bacteriol. 183: in commercial egg and poultry production. Appl Environ Microbiol. 77: 2348–2358. 4273–4279. NARMS. 2010. National Antimicrobial Resistance Monitoring System- Folster JP, et al. 2012. Characterization of multidrug-resistant Salmonella Enteric Bacteria (NARMS). Bethesda (MD): Executive Report FDA. enterica serovar Heidelberg from a ground turkey-associated outbreak Norman A, Hansen LH, She Q, Sorensen SJ. 2008. Nucleotide sequence of in the United States in 2011. Antimicrob Agents Chemother. 56: pOLA52: a conjugative IncX1 plasmid from Escherichia coli which en- 3465–3466. ables biofilm formation and multidrug efflux. Plasmid 60:59–74. Francia MV, et al. 2004. A classification scheme for mobilization regions of Paradis E, Claude J, Strimmer K. 2004. APE: analyses of phylogenetics and bacterial plasmids. FEMS Microbiol Rev. 28:79–100. evolution in R language. Bioinformatics 20:289–290. Fricke WF, et al. 2011. Comparative genomics of 28 Salmonella enterica Patchanee P, Zewde BM, Tadesse DA, Hoet A, Gebreyes WA. 2008. isolates: evidence for CRISPR-mediated adaptive sublineage evolution. Characterization of multidrug-resistant Salmonella enterica serovar J Bacteriol. 193:3556–3568. Heidelberg isolated from humans and animals. Foodborne Pathog Gardner SN, Slezak T. 2010. Scalable SNP analyses of 100+ bacterial or Dis. 5:839–851. viral genomes. J Forensic Res. 01:107. Pelludat C, Mirold S, Hardt WD. 2003. The SopEPhi phage integrates Gokulan K, et al. 2013. Impact of plasmids, including those into the ssrA gene of Salmonella enterica serovar Typhimurium A36 encodingVirB4/D4 type IV secretion systems, on Salmonella enterica and is closely related to the Fels-2 prophage. J Bacteriol. 185: serovar Heidelberg virulence in macrophages and epithelial cells. PLoS 5182–5191. One 8:e77866. Pritchard JK, Stephens M, Donnelly P. 2000. Inference of population struc- Gunn JS, Alpuche-Aranda CM, Loomis WP, Belden WJ, Miller SI. 1995. ture using multilocus genotype data. Genetics 155:945–959. Characterization of the Salmonella typhimurium pagC/pagD chromo- Rosenberg NA. 2004. Distruct: a program for the graphical display of somal region. J Bacteriol. 177:5040–5047. population structure. Mol Ecol. 4:137–138. Han J, et al. 2011. Comparison of Salmonella enterica serovar Heidelberg Sahu SN, et al. 2012. Genomic analysis of immune response against Vibrio isolates from human patients with those from animal and food cholerae hemolysin in Caenorhabditis elegans. PLoS One 7:e38200. sources. J Clin Microbiol. 49:1130–1133. Scallan E, et al. 2011. Foodborne illness acquired in the United States— Han J, et al. 2012. DNA sequence analysis of plasmids from multidrug major pathogens. Emerg Infect Dis. 17:7–15. resistant Salmonella enterica serotype Heidelberg isolates. PLoS One 7: Sprando RL, Olejnik N, Cinar HN, Ferguson M. 2009. A method to e51160. rank order water soluble compounds according to their toxicity Hoffmann M, et al. 2012. Genome sequences of five Salmonella enterica using Caenorhabditis elegans, a complex object parametric ana- serovar Heidelberg isolates associated with a 2011 multistate outbreak lyzer and sorter, and axenic liquid media. Food Chem Toxicol. 47: in the United States. J Bacteriol. 194:3274–3275. 722–728. Huehn S, Bunge C, Junker E, Helmuth R, Malorny B. 2009. Poultry-asso- Sunde M, Nortstro ¨ m M. 2005. The genetic backround for streptomycin ciated Salmonella enterica subsp. enterica serovar 4,12:d:- reveals high resistance in Escherichia coli infuences the distribution of MICs. clonality and a distinct pathogenicity gene repertoire. Appl Environ Antimicrob Agents Chemother. 56:87–90. Microbiol. 75:1011–1020. Tamura K, Dudley J, Nei M, Kumar S. 2007. MEGA4: molecular evolution- Huehn S, et al. 2010. Virulotyping and antimicrobial resistance typing of ary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 24: Salmonella enterica serovars relevant to human health in Europe. 1596–1599. Foodborne Pathog Dis. 7:523–535. Timme RE, et al. 2013. Phylogenetic diversity of enteric pathogen Johnson TJ, et al. 2012. Expansion of the IncX plasmid family for improved Salmonella enterica subsp. enterica inferred from genome-wide refer- identification and typing of novel plasmids in drug-resistant ence-free SNP characters. Genome Biol Evol. 5(11):2109–2123. Enterobacteriaceae. Plasmid 68:43–50. Unterholzner SJ, Hailer B, Poppenberger B, Rozhon W. 2013. Johnson TJ, Nolan LK. 2009. Plasmid replicon typing. Methods Mol Biol. Characterisation of the stbD/E toxin-antitoxin system of pEP36, a 551:27–35. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1067 Hoffmann et al. GBE plasmid of the plant pathogen Erwinia pyrifoliae. Plasmid 70: Zhou Y, Liang Y, Lynch KH, Dennis JJ, Wishart DS. 2011. PHAST: a fast 216–225. phage search tool. Nucleic Acids Res. 39:W347–W352. WHO. 2013. Food safety: general information related to microbiological Zwickl D. 2006a. Genetic algorithm approaches for the phylogenetic anal- risks in food. Available from: http://www.who.int/foodsafety/en/. ysis of large biological sequence dataset under the maximum likeli- Wilmshurst P, Sutcliffe H. 1995. Splenic abscess due to Salmonella hood criterion. Austin (TX): University of Texas. Heidelberg. Clin Infect Dis. 21:1065. Zwickl DJ. 2006b. Genetic algorithm approaches for the phylogenetic Yamaguchi Y, Inouye M. 2011. Regulation of growth and death in Escher- analysis of large biological sequence datasets under the maximum ichia coli by toxin-antitoxin systems. Nat Rev Microbiol. 9:779–790. likelihood criterion [PhD dissertation]. [Austin (TX)]: The University of Zhao S, et al. 2008. Antimicrobial resistance in Salmonella enterica serovar Texas. Heidelberg isolates from retail meats, including poultry, from 2002 to Associate editor: Takashi Gojobori 2006. Appl Environ Microbiol. 74:6656–6662. 1068 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Comparative Genomic Analysis and Virulence Differences in Closely Related Salmonella enterica Serotype Heidelberg Isolates from Humans, Retail Meats, and Animals

Download PDF
23 pages

Loading next page...
 
/lp/oxford-university-press/comparative-genomic-analysis-and-virulence-differences-in-closely-IyfeMfGQal

References (74)

Publisher
Oxford University Press
Copyright
Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution 2014. This work is written by US Government employees and is in the public domain in the US.
ISSN
1759-6653
eISSN
1759-6653
DOI
10.1093/gbe/evu079
pmid
24732280
Publisher site
See Article on Publisher Site

Abstract

Salmonella enterica subsp. enterica serovar Heidelberg (S. Heidelberg) is one of the top serovars causing human salmonellosis. Recently, an antibiotic-resistant strain of this serovar was implicated in a large 2011 multistate outbreak resulting from consumption of contaminated ground turkey that involved 136 confirmed cases, with one death. In this study, we assessed the evolutionary diversity of 44 S. Heidelberg isolates using whole-genome sequencing (WGS) generated by the 454 GS FLX (Roche) platform. The isolates, including 30 with nearly indistinguishable (one band difference) Xbal pulsed-field gel electrophoresis patterns (JF6X01.0032, JF6X01.0058), were collected from various sources between 1982 and 2011 and included nine isolates associated with the 2011 outbreak. Additionally, we determined the complete sequence for the chromosome and three plasmids from a clinical isolate associated with the 2011 outbreak using the Pacific Biosciences (PacBio) system. Using single-nucleotide polymorphism (SNP) analyses, we were able to distinguish highly clonal isolates, including strains isolated at different times in the same year. The isolates from the recent 2011 outbreak clustered together with a mean SNP variation of only 17 SNPs. The S. Heidelberg isolates carried a variety of phages, such as prophage P22, P4, lambda-like prophage Gifsy-2, and the P2-like phage which carries the sopE1 gene, virulence genes including 62 pathogenicity, and 13 fimbrial markers and resistance plasmids of the incompatibility (Inc)I1, IncA/C, and IncHI2 groups. Twenty-one strains contained an IncX plasmid carrying a type IV secretion system. On the basis of the recent and historical isolates used in this study, our results demonstrated that, in addition to providing detailed genetic information for the isolates, WGS can identify SNP targets that can be utilized for differentiating highly clonal S. Heidelberg isolates. Key words: outbreak, antimicrobial resistance, plasmid, SNP analysis, trace-back. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution 2014. This work is written by US Government employees and is in the public domain in the US. 1046 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE to determine whether the outbreak isolates belonged to a Introduction new strain with potentially higher pathogenicity or whether Disease caused by nontyphoidal serotypes of Salmonella these isolates were typical members of a common PFGE type enterica is the leading cause of food-related death in the present in retail poultry. In addition, the data show the role of United States (Scallan et al. 2011). Salmonella enterica sup- transmissible mobile genetic elements in the evolution of vir- species (supsp.) enterica serovar Heidelberg (S. Heidelberg) has ulence and resistance among S. Heidelberg. been among the most frequently isolated serovars in clinical cases of salmonellosis, causing an estimated 84,000 illnesses Materials and Methods in the United States annually (Foley et al. 2011; Han et al. 2011). It is frequently isolated from poultry and poultry Bacterial Strains, Growth Condition, and Characterization meat (Zhao et al. 2008). S. Heidelberg was the seventh Isolates were chosen for the study based on similarity or dis- most common serovar isolated from humans in 2011 (CDC similarity by PFGE to S. Heidelberg isolated from a large 2013) and was found in chicken breast, ground turkey, and ground turkey-associated outbreak in 2011 (supplementary eggs, the main sources of S. Heidelberg infections (Chittick fig. S1, Supplementary Material online) (CDC 2011; Folster et al. 2006; NARMS 2010). S. Heidelberg also tends to be et al. 2012). Forty-four isolates of S. Heidelberg were included associated with invasive diseases such as septicemia and myo- in the study (table 1). Isolates were obtained from animal carditis (Wilmshurst and Sutcliffe 1995). After S. enterica (n¼ 9), retail meat (n¼ 27), and human clinical (n¼ 7) sources supsp. enterica serovar Typhimurium (S. Typhimurium), in the United States, and one of the isolates was obtained S. Heidelberg is the serovar of Salmonella most often associ- from an unknown sample source in Brazil. Four isolates ated with Salmonella-related deaths in the United States were collected from 1982 to 1987; 38 of the isolates were (Kennedy et al. 2004; Patchanee et al. 2008; Crump et al. collected from 2002 to 2011, and collection dates for two of 2011). The National Antimicrobial Resistance Monitoring the isolates were unknown. Five of the isolates were from the System (NARMS), which is responsible for monitoring antimi- Salmonella Reference Collection A (SARA) (Beltran et al. crobial resistance in Salmonella estimated that 65% of the 1991), and nine of the isolates were collected in the course S. Heidelberg isolates from ground turkey in 2010 were resis- of the investigation of a ground turkey-associated outbreak in tant to 3 antimicrobial classes. Presently, the antimicrobial 2011 (Folsteretal. 2012). agents to which this serovar is most commonly resistant are Salmonella isolates were cultured on trypticase soy agar ceftriaxone (a drug of choice for treatment), along with resis- (TSA; Becton, Dickinson, NJ) and in trypticase soy broth tance to streptomycin, tetracycline, sulfamethoxazole, chlor- (TSA; Becton) overnight at 37 C. All isolates used for WGS amphenicol, and trimethoprim-sulfamethoxazole (NARMS were serotyped by conventional methods and tested for anti- 2010). microbial susceptibility according to the NARMS standard pro- The Centers for Disease Control and Prevention (CDC) in- tocol as previously described (Zhao et al. 2008). Antimicrobial vestigated a multistate (34 states) outbreak of antimicrobial- susceptibility testing, using a panel consisting of 15 antimicro- resistant S. Heidelberg infections comprised of 136 confirmed bial agents (amikacin, ampicillin, amoxicillin-clavulanic acid, cases between February 27 and September 13, 2011. Among cefoxitin, ceftiofur, ceftriaxone, chloramphenicol, ciprofloxa- the 94 case patients for which there was available informa- cin, gentamicin, kanamycin, nalidixic acid, streptomycin, sulfi- tion, 37 (39%) had been hospitalized and one patient died. soxazole, tetracycline, and trimethoprim-sulfamethoxazole) Collaborative investigative efforts by state and federal officials was performed according to the NARMS methodology with implicated ground turkey as the source of this outbreak, and the Sensititre automated antimicrobial susceptibility system as a result, 36 million pounds of ground turkey meat were (Trek Diagnostic Systems, Westlake, OH). Isolate antimicrobial recalled (CDC 2011; Folsteretal. 2012). resistance was determined by comparison of MICs to values Recently, our investigative capabilities have been greatly established MICs by the Clinical and Laboratory Standards enhanced with the development and increasing feasibility of Institute (CLSI). PFGE was performed according to the CDC whole-genome sequencing (WGS) as a molecular epidemio- PulseNet protocol (http://www.cdc.gov/pulsenet/pathogens/ logical tool to complement current foodborne outbreak inves- index.html, last accessed April 23, 2014). Genomic DNA of tigation techniques. WGS is of particular interest because it each strain was isolated from overnight cultures using DNeasy provides definitive data for distinguishing outbreak isolates Blood and Tissue Kit (Qiagen, CA). All S. Heidelberg isolates were stored in TSB containing 15% glycerol at 80 C. from nonoutbreak isolates in common and highly clonal pop- ulations (Allard et al. 2012, 2013). In this study, we sought to Genome Sequencing, Assembly, and Annotation determine how effectively WGS and single nucleotide poly- morphisms (SNPs) analysis would differentiate outbreak iso- We performed shotgun sequencing of the 44 S. Heidelberg lates of S. Heidelberg from nonoutbreak isolates that share the using the Genome Sequencer FLX 454 (Roche, Branford, CT) same Xbaland BlnI pulsed-field gel electrophoresis (PFGE) pat- and the GS FLX Titanium Sequencing Kit XLR70 according to terns. Using WGS data and virulence assays, we also wanted the manufacturer’s protocol to generate an average genome Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1047 Hoffmann et al. GBE 1048 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Table 1 Metadata of S. Heidelberg Isolates Used in This Study Salmonella enterica susp. Source Source Source Isolated Accession Closed Plasmid BioProject Plasmid Profile Phage Profile enterica serovar and strain Country State Type Date Heidelberg str. 20752 USA NC Animal 2002 AMNR00000000 78501 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 21381 USA CT Food 2002 AMMZ00000000 80505 mob1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24359 USA MO Animal 2003 AMNP00000000 78497 mob1, IncI1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24388 USA MO Animal 2003 AMNQ00000000 78499 mob1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24390 USA MO Animal 2003 AMMY00000000 78475 mob1, IncI P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24391 USA MO Animal 2003 AMNO00000000 78495 mob1, IncI P22, Gifsy-2 like, P4, P2-like Heidelberg str. 24393 USA MO Animal 2003 AMNM00000000 78491 IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 29169 USA CT Food 2003 APIX00000000 80527 mob1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 32507 USA GA Food 2003 AMNL00000000 80525 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41563 USA OR Human 2011 AJGX00000000 pCFSAN000312_01 78487 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41565 USA WA Food 2011 AJHA00000000 80521 mob1, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41566 USA WA Food 2011 AJGZ00000000 78489 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41567 USA WA Food 2011 AMNG00000000 80511 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41573 USA OH Human 2011 AJGY00000000 pCFSAN000316_01 80513 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41576 USA OH Human 2011 AMNH00000000 80515 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41578 USA OH Human 2011 CP004086 pSEEH1578_01,_02,_03 78477 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41579 USA OH Food 2011 AJGW00000000 78479 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 41584 USA MI Human 2011 AMNI00000000 80517 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. 82-2052 USA ME Human 1982 AMMX00000000 78473 IncI1, VirB/D4 P22, Gifsy-2 like, Fels-2, ST64B Heidelberg str. N1514 USA GA Food 2004 AKYM00000000 69871 VirB/D4 P22, Gifsy-2 like Heidelberg str. N1536 USA GA Food 2004 AKYN00000000 69869 IncI, VirB/D4 P22, Gifsy-2 like Heidelberg str. N15757 USA GA Food 2007 AMND00000000 80509 mob1, IncI1, VirB/D4, IncHI2 P22, Gifsy-2 like Heidelberg str. N18393 USA CT Food 2008 AMNF00000000 78485 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18413 USA GA Food 2008 AMNS00000000 80529 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18440 USA MD Food 2008 AMNT00000000 80531 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N18453 USA MD Food 2008 AMNU00000000 80533 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N189 USA CT Food 2004 AMMP00000000 69873 mob2, IncI1, VB/D4 P22, Gifsy-2 like Heidelberg str. N19871 USA NY Food 2008 AMNC00000000 78481 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N19992 USA CA Food 2009 AMMC00000000 66709 ColE1 P22, Gifsy-2 like, P2-like Heidelberg str. N20134 USA GA Food 2009 AMNN00000000 78493 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N26457 USA CO Food 2010 AMNE00000000 78483 mob1, IncI1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. N29341 USA MN Food 2011 AMNJ00000000 80519 mob1, IncI1, VirB/D4 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N30678 USA NM Food 2011 AMNK00000000 80523 mob1, IncI1, VirB/D4, IncHI2 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N418 USA NM Food 2004 AMMB00000000 66707 mob1, IncA/C P22, Gifsy-2 like, P4, P2-like Heidelberg str. N4403 USA CA Food 2005 AMMT00000000 ‘ 69861 mob2, IncI1, VirB/D4 P22, Gifsy-2 like Heidelberg str. N4496 USA CO Food 2005 AMMQ00000000 69867 VirB/D4 P22, Gifsy-2 like Heidelberg str. N4541 USA CT Food 2005 AMMV00000000 69857 VirB/D4 P22, Gifsy-2 like Heidelberg str. N4630 USA GA Food 2005 AMNA00000000 80507 mob1, IncI1 P22, Gifsy-2 like, P4, P2-like Heidelberg str. N653 USA OR Food 2004 AMMW00000000 pCFSAN000412_01, _02 69855 mob2, IncI1, VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA31 USA MD Animal 1987 AMLV00000000 pCFSAN000441_01 66679 VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA32 USA TX Animal 1986 AMLU00000000 pCFSAN000442_01 66677 VirB/D4 P22, Gifsy-2 like Heidelberg str. SARA35 Brazil n/a Unknown n/a AMLT00000000 pCFSAN000443_01 66675 ColE1, 8 kb resistance plasmid P22, Gifsy-2 like Heidelberg str. SARA37 USA CO Animal 1987 AMLR00000000 66673 IncI1, VirB/D4 P22, Gifsy-2 like, Fels-2, ST64B Heidelberg str. SARA 39 USA NC Human n/a AMMJ00000000 69887 IncHI2 P22, Gifsy-2 like, Fels-2 First released (Hoffmann et al. 2012). First released (Timme et al. 2013). Genomic Analysis and Virulence Differences GBE coverage of 23. De novo assemblies were performed using SNP matrix. We determined branch support by performing Roche’s Newbler software (v.2.6) with the resulting contigs 1,000 bootstrap replicates. being annotated using the NCBIs Prokaryotic Genomes We elucidated the relationships among isolates using the Automatic Annotation Pipeline (Klimke et al. 2009). Bayesian clustering method implemented in STRUCTURE Using DNA from the 2011 clinical outbreak isolate 41578, v2.3.2 (Pritchard et al. 2000; Falush et al. 2003). The analyses we also sought to close the genome using the Pacific were based on the kSNP matrix, and we ran ten replicate Biosciences (PacBio) sequencing platform. Specifically, we pre- analyses at K¼ 2–9 under the admixture model with corre- pared a single 10-kb library that was sequenced using the C lated allele frequencies being performed and visualized with chemistry on eight single-molecule real-time cells with a 90- DISTRUCT v1.1. (Rosenberg 2004). Each independent run min collection protocol on the PacBio RS. The 10-kb continu- consisted of 50,000 generations serving as burnin followed ous-long-read data were de novo assembled using the PacBio by 100,000 generations. The K statistic (Evanno et al. 2005) hierarchical genome assembly process/Quiver software pack- was used to identify the k value that best fits the data. Pairwise age, followed by Minimus 2, and they were polished with genetic distances, calculated as the number of nucleotide dif- Quiver. ferences, 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. Comparative and Phylogenetic Genome Analysis Alignments of core genes were accomplished using For comparative analysis, in addition to the 44 isolates se- MUSCLE with default settings; these alignments were then quenced in this study, we included sequence data available used to identify variable core genes. 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 Prediction of Prophages, Plasmids, Resistance, and (ABEL00000000) at the NCBI genome database. Virulence Genes Phylogenetic informative SNP sites (i.e., SNPs shared by two Prophages were identified using PHAST (Zhou et al. 2011). or more strains in the alignment) were identified by two dif- Sequences for intact phages were extracted from the original ferent methods. The first involved mapping the 454 reads to contig, in which they were found and mapped against the the complete reference genome of S. Heidelberg strain SL476 entire assembled data to determine whether the phage was using Roche Newbler software GS Reference Mapper (v.2.8.). present in an isolate at multiple sites on different contigs. SNP sites were then found using the SNP calling function of Contigs that could not be aligned to the reference genome GS Reference Mapper. SNP positions were defined, where were evaluated using BLAST to identify plasmids. The plasmids one or more isolates differed from strain SL476 with coverage were analyzed by comparative analysis using MAUVE algo- 10 and with 95% of the reads containing the SNP after rithm (Darling et al. 2010) and with the comparative analysis having filtered out the SNPs from homopolymer artifacts, fol- tools of RAST (Aziz et al. 2008). For plasmid closure, a con- lowed by a custom pipeline to construct a SNP matrix (Allard catenated sequence was generated from a 500-bp sequence et al. 2013). S. Heidelberg SL486 could not be included in this from each end of the contig. This artificial sequence from the analysis as its raw sequence data were not available. The single contig of the plasmid in study was used as a reference second SNP detection method was a reference free k-mer- to map all the 454 raw reads, using runMapping from based approach implemented in the program kSNP Newbler. If there were mapped reads that covered the con- (Gardner and Slezak 2010). This is a collection of Perl scripts junction point in the reference, the contig was a closed that aggregate the results from Jellyfish v1.1.3 (used for k-mer plasmid. counting) (Marcais and Kingsford 2011) and MUMmer v3.22 To determine the incompatibility (Inc) groups for plasmids, (used to align k-mers and detect variable positions) (Kurtz we used BLAST to find sequences described by Johnson and et al. 2004). Analyses with kSNP were based on a k-mer Nolan (2009) for specific Inc groups that would produce the- length of 25. oretical PCR amplicons for known Inc group sequences. To construct evolutionary relationships among the isolates, Resistance and virulence genes were identified by mapping the maximum-likelihood (ML) method was implemented in sequence data available at an in-house database, consisting the Genetic Algorithm for Rapid Likelihood Inference of 1,379 resistance genes, and 107 previously characterized (GARLI) software (Zwickl 2006b). GARLI analyses were per- virulence genes (84 pathogenicity and 23 fimbrial markers; formed using a web service (Bazinet and Cummings 2011) Huehn et al. 2009). Using a presence/absence matrix of that uses a special programming library and associated tools (Bazinet et al. 2007). All ML trees were constructed with genes conferring resistance to antibiotics and disinfectant GTR+I+G nucleotide substitution model. The 100 replicate agents, we constructed a similarity tree based on binary dis- runs of the nonbootstrapped data set were conducted to tances under neighbor-joining algorithm for tree construction; identify the most probable phylogeny based on our observed topological support was assessed based on 100 bootstrap Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1049 Hoffmann et al. GBE replicates. The analyses were performed using the R package whereas 20 isolates had composite XbaI/BlnI PFGE patterns ape (Paradis et al. 2004; Liu et al. 2011). 98% similar to isolates from the 2011 outbreak (supplemen- tary fig. S1, Supplementary Material online). Nine reference isolates showed the same PFGE pattern, as well as the same Caenorhabditis elegans Assay antibiotic resistance profile, noted for most of the outbreak Pathogenicity was evaluated using the Caenorhabditis elegans isolates (supplementary fig. S1, Supplementary Material survival assay. Caenorhabditis elegans strain SS104 glp-4 (bn2) online). was acquired from the Caenorhabditis Genetics Center. Worm cultures were maintained at 16 C, which is the per- Draft Genome Size missive temperature for this temperature-sensitive sterile mutant of C. elegans. Strains were cultured in C. elegans hab- Shotgun sequencing produced 44 draft genomes with an av- itation media (CeHM) in tissue culture flasks on a platform erage coverage of 23 and a minimum of 4,700 genes. The shaker (Sprando et al. 2009). Nematodes were bleached genomic size of the 44 S. Heidelberg isolates (including (0.5 M NaOH, 1% hypochlorite) to collect eggs, which were extrachromosomal DNA) ranged from 4.671 to 5.130 Mb incubatedinM9media for24 h to bringthem to synchronized (fig. 1). Differences in genomic size were due primarily to L1 stage and then transferred to CeHM. To produce synchro- the presence or absence of mobile genetic elements, including nized L4 stage worms, L1 worms were grown for an addi- phages and plasmids, especially the carriage of IncHI2 or IncA/ tional 72 h in CeHM. C plasmids (fig. 1). Next to thegenome sizevariation, figure 1 Pathogen lawns for survival assays with test strains, as well shows the estimated N50 sizes within S. Heidelberg draft as the food bacteria OP50, were prepared by inoculating genome sequences. The estimated N50 value is a rough esti- Nematode Growth Media (in 6-cm Petri plates) with 50 mlof mate of the quality and coverage of the draft genomes, and it an overnight bacterial culture. Plates were incubated over- represents the average contig size after assembly with the night at room temperature before worms were added. We Newbler software. used 60–80 synchronized L4 worms for each treatment group. Worms were scored every 24 h for survival (Aballay Phylogenetic Relationships among S. Heidelberg Using et al. 2000). SNP Analysis Animal survival was plotted using Kaplan–Meier survival Phylogenetically informative SNP sites were identified using curves and analyzed by log rank test using GraphPad Prism two independent methods. Figure 2 shows a ML tree based (GraphPad Software, Inc., La Jolla, CA). Survival assays were on SNP analysis of 44 S. Heidelberg isolates mapped to repeated at least two times. Survival curves resulting in strain SL476. Among the collective 40,716 variable SNPs P values< 0.05 relative to control were considered signifi- identified, 12,187 were determined to be “informative” cantly different. (i.e., SNPs shared by at least two isolates). The ML tree shown in figure 2 clearly demonstrates that S. Heidelberg Results formed a monophyletic group distinct from the outgroup comprised Salmonella enterica supsp. enterica serovar Antimicrobial Resistance and PFGE Two representative isolates from the 2011 ground turkey-as- sociated outbreak were previously reported to be resistant to ampicillin, gentamicin, streptomycin, and tetracycline (Folster et al. 2012). Antimicrobial susceptibility testing of seven addi- tional 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 FIG.1.—The number of assembled bases (Mb) and N50 contig size associated isolates exhibited indistinguishable, or nearly indis- (kb) for each sequenced S. Heidelberg isolate. Samples are colored accord- tinguishable (one band difference), XbalPFGEpatterns ing to the presence of antimicrobial resistance plasmids. No antimicrobial JF6X01.0032 or JF6X01.0058, and BlnI pattern JF6A26.0076 resistance plasmid, filled triangles; antimicrobial resistance Inc-I plasmid, and one food isolate had Xbal PFGE pattern JF6X01.0058 and filled diamonds; antimicrobial resistance Inc-H1/2 plasmid, filled circles; BlnI pattern JF6A26.0017. Of the 35 isolates not associated antimicrobial resistance plasmid Inc-I and H1/2, filled circles with borders; with the outbreak, 15 had composite Xbal/BlnI PFGE patterns antimicrobial resistance plasmid Inc-AC, filled squares; antimicrobial resis- tance plasmid IncI and Inc-AC, filled squares with borders. indistinguishable from ground turkey outbreak isolates 1050 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE 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. Newport (S. Newport) and S. Typhimurium. The SNP diver- trend for temporally isolated strains to cluster closely together, gence between S. Heidelberg and the two outgroups, one isolate, N29341, that was isolated in Minnesota in 2011 S. Newport and S. Typhimurium, is 24,724 and 24,305 from ground turkey by NARMS, clustered tightly together SNPs, respectively. The resultant tree shows two more inter- with the outbreak isolates received from CDC, whereas all esting points. First, isolates having a similar XbalPFGEpattern other reference isolates having the same PFGE pattern not (supplementary fig. S1, Supplementary Material online) clus- associated with the ground turkey outbreak are distinct tered together. Second, there is a tendency for isolates to from the isolates of the 2011 outbreak (fig. 2). These results cluster based on date of collection, such as noted for isolates suggest that N29341 might be epidemiologically related to from 2003, 2008, and 2011. In addition, the 2011 clinical and the 2011 outbreaks strains, as Minnesota was one of the food isolates that were originally thought to be associated states involved in that outbreak. with the ground turkey outbreak isolated in different states The SNP matrix, including the same S. Newport and clustered together with 80% bootstrap support (fig. 2), which S. Typhimurium outgroups, built using the k-mer strategy of could indicate that they are all members of a clonal commu- kSNP consisted of 41,432 variable SNP positions with 10,873 nity likely derived from the same source. Consistent with the being informative. Using the k-mer strategy, we identified 716 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1051 Hoffmann et al. GBE more SNPs than were noted within the in-house pipeline, present in genes shared by all 46 S. Heidelberg isolates and the where 40,716 variable SNPs were observed. Notably, the two outgroups. Only 718 SNPs, with 284 being informative, SNP diversity within the highly clonal S. Heidelberg isolates were found in genes shared only by S. Heidelberg isolates. The was lower, only 4,053 SNPs were found with 1,394 being SNP matrix of the core genes was used to construct a ML informative. The resulting ML tree (fig. 3) was well resolved phylogenetic tree (supplementary fig. S2, Supplementary and congruent with the tree topology in figure 2.Toshow Material online), which also confirms the tree topology in that the tree topology is not an artifact of mobile elements figures 2 and 3. that are not present in all isolates, we constructed a kSNP In this study, we determined that all 46 S. Heidelberg iso- matrix using SNPs present only on core genes. This produced lates, including S. Heidelberg SL476 and S. Heidelberg SL486, a matrix of 24,447 SNPs, of which 7,630 were informative and partitioned into three paraphyletic groups, each having 100% FIG.3.—ML tree based for the 44 S. Heidelberg isolates and two previously reported S. Heidelberg genome sequences SL476 and SL486 (Fricke et al. 2011). A total of 4,053 SNPs with 1,394 being informative were found based on k-mer analysis using kSNP. ML trees were generated as described in figure 2. Bootstrap values (1,000 iterations) are reported below each node. The numbers of unambiguous substitutions that mapped to the tree only once and are greater than zero are given above each node in blue. The numbers in parenthesis represent the nodes in table 3. To the right of the tree, two Distruct plots were reconstructed with the same SNP matrix—one including all 46 S. Heidelberg isolates and, adjacent to that, another with only those isolates from group 3—to present a fine-scale structure is shown. The Distruct plot was generated using a model-based Bayesian clustering method implemented in Structure v2.3.2 and visualized with DISTRUCT v1.1. 10 replicate analysis at K¼ 2–9 under the admixture model with correlated allele frequencies were performed. Each independent run consisted of 50,000 generations serving as burnin followed by 100,000 generations. Different colors represent the different clusters and each bar represents an individual isolate. The fraction of the bar that is a given color represents the coefficient of membership to that cluster (e.g., multicolored bars indicate membership to multiple groups indicative of admixture). 1052 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE bootstrap support (fig. 3). Group 1 contains three historical 41578 and 41579 are quite distinct from isolates N19871, N18413, N18440, and N18453, which have 44 SNPs that isolates (two clinical and one turkey isolate) with a PFGE pat- allow complete discrimination between the two groups tern, for which there was no assigned Pulse-Net designation (fig. 3). Furthermore, isolate 41565, which has the same (supplementary fig. S1, Supplementary Material online). PFGE and antibiotic resistance profile as the 2003 isolates Group 2 is comprised of 11 isolates, including clinical, food, 32507, 24393, and 20752, is more closely related to the and animal isolates, with XbaI PFGE pattern JF6X01.0022, as other isolates associated with the 2011 outbreak, even with well as the isolates SL476 and N19992 having different XbaI having a different BlnI pattern and antibiotic resistance profile patterns JF6X01.0133 and JF6X01.0326, respectively. Group (fig. 3). 3 is composed of 30 isolates from food, animal, and clinical To provide an additional picture as to the number of distinct sources, including the outbreak isolates, all having a very sim- groups the sampled serovars segregate into, we constructed ilar XbaI pattern (JF6X01.0032, JF6X01.0034, JF6X01.0058). groups using the program STRUCTURE, the results of which The distribution of PFGE patterns on the tree topology were visualized using DISTRUCT (fig. 3). The results are con- based on SNPs differs between the two enzymes XbaIand sistent with the groupings formed through the phylogenetic BlnI. For example, N19992 (group 2) and N18393 (group 3) analyses in that there are three distinct groups (each repre- are quite distinct on the ML tree, even though both have the sented by a different color [fig. 3]). same BlnI pattern JF6A26.0015 and are adjacent to each other on the PFGE dendrogram (supplementary fig. S1, Differences between Outbreak Isolates and Nonoutbreak Supplementary Material online). Pairwise SNP variation be- Isolates in Group 3 tween the three distinct S. Heidelberg lineages seen in figure 3 are presented in table 2 and support the phylogenetic To further investigate the relationships among outbreak iso- partitions revealed above. Mean divergence among the three lates, we constructed the Distruct plot for group 3. Within that clades ranged from 200 to 712 nt differences, whereas the group, there is support for four subgroups; these four groups mean intragroup nucleotide differences ranged from 47 to also are found in the phylogenetic analysis, each with a75% 155 (table 2). bootstrap support. Pairwise SNP variation between the four Notably, the mean SNP variation within the outbreak iso- sublineages also support the phylogenetic partitions described lates is only 17 SNPs. The phylogenetic analysis clustered the above (table 2 and fig. 3). Mean divergence among the four outbreak isolates together even though the PFGE profile (sup- subclades ranged from 36 to 153 nt differences, whereas the plementary fig. S1, Supplementary Material online) is slightly mean intragroup nucleotide differences ranged from 15 to 30 different (XbaI pattern JF6X01.0032, JF6X01.0058; BlnIpat- (table 2). The 2008 strains seem to have quite a unique SNP tern JF6A26.0017 and JF6A26.0076). For example, isolates profile when compared with the remaining isolates of group 41578 and 41579, sharing the same PFGE and antibiotic re- 3. The ML tree, the Distruct plot, and the pairwise SNP varia- sistance profile with isolates N19871, N18413, N18440, and tion show the evolutionary changes occurring among isolates N18453 (isolated in 2008), are more closely related to the from 2003 to 2009 and 2011, revealing the apparent emer- other outbreak isolates, even though they have a different gence of the genetically unique S. Heidelberg strain that was XbaI pattern. We also found that the outbreak isolates responsible for the 2011 outbreak. Table 2 (Mean Pairwise Distance (Number of Nucleotide Differences) between the (A)Three Major S. Heidelberg Groups in figure 3 and the (B) Four S. Heidelberg Subgroups from Group 3 in figure 3 S. Heidelberg Gr. 1 S. Heidelberg Gr. 2 S. Heidelberg Gr. 3 (A) S. Heidelberg Gr. 1 155 S. Heidelberg Gr. 2 535 101 S. Heidelberg Gr. 3 712 200 47 Gr. 3/I Gr. 3/II Gr. 3/III Gr. 3/IV (B) S. Heidelberg Gr. 3/I n/c S. Heidelberg Gr. 3/II 144 30 S. Heidelberg Gr. 3/III 139 42 15 S. Heidelberg Gr. 3/IV 153 46 51 27 Outbreak isolates 151 46 65 36 Intragroup mean pairwise distance (number of nucleotide differences). Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1053 Hoffmann et al. GBE Table 3 lists variable genes having unique nucleotide sub- protein (15 SNPs), terminase small subunit (14 SNPs), outer stitutions which defined specific sublineages that could be membrane lipoprotein Blc (7 SNPs), and alpha amylase used to separate them from other highly clonal (6 SNPs), were found to have a disproportionate number of S. Heidelberg isolates in group 3. Table 3 also lists their SNP SNPs that could be useful for rapid subtyping the serovar change and, if applicable, the gene carrying the SNPs and the Heidelberg and/or as targets of resequencing. Interestingly, resulting amino acid change. Twenty-three unique SNPs, with although these genes have several SNPs, the haplotype diver- 14 being nonsynonymous (i.e., causing an amino acid change) sity is fairly low, with only a few haplotypes (2–4) among 46 SNPs, were found to be unique to group 3 and clustering isolates being identified. This fact suggests that, even though together isolates with XbaI pattern (JF6X01.0032, thereare geneswithup to35SNPs, thereare only twotofour JF6X01.0034, JF6X01.0058). These nucleotide substitutions unique sequences among S. Heidelberg that are congruent were found on 20 different genes that could be useful for with our tree topology for three different groups of subtyping of strains belonging to group 3 (table 3 and fig. 3). S. Heidelberg. Furthermore, 24 unique SNPs (11 being nonsynonymous) were found at the node that characterizes the 2008 ground Complete Genome from Clinical Isolate 41578 turkey isolates, which seemed to have their own unique SNP profile distinct from that of the remaining group 3 isolates. In the course of performing this study, we demonstrated that Twelve of the 24 SNPs are located on the insertion sequence it was possible to rapidly determine the complete sequence for (IS) 2 element with five SNPs being within the IS2 TnpA trans- the chromosome and three plasmids from a clinical isolate posase and seven SNPs on the IS2 TnpB transposase genes associated with the 2011 outbreak using the Pacific (table 3 and fig. 3). Biosciences (PacBio) system. A single contig of 4,793,479 bp Importantly, two nonsynonymous SNPs were found to be (GC content 52.2%) representing the complete chromosome only present in all isolates associated with the 2011 outbreak. and three contigs of 4,773 bp, 35,297 bp, and 117,929 bp These two SNPs were located within two different variable representing three plasmids were generated. The chromo- core genes: a putative hydrolase with amino acid change some contains 12 different Salmonella pathogenicity islands P653Q, and the gluconate transporter GntP gene with and four prophages (prophage P22, lambda-like prophage amino acid change L187V (table 3 and fig. 3). Five SNPs, Gifsy-2, prophage P4, and prophage P2-like; fig. 5). The larg- with four being nonsynonymous, were unique to isolates est plasmid pSEEH1578_01 is a resistance plasmid with incom- 41563, 41566, 41567, 41573, 42578, 41579, and N29341 patibility group IncI1 encoding the resistance to gentamicin (associated to the outbreak) and not observed in 41565 (aacC), streptomycin (aadA1), tetracycline (tetA, tetR(A)), (table 3 and fig. 3). From these five SNPs, three were in plas- and ampicillin (bla )(fig. 5). The plasmid pSEEH1578_02 tem1 mid sequences, whereas two were found on variable core is a VirB/D4 plasmid that carries the type IV secretion system genes: the flagellin CDS resulting in an amino acid change (T4SS; fig. 5), and plasmid pSEEH1578_03 is a mobilization from threonine to serine at position 282 (T282S), and plasmid that carries genes capable of promoting plasmid mo- within a putative hypothetical protein (possibly functioning bilization such as mbeA, mbeB, mbeC,and mbeD. as H+ gluconate symporter related to permeases) producing an amino acid change from glycine to serine at position 358 (G358S) (table 3). Identification of Prophages PHAST analysis identified six different phages among the Genetic Variation among S. Heidelberg S. Heidelberg isolates. Only those phages that PHAST deter- mined as being intact were further analyzed. All 44 isolates At a minimum, 273 genes involved in a variety of different contained prophage P22 and lambda-like prophage Gifsy-2. functions, such as DNA replication and repair, cell division, Group 1, which was composed of three historical isolates transcription, metabolism, and virulence, vary among the (SARA 39, SARA 37, and 82-2052), also contained prophage S. Heidelberg isolates observed in clinical, animals, and/or Fels-2, not observed in any of the group 2 or group 3 isolates. food samples, including the isolates associated with the Additionally, SARA 37 and 82-2053 also contained prophage 2011 outbreak. The variable core genes, including the ST64B (table 1 and fig. 2). number of SNPs and haplotypes are listed in supplementary Heidelberg isolates belonging to group 3, which includes table S1, Supplementary Material online. Figure 4 shows a the 2011 outbreak isolates, contained prophage P4 and P2- histogram comparing strains based on: 1) the number of like, neither of which was seen in any of the isolates from SNPs according to the number of core genes and 2) the hap- lotype diversity according to the number of core genes. Only a group 1 and 2, excepting isolate N19992, a group 2 isolate few mutation “hotspot” genes, such as tailspike (35 SNPs), that was confirmed to contain the P2-like prophage (table 1 scaffolding protein (22 SNPs), DNA polymerase V subunit and fig. 2). The phage presence and absence matrix shows umuD (20 SNPs), phage lysozyme (18 SNPs), gifsy-2 prophage correlation with the tree phylogeny shown in figure 2. 1054 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Table 3 Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study Node on Tree SNP Location NT Change AA Change NT Position Locus Tag 1 D-Mannonate oxidoreductase A!GD!G 839 SEEHRA37_03221 1 Translocation protein TolB T!A P 1152 SEEHRA37_24108 1 Glutathione S-transferase C!AW!Stopp 233 SEEHRA37_20055 1 Nonannotated C!TR!C n/a SEEHRA37_contig25 1 Hypothetical protein/putative transcriptional regulator G!TV!F 376 SEEHRA37_19825 1 Ureidoglycolate dehydrogenase AllD C!TT!I 20 SEEHRA37_20355 1 Starvation-inducible outer membrane lipoprotein C!T S 63 SEEHRA37_07830 1 ATPase recombination and repair protein G!T L 966 SEEHRA37_05779 1 DNA-specific endonuclease I A!CK!T 341 SEEHRA37_19810 1 Scaffolding protein C!A V 135 SEEHRA37_21903 1 Peptidase T G!TV!L 976 SEEHRA37_13506 1 Predicted metal-dependent membrane protease C!A A 108 SEEHRA37_15707 1 Hypothetical protein A!GN!D 418 SEEHRA37_04201 1 Glycine/serine hydroxymethyltransferase GlyA T!GS!A 145 SEEHRA37_00365 1 Preprotein translocase subunit SecF G!A A 534 SEEHRA37_22413 1 Hypothetical protein A!GT!A 1879 SEEHRA37_01714 1 Ribosomal protein L21 RplU A!CI!L 220 SEEHRA37_20845 1 Nonannotated C!GQ!E n/a SEEHRA37_contig19 1 Isoaspartyl peptidase C!TA!V 380 SEEHRA37_16329 1 Lysine-N-methylase C!T G 57 SEEHRA37_07110 1 Hypothetical protein T!C L 162 SEEHRA37_23678 1 Hypothetical protein A!CK!T 227 SEEHRA37_07475 1 Nonannotated T!CI!V n/a SEEHRA37_contig14 2 Hypothetical protein A!GQ!R 851 SEEHRA37_20355 2SensorkinaseDpiB A!TI!F 67 SEEHRA37_10675 2 Notannotated C!T P 63 SEEHRA37_contig6 2 Hypothetical protein G!A F n/a SEEHRA37_contig6 2 Nonannotated G!AG!S n/a SEEHRA37_contig5 2 Proline dehydrogenase PutA C!A G 975 SEEHRA37_04036 2 Electron transport complex protein RnfC T!C I 276 SEEHRA37_02971 2 Protein involved in chromosome partitioning MukB G!CG!A 1514 SEEHRA37_14706 2 N-ethylmaleimide reductase A!T G 351 SEEHRA37_02866 2 Hypothetical protein C!T L 1039 SEEHRA37_16509 2 Multidrug efflux system subunit MdtC A!TQ!H 2067 SEEHRA37_00060 3 Putrescine/spermidine ABC transporter ATPase G!TP!Q 47 SEEHRA37_13511 protein PotA 3 Pilus assembly protein, porin PapC G!TR!S 4 SEEHRA37_18324 3 IS2 repressor TnpA G!A T 93 SEEH8393_18082 3 IS2 transposase TnpB A!G T 195 SEEH8393_18077 3 IS2 transposase TnpB A!GN!D 502 SEEH8393_18077 3 IS2 transposase TnpB G!A A 483 SEEH8393_18077 3 Not annotated (plasmid sequence) T!CA!V n/a SEEH8393_contig56 3 IS2 repressor TnpA C!T Y 309 SEEH8393_18082 3 IS2 repressor TnpA A!G E 195 SEEH8393_18082 3 IS2 transposase TnpB C!A V 525 SEEH8393_18077 3 IS2 transposase TnpB C!A P 192 SEEH8393_18077 3 IS2 transposase TnpB A!G A 789 SEEH8393_18077 3YacA T!G V 123 SEEH8393_01004 3 Hypothetical protein G!AA!T 232 SEEHRA37_09669 3 Multifunctional fatty acid oxidation complex G!TM!I 2094 SEEHRA37_15762 subunit alpha FadJ (continued) Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1055 Hoffmann et al. GBE Table 3 Continued Node on Tree SNP Location NT Change AA Change NT Position Locus Tag 3 tRNA uridine 5-carboxymethylaminomethyl G!AG!S 1199 SEEHRA37_11040 modification enzyme GidA 3 Permease DsdX C!AA!D 1166 SEEHRA37_02334 3 IS2 transposase TnpB T!C R 549 SEEH8393_18077 3 IS2 repressor TnpA C!T A 225 SEEH8393_18082 3 IS2 transposase TnpB C!AP!T 319 SEEH8393_18077 3 IS2 repressor TnpA C!T L 102 SEEH8393_18082 3 Arginyl-tRNA synthetase ArgS A!T G 237 SEEHRA37_07325 3 Nonannotated C!AA!N n/a SEEHRA37_contig18 3 Hypothetical protein G!AD!N 103 SEEHRA_37_05617 4 Not annotated G!AV!L n/a SEERA37_contig2 4 Sodium/glucose cotransporter G!AS!N 1037 SEERA37_14011 5 Phosphoserine phosphatase C!AP!Q 653 SEEHRA37_11265 5 Gluconate transporter GntP C!GL!V 187 SEEHRA37_18800 6 Flagellin and related hook-associated proteins G!CT!S 282 SEEHRA37_18390 6 Electron transfer flavoprotein, beta subunit G!AG!S 358 SEEHRA37_02456 6 Hypothetical protein (plasmid sequence) C!AL!M 562 SEEH8393_17461 6 Not annotated (plasmid sequence) G!AG!R n/a SEEH8393_contig56 6 Not annotated (plasmid sequence) C!TG!S n/a SEEHRA39_contig149 7 Type III secretion protein SopE T!AH!Q 560 SEEHRA37_14371 FIG.4.—(A) Histogram showing the number of SNPs per core genes. (B) Histogram showing haplotype diversity for all variable, core genes. Identification of Virulence Genes among S. Heidelberg S. Heidelberg isolates (supplementary table S2, Supplementary From the 107 known and recognized Salmonella-associated Material online). One outer membrane fimbrial usher gene, virulence genes, we identified 62 pathogenicity and 13 fim- safC, was only found in the three historical isolates (SARA 39, brial markers that were highly conserved among the 44 SARA 37, and 82-2052), all belonging to group 1. All other 1056 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE A B SPI-1 SPI-9 SPI-13 SPI-12 P-4 SPI-3 S. Heidelberg 41578 P-2 CP004086 tetA pSEEH1578_01 tetR(A) SPI-4 SPI-2 aadA1 SPI-11 SPI-5 blaTEM1 SPI-14 P-22 SPI-6 Gifsy-2 aacC SPI-16 virB3-4 virB5 virB2 virB6 virB1 virB8 virB9 virD2 virB10 hicB hicA virB11 virD4 pSEEH1578_02 stbD repA 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 areshownindarkblueand 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. detected virulence genes, observed to be located on the chromosomal msgA gene, which is essential to Salmonella Salmonella pathogenicity islands (SPI-1, SPI-2, SPI-3, SPI-4, mouse virulence (Gunn et al. 1995), and is common in S. and SPI-5), the lambda-like prophage Gifsy-2, and/or islets, Typhimurium (Huehn et al. 2010). Other genes present were present in all S. Heidelberg isolates, suggesting that all among all isolates encoded proteins involved in type III secre- isolates are pathogenic (supplementary table S2, Supplemen- tion system (T3SS) and adhesins. Remarkably, the sopE1 viru- tary Material online). For example, all isolates carried a lence gene, involved in the translocation of effector proteins Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1057 Hoffmann et al. GBE into host cells, was found in all of our strains. Moreover, al- isolates belonging to group 3 in figure 4, while subclade 2 though sopE1 is known to be carried by P2-like phages, in this contains one clinical isolate and isolate N653, the unique iso- study we also identified this gene on lambda-like phage Gifsy- late found to carry two VB/D4 plasmids (the other plasmid partitioned into clade III). Moreover, using the primers de- 2, which is present in all S. Heidelberg isolates. Consequently, signed by Johnson et al. (2012) for the relaxase gene (taxC) all S. Heidelberg isolates have the sopE1 virulence gene. Group (VirD2 component) to screen for the four different IncX sub- 3 isolates, which also contain the P2-like phage, have the groups, we determined that the VB/D4 plasmids from clade I, sopE1 gene twice in their chromosome; one gene is carried subclade 1 from clade II, and clade III belonged to the IncX1 by lambda-like phage Gifsy-2 and the other one is located on plasmid subgroup, whereas the two plasmids from clade II the P2-like phage. Analysis confirmed that these two genes subclade 2 were deemed to be members of the IncX4 shared a high degree of similarity having no more than eight subgroup. mutations between the two different sopE1 alleles. In con- Interestingly, isolate N653 carried two plasmids, one being trast, virulence determinants known to be associated with in the IncX1 and another in the IncX4 incompatibility groups. It SPI-7, prophage Gifsy-1, Gifsy-3, Fels-1, and Salmonella plas- appears that there are significant differences among these mid virulence (spv) were not identified in any of the 44 isolates groups that allow simultaneous maintenance of two plasmids examined in this study, confirming sequencing results that belonging to the same incompatibility group, a heretofore none of the strains carried any of these mobile genetic believed unachievable occurrence. Given the notion that elements. two plasmids in the same Inc group cannot simultaneously coexist within the same bacterium, another plausible explana- Description of VirB/D4 Plasmids tion could be that IncX4 might actually belong to a different Several plasmids were found among the 44 S. Heidelberg incompatibility group other than that of IncX1, which was isolates sequenced in this study (table 1). Electrophoresis of only previously experimentally confirmed to belong to IncX samples obtained using a basic alkaline plasmid purification (Norman et al. 2008). Consistent with this hypothesis was procedure demonstrated the presence of endogenous plas- the fact that a BLAST search with the IncX4 repA gene (pre- mids that had not integrated into the host genome (data not sent only in both putative IncX4 plasmids) sequence identified shown). The plasmid sizes varied from 3.7 to over 200 kb. significant sequence similarity to only the Citrobacter roden- Twenty-one S. Heidelberg isolates, including all five human tium ICC168 plasmid pCROD2 (previously grouped to IncX4 isolates and two ground turkey isolates associated with the [Johnson et al. 2012]) and two plasmids from S. Heidelberg 2011 outbreak, contained a plasmid ranging in size from 33 to (pSH163_34, pSH696_34), suggesting that the IncX4 group is 40 kb, which carried genes (virB1,virB2,virB3-4,virB5,virB6, a new incompatibility group distinct from that of IncX1. If this virB7,virB8,virB9,virB10,virB11,virD2, and virD4) associated is confirmed to the case, coexistence of these two “IncX” withthe VirB/D4T4SS(table 1). Three specific roles for the plasmids would be explained, although additional research bacterial T4SS have been identified: 1) to facilitate transloca- needs to be conducted to confirm this possibility. tion of DNA via a conjugative mechanism to recipient cells; 2) In this study, eight of the plasmids carrying the VirB/D4 to facilitate translocation of effector molecules, such as pro- T4SS were completely sequenced including the two plasmids teins, to eukaryotic target cells; and 3) to function in DNA from isolate, N653 (table 1). The annotated sequence dem- uptake or dissemination from or into the environment milieu onstrated that the eight plasmids carry, in addition to the VirB/ (Christie et al. 2005; Alvarez-Martinez and Christie 2009). D4 T4SS systems, genes encoding other virulence-associated To further study the diversity among all 22 identified VB/D4 determinants and several proteins involved in a variety of met- plasmids, we identified SNPs that might be useful in differen- abolic processes, for example, hemolysin expression-modulat- tiating these plasmids. The SNP matrix from the 22 plasmids ing gene ymoA, DNA-binding protein genes (hns), (stpA), the and the reference plasmid built using the k-mer strategy of toxin/antitoxin stability genes (stbE), (stbD), IncN plasmid Kika kSNP consisted of 338 informative SNP positions. Figure 6 gene (kika), and DNA topoisomerase III (topB). shows a ML tree based on SNP analysis. The phylogeny A progressive Mauve alignment using the data from the placed the VB/D4 plasmids into three well-supported clades eight completely sequenced plasmids and the reference (95%, 100%, and 100% bootstrap support). Clade I con- S. Heidelberg plasmid pSARA30 (supplementary fig. S3, tained two plasmids isolated from the historical isolates, Supplementary Material online) separated the plasmids into SARA 37 and 82-2052, belonging to group 1 on the ML three distinct groups. As expected, the plasmids partitioned in tree shown in figure 3. Clade III is composed of ten plasmids, different clades in the plasmid VirB/D4 ML tree (fig. 6)were including the reference plasmid pSARA30 and 9 isolates be- also found by Mauve analysis to belong to different groups. longing to group 2 in figure 4. Clade II contains 11 plasmids, Group 1 included plasmids isolated from S. Heidelberg SARA including 7 plasmids from 7 isolates associated with the 2011 31, SARA 32, N4403, and N653, all having a size of approx- outbreak. Clade II can be further separated into two sub- imately ~38.0 kb and showing the highest degree of sequence clades, one of which contains nine plasmids derived from similarity to the S. Heidelberg plasmid pSARA30. Unique to 1058 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE 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. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1059 Hoffmann et al. GBE this group are several hypothetical proteins: truncated trans- mobilization plasmid, the VirB/D4 T4SS plasmid, the IncI1, and posase, DEAD/DEAH box helicase domain protein, and puta- the IncHI1. Except SARA 31, in which phage P4 and P2-like tive N-acetyltransferase. The second group is composed of the were not present, all other isolates had the same phage plasmids isolated from S. Heidelberg 41573 (~38 kb) and profile. 41563 (~33 kb)—both involved in the 2011 outbreak. The The results, which showed that S. Heidelberg isolates are best BLAST hit for them is Escherichia coli plasmid pDKX1- significantly (P< 0.0001) more pathogenic than the E. coli OP TEM-52 (GenBank accession number: JQ269336). Attributes 50 control strain, failed to demonstrate any appreciable dif- unique to this group are several hypothetical proteins, a ference in pathogenicity between strains carrying the VirB/D4 replicon initiation protein, and an ATPase central domain- T4SS plasmid and those strains that did not (fig. 7). However, containing protein. The third group contains plasmids isolated we also identified some VirB/D4 T4SS components on the from S. Heidelberg 41578, a 2011 outbreak clinical isolate, resistance plasmids IncA/C, IncHI2, and IncI1. Interestingly, and the chicken breast isolate N653. Both plasmids are ~35 kb isolate 29169, which did not carry any large plasmid and in size and members of the IncX4 subgroup. They share a high any of the T4SS components, was found to be significantly degree of sequence similarity to C. rodentium ICC168 plasmid (P< 0.0001) less pathogenic using this assay than any of the pCROD2 (FN543504). Unique to this group are some meta- other isolates evaluated. Isolate 41565, which carried the bolic genes, such as putative toxin–antitoxin system genes small mobility plasmid as well as additional IncHI2 plasmid, (hicA), (hicB), DnaJ-class molecular chaperone (dnaJ) (heat demonstrated the second weakest pathogenic potential. shock protein), plasmid replication gene (repA), site-specific Isolate SARA 31, which carried the VB/D4 T4SS plasmid but recombinase (xerD), and acetyl-CoA carboxylase beta subunit no other resistance plasmid, showed a survival curve profile (accD). With regards to the VirB/D4 gene clusters found on the similar to that of isolate N418, which did not have the VB/D4 IncX or IncX-like plasmids in this study, a considerable variation T4SS plasmid but carried the IncA/C plasmid. was observed, sometimes having <55% sequence identity The highest degree of pathogenicity was seen in N30678, between homologous genes (supplementary figs. S3 and which carries four plasmids including resistance and VB/D4 S4, Supplementary Material online), a phenomenon that has plasmids, and it is significantly (P< 0.05) different from all been observed elsewhere among the Enterobacteriaceae other isolates except the clinical outbreak isolate 41578, (Johnson et al. 2012). Comparison of the plasmids from the which also consist resistance and VB/D4 plasmids (fig. 7). three groups demonstrated that most of the T4SS-associated The data show that isolates not carrying any T4SS components genes divergence was seen in plasmids derived from group 3 tend to be significant less (P< 0.0001) pathogenic than those isolates. isolates that do carry the VB/D4 T4SS plasmid. Moreover, Noteworthy was the observation that three of these eight those strains that carry both the VB/D4 T4SS plasmid and an plasmids, each isolated from different 2011 clinical outbreak antibiotic resistance plasmid, which also carries the T4SS strains (41563 OR Blood, 41578 OH stool, 41573 OH Urine), genes, demonstrated a pathogenic potential greater than iso- differed in DNA sequence and gene carriage, suggestive of lates carrying only one of these plasmids. extensive horizontal gene transfer (HGT) or a mixed subpop- ulation of S. Heidelberg isolates. Caenorhabditis elegans System 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 po- tential of our strains. This system has been used successfully as an invertebrate host model to assess the virulence determi- nants 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 FIG.7.—Caenorhabditis elegans survival data from six S. Heidelberg small mobilization plasmid; strain SARA 31 carries only the isolates. The figure shows that the six S. Heidelberg isolates (29169, VirB/D4 T4SS plasmid; strain 418 carries a small mobilization SARA 31, 418, 41565, 41578, and N30678) are significantly plasmid and an IncA/C plasmid; strain 41565 carries a small (P< 0.0001) more pathogenic than the Escherichia coli OP 50 control mobilization plasmid and an IncHI2 plasmid; strain 41578 has strain. Further the figures show that isolates not carrying any T4SS com- a small mobilization plasmid, the VirB/D4 T4SS plasmid, and ponents tend to be significant less pathogenic than those isolates that do carry them. an IncI1 plasmid; and N30678 has four plasmids: the small 1060 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Identification of Antimicrobial Resistance Genes and Figure 8 shows a neighbor-joining tree generated from a Plasmids presence/absence matrix of resistance genes for the 30 group 3 Heidelberg isolates having a common PFGE profile. As ex- Analysis of WGS data showed 27 plasmid-associated antibi- pected, resistance tended to segregate according to the par- otic resistance genes among the 44 isolates and included ticular Inc group plasmid they carried and not by clonality or genes expected to produce resistance to aminoglycosides, evolutionary relationship. This observation indicated that the beta-lactams, phenicols, folate pathway inhibitors, tetracy- presence of specific resistance genes correlates with the par- clines, and cephems. These genes were not detected in the ticular Inc group of the plasmids. For example, the resistance eight pan susceptible isolates that also did not carry any of the genes (aphA1, bla , sul2, cmlA1, floR) were only detected CMY2 antibiotic resistance plasmids. Resistance phenotypes corre- in isolates that carry an IncA/C plasmid, whereas only isolates lated with genotypes in all strains examined (table 4). carrying an IncHI2 plasmid were found to contain the resis- However, some strains with the same genotype exhibited dif- tance genes aph, sph, tetC, tetD,and ble. Similarly, the aacC ferent antimicrobial resistance patterns. For example, both gene only appears to be present in those isolates identified as strains N18440 and N18413 contain the aadA1genebut carrying an IncI1 plasmid. The association of different resis- only isolate N18440 showed resistance to streptomycin. In tance genes with specific Inc groups is well known. this study, streptomycin resistance was defined using the MIC value >32 mg/l determined by CLSI. However, studies Mobilization Plasmids have shown that strains harboring an aadA gene cassette Plasmids below 8 kb did not carry any antimicrobial resistance can have MICs of 16 mg/l, which would classify these strains or virulence-associated genes. Instead, these smaller plasmids as streptomycin-susceptible (Sunde and Norstro ¨ m2005). tended to carry genes capable of promoting plasmid mobili- Among the 36 S. Heidelberg isolates resistant to one or zation, such as the mbeA, mbeB, mbeC, mbeD genes. more antibiotics, 24 contained sequences indicative of pres- S. Heidelberg isolates N189, N653, and N4403, which were ence of an IncI1 plasmid, 5 had an IncHI2 plasmid, and 1 originally isolated from chicken breast, carry a small, 3,772 bp carried an IncA/C plasmid. Three isolates contained two plas- plasmid that is identical to the S. Heidelberg SL476 plasmid mids, one with IncI1+ IncHI2, and two isolates each having pSL476_3 (GenBank accession number: CP001119) (Fricke IncI1+ IncA/C replicons (table 4). Besides isolate 41565, which et al. 2011). Allgroup3isolates (fig. 3), except isolate carries an IncHI2 plasmid, all Heidelberg isolates associated 24393, have similar XbaI patterns (JF6X01.0032, with 2011 outbreak carry an IncI1 plasmid, which was de- JF6X01.0033, JF6X01.0034, JF6X01.0058), and carry a tected the most overall among our strains. The IncI1 plasmid 4,473 bp plasmid that shares high sequence similarity with carries several antibiotic resistance genes and components of the S. Heidelberg plasmid pSH1148_4.8 (GenBank accession the T4SS. The IncHI2 plasmid was found in eight isolates and number: JX494965). SARA 35, which was isolated in Brazil, carries several antibiotic resistance genes and genes conferring was found to carry a 6,647 bp ColE1 plasmid resistance to quaternary ammonium compounds, as well as (pCFSAN000443_01) that is similar to E. coli plasmid ColE1 genes encoding the conjugative T4SS. Furthermore, unique to (GenBank accession number: J01566). the IncHI2 plasmid is the presence of genes encoding resis- tance to heavy metals, including tellurium, copper, cadmium, Discussion mercury, and silver. The three isolates resistant to the most (ten or more) anti- Within the last several years, S. Heidelberg has been identified biotics carried the IncA/C multidrug resistance plasmid. In ad- as one of the top serovars responsible for human illness, in- dition to carrying the genes encoding resistance to different cluding a recent multistate outbreak of an antibiotic-resistant antibiotics and quaternary ammonium compounds, this plas- strain associated with consumption of contaminated ground mid also carries the conjugative T4SS transfer system. All three turkey. Using whole-genome sequences, we characterized the IncA/C plasmids carried the resistance element bla -blc- CMY-2 genetic diversity of 46 S. Heidelberg isolated over ~30 years. sugE2, which was only found in one other isolate, N189, on This study illustrates the novel ways in which WGS, combined an IncI1 plasmid. with phylogenetic analysis, can be used to investigate and SARA 35, which was only resistant to ampicillin, did not characterize the diversity of S. Heidelberg. Using this strategy carry a multidrug resistance plasmid but, instead, carried a with outbreak isolates, we were also able to comprehensively smaller plasmid (~8.0 kb) having only a single determinant characterize and differentiate among highly clonal that correlated with resistance to b-lactams (tnpA-bla ). S. Heidelberg isolates. For example, SNP analysis in combina- TEM-1 This plasmid shared a high degree of sequence similarity to tion with phylogenetic analysis is able to determine the phy- S. Typhimurium pAnkS (accession number: DQ916413), logenetic relationships among S. Heidelberg isolates, easily which also carried the transposon resistance element (tnpA- separating them from other S. Heidelberg strains, including bla ). The Inc group for this plasmid could not be deter- isolates having the same XbaIand BlnI PFGE patterns and TEM-1 mined from the sequence. thought to be highly clonal. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1061 Hoffmann et al. GBE 1062 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Table 4 Antibiotic Resistance and Plasmid Characteristics of S. Heidelberg Isolates Used in This Study Isolates Antimicrobial Resistance Profile MDR Plasmids Resistance Genes on Plasmids SARA37 STR, TET IncI1 aph, sph, strA, strB, tetC, ble SARA35 AMP ca 8 kb bla TEM1 SARA32 Susceptible None None SARA31 Susceptible None None N418 AMC,AMP,FOX,AXO,CHL,GEN,KAN,STR,SUL,TET,TIO IncA/C aac3-VI, aadA1, aadA2, aadB, aphA1, strA, strB, bla ,bla ,sul1, sul2,cmlA1, CMY2 TEM1 floR, tetA, tetR(A), qacDE, sugE2 N19992 Susceptible None None SARA39 SUL IncHI2 aadA1, sul1, tetG, tetR(G), qacDE N189 AMC,AMP,FOX,AXO,TIO IncI1 bla , sugE2 CMY2 N1514 Susceptible None None N1536 AMP, GEN,STR,SUL IncI1 aac3-VI, aadA1, bla , sul1, qacDE TEM1 N4496 Susceptible None None N4403 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE N4541 Susceptible None None N653 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 82-2052 KAN, STR, TET IncI1 aph, sph, strA, strB, tetC, ble 24390 GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 21381 AMC,AMP,FOX,AXO,GEN,KAN,STR,SUL,TET,TIO IncI1; IncA/C aac3-VI, aadA1, aphA1, strA, strB, bla ,sul1, sul2,tetA, tetR(A), qacDE, sugE2 CMY2 N4630 GEN,STR,SUL IncI1 aac3-VI, aadA1,sul1, qacDE N19871 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N15757 AMP,KAN,STR,TET IncI1; IncHI2 aph, sph, strA, strB, bla , tetC, tetD, ble TEM1 N26457 AMC,AMP,FOX,AXO,GEN,KAN,STR,SUL,TET,TIO IncI1; IncA/C aac3-VI, aadA1, aphA1, strA, strB, bla ,sul1, sul2,tetA, tetR(A), qacDE, sugE2 CMY2 N18393 AMP,GEN,TET IncI1 aacC, aadA1, blaTEM1, tetA, tetR(A) N29341 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41578 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 41563 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41567 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41573 AMP,GEN,TET IncI1 aacC, aadA1,bla , tetA, tetR(A) TEM1 41576 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41584 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB, bla , tetA, tetR(A) TEM1 41579 AMP,GEN,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 41566 AMP,GEN,STR,TET IncI1 aacC, aadA1, strA, strB,bla , tetA, tetR(A) TEM1 41565 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N30678 AMP,GEN,KAN,STR,SUL,TET IncI1; IncHI2 aac3-VI, aadA1, aph, sph, strA, strB, bla , sul1, tetC, tetD, ble, qacDE TEM1 32507 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble 24393 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N20134 AMP,GEN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 29169 Susceptible None None 24391 KAN,STR,SUL IncI1 aac3-VI, aadA1, sul1, qacDE 24359 GEN,KAN,STR,SUL,TET IncI1; IncHI2 aac3-VI, aadA1, aph, sph, strA, strB, sul1, tetC, tetD, ble, qacDE 24388 Susceptible None None 20752 KAN,STR,TET IncHI2 aph, sph, strA, strB, tetC, tetD, ble N18413 AMP,GEN,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N18440 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 N18453 AMP,GEN,STR,TET IncI1 aacC,aadA1,bla , tetA, tetR(A) TEM1 Genomic Analysis and Virulence Differences GBE 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 plasmidare mapped to the tree. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1063 Hoffmann et al. GBE We also found that the SNP analysis clustered isolates to carry cassette genes (morons) encoding factors that en- together that did not have the same PFGE profiles, such as hance the proliferation and dissemination of the prophage the 2011 outbreak isolates, which were distinct on the PFGE by improving the fitness and /or virulence of the lysogen two-enzyme dendrogram. SNP analysis tightly clustered (Boyd and Brussow 2002; Pelludat et al. 2003). This possibility the clinical and ground turkey outbreak isolates together has led to the assumption that phage-mediated distribution of and distinguished these from other apparently clonal iso- virulence factors and fitness factors is a key driving force in the lates that shared the same PFGE and antibiotic resistance pro- optimization of Salmonella–host interactions and the emer- files. This result confirms the epidemiological data that the gence of new epidemic clones (Brussow et al. 2004). For ex- clinical and ground turkey isolates belong to the same out- ample, lambda-like phage Gifsy-2, which is present among all break, and are separate from other strains grouped with them S. Heidelberg isolates examined in the study, carries the viru- by PFGE. lence factor periplasmic CuZn-superoxide dismutase (SodCI), We found two nonsynonymous SNPs that were present an enzyme that has the potential to enhance the fitness of a only among the outbreak isolates. The first of these SNPs Gifsy-2 lysogen (Brussow et al. 2004). Moreover, the pro- was found to be located on the gluconate transporter gntP phage P4 and prophage P2-like were associated with isolates gene. The biological role of GntP might be to permit the host from group 3 which included the outbreak isolates. The P2- to take up gluconate when only tiny amounts are present, like phage carries the sopE1 gene that encodes part of the thereby conferring an advantage to the host in a competitive T3SS that translocates bacterial effector proteins directly into environment (Klemm et al. 1996). The second of these SNPs the cytosol of host cells. The phage P2-like has also been was located within a gene encoding a putative hydrolase identified in S. Typhimurium strain DT49/DT204 that caused (phosphoserine phosphatase). This gene is involved in glycine a major outbreaks in Europe in 1970s and 1980s. It has been and serine metabolism, and resulted in a proline to glutamine speculated that lysogenic conversion with sopE1 was an im- amino acid change that is known to change the secondary portant step in the emergence of this epidemic strain (Mirold structure of this enzyme. et al. 2001). Interestingly, Mirold et al. reported that sopE1 is SNP-based phylogeny clearly has a distinct advantage in not restricted to a certain bacteriophage as a “vehicle,” iden- subtyping clonal isolates because the method offers nucleo- tifying sopE1 in other serovars present on a Gifsy-like phage tide base resolution. We detected 4,053 SNPs within highly (Mirold et al. 2001). We also saw the same phenomenon clonal S. Heidelberg isolates. Compared with the PFGE analy- among the S. Heidelberg of this study where the sopE1 ses that grouped isolates together having the same Xbaland gene, which was found to be present among all isolates, BlnI patterns, the resolution of the SNP phylogeny conclusively was carried by both the P2-like and the lambda-like Gifsy-2 demonstrated that genome sequencing distinguishes be- phages. tween outbreak and nonoutbreak isolates that shared the The exchanges of genetic cassettes between unrelated same XbaIand BlnI patterns. The results from the ML tree phages further increase the degree of HGT noted among generated from the SNP matrix constructed with kSNP are S. Heidelberg strains, clearly demonstrating the importance concordant with those observed based on the SNP matrix of phages in the emergence and evolution of Salmonella path- developed with our custom pipeline. The study showed that ogens. To date, this is the first time that phages have been SNP analysis is very consistent and reliable, confirming that characterized by name among S. Heidelberg isolates and that even when different methods were used to identify SNPs, a P2-like phage was identified in S. Heidelberg that carries the there are no appreciable changes in the overall phylogeny of sopE1 gene. This phage and, subsequently, the sopE1 deter- the strains included in the analyses. minant, was present in the outbreak isolates, possibly serving Only a few mutation “hotspot” genes were found to have as one of the virulence factors giving rise to these pathogens. a disproportional number of SNPs and that might serve as Furthermore, we identified the presence of 74 virulence de- targets useful for rapidly subtyping serovar Heidelberg. terminants in all S. Heidelberg isolates. Although these genes have several SNPs, the haplotype diver- In the course of performing this study, both newly se- sity is fairly low, indicating there are only few unique quenced plasmids and plasmids sharing a high sequence sim- sequences among S. Heidelberg. This information supports ilarity to previously described plasmids (Han et al. 2012)were the value of SNP analysis as a useful tool for subtyping this identified. Similar to the phage presence and absence tree serovar. topology previously described, the mobilization plasmids ab- WGS data and comparative genomics not only offer a sence/presence seemed to correlate with the tree topology useful subtyping method to differentiate closely related bac- according to the plasmid(s) they carried. For instance, the teria, but they also provide a better understanding of patho- 3,772 bp mobilization plasmid was only identified in isolates genicity and evolution. Our results show substantial loss and belonging to paraphyletic group 2, whereas the 4,473 bp plas- gain of plasmids and phages between these S. Heidelberg mid was found only in isolates which belong to the paraphy- isolates that were consistent with previously generated tree letic group 3 on the ML tree. The mobilization plasmids carry topologies. The phages found among these strains are known the genes necessary to encode the proteins involved with 1064 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE relaxosome formation and processing. Moreover, these plas- stability to assure that the plasmid as well as its associated mids are efficiently transferred by other plasmids, such as IncI1 virulence determinants are inherited by all bacterial progeny and IncX (Francia et al. 2004), two plasmids frequently found (Yamaguchi and Inouye 2011; Unterholzner et al. 2013). As among the S. Heidelberg in this and other studies (Folsteretal. such, these genes serve to stabilize the “pathogenic poten- 2012; Han et al. 2012). The mobilization plasmids tend to be tial” of these isolates. retained and are, continually carried by specific groups of Finding that the T4SS system was present in the majority of S. Heidelberg isolates, suggesting that they must play an es- the pathogenic isolates and carried on plasmids also bearing sential role in this serovar survival, although the precise role is determinants associated with assuring vertical transmission of unclear. Plasmid sequence analysis detected additional plas- the plasmid to progeny bacteria, we sought to assess the role mids that encode several S. Heidelberg virulence determi- of the T4SS in overall pathogenicity. Using the C. elegans vir- nants, including T4SS and antibiotic resistance genes. Unlike ulence assay, we confirmed that S. Heidelberg isolates that the smaller plasmids, the loss and gain of the larger plasmids carried components of the T4SS were more pathogenic in this (ranging from 33 to >200 kb) did not correlate with the tree model than isolates that did not carry the system, suggesting topology. The plasmid analyses performed in this and previous the T4SS plays a role in the pathogenicity of S. Heidelberg. studies (Folsteretal. 2012; Han et al. 2012) demonstrated that Much more inclusive investigations must be performed to fully S. Heidelberg isolates have at least four transmissible plasmids, understand the role of plasmid-encoded genes in pathogen- including an IncA/C, IncI, and IncHI2 and IncX plasmid. The esis, with special emphasis on the role of virulence-associated IncA/C, IncI, and IncHI2 plasmids carry genes that are impor- VirB/D4 T4SS secretory mechanism. tant for virulence, colonization, and persistence, as well as In summary, the data presented in this study support the those encoding resistance determinants to heavy metals, dis- role of transmissible mobile genetic elements in the evolution infection, and antimicrobial agents. We determined that most of virulence and resistance among S. Heidelberg. Understand- of the antimicrobial-resistance phenotypes could be ac- ing the results and the role of these elements is very important counted for by resistance genes carried on plasmids. We did since S. Heidelberg is the second most common serovar, fol- not identify an instance of resistance occurring in the absence lowing S. Typhimurium, responsible for Salmonella-related of a corresponding gene. This shows that genetic prediction of deaths in the United States (Kennedy et al. 2004; Patchanee phenotypic resistance may be possible and may become et al. 2008). Although the mechanisms responsible for necessary as clinical diagnostics moves toward culture- S. Typhimurium pathogenicity and colonization are well stud- independent technologies. More studies are necessary to de- ied, fewer studies are available for S. Heidelberg (Han et al. termine how well specific DNA sequences, and combinations 2011). Presumably, the transmissible plasmids identified in our of sequences, predict resistance or susceptibility to various strain set play a major role in the spread of virulence and antimicrobial agents. resistance genes among S. Heidelberg and, most likely, Besides the IncI1 and the IncHI2 antibiotic resistance plas- other Salmonella serovars. mids, an IncX plasmid that encodes the structural components The results from this and other studies suggest that plas- of the VirB/D4 T4SS apparatus was detected in most of the mids are being exchanged quickly. For example, it was noted outbreak isolates. The VirB/D4 T4SS apparatus is important in that some strains associated with the 2011 outbreak isolated that it assists with “rapid dissemination of antibiotic resistance from the same state had different plasmid profiles. In addition and virulence determinants” (Bhatty et al. 2013). Unlike other to the ecological benefits conferred as a consequence of secretion systems, the T4SS is not only capable of acting in a this rapid plasmid exchange, large plasmids can change conjugative fashion disseminating DNA to other bacteria, PFGE patterns. We noted that, while most of the se- thereby contributing to genome plasticity and the evolution quenced isolates associated with the 2011 outbreak had of infectious pathogens through dissemination of antibiotic BlnI PFGE pattern JF6A26.0076, a single isolate, 41565, had resistance and virulence genes, but also capable of transfer- a different, JF6A26.0017 BlnI pattern, most likely due to ring protein effectors across both the bacterial and eukaryotic the presence of an IncHI2 plasmid rather than IncI1. Other plasma membrane directly into the cytosol of the target cells, closely related isolates not associated with the outbreak car- thereby contributing directly to bacterial pathogenicity (Juhas rying the IncHI2 plasmid also exhibited BlnIPFGE pattern et al. 2008; Gokulan et al. 2013). The T4SS was first described JF6A26.0017. for Agrobacterium tumefaciens and consists of 11 proteins The results determined in this study have convincingly encoded by the virB gene complex. Of particular interest shown that WGS analysis could be a useful tool for identifying was the fact that some plasmids carrying the T4SS, such as the source of contamination and studying the short-term evo- those of the IncX group, also carry systems that confer plasmid lution of these epidemic clones. In the context of the 2011 stability among the strains possessing the plasmid. For in- ground turkey-associated outbreak, SNP analysis provided an stance, the IncX plasmids carry the toxin–antitoxin system excellent tool for subtyping S. Heidelberg and clustering clo- genes (hicA), (hicB), as well as the toxin/antitoxin stability sely related isolates together. Using these data, we deter- genes (stbE), (stbD), both of which contribute to the plasmid mined that the antibiotic-resistant outbreak isolates diverged Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1065 Hoffmann et al. GBE from environmental isolates through the gain of several viru- Disclosure forms provided by the authors will be available lence factors, such as sopE1 gene carried on Phage P2 or the with the full text of this article. This work was supported by T4SS system carried on the IncX plasmid. Of particular interest, the Center for Veterinary Medicine and the Center for Food other isolates in paraphyletic group 3 on the ML tree carried Safety and Applied Nutrition at the US Food and Drug the same virulence factors as the 2011 outbreak strain, sug- Administration and by an appointment of MH by the Joint gesting that this strain has been seen earlier in the environ- Institute for Food Safety and Applied Nutrition (JIFSAN), ment and that the outbreak occurred as a consequence of University of Maryland, College Park, MD. changes in outside influences such as the dose or the handling of the food. Literature Cited CDC projected that for every case of Salmonella reported Aballay A, Yorgey P, Ausubel FM. 2000. Salmonella typhimurium prolifer- that 29.3 cases occur in the community and 28.3 cases goes ates and establishes a persistent infection in the intestine of undiagnosed or unreported (Scallan et al. 2011), suggesting Caenorhabditis elegans. Curr Biol. 10:1539–1542. that the 2011 outbreak actually may have sickened over Allard MW, et al. 2012. High resolution clustering of Salmonella enterica 4,000 people in the United States. The World Health serovar Montevideo strains using a next-generation sequencing ap- Organization has declared that most countries have docu- proach. BMC Genomics 13:32. Allard MW, et al. 2013. On the evolutionary history, population genetics mented a significant increase of microorganism-mediated and diversity among isolates of Salmonella enteritidis PFGE pattern foodborne illnesses (WHO 2013). Therefore, it is paramount JEGX01.0004. PLoS One 8:e55254. to have a robust identification system that can distinguish Alvarez-Martinez CE, Christie PJ. 2009. Biological diversity of prokaryotic closely related strains in an outbreak situation. To date, type IV secretion systems. Microbiol Mol Biol Rev. 73:775–808. DNA “fingerprinting” facilitated by PFGE at PulseNet has Aziz RK, et al. 2008. The RAST Server: rapid annotations using subsystems technology. BMC Genomics 9:75. often been used to establish the source of clinical isolates. Bazinet AL, Cummings MP. 2011. Computing the tree of life—leveraging However, several examples in the past have shown that PFGE the power of desktop and service grids. In: Proceedings of the Fifth may not be as useful for some serotypes that show less di- Workshop on Desktop Grids and Volunteer Computing Systems versitybyPFGE(CDC 2010a, 2010b). Therefore, it is essential PCGrid 2011. Available from: http://lattice.umiacs.umd.edu/lattice- to develop and evaluate other methods for their ability to files/publications/lattice/computing_tree_of_life.pdf. Bazinet AL, Myers DS, Fuetsch J, Cummings MP. 2007. Grid services base effectively and accurately discriminate outbreak isolates library: a high-level, procedural application program interface for writ- from nonoutbreak isolates, like the combination of WGS ing Globus-based Grid services. Future Gener Comput Syst. 23: and computationally efficient methods of analysis. Such a 517–522. combination of approaches also will help to genetically link Beltran P, et al. 1991. Reference collection of strains of the Salmonella the implicated food sources of contamination with farm or typhimurium complex from natural populations. J Gen Microbiol. 137: factories. As more WGS data are gathered and analytical 601–606. Bhatty M, Laverde Gomez JA, Christie PJ. 2013. The expanding bacterial tools developed in combination with epidemiology investiga- type IV secretion lexicon. Res Microbiol. 164(6):620–639. tions, we anticipate significant improvements in timeliness of Boyd EF, Brussow H. 2002. Common themes among bacteriophage- outbreak investigations and attribution of human infections encoded virulence factors and diversity among the bacteriophages to specific sources. involved. Trends Microbiol. 10:521–529. Brussow H, Canchaya C, Hardt WD. 2004. Phages and the evolu- tion of bacterial pathogens: from genomic rearrangements to lyso- Supplementary Material genic conversion. Microbiol Mol Biol Rev. 68:560–602, table of contents. Supplementary tables S1 and S2 and figures S1–S4 are avail- CDC. 2010a. Centers for Disease Control and Prevention. 2010. Multistate able at Genome Biology and Evolution online (http://www. outbreak of human Salmonella Montevideo infections (final ubdate). gbe.oxfordjournals.org/). Atlanta (GA): Centers for Disease Control and Prevention. Available from: http://www.cdc.gov/salmonella/montevideo. CDC. 2010b. Centers for Disease Control and Prevention. 2011. Multistate Acknowledgments outbreak of human Salmonella Enteritidis infections associated with shell eggs (final ubdate). Atlanta (GA): Centers for Disease The authors thank Ruth Timme for submission of the draft Control and Prevention. Available from: http://www.cdc.gov/salmo genomes to NCBI and the NCBI rapid annotation pipeline nella/enteritidis. team, Bill Klimke, Dmitry Dernovoy, Stacy Ciufo, Kathleen CDC. 2011. Centers for Disease Control and Prevention. 2011. Multistate O’Neill, Azat Badretdin, and Tatiana Tatusova, for key outbreak of human Salmonella Heidelberg infections linked to ground turkey. Atlanta (GA): Centers for Disease Control and Prevention. genome annotation service. Further they thank Yuansha Available from: http://www.cdc.gov/salmonella/heidelberg/. Chen for sharing the CVM/DAFM database of antimicrobial CDC. 2013. Centers for Disease Control and Prevention. 2013. resistance genes. No human subjects or animals were used in National enteric disease surveillance: Salmonella annual report, this study. All authors have read the manuscript and agreed to 2011. Atlanta (GA): Centers for Disease Control and Prevention. its contents, subject matter, and author line order. These data Available from: http://www.cdc.gov/ncezid/dfwed/PDFs/salmonella- are novel and have not been previously published elsewhere. annual-report-2011-508c.pdf. 1066 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 Genomic Analysis and Virulence Differences GBE Chittick P, Sulka A, Tauxe RV, Fry AM. 2006. A summary of national Juhas M, Crook DW, Hood DW. 2008. Type IV secretion systems: tools of reports of foodborne outbreaks of Salmonella Heidelberg infections bacterial horizontal gene transfer and virulence. Cell Microbiol. 10: in the United States: clues for disease prevention. J Food Prot. 69: 2377–2386. 1150–1153. Kennedy M, et al. 2004. Hospitalizations and deaths due to Salmonella Christie PJ, Atmakuri K, Krishnamoorthy V, Jakubowski S, Cascales E. infections, FoodNet, 1996–1999. Clin Infect Dis. 38(Suppl. 3): 2005. Biogenesis, architecture, and function of bacterial type IV secre- S142–S148. tion systems. Annu Rev Microbiol. 59:451–485. Klemm P, Tong S, Nielsen H, Conway T. 1996. The gntP gene of Crump JA, et al. 2011. Antimicrobial resistance among invasive nontyphoi- Escherichia coli involved in gluconate uptake. J Bacteriol. 178:61–67. dal Salmonella enterica isolates in the United States: National Klimke W, et al. 2009. The National Center for Biotechnology Antimicrobial Resistance Monitoring System, 1996 to 2007. Information’s Protein Clusters Database. Nucleic Acids Res. 37: Antimicrob Agents Chemother. 55:1148–1154. D216–D223. Darling AE, Mau B, Perna NT. 2010. progressiveMauve: multiple genome Kurtz S, et al. 2004. Versatile and open software for comparing large alignment with gene gain, loss and rearrangement. PLoS One 5: genomes. Genome Biol. 5:R12. e11147. Liu F, et al. 2011. Novel virulence gene and clustered regularly interspaced Edgar RC. 2010. Search and clustering orders of magnitude faster than short palindromic repeat (CRISPR) multilocus sequence typing scheme BLAST. Bioinformatics 26:2460–2461. forsubtypingofthe majorserovarsof Salmonella enterica subsp. Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters of enterica. Appl Environ Microbiol. 77:1946–1956. individuals using the software STRUCTURE: a simulation study. Mol Marcais G, Kingsford C. 2011. A fast, lock-free approach for efficient Ecol. 14:2611–2620. parallel counting of occurrences of k-mers. Bioinformatics 27: Falush D, Stephens M, Pritchard JK. 2003. Inference of population struc- 764–770. ture using multilocus genotype data: linked loci and correlated allele Mirold S, et al. 2001. Salmonella host cell invasion emerged by acquisition frequencies. Genetics 164:1567–1587. of a mosaic of separate genetic elements, including Salmonella Foley SL, et al. 2011. Population dynamics of Salmonella enterica serotypes pathogenicity island 1 (SPI1), SPI5, and sopE2. J Bacteriol. 183: in commercial egg and poultry production. Appl Environ Microbiol. 77: 2348–2358. 4273–4279. NARMS. 2010. National Antimicrobial Resistance Monitoring System- Folster JP, et al. 2012. Characterization of multidrug-resistant Salmonella Enteric Bacteria (NARMS). Bethesda (MD): Executive Report FDA. enterica serovar Heidelberg from a ground turkey-associated outbreak Norman A, Hansen LH, She Q, Sorensen SJ. 2008. Nucleotide sequence of in the United States in 2011. Antimicrob Agents Chemother. 56: pOLA52: a conjugative IncX1 plasmid from Escherichia coli which en- 3465–3466. ables biofilm formation and multidrug efflux. Plasmid 60:59–74. Francia MV, et al. 2004. A classification scheme for mobilization regions of Paradis E, Claude J, Strimmer K. 2004. APE: analyses of phylogenetics and bacterial plasmids. FEMS Microbiol Rev. 28:79–100. evolution in R language. Bioinformatics 20:289–290. Fricke WF, et al. 2011. Comparative genomics of 28 Salmonella enterica Patchanee P, Zewde BM, Tadesse DA, Hoet A, Gebreyes WA. 2008. isolates: evidence for CRISPR-mediated adaptive sublineage evolution. Characterization of multidrug-resistant Salmonella enterica serovar J Bacteriol. 193:3556–3568. Heidelberg isolated from humans and animals. Foodborne Pathog Gardner SN, Slezak T. 2010. Scalable SNP analyses of 100+ bacterial or Dis. 5:839–851. viral genomes. J Forensic Res. 01:107. Pelludat C, Mirold S, Hardt WD. 2003. The SopEPhi phage integrates Gokulan K, et al. 2013. Impact of plasmids, including those into the ssrA gene of Salmonella enterica serovar Typhimurium A36 encodingVirB4/D4 type IV secretion systems, on Salmonella enterica and is closely related to the Fels-2 prophage. J Bacteriol. 185: serovar Heidelberg virulence in macrophages and epithelial cells. PLoS 5182–5191. One 8:e77866. Pritchard JK, Stephens M, Donnelly P. 2000. Inference of population struc- Gunn JS, Alpuche-Aranda CM, Loomis WP, Belden WJ, Miller SI. 1995. ture using multilocus genotype data. Genetics 155:945–959. Characterization of the Salmonella typhimurium pagC/pagD chromo- Rosenberg NA. 2004. Distruct: a program for the graphical display of somal region. J Bacteriol. 177:5040–5047. population structure. Mol Ecol. 4:137–138. Han J, et al. 2011. Comparison of Salmonella enterica serovar Heidelberg Sahu SN, et al. 2012. Genomic analysis of immune response against Vibrio isolates from human patients with those from animal and food cholerae hemolysin in Caenorhabditis elegans. PLoS One 7:e38200. sources. J Clin Microbiol. 49:1130–1133. Scallan E, et al. 2011. Foodborne illness acquired in the United States— Han J, et al. 2012. DNA sequence analysis of plasmids from multidrug major pathogens. Emerg Infect Dis. 17:7–15. resistant Salmonella enterica serotype Heidelberg isolates. PLoS One 7: Sprando RL, Olejnik N, Cinar HN, Ferguson M. 2009. A method to e51160. rank order water soluble compounds according to their toxicity Hoffmann M, et al. 2012. Genome sequences of five Salmonella enterica using Caenorhabditis elegans, a complex object parametric ana- serovar Heidelberg isolates associated with a 2011 multistate outbreak lyzer and sorter, and axenic liquid media. Food Chem Toxicol. 47: in the United States. J Bacteriol. 194:3274–3275. 722–728. Huehn S, Bunge C, Junker E, Helmuth R, Malorny B. 2009. Poultry-asso- Sunde M, Nortstro ¨ m M. 2005. The genetic backround for streptomycin ciated Salmonella enterica subsp. enterica serovar 4,12:d:- reveals high resistance in Escherichia coli infuences the distribution of MICs. clonality and a distinct pathogenicity gene repertoire. Appl Environ Antimicrob Agents Chemother. 56:87–90. Microbiol. 75:1011–1020. Tamura K, Dudley J, Nei M, Kumar S. 2007. MEGA4: molecular evolution- Huehn S, et al. 2010. Virulotyping and antimicrobial resistance typing of ary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 24: Salmonella enterica serovars relevant to human health in Europe. 1596–1599. Foodborne Pathog Dis. 7:523–535. Timme RE, et al. 2013. Phylogenetic diversity of enteric pathogen Johnson TJ, et al. 2012. Expansion of the IncX plasmid family for improved Salmonella enterica subsp. enterica inferred from genome-wide refer- identification and typing of novel plasmids in drug-resistant ence-free SNP characters. Genome Biol Evol. 5(11):2109–2123. Enterobacteriaceae. Plasmid 68:43–50. Unterholzner SJ, Hailer B, Poppenberger B, Rozhon W. 2013. Johnson TJ, Nolan LK. 2009. Plasmid replicon typing. Methods Mol Biol. Characterisation of the stbD/E toxin-antitoxin system of pEP36, a 551:27–35. Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014 1067 Hoffmann et al. GBE plasmid of the plant pathogen Erwinia pyrifoliae. Plasmid 70: Zhou Y, Liang Y, Lynch KH, Dennis JJ, Wishart DS. 2011. PHAST: a fast 216–225. phage search tool. Nucleic Acids Res. 39:W347–W352. WHO. 2013. Food safety: general information related to microbiological Zwickl D. 2006a. Genetic algorithm approaches for the phylogenetic anal- risks in food. Available from: http://www.who.int/foodsafety/en/. ysis of large biological sequence dataset under the maximum likeli- Wilmshurst P, Sutcliffe H. 1995. Splenic abscess due to Salmonella hood criterion. Austin (TX): University of Texas. Heidelberg. Clin Infect Dis. 21:1065. Zwickl DJ. 2006b. Genetic algorithm approaches for the phylogenetic Yamaguchi Y, Inouye M. 2011. Regulation of growth and death in Escher- analysis of large biological sequence datasets under the maximum ichia coli by toxin-antitoxin systems. Nat Rev Microbiol. 9:779–790. likelihood criterion [PhD dissertation]. [Austin (TX)]: The University of Zhao S, et al. 2008. Antimicrobial resistance in Salmonella enterica serovar Texas. Heidelberg isolates from retail meats, including poultry, from 2002 to Associate editor: Takashi Gojobori 2006. Appl Environ Microbiol. 74:6656–6662. 1068 Genome Biol. Evol. 6(5):1046–1068. doi:10.1093/gbe/evu079 Advance Access publication April 14, 2014

Journal

Genome Biology and Evolution – Oxford University Press

Published: May 1, 2014

There are no references for this article.