Enrichment dynamics of Listeria monocytogenes and the associated microbiome from naturally contaminated ice cream linked to a listeriosis outbreak

1Optimization of enrichment methods to culture target pathogens from complex environmental, food and clinical samples is an ongoing challenge. Traditionally, samples are incubated in nonselective and/or selective enrichment broths and then plated onto selective media. Enrichment methods for specific foodborne pathogens will benefit from an improved understanding of the taxonomic diversity and relative abundance of microbiota that co-culture during enrichment. Here, we use culture independent next generation sequencing (NGS) to characterize the microbiome at four hour intervals using three different enrichment methods used for recovery of Listeria monocytogenes from naturally contaminated ice cream.

2L. monocytogenes was first reported in 1926 by Murray, Webb and Swann as the causative agent of illness in >rabbits and guinea pigs in a laboratory breeding unit. Although it was long suspected that food might be a mode of transmission for human listeriosis, it was not until after 1980 that several outbreaks conclusively linked L. monocytogenes to foods including, coleslaw, milk, cheese, meat, pâté, and jellied pork tongues. Listeriosis outbreaks in the United States over the past five years have been associated with contaminated cheeses, stone fruits, ice cream, cantaloupes and caramel apples. Results of the study presented here were obtained from bacteriological analysis of samples from the 2010 to 2015 listeriosis outbreak with several case-patients linked to milkshakes made from contaminated ice cream. Analysis of L. monocytogenes in ice cream samples manufactured in the implicated production line provided information about the prevalence and level of L. monocytogenes and activity of L. monocytogenes in milkshakes prepared from the ice cream. These analyses did not identify Listeria species other than L. monocytogenes in the ice cream samples.

3A variety of enrichment media and methods have been developed for detection of L. monocytogenes. The three commonly employed methods examined in this study are as follows:

41) BLEB as described in the U.S. Food and Drug Administration (FDA) Bacteriological Analytical Manual (BAM): 4 h incubation at 30 °C in buffered Listeria enrichment broth (BLEB) without antibiotics, followed by 44 h incubation at 30 °C in BLEB with antibiotics; 2) HFB-FB as described in the International Organization of Standardization 11290–1: 24 h enrichment in Half-Fraser broth (HFB) at 30 °C followed by 24 h in Fraser broth (FB) at 37 °C; and 3) UVM-FB as described in the U.S. Department of Agriculture (USDA) Microbiological Laboratory Guidebook (MLG): 24 h incubation at 30 °C in University of Vermont modified broth (UVM) for 24 h, followed by 24 h in FB at 37 °C.

5Once L. monocytogenes is detected and isolated from a food or environmental source, pulsed-field gel electrophoresis (PFGE) has been a widely applied method for subtyping isolates. PFGE has been the gold standard for the FDA and the Centers for Disease Control and Prevention (CDC) for foodborne outbreak investigations for over 15 years. Recently however, whole genome sequencing (WGS) has been employed to improve resolution of PFGE for identifying closely related strains. Currently, both PFGE and WGS require isolation of confirmed cultures, which typically requires at least 5 to 7 days for L. monocytogenes.

6L. monocytogenes from the same lot of ice-cream examined here was previously enumerated using an Most Probable Number (MPN) method. All samples from the lot tested positive for L. monocytogenes with a geometric mean of 3.35 MPN −1g. To complement enumeration data, we used a metagenomic approach to examine how L. monocytogenes and other members of the microbiota of the naturally contaminated ice cream responded to three commonly used enrichment methods (BLEB (FDA), HFB-FB (USDA), UVM-FB (ISO)). Ribosomal RNA amplicons and metagenomic sequence data from time-points every 4 h during a 48 h period were used to describe the taxonomic composition of co-enriching bacterial taxa and provide data to explore whether a hybrid culture/metagenomic approach can be used to source track target pathogens before they are fully isolated.

716S rRNA amplicon sequencing revealed that proportional abundances of L. monocytogenes remained low (0–10%) until 24 h of enrichment, even though each enrichment method employed selective antimicrobials (Fig. 1). After 24 h of selective enrichment, relative abundances of L. monocytogenes increased at each successive time point until 40 h, at which time relative abundances ranged from 90% to near 100% of the microbial community. Interestingly, L. monocytogenes was found to be more abundant at hour 0 than at the three subsequent time points. BLEB and HFB-FB enrichment resulted in higher proportional abundances (although not statistically significant) of L. monocytogenes at 24 to 36 h than UVM-FB, while all enrichment yielded equivalent results at later time points.

8An examination of other community members, using 16S rRNA gene amplicons, revealed a predominance of Anoxybacillus, followed by Serratia, Geobacillus, and Streptococcus species (Fig. 1). During enrichment for all three methods (BLEB, UVM-FB, HFB-FB), Bacillaceae genera, Anoxybacillus and Geobacillus increased rapidly from a combined relative abundance of approximately 45% at hour 0 to almost 90% at hours 4, 8 and 12 (Fig. 1b). Taxonomy based on shotgun sequencing also supported the presence of Anoxybacillus and Geobacillus species (Fig. 2). Species of both genera are reported as moderately thermophilic and in this study appeared to have an advantage over L. monocytogenes during early incubation at 30 °C. Additionally, shotgun metagenomic data suggested the presence of two other thermophiles, Thermus parvatiensis and T. thermophilus (Fig. 2). Anoxybacillus spp. have an optimum growth temperature (OGT) ranging from 50 to 62°C and their close relatives, Geobacillus spp., have a slightly higher OGT of 55 to 65°C; Thermus spp. have an OGT ranging from 50 to 82°C.

9Another interesting observation of microbial dynamics during enrichment was the change in relative abundance of Serratia during the first 24 h of enrichment which appeared to mirror that of Listeria (Fig. 1). However, after 24 h of enrichment Serratia was outcompeted by other members of the microbial community, mainly Listeria monocytogenes in all three enrichment broths. Assuming that hour 0 (initiation of the experiment) resembled the mixed culture microbiota in ice cream, Serratia was a well-represented constituent of that microbiota, comprising approximately 20–30%. As observed for L. monocytogenes, Anoxybacillus and Geobacillus spp. outcompeted Serratia spp. during the 4 to 12 h time points. In terms of other taxa present in the ice cream microbiome, relative abundances of Streptocococcus and Erwinia remained consistent during the first 24 h of enrichment until L. monocytogenes began to outcompete the rest of the community, while the relative abundances of Pseudomonas species increased after 24 h and decreased after 40 h of enrichment. In the absence of selective antimicrobial agents, we observed a completely different microbial community dynamics during the 48 h of non-selective enrichment (Additional file 2: Figure S2). At 0 to 8 h, we observed a very similar pattern to that in selective enrichments; however, after 16 h of non-selective enrichment, Bacillus species dominated the ice cream microbiome, comprising ~80% of the total population. After 24 h, Lactococcus emerged and became the slightly dominant species over Bacillus for the remainder of the non-selective enrichment.

10Though 16S rRNA gene sequencing data provided a detailed description of the ice cream microbiome and resulted in the detection of L. monocytogenes, a deeper sequence analysis utilizing shotgun metagenomics was necessary to characterize strains of L. monocytogenes. Analysis of these samples revealed that near 100% genome coverage can be achieved as early as 24 h of enrichment (Table 1). L. monocytogenes-specific sequence reads constituted 0.20 to 0.76% of the metagenome, varying according to enrichment broths employed. At 24 h, we achieved 23 to 97% coverage of the Listeria genome at a 7.5× to 12× depth of coverage with 15 to 35 M total metagenomic sequence reads. A similar trend was also observed after 28 h, however, shotgun metagenomic reads from this time point yielded slightly higher genome coverages (44 to 98%) and the depth of coverage ranged from 14.5× to 95×. The highest coverage among these two time points and three enrichment schemes was achieved in BLEB ice cream enrichments at hour 28, which had 45 M total sequence reads, (>95× depth) with near complete genome coverage (>97%), even though the UVM-FB ice cream enrichments contained 82 M total reads (15× depth, 44% genome coverage) (Table 1).

11Interestingly three variants of 16S rRNA gene sequences were observed in L. monocytogenes populations from the selective enrichments. Variable nucleotides occurred in the V2 region of the 16S rRNA gene (at positions 159 and 174, E. coli) with one variant comprised of AC nucleotides at those positions, a second with GT, and a third, the most abundant, possessing GC (Fig. 3, Table 2). Analysis of 68 closed genomes, available from the PATRIC database, revealed that the majority of L. monocytogenes genomes encoded six copies of the 16S rRNA gene (85.3%), with the remaining genomes encoding five (11.8%), or four (5.9%) copies. Six 16S rRNA gene variant sequence motifs were identified, with 96% of the 68 genomes harboring at least one of the variants identified in this study (GC, GT, or AC). Interestingly, among the 68 closed genomes, two sub-variants of type GC were identified (type GC.2 and GC.3), characterized by changes to nucleotides in the G and C allele positions. These minor variants were found in five L. monocytogenes lineage I genomes, with type GC.3 comprising the sole 16S rRNA gene sequence variant present in three of these five genomes. Type GC was present in the four lineage III genomes analyzed. Overall, the distribution of major 16S rRNA gene sequence variants was not significantly associated with lineage.

12Analysis of intra-genomic distribution of 16S rRNA gene variants (i.e. variant types representing the complete complement of 16S rRNA gene copies within a single genome) revealed the majority of L. monocytogenes genomes encoded a single 16S rRNA gene variant type (79.1%), though genomes containing multiple variant types (20.6%) were common. Overall, types GC or AC were most prevalent in genomes encoding a single variant, 51.5 and 20.6%, respectively; while types GC and GT represented the dominant type in mixed 16S rRNA gene variant genomes, 19.1 and 17.6%, respectively. Interestingly, type GT was predominantly found in mixed variant genomes, and accounted for the sole variant in only two genomes, whereas type AC was present in a single mixed variant genome (Table 2).

13It is noteworthy that the AC 16S rRNA gene variant, which was the least abundant of the three types up to hour 40, and less common among the reference set of 68 closed genomes, overtook type GT at hours 44 and 48 across all of the media employed in the study. Type GC remained the most abundant throughout all time points (Fig. 3). This observation was corroborated by analysis of shotgun datasets against a Listeria specialty database, suggesting the possibility that of three potentially distinct L. monocytogenes strains were present. CosmosID identified strain- specific biomarkers (AACABABA, AACABB, and AACABD) in the metagenomes, which was evidence for three putative L. monocytogenes variants (Fig. 4). Furthermore, we observed an interesting strain interplay in terms of their abundance over time where strain AACABD was most abundant during the first 20 h of enrichment and became the least abundant at later time points, due to the rapid upsurge of strain AACABABA, beginning at 28 h (Fig. 4). This same phenomenon was observed in the 16S data which suggests we may be documenting something biologically relevant.

14Every food has its own innate or imparted microbiome that will respond to enrichment conditions according to complex eco-physiology. Frequently, antagonistic microorganisms are co-enriched along with target pathogens. For example, Paenibacillus spp.,which are capable of inhibiting and killing Salmonella, were also enriched using protocols outlined in the FDA BAM to recover Salmonella from tomatoes and cilantro. Thus, continued optimization of reference enrichment protocols is still needed. In the case of L. monocytogenes, additional selective agents as well as changes in media formulations may improve efficiency of recovery; however, the situation may be more challenging in this case due to relatedness of L. monocytogenes to its co-enriching Bacilli relatives. The increase (almost doubling) of Anoxybacillus and Geobacillus spp. during the first 8 h of enrichment of the ice cream microbiome suggested these bacteria were able to adapt more readily to the environmental changes inherent to this study (Fig. 1).

15Length of lag phase and growth rate typically depend on specific environmental parameters, as well as fitness of the bacterial cells. Many factors have been reported to play a role in length of lag time and/or growth rate for a given bacterial species, including nutritional content, pH, physical environment, inorganic nutrients, temperature, rate of temperature change, water activity, gas atmosphere, inhibitors, spore germination, and initial cell levels, as well as fitness, age, size, and health of individual cells. Although lag time and/or growth rate modeling has been proposed to describe what should occur within a specific set of parameters, such results were variable. Hence, it is evident that relationships between growth environments and lag time and growth rate are complex, especially for testing foodborne pathogens in food and environmental samples. Additionally, the effect of co-occurring bacterial species in enrichments has not been extensively described or considered in the literature since the analytic tools to do so have not been available and it is not always possible to obtain a sufficient number of naturally contaminated samples that are 100% positive for the target pathogen with relatively homogeneous levels of contamination.

16The sequence depth and genome coverage achieved in the ice cream metagenomes was sufficient to generate draft L. monocytogenes genomes after 24–28 h of enrichment, which is considerably shorter than the 96 h typically required to enrich and isolate viable bacteria from food samples. Near complete genome coverage would facilitate high-resolution phylogenetic analysis, identification of virulence and antibiotic resistance factors, as well as subtyping. Thus this approach might be applicable for rapid and accurate microbial forensics in future work but remains to be more extensively validated.

17Type GC represented the predominant allele and was present among the majority of strains belonging to lineage I (808/914 total lineage I strains examined), and was the representative pattern in 30 of the outbreak isolates sequenced in FDA labs. Type GC was also heterogeneously distributed within and between multilocus sequence typing (MLST) clonal groups, designated as clonal complex (CC), for strains within lineage II (361/914 total lineage II strains) Type GC sequences were also present in three lineage III genomes. Type GT and type AC were rare among lineage I (19/914) and lineage II (1/914) strains, respectively. The GC type was significantly prevalent among the human clinical strains from both lineages. The two variable positions (159 and 174) have been shown to contribute to the structure of helix 8, located in the 5′ major domain of the 16S rRNA gene (E. coli), which forms direct contact with nearby helices (e.g., helix 6) in the small ribosomal subunit. As correct 16S rRNA gene folding is important for proper ribosome assembly, translational kinetics, and fidelity it is intriguing to consider how specific allelic changes within this region may affect L. monocytogenes growth under infection-associated conditions.

18The ability to rapidly and accurately identify the etiologic agent of a foodborne outbreak from a variety of food matrices is critical for preserving public health. Culture-based recovery methods have been used for over 100 years and many have been optimized for recovery of major bacterial and fungal pathogens. We demonstrated that shotgun metagenomic sequencing facilitated the characterization of enrichment microbiota dynamics and identified the presence of competitive species co-occurring during enrichment. As L. monocytogenes was enriched, so were other bacterial species. The taxonomy and dynamics of co-enriching species are likely unique to each food commodity and may play a significant, as yet unidentified, role in the recovery and growth of target pathogens. Lag times and growth rates for certain pathogens may be significantly influenced by those co-enriching species “native” to the food commodity or its production environment. Thus, the approach taken in this study, namely to use a culture independent method to describe culture dependent dynamics of L. monocytogenes recovery and growth, provides insight and makes it possible to improve both sequence- and culture-based L. monocytogenes detection methods. The validation that Resphera Insight can accurately identify L. monocytogenes using 16S rRNA amplicons is quite useful. As is the likelihood that draft assemblies of L. monocytogenes from shotgun metagenomic data may provide equivalent trace-back utility to WGS genomes - although more work will be needed to confirm this utility.

19Ice cream samples (80–85 g /scoop) were aseptically transferred to Whirlpak bags, allowed to stand at room temperature for 20 to 30 min to fully melt. Analytical units were set up as follows: 25 g portions were added to 225 ml of enrichment broth, and stomached for 1 min. Four scoops were used for each of three enrichment protocols.

20Ice cream samples were incubated at 30 °C for 48 h in buffered Listeria enrichment broth (BLEB) (Cat. CM0897, Thermo Scientific Inc., Waltham, MA). Selective supplements (acriflavin hydrochloride 10 mg/l, nalidixic acid 40 mg/l and cycloheximide 50 mg/l, Cat. SR0149, Oxoid, UK) were added after 4 h of incubation.

21Ice cream samples were incubated at 30 °C for 24 h in University of Vermont modified broth (UVM) after which, 0.1 ml of each culture was transferred to Fraser broth (FB) and incubated at 37 °C for 24 h. Selective agents were added to UVM and FB prior to enrichment.

22Ice cream samples were incubated at 30 °C for 24 h in Half-Fraser broth (HFB) after which, 0.1 ml of culture was transferred to FB and incubated at 37 °C for 24 h. Selective agents were added to HFB and FB prior to the enrichment.

23For each of the enrichment protocols, 4 ml was taken from each sample at hour 0. After that, every 4 h during the 48 h incubation, enrichments were stomached for 1 min and 4 ml samples taken from each of the four replicates of the three enrichment mixes. Samples were immediately frozen at −20 °C for subsequent DNA analysis.

24Genomic DNA was extracted using DNeasy Blood and Tissue kit (Cat No. 69506, Qiagen, Germantown, MD, USA) following the protocol for Gram-positive bacteria with minor modification: 1.5 ml of the culture was pelleted (5000 × g, 15 min) and the pellet resuspended in 200 ml of enzymatic lysis buffer containing 20 mM Tris-HCl (pH-8.0), 2 mM Sodium EDTA, 1.2% Triton X- 100, 20 mg/ml of lysozyme. The samples were subsequently incubated for 60 min at 37 °C. DNA was extracted following manufacturer’s protocol.

25Sequencing libraries were prepared using the Truseq Nano prep kit (Illumina, SanDiego, CA, USA), according to the manufacturer’s specifications.

26PCR for 16S rRNA gene sequencing targeted the first 330 bases of the V1-V2 region, using the following PCR primers 27 F1 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGAGTTTGATCMTGGCTCAG and 336R1 GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGACTGCTGCSYCCCGTAGGAGTCT.

27PCR cycling conditions consisted of an initial denaturation at 94 °C for 2 min, followed by 25 cycles of 94 °C for 40 s, 56 °C for 15 s, 68 °C for 40 s, and a final extension at 68 °C for 5 min. The amplicons were indexed using Illumina Nextera indexing primers following the manufacturer’s instructions (Illumina, San Diego, CA, USA).

2816S rRNA gene amplicon sequencing was performed on an Illumina MiSeq using the 500 cycle V2 chemistry and cartridges. Shotgun data was obtained by sequencing on an Illumina NextSeq, using V2 chemistry (2 by 150).

29Raw paired-end reads from the MiSeq platform were merged into consensus fragments by FLASH and subsequently filtered for quality (max error rate 1%) and length (minimum 300 bp) using Trimmomatic and QIIME. Spurious hits to the PhiX control genome were identified using BLASTN and removed. Passing sequences were trimmed of primers, evaluated for chimeras with UCLUST (de novo mode), and screened for chloroplast and mitochondrial contaminants using the RDP classifier with a threshold of 0.5. Sequences were further evaluated for unknown contaminants using a sensitive BLASTN search against the GreenGenes 16S database. High-quality 16S sequences were submitted for high-resolution taxonomic profiling using Resphera Insight (Baltimore, MD, www.respherabio.com). Sequence counts were rarefied to 5000 sequences for each independent replicate for downstream analyses.

30To evaluate the potential for multiple cultured strains in the enrichments, all 16S rRNA gene sequences assigned unambiguously to L. monocytogenes by Resphera Insight were aligned by PYNAST to a smaller template multiple sequence alignment (MSA) of 5000 randomly selected sequences from the same set (generated by MUSCLE). The full PYNAST MSA was filtered for positions with > 10% gaps and passing positions submitted for entropy calculation. MSA Positions with measured entropy > 0.7 were utilized to assign putative strain membership.

31To perform an external validation of the species level accuracy of Resphera Insight, Resphera Biosciences and Center for Food Safety and Applied Nutrition (CFSAN) scientists collaborated to interrogate 1695 whole-genome shotgun datasets from the GenomeTrakr Project (NCBI Project ID PRJNA183844) designated as L. monocytogenes isolates. Raw paired-end sequences were filtered for quality and length, followed by merging of overlapping sequences using FLASH. Merged reads were screened for 16S rRNA fragments using Bowtie2 against a broad database of 16S rRNA gene sequences, with additional BLAST-based filtering to confirm location specific query matches to a reference L. monocytogenes 16S rRNA gene (L. monocytogenes strain 07PF0776; NCBI accession NR_102780.1).

32Passing sequences were submitted to Resphera Insight for high-resolution taxonomic identification. The primary measure of performance was the Diagnostic True Positive rate (DTP), defined as the percentage of reads with an unambiguous assignment to L. monocytogenes and differences in accuracy associated with changes in read length and gene position were evaluated. We also computed Sensitivity (SN), defined as percentage of reads with an unambiguous or ambiguous assignment to L. monocytogenes.

33Overall, across all 1695 isolates, for fragments ≥200 bp originating within 16S reference gene positions 1–36, Resphera Insight achieved a mean diagnostic true positive rate of 99.43% and a mean sensitivity of 99.94%, with a misassignment rate of 0.06%. We observed increased DTP rates associated with increasing fragment lengths and a loss of species level resolution 3′ to 16S reference gene position 200 (Additional file 1: Figure S1). For fragments <200 bp, we also found a loss of DTP resolution, as more assignments became ambiguous, but sensitivity was not reduced.

34To establish benchmark false positive rates of L. monocytogenes, we simulated amplicon fragments from our primer region for eight closely related Listeria species (L. fleischmannii, L. grayi, L. innocua, L. ivanovii, L. rocourtiae, L. seeligeri, L. weihenstephanensis, L. welshimeri). A total of 10,000 sequences per species were simulated, lengths 250–500 bp, with a random nucleotide error rate of 0.5%. Overall, the average sensitivity of detection of these organisms was 99.9%, with a misclassification rate (assignments to L. monocytogenes) of 0.07%.

35Unassembled metagenomic sequencing reads were directly analyzed by Genius bioinformatics software package (CosmosID Inc., Rockville, MD), described elsewhere to achieve identification at the species, subspecies, and/or strain level and quantification of relative abundance. Briefly, the system utilizes curated genome databases and a high performance data-mining algorithm to rapidly disambiguate millions of metagenomic sequence reads into discrete microbial taxa. The pipeline has two separable comparators. The first consists of a pre-computation phase and a per-sample computation. The input to the pre-computation phase is a curated reference microbial database and its output is a whole genome phylogeny tree, together with sets of fixed length n-mer fingerprints (biomarkers) that are uniquely identified to distinct nodes of the tree. The second per-sample, computational phase searches the millions of sequence reads against the fingerprint sets. The resulting statistics are analyzed to give fine-grain composition and relative abundance estimates at all nodes of the tree. Overall classification precision is maintained through aggregation statistics.

36Furthermore, CosmosID and CFSAN collaboratively developed a specialized Listeria database including genome sequences of all L. monocytogenes available at the time of this analysis. As for the pre-computational phase describe above, the L. monocytogenes database was organized as a phylogenetic tree with thousands of unique and shared biomarkers specific to distinct clades, branches, nodes and leaves of the tree along with specific alphabetical fingerprints reflective of their phylogenetic hierarchies. All metagenomic datasets were analyzed against this L. monocytogenes database to investigate the potential presence of multiple strains of L. monocytogenes in the ice cream samples. Predicted % genome coverages of L. monocytogenes were estimated based on the mean of % total match – a statistics derived from CosmosID algorithm that approximates genome coverage. Predicted depths of coverage were calculated as X = (# of Listeria reads X Read length) / (L. monocytogenes genome size X Predicted % genome coverage).