High-throughput genome sequencing of two Listeria monocytogenes clinical isolates during a large foodborne outbreak

Primary outbreak isolate harbors unique plasmid pLM5578 and prophage ϕLMC1

Public health responses to foodborne outbreaks require bacterial subtyping data to definitively link affected individuals with contaminated food sources and other suspected patients. The genetic relatedness between clinical and food isolates is routinely determined using methods such as PFGE, however, with the latest generation of DNA sequencing technologies it is theoretically possible to now determine the complete genetic blueprint of bacterial isolates within the scope of an outbreak. In this study, draft genome sequences suitable for comparative analysis between outbreak isolates and to reference genomes were available within three days. These data allowed for discovery of novel prophage and genomic islands, and the subsequent sequence closure and the confirmation of sequence variants allowed for annotation of polymorphisms and other traits associated with micro-diversity.

Clinical isolate 08-5578 was selected as the reference outbreak strain (serotype 1/2a; PFGE patterns LMACI.0040 and LMAAI.0001; Table 1) and clinical isolate 08-5923 differed only in the Asc I restriction pattern (LMACI.0001; Fig. 1). High-throughput pyrosequencing was performed on both isolates and closed circular chromosomal sequences were obtained for each isolate (Table 2; Fig. 2). A 77 kbp contig present in the 08-5578 sequence assembly but absent in 08-5923 was identified as a circular plasmid and designated pLM5578 (Table 2). Comparison of the chromosomal sequences for isolates 08-5578 and 08-5923 with publicly available Listeria spp. genomes (Fig. 3A) and with select loci used for MLST (Fig. 3B) indicated that these isolates were members of evolutionary lineage II and clonal complex 8 (discussed further below).

The Asc I restriction patterns of 08-5578 and 08-5923 were indistinguishable, except for a single band shift of approximately 32 kbp (from 476 kbp in 08-5578 to 444 kbp in 08-5923; Fig. 1). Accordingly, comparative sequence alignment revealed a 33 kbp contiguous region integrated adjacent to (but not disrupting) the tRNA-Ser gene that was unique to 08-5578 and encoding several putative bacteriophage-related coding sequences (designated ϕLMC1) (Fig. 4). BLAST analyses indicated that this prophage was unique among sequenced Listeria genomes, although it was comprised of coding sequences similar to determinants from several known phages associated with Gram-positive microbes, including L. monocytogenes phages A006, A118 and B025 (see Additional file 1). The presence of prophage ϕLMC1 accounted for the Asc I restriction pattern difference between isolates 08-5578 and 08-5923.

The extrachromosomal circular plasmid (pLM5578) harbored in isolate 08-5578 encoded putative determinants for replication, partitioning, heavy metal resistance and DNA translocation (Fig. 5a, and Additional file 2). Many of these features were homologous to sequences encoded by plasmids from L. innocua (pLI100), L. monocytogenes serotype 4b isolate H7858 (pLM80) and Bacillus anthracis (pXO2), including a replication and regulation control center (pLM5578_17 - 20). Two regions of pLM5578 were related to heavy metal resistance, including a transposon-associated cadA that exhibited high sequence identity to pLI100. A second, distinct region unique to pLM5578 was comprised of cadA (second copy) and cadC genes, with the latter encoding a cadmium-efflux ArsR-family regulatory accessory protein (repressor) homolog. Resistance to the heavy metal cadmium has been recently discovered to be associated with resistance to quaternary ammonium compounds (sanitizers), which may contribute to L. monocytogenes persistence in food production facilities [35]. Notably, several coding sequences exhibited sequence similarity to type IV secretion systems (T4SS), including DNA processing (virD2), DNA coupling (virD4) and membrane channel formation (virB1, -B4, -B5, -B6, -B11) (see Additional file 2). Additional predicted membrane-associated proteins were present within this region, and all coding sequences were oriented in the same direction (Fig. 5A), which is consistent with the operon-style organization of most known T4SS. Varying degrees of similarity between pLI100, pLM80 and pXO2 plasmids have already been established [21], and our data further support a chimeric nature for Listeria plasmids.

pLM5578 and ϕLMC1 exhibit non-homogeneous distribution within a wider panel of outbreak-associated isolates

During the outbreak, serotype 1/2a L. monocytogenes isolates were recovered from clinical samples, contaminated food products and the food processing environment that shared the two Asc I PFGE patterns observed for the sequenced isolates (Table 1). To determine the distribution of ϕLMC1 and pLM5578 amongst outbreak-associated isolates, PCR-based screening was conducted by targeting ϕLMC1 loci encoding a putative phage terminase large subunit and phage tail tape measure protein, and pLM5578 loci encoding virD4 and fic (Table 1). All isolates with Asc I pattern LMACI.0001 lacked the ϕLMC1 determinants by PCR, confirming that this novel sequence insertion accounts for the restriction enzyme pattern shift relative to PFGE pattern LMACI.0040. Bacteriophage insertions have previously been recognized as contributing to macromolecular genetic diversity between related outbreak-associated L. monocytogenes strains [12, 36].

The PCR screening assays for plasmid determinants confirmed that pLM5578 was absent in isolate 08-5923, and this plasmid was also absent in other outbreak-related isolates (independent of PFGE profile), including one clinical isolate (Table 1). The sequence of pLM5578 revealed an absence of Asc I restriction sites and presence of a single Apa I restriction site, but this element did not contribute to the PFGE patterns resulting from these enzymes. In contrast, plasmid pLM80 (containing a single Apa I recognition site) had been reported as visible in Apa I PFGE digests [36]. Just as a variable presence of plasmid pLM5578 was observed in this Canadian outbreak, plasmid pLM80 was sporadically present in isolates during the 1998-99 multi-state serotype 4b frankfurter USA outbreak [36]. Lack of stable vertical transmission of pLM5578 during the outbreak or laboratory culture (despite the presence of a partitioning determinant, pLM5578_15) might account for plasmid absence. Alternatively, pLM5578 may have been carried only in a subset of L. monocytogenes strains that subsequently caused the outbreak. Presence of determinants related to bacterial conjugation (Fig. 5A) further confounds any speculation on the transmission of this plasmid.

Whole genome comparisons enable construction of an evolutionary model of the outbreak isolates

The genomes of 08-5578 and 08-5923 were collinear, and with the exception of the prophage and plasmid unique to 08-5578, all genetic diversity was accounted for by 2 short indels and 28 confirmed SNPs. The two indels included a single C/T base pair insertion in 08-5578 that introduced a frameshift truncation into coding sequence LM5578_1509 (buk, encoding a putative butyrate kinase). Secondly, there was a 21 bp (7 codon) in-frame deletion in the gltX gene (LM5578_0279) of 08-5578 encoding a putative glutamyl-tRNA synthetase. DNA sequencing confirmed that this gltX deletion was present in all examined LMACI.0040 isolates but that a wild type allele was present in all isolates typed as LMACI.0001 (Table 1). Alternatively, a wild type buk gene was only observed in isolate 08-5923 and a truncated butyrate kinase was encoded in all other isolates, independent of PFGE pattern (Table 1). Directed PCR and DNA sequencing of the SNP sites within the panel of clinical, food and environmental isolates was also performed (Table 1). Of the 28 SNPs, 27 were only present in isolate 08-5923. The singular remaining non-coding SNP was observed in each of LMACI.0040 isolates (at intergenic coordinate 2691224 of 08-5578) but was absent in all of the LMACI.0001 isolates, including 08-5923.

We propose a model to describe the chromosomal evolution of strains involved in this nation-wide outbreak using the distribution and segregation of genetic traits such as SNPs, indels and prophage (Fig. 6). This model assumes that the last common ancestor (LCA) of all strains encoded wild type buk and gltX loci, ϕLMC1 was absent, and the wild type state of all SNP positions was represented by the sequences in 08-5578 with the exception of intergenic SNP at 2691224. Isolate 08-5923 is predicted to be a direct descendent of the LCA, as it is the only isolate encoding a wild type buk gene. The remainder of the LMACI.0001 isolates differed from the LCA by only the single frameshift in buk, generating ancestor 2 (Fig. 6). From ancestor 2, acquisition of three independent traits (ϕLMC1, the 21 bp deletion in gltX and the intergenic SNP at coordinate 2691224) occurred, resulting in a third lineage defined by the LMACI.0040 PFGE pattern. In light of the high numbers of non-synonymous SNPs in 08-5923 and macromolecular genetic changes such as bacteriophage integration in the lineage represented by 08-5578, we speculate that three distinct, but highly related strains were circulating in this outbreak (Fig. 6). The diversifying SNPs in 08-5923 may have occurred in response to adverse environmental conditions prior to infection or host-specific conditions during infection, and ϕLMC1 likely arose by acquisition from other natural populations in the affected production facility.

Genome sequence alignments reveal that the outbreak isolates belong to evolutionary lineage II - clonal complex 8

Phylogenetic comparison to Listeria genome sequences representing each of the three L. monocytogenes evolutionary lineages, L. welshimeri and L. innocua indicated that 08-5578 and 08-5923 are both lineage II strains (Fig. 3A). The most closely related genome was strain EGDe, a serotype 1/2a lab-passaged animal isolate from 1924 [22]. For detection of large-scale genomic rearrangements, a dot plot comparison of 08-5578 versus the EGDe genome sequence was performed (see Additional file 3). This revealed that the common chromosomal sequences were collinear, with the exception of a large 50 kbp inverted segment around the putative origin of replication (coordinate 1 of EGDe) (see Additional file 3). Symmetric inversions around the origin of replication have previously been observed between other bacterial genome sequences [38]. In contrast to only 28 SNPs being observed between 08-5578 and 08-5923, 24660 high quality SNPs and 941 estimated indels were revealed in a comparison of 08-5578 and EGDe (data not shown). These included 17874 synonymous SNPs, 4692 non-synonymous SNPs and 475 indels occurring within 1214 predicted coding sequences, with 1 to 167 non-synonymous SNPs and/or indels occurring per gene. Moreover, two large contiguous regions present in both 08-5578 and 08-5923 represent probable prophages (designated ϕLMC2 and ϕLMC3) that are absent in EGDe (Fig. 2, see Additional file 3). ϕLMC2 was highly related to phage B025, and ϕLMC3 was inserted at comK and related to Listeria phage A006.

Comparative analysis of the regions used for multi-locus sequence typing (MLST) indicated that isolates 08-5578 and 08-5923 are more closely related to lineage II clonal complex 8 (CC8) than to EGDe (Fig. 3B). The MLST determinants used by Ragon et. al.[24] were extracted in silico from the genome sequences and 08-5923 was identical to a previously characterized CC8, sequence type 120 (ST120) strain (encoding abcZ allele 5, bglA allele 6, cat allele 2, dapE allele 29, dat allele 5, ldh allele 3, lhkA allele 1). The extracted sequences for 08-5578 were also identical to ST120 except for one SNP within the abcZ locus not currently listed in the Institute Pasteur's L. monocytogenes MLST database. The sole previous member of ST120 was identified as a clinical serotype 1/2a isolate derived from a CNS infection in New Zealand in 1995 [24]. Phylogenetic trees based on MLST and whole genome sequences were relatively congruent when considering evolutionary lineages, however genome sequences provide a more precise placement relative to other examined isolates. The pyrosequencing data also provided the opportunity to analyze the entire genome for micro and macro-diversity, including the full repertoire of phylogenetically relevant loci.

Virulence-associated determinants of the outbreak isolates informs lineage- and strain-specific virulence potential

The surface protein internalin, encoded by the inlA locus, is a principal virulence determinant promoting mammalian host cell invasion via interaction with E-cadherin on epithelial cells [39]. Truncation variants of inlA have been associated with decreased pathogenicity and significantly reduced invasiveness for human intestinal epithelial cells [40, 41]. Using nomenclature assigned by Ragon et al.[24], isolates 08-5578 and 08-5923 both encode a full length inlA allele type 2 that has previously been observed in isolates related to human illness [24]. This singular determinant however cannot account for differences in pathogenicity since foodborne isolates have been observed with intact inlA loci [42] and inlA allele 2 is also encoded by environmental isolates [43]. L. monocytogenes does encode multiple additional internalin and internalin-like coding sequences, although the function of each paralog is currently unknown [21, 44]. Isolates 08-5578 and 08-5923 both encode inlB, inlC, inlC2, inlD, inlE, inlF, inlG, inlI, inlJ and 11 other leucine-rich internalin-like coding sequences (see Additional file 4). The presence of an intact inlA locus and this compilation of other internalin-like loci may account in part for the pathogenicity of the serotype 1/2a outbreak strains sequenced in this study, but this requires further examination.

The Listeria pathogenicity island 1 (LIPI-1) encodes six significant virulence-associated loci: prfA (pleiotropic virulence transcriptional regulator), plcA, plcB (both encoding phospholipases), hly (listeriolysin O), mpl (metalloprotease) and actA (involved in actin-mediated motility). Previous phylogenetic analysis of this region from food, clinical, animal and environmental isolates indicated that this gene cluster was genetically diverse and clustered based on the three L. monocytogenes evolutionary lineages [25, 45]. BLASTn comparative analyses of this ~9 kbp region from 08-5578 and 08-5923 identified a very high sequence identity with other lineage II isolates, including 99.5% identity to the EDGe reference genome; ≥ 98% identity to all 23 lineage II isolates characterized by Ward et al.[25], and ≥ 98% nucleotide identity to all 21 lineage II isolates characterized by Orsi et al.[45]. In comparison, all LIPI-1 regions encoded by evolutionary lineage I and III strains in these previous studies were ≤ 95% identical to our lineage II strains.

PrfA regulates additional virulence determinants encoded elsewhere on the chromosome, such as bsh encoding a bile salt hydrolase that promotes survival within the gut [46] and uhpT encoding for a hexose phosphate permease for utilization of host carbon sources [47]. These latter two determinants were also encoded within 08-5578 and 08-5923. The presence of these determinants along with LIPI-1 is well conserved across L. monocytogenes strains independent of evolutionary lineage, so are unlikely to enhance virulence potential unless the allelic diversity observed between lineages results in a phenotypic change. In general, the contribution of lineage-specific genome differences to virulence and pathogenicity of L. monocytogenes are not well understood. For example, recent in silico analyses have indicated that while there are lineage II-specific genomic regions not present in lineage I, their role in virulence is not clear [30]. Furthermore, the presence of the Listeria pathogenicity island LIPI-3 operon encoding Listeriolysin S has been associated with evolutionary lineage I strains causing foodborne outbreaks [48]. This region was absent in the genomes of both 08-5578 and 08-5923.

Outbreak isolates harbor a putative mobile genetic island encoding translocation functions

A large 49.8 kbp contiguous region (coordinates 1836435-1886209 of 08-5578; coding sequences LM5578_1850 to LM5578_1903) was present in both sequenced genomes yet absent in all publicly available Listeria genome sequences to date, including EGDe (Fig. 2). The bordering coding sequences LM5578_1849 and LM5578_1904 were each 98% identical to contiguous EGDe coding sequences lmo1702 and lmo1703, respectively, implying that the 50 kbp operon-like structured region represents a genomic insertion within the ancestral chromosome of these isolates. Accordingly, putative serine recombinases are encoded in this region (loci LM5578_1855-58) and 16 bp imperfect inverted repeats are present at the borders in the intergenic regions between loci 1849/50 and 1903/04 (Fig. 5B), indicating that this is a horizontally acquired genetic island, hereafter designated Listeria genomic island 1 (LGI1).

Coding sequences within LGI1 exhibited sequence homology and were similarly organized as contiguous regions present within several environmental firmicutes, including Clostridium kluyveri, C. bolteae and Desulfitobacterium hafniense. These LGI1-like regions also appeared to be horizontally-acquired, based on skewed G/C content relative to neighboring sequences [49]. The genetic organization and predicted functions of several LGI1 loci resembled the proposed B. anthracis pXO1 plasmid-encoded secretion system [50], including several coding sequences that resembled putative type II and type IV secretion systems (T4SS) (Fig. 5B; and Additional file 2). Prototypical members of these systems are involved in pilus biogenesis and translocation of DNA-protein complexes or virulence effectors, respectively. Canonical T4SS genes predicted within LGI1 include virB4, virD4, and virB11, which encode ATPases involved in substrate recruitment to the transfer complex and substrate translocation, and virB5 and virB6 subunit genes, which form the core membrane-spanning transfer complex (Fig. 5B; and Additional file 2). Genes that likely contribute to pilus-like surface appendages were also detected (cpa, tad), and the presence of a dnaG gene encoding a putative primase neighboring the virD4 coupling protein homolog suggests that this genetic island may be mobilizable (see Additional file 2). Besides pXO1, Type IV secretion-like systems (T4SLS) features have also been described in Gram-positive plasmids such as B. anthracis plasmid pXO2 [51, 52] and the pheromone-inducible conjugative plasmid pCF10 of Enterococcus faecalis[53].

Gram-positive T4SLS also encode a putative peptidoglycan hydrolase and adhesin, both of which are predicted functions in LGI1 (LM5788_1873 and LM5578_1866, respectively). LM5788_1873 encodes a C-terminal NlpC/P60 domain, and is predicted to be involved in assembly of the core substrate transfer complex as a VirB1-like peptidoglycan hydrolase. A P60 domain is also present in the autolytic virulence determinant designated Listeria invasion associated protein (Iap) [54]. LM5578_1866 is predicted to encode a Sel1-like repeat (SLR) family protein that contains helical domains mediating protein-protein interactions [55]. Although protein adhesins are often used in Gram-positive bacteria to mediate conjugative attachment to target cells, SLR proteins also are involved for eukaryotic cell entry by Legionella pneumophilia[56, 57]. Whether this T4SLS is involved in translocation of pathogenicity effector molecule(s) (in addition to probable DNA-protein translocation) is tempting to speculate, but further study is required.

LGI1 also encoded a homolog (LM5578_1862) to the multidrug efflux proton:drug antiporter EmrE implicated in resistance to toxic cationic hydrophobic compounds such as quaternary ammonium compounds and tetracycline [58]. This gene was flanked by a MarR-family transcriptional regulator loci (1425) and a putative DNA-directed RNA polymerase sigma-24 subunit (rpoE, 1861) encoding a specialized ECF (extracytoplasmic function) family sigma factor, which is part of the bacterial stress response regulon [59]. A two component signal transduction system (sensor histidine kinase and response regulator; LM5578_1852 and 1851) and restriction modification components (LM5578_1850 and 1853) also were present. Cumulatively, LGI1 is unique compared to all currently sequenced L. monocytogenes but the contribution to pathogenicity or environmental persistence is unconfirmed.

Bacterial isolates

Clinical, food and environmental isolates (Table 1) were collected by the Ontario Central Public Health Laboratory and the Canadian Food Inspection Agency during an outbreak investigation of L. monocytogenes and transferred to the National Microbiology Laboratory (NML) for additional subtyping. The clinical isolates examined in this study were collected from individual outbreak-associated cases.

Serotyping and Pulsed-field Gel Electrophoresis

Serotyping was performed by slide agglutination with antisera prepared at the NML according to Seeliger and Höhne [60]. Molecular serotyping was completed using a multiplex PCR scheme [61]. Pulsed-field gel electrophoresis (PFGE) was performed according to the PulseNet standardized protocol using restriction enzymes Asc I and Apa I. PFGE patterns were designated using BioNumerics after comparison to the PulseNet Canada database.

Genome sequencing and bioinformatics

Genome sequencing was performed on the Roche GS FLX™ standard platform, as per the manufacturer's recommendations. The relative genome coverage achieved was approximately 40× for 08-5578 (Q40+ bases = 99.88%) and 36× for 08-5923 (Q40+ bases = 99.91%). Within the 08-5578 project the coverage of pLM5578 was 68×, indicating that the plasmid copy number was between 1-2 copies per cell. De novo sequence assembly was completed using Roche's Newbler assembler and the Staden software package [62] after a fosmid mate-pair end library was constructed. Gap closure was completed by Sanger-based sequencing of select regions to bridge contigs. Presumptive insertion/deletions (indels) associated with homopolymeric regions were verified or corrected after Sanger-based sequencing of directed PCR amplicon templates generated from the corresponding regions in both 08-5578 and 08-5923 by proof-reading Taq polymerase. Sequences were deposited into GenBank under the accession numbers CP001602 (08-5578), CP001603 (pLM5578) and CP001604 (08-5923).

Pyrosequencing was completed within three days to provide greater than 99% of the genome coverage. This was sufficient for preliminary analysis such as identification of plasmid and phage determinants and candidate SNP's. Closed circular chromosomal sequences were obtained within an additional 5 weeks after fosmid libraries were constructed, contigs ordered and all gaps closed. Three additional weeks were required for resequencing of SNP and homopolymeric regions to confirm polymorphisms and indels between the two sequenced isolates. The costs associated with this project (in US dollars) were $5000 for the GS-FLX run covering both isolates, $500/isolate for the fosmid libraries and related consumables, and $2000/isolate for primer synthesis and Sanger-based sequencing for gap closure and resequencing.

Annotation of coding sequences for isolate 08-5578 was completed using GenDB v2.2 [63], BASys [64], Glimmer3 [65] and by comparison to RefSeq annotations completed for L. monocytogenes strain EGDe [22]. Sequences encoding rRNA and tRNA were identified by RNAmmer version 1.2 [66] and tRNAscan-SE version 1.23 [67], respectively. Artemis was used for manual annotation and additional manipulations [68]. Multiple sequence alignments of whole genomes and select loci were generated using MAUVE [69] and ClustalW [70], respectively, and maximum likelihood phylogenetic trees were generated using the Phylogeny Inference Package (Phylip version 3.68). Distant homologies between hypothetical coding sequences and entries at the PDB, COG, PFAM and SCOP databases were detected using the FFAS03 server [71]. Single nucleotide polymorphisms (SNPs) were identified between the sequenced isolates using NUCMER and dot plots were generated using MUMmerplot [72]. The circular representation of the genome of 08-5578 was created using DNAplotter [73] and feature tables generated using Artemis.

Reference Listeria genome sequences were obtained from GenBank for isolates EGDe (accession number NC_003210), F2365 (NC_002973), HCC23 (NC_011660), Clip11262 (NC_003212; L. innocua), SLCC5334 (NC_008555; L. welshimeri) and L. monocytogenes phage B025 (NC_009812). Accession numbers for plasmid DNA sequences of pLI100 and pXO2 were NC_003383 and NC_001496, respectively. In addition, the DNA sequences of L. monocytogenes strains F6854 and H7858 [21] and plasmid pLM80 from strain H7858 were obtained from the J. Craig Venter Institute Comprehensive Microbial Resource http://cmr.jcvi.org/cgi-bin/CMR/CmrHomePage.cgi. A cross match table (Additional file 5) of the coding sequences between all of these L. monocytogenes genomes was generated using the publically available annotations as of September 18^th, 2009 with a reciprocal BLASTp alignment threshold of 80% identity and 80% hsp. Draft genomes sequences at the Broad Institute were also screened for select regions using BLAST available at their database http://www.broad.mit.edu/annotation/genome/listeria_group/MultiHome.html. MLST loci data were downloaded from the Institute Pasteur Listeria monocytogenes MLST Database http://www.pasteur.fr/recherche/genopole/PF8/mlst/Lmono.html.

PCR screening and SNP confirmation

All polymerase chain reactions were performed with Invitrogen HiFi Platinum proof-reading Taq polymerase following the manufacturer's directions, and with 1 uM of each oligonucleotide (see Additional file 6) using the following thermocycling conditions: 94°C for 5 min.; 35 cycles of 94°C for 30 sec., 50°C for 30 sec. and 68°C for 45 sec.; followed by 68°C for 5 min. Amplicons were visualized after electrophoresis in 1% agarose-TBE gels. Confirmation of SNP sites was achieved by Sanger-based sequencing of targeted amplicons using the same oligonucleotides as used for PCR amplification.