ABSTRACT
Spoilage and pathogenic spore-forming bacteria are a major cause of concern for producers of dairy products. Traditional agar-based detection methods employed by the dairy industry have limitations with respect to their sensitivity and specificity. The aim of this study was to identify low-abundance sporeformers in samples of a powdered dairy product, whey powder, produced monthly over 1 year, using novel culture-independent shotgun metagenomics-based approaches. Although mesophilic sporeformers were the main target of this study, in one instance thermophilic sporeformers were also targeted using this culture-independent approach. For comparative purposes, mesophilic and thermophilic sporeformers were also tested for within the same sample using culture-based approaches. Ultimately, the approaches taken highlighted differences in the taxa identified due to treatment and isolation methods. Despite this, low levels of transient, mesophilic, and in some cases potentially pathogenic sporeformers were consistently detected in powder samples. Although the specific sporeformers changed from one month to the next, it was apparent that 3 groups of mesophilic sporeformers, namely, Bacillus cereus, Bacillus licheniformis/Bacillus paralicheniformis, and a third, more heterogeneous group containing Brevibacillus brevis, dominated across the 12 samples. Total thermophilic sporeformer taxonomy was considerably different from mesophilic taxonomy, as well as from the culturable thermophilic taxonomy, in the one sample analyzed by all four approaches. Ultimately, through the application of shotgun metagenomic sequencing to dairy powders, the potential for this technology to facilitate the detection of undesirable bacteria present in these food ingredients is highlighted.
IMPORTANCE The ability of sporeformers to remain dormant in a desiccated state is of concern from a safety and spoilage perspective in dairy powder. Traditional culturing techniques are slow and provide little information without further investigation. We describe the identification of mesophilic sporeformers present in powders produced over 1 year, using novel shotgun metagenomic sequencing. This method allows detection and identification of possible pathogens and spoilage bacteria in parallel. Strain-level analysis and functional gene analysis, such as identification of toxin genes, were also performed. This approach has the potential to be of great value with respect to the detection of spore-forming bacteria and could allow a processor to make an informed decision surrounding process changes to reduce the risk of spore contamination.
INTRODUCTION
Milk and resultant dairy products can become contaminated by bacteria from a number of sources, including production and processing facility contaminants. Soil, bedding, feed, feces, and teat surfaces all harbor bacteria that can transfer to raw milk (1, 2), and milking and housing practices are a potential contributor to raw milk contamination (3). Although many of the bacteria present in milk are killed by traditional processing techniques, bacterial spores can survive heat treatment and desiccation (1, 4). Furthermore, within processing facilities, microbial biofilms formed on equipment surfaces can be persistent. These are frequently resistant to cleaning, and cells, including spore-forming bacteria, from these biofilms can slough off to contaminate products (5, 6). Furthermore, spores, regardless of their origin, can withstand further processing and remain in a dormant form in powdered dairy products thereafter (7, 8).
On the basis of culture-based analyses, sporeformers identified frequently in dairy powders include representatives of the genera Bacillus, Geobacillus, Anoxybacillus, Lysinibacillus, Brevibacillus, and Paenibacillus (9, 10), as well as Clostridium (11, 12). Some of these sporeformers, including Bacillus licheniformis, Bacillus subtilis, Bacillus amyloliquefaciens, members of Bacillus cereus sensu lato (13), and Geobacillus species (14), are associated with spoilage of dairy products. Also, some members of B. cereus sensu lato have the potential to cause food poisoning. More specifically, B. cereus can cause diarrheal or emetic food poisoning. Diarrheal food poisoning is caused by enterotoxins, i.e., nonhemolytic enterotoxin (Nhe), cytotoxin K (CytK) or hemolysin BL (HBL), with symptom onset occurring 8 to 16 h after ingestion (15, 16). Emetic food poisoning is caused by the emetic toxin cereulide (Ces), which is produced by a nonribosomal peptide synthetase and results in vomiting within a few hours of ingestion (16). Toxin-producing Clostridium species, such as Clostridium botulinum and Clostridium perfringens (7), are also a concern but are less common in dairy powders than bacteria of the class Bacilli (17).
Despite safety and economic concerns, there is little legislation governing the total numbers of permissible sporeformers in dairy powders. However, processors and customers generally set strict guidelines to ensure high standards. U.S. powders destined for international customers have limits of <500 CFU · g−1 thermophilic sporeformers and <1,000 CFU · g−1 mesophilic sporeformers (18). These limits reflect the fact that traditional culture-based detection and enumeration methods for spore-forming bacteria rely on a variety of temperature treatments to differentiate between heat-resistant and highly heat-resistant sporeformers, incubation temperatures (to differentiate among psychrotolerant, mesophilic, and thermophilic sporeformers), and incubation conditions (to differentiate between aerobic and anaerobic sporeformers) (17, 18). A number of selective and chromogenic agars may also be employed to select for and/or identify pathogenic groups, such as Clostridium species and B. cereus group (19, 20). Further biochemical or molecular methods can be subsequently utilized to identify the species isolated (1).
The heavy reliance on culture-based assays has been highlighted as a concern in recent years. These methods are labor intensive, requiring many medium types, sample treatments, incubation temperatures, and conditions per sample to get an overview of the total numbers of different groups of bacteria present. These need to be followed by additional testing to identify or confirm identity of the species present and, where necessary, further analysis to determine if species present contain toxin genes (21).
Here, shotgun metagenomic DNA sequencing (22) is employed to identify bacteria present in dairy powders, as an extension of the recent application of the technology to identify pathogenic strains in a fermented dairy beverage (23) or spoilage microorganisms in cheese (24). Shotgun metagenomics can generate a considerable volume of data, enabling the detection of culturable and nonculturable microorganisms, with accurate identification to species level and sometimes to strain level. It also allows an investigation of specific gene sequences of interest, such as, for example, antimicrobial resistance or toxin-encoding genes (25, 26). However, this analysis remains expensive and poses data analysis-related difficulties due to the frequently large quantity of data generated (22). With a view to beginning to bridge the gap toward the application of this powerful technology to food quality and safety testing, this study employed shotgun metagenomics to test 12 powdered dairy samples produced monthly at the same location over one calendar year to determine the mesophilic spore content of these powders. For the purpose of comparison, one sample (from August [A]) was examined in greater depth to assess the mesophilic spore content (AM) and thermophilic spore content (AT), as well as the culturable mesophilic spore content (AMP) and culturable thermophilic spore content (ATP).
RESULTS
Shotgun metagenomics can be used to identify and determine the functional potential of low-abundance sporeformers present in dairy powders.A total of 12 whey powder samples, produced monthly at the same location over one calendar year, were sampled upon exit from the spray drier prior to packing and collected for spore analysis. Following spore pasteurization of reconstituted powders, in which samples were incubated at 80°C for 12 min to select for spore-forming bacteria, a total of <400 CFU · g−1 mesophilic and thermophilic sporeformers was enumerated by plating whey solutions on brain heart infusion (BHI) agar and incubating at 30°C and 55°C, respectively, for 24 h. These are within previously described consumer specified limits of <500 CFU · g−1 thermophilic sporeformers and <1,000 CFU · g−1 mesophilic sporeformers (18).
To facilitate more in-depth characterization, sporeformer enrichment was performed and DNA extraction and shotgun metagenomic sequencing were completed. The enrichment approach taken addressed the risk of reads arising from sequencing of DNA from dead cells, but it had the drawback of underrepresenting spores that did not germinate under the conditions employed. The average number of reads after quality filtering and trimming per sample was 1,106,747. Kaiju, Kraken, and MetaPhlAn2 were all used to assess each package's relative ability to taxonomically assign reads, and results from all three tools are presented. Kraken assigned the greatest percentage of reads at the species level. This, coupled with the previously reported possibility of high levels of false positives resulting from Kaiju assignment (27) and the fact that MetaPhlAn2 works off only a subset of marker genes per species (28), is why Kraken was preferentially employed, with a filter threshold of 0.2 to increase precision without detrimentally impacting sensitivity. Furthermore, to reduce the possibility of false positives (27), taxa were included only if present at a minimum of 1% relative abundance in at least one sample; otherwise, reads were assigned as “others.”
B. cereus was found to be the dominant species in 7 of the 12 monthly mesophilic sporeformer-enriched samples, i.e., those from January, February, March, May, July, October, and November. Bacillus thuringiensis, a member of B. cereus sensu lato, was the next most abundant species in these samples (Fig. 1). Among the other mesophilic sporeformer-enriched samples, B. licheniformis was present as the dominant species in April and August. It was also present in the September and December samples but was not dominant. B. subtilis was dominant in the September sample but was also detected at lower relative abundances in August and December. B. paralicheniformis was the dominant species in the December sample and was also present, but not dominant, in April, August, and September (Fig. 1). The June sample was dominated by Brevibacillus brevis and Streptococcus thermophilus. S. thermophilus was also detected at lower levels in other samples, despite not being a sporeformer. All species identified as dominant were confirmed by Kaiju- and MetaPhlAn2-derived results (see Fig. S1 and S2 in the supplemental material).
Relative abundance of the most abundant species (shotgun metagenomic reads assigned by Kraken with a 0.2 filter threshold and minimum of 1% relative abundance in at least one sample) in samples enriched for mesophilic sporeformers.
Spearman correlations were utilized to evaluate the population relationships in samples with taxonomic analysis and functional gene analysis. Spearman correlations are used for nonnormally distributed data that are either skewed or ordinal. Statistical analysis was performed on correlations, and asymptotic P values were generated using rcorr in the R package Hmisc and corrected with Benjamini-Hochberg corrections for multiple comparisons. Significant negative correlations (P < 0.05) between B. cereus and B. licheniformis and between B. cereus and B. paralicheniformis were noted. Significant positive correlations (P < 0.05) were seen between B. licheniformis and B. paralicheniformis, Bacteroides vulgatus and Bacteroides fragilis, B. vulgatus and others, and B. fragilis and others. As expected, many sporulation-associated functional gene groups were significantly positively correlated (P < 0.05), including those for dormancy and sporulation with those for virulence disease/defense, regulation/cell signaling, iron acquisition/metabolism, and cell wall/capsule formation (see Fig. S3 in the supplemental material).
FASTA reads were aligned to the Bagel 3 bacteriocin database using DIAMOND to determine if bacteriocin production potential could be influencing current and future relationship dynamics. This showed a number of potential type I, type II, and type III bacteriocin genes in each sample, with the highest relative abundance per sample going to type I bacteriocin genes (Fig. 2). Bacteriocin genes of note included two lichenicidin peptide genes present in all samples containing B. licheniformis, genes for sporulation killing factor (skfA) and subtilosin seen in samples containing B. subtilis, and thuricin genes in samples containing B. thuringiensis (Fig. 2A).
Percentage of total reads per sample attributed to bacteriocin genes found in the Bagel 3 database. (A) Percentage of reads attributed to type I bacteriocin genes. (B) Percentage of reads attributed to type II bacteriocin genes. (C) Percentage of reads attributed to type III bacteriocin genes.
The August sample was selected for a parallel culture-based analysis and was enriched for culturable mesophilic sporeformers by culturing on BHI agar and incubating at 30°C for 48 h, following 48 h enrichment of a spore-pasteurized sample at 30°C. This was labeled August plate-cultured mesophilic sporeformer-enriched sample (AMP). An additional aliquot of the same spore-pasteurized sample was enriched for thermophilic sporeformers by incubating at 55°C for 48 h; this was labeled August thermophilic sporeformer-enriched sample (AT). Subsequently, this sample was plated on BHI agar and incubated at 55°C for 48 h. This was labeled August plate-cultured thermophilic sporeformer-enriched sample (ATP). The original August mesophilic sporeformer-enriched sample was relabeled AM for the purpose of comparison. In the case of the plate-cultured samples (AMP and ATP), colonies were scraped off, and DNA was extracted in all instances and sequenced as previously described. Comparative analysis revealed that the AMP and AM samples had very similar profiles, highlighting that plating did not bias results (Fig. 3). However, the taxonomic profiles of the AT and ATP samples were very different; i.e., the dominant species present in the AT sample was Thermoanaerobacterium thermosaccharolyticum, whereas the equivalent postplating sample (ATP) was comprised of B. brevis and Bacteroides species, as well as a variety of other species at low abundance (Fig. 3).
Relative abundance of the most abundant species (assigned by Kraken with a 0.2 filter threshold and minimum 1% relative abundance in at least one sample) in the August sample subjected to different enrichment temperatures and conditions.
Beta diversity analysis highlights that samples clustered according to the dominant species present.Bray-Curtis beta diversity analysis of mesophilic sporeformer reads showed that the 12 samples clustered into 3 distinct groups (individual data points within each group being colocated; Fig. 4A). One cluster consisted of samples that contained B. cereus (i.e., January, February, March, May, July, October, and November); a second cluster consisted of samples from the months April, August, September, and December, which all contained a high relative abundance of B. licheniformis; and a third cluster corresponded to the June sample, which had a high relative abundance of B. brevis (Fig. 4A).
Alpha and beta diversity metrics for mesophilic sporeformer-enriched samples and August samples subject to different incubation temperatures and conditions. (A) Bray-Curtis weighted classical multidimensional scaling (MDS) of mesophilic sporeformer-enriched samples. (B) Simpson alpha diversity of mesophilic sporeformer-enriched samples. (C) Bray-Curtis weighted classical multidimensional scaling (MDS) of August samples. (D) Simpson alpha diversity of August samples under different conditions.
Alpha diversity analysis did not reveal any notable pattern other than the observation that the June, September, and December samples had the highest diversities (Fig. 4B). As might be expected from taxonomic results, beta diversity analysis of the August sample, which was enriched in a variety of ways, showed that the AM and AMP samples were least dissimilar, whereas the AT and ATP samples differed from each other and from the AM and AMP samples (Fig. 4C). Among these samples, alpha diversity in the ATP sample was higher than those of all others (Fig. 4D).
Toxin gene analysis revealed the presence of potentially toxigenic B. cereus.Further analysis was performed to determine if the B. cereus strains detected in some samples might be capable of causing emetic and enterotoxic food poisoning. For this purpose, toxin genes were aligned to sequence reads and verified by aligning them to assembled contigs. The genes used and the genome sequences from which they originated are shown in Table 1. Of the 11 toxin genes screened for among all 12 mesophilic enriched samples, the two nonhemolytic enterotoxin (Nhe) genes nhe L1 and nhe L2 were detected in all 7 samples previously shown to contain high relative abundances of B. cereus. The cytotoxin K-encoding gene was detected in 3 samples (from March, October, and November). The presence of toxin genes in the sequence reads was verified by alignment and visualization using Mauve. Toxin gene sequences were aligned to sample contigs, which were assembled using the IDBA-UD algorithm (see Fig. S4 in the supplemental material). Among samples in which all 3 toxin genes were present, toxin genes accounted for close to 0.1% of reads, and when just the two nhe genes were present, they accounted for close to 0.06% of reads (Fig. 5).
Toxin genes aligned to sequence reads
Percentage of total reads per sample attributed to B. cereus toxin genes.
Strain-level analysis revealed the absence of evidence for persistent contamination.As the samples that contained high relative abundances of B. cereus appeared to have different toxin profiles, it was decided to carry out PanPhlAn-based strain analysis to more accurately determine if one strain was dominant across these samples, potentially indicating strain persistence in the processing facility or transient colonization by different strains. Sequence reads for all samples were aligned to a pangenome created from 32 complete B. cereus genomes downloaded from the NCBI database. This analysis established that the presence/absence of B. cereus was consistent with that in previous Kraken analysis. PanPhlAn showed that the 7 samples containing B. cereus contained reads that clustered with 5 strains of B. cereus by Euclidean distance, as visualized using GraPhlAn (Fig. 6). The percentage match from sequence reads to the 32 genomes used in the pangenome was also examined (Table 2). Overall, as shown in Table 2, none of the strains identified in the 7 samples are a complete match to any of the 32 strains with which they were compared. Similarly, none of the strains identified in the 7 samples appear to be exactly the same, with the January and February samples being most alike (Fig. 6).
GraPhlAn visualization of PanPhlAn clustering of 7 samples containing B. cereus with representative strains of the pangenome. The 7 samples contain reads that cluster with at least 5 distinct strains.
Percentage match of B. cereus strains present in samples of 32 reference strains
DISCUSSION
Traditionally, the detection and identification of bacterial sporeformers have involved culturing under different temperatures and conditions (17, 18), together with the use of selective agars to identify pathogenic species (19, 20), followed by further biochemical or molecular approaches for confirmation (1). This study highlights the ability of next-generation shotgun sequencing of enriched samples to identify low-level persistent or transient spore contamination in a dairy powder. Although many of the species identified are similar to those identified in previous studies using traditional detection methods (17), the approach taken here has the potential to rapidly identify these species while allowing strain-level analysis and, as a result, the tracking of persistent microbes in products. It highlights that extremely low levels of potentially pathogenic bacteria can be present, although in this instance, these bacteria are unlikely to be from persistent contamination in the processing facility but rather are thought to represent a transient low-level contamination. Due to the sensitivity of this approach, new guidelines and standards may need to be introduced to ensure that the risks associated with detection of low-level contamination of B. cereus in powders are adequately reflected.
The combined approach of high-temperature treatment coupled with shotgun sequencing employed in this study was selected for a number of reasons. Prior to high-temperature activation, the spore-forming bacteria present in dairy powders are likely to be in a spore form, from which DNA extraction can be difficult (29). Furthermore, in the absence of enrichment, the possibly high level of DNA remaining from deceased bacteria and the presence of host bovine cells in the powders could also be an issue, as DNA from these sources will also be sequenced when untargeted shotgun metagenomic sequencing is employed (22), thereby necessitating the use of additional steps and expensive kits and reagents to deplete DNA from these other sources (30–34). DNA amplification would likely be necessary thereafter due to the low yield of DNA from the low level of sporeformers present in samples. Amplification methods add extra bias, costs, and potential for contamination (34–36). A more targeted, economical approach was used to focus on sporeformers that could potentially germinate in rehydrated powder if given favorable temperatures, without the extra bias and cost of kits. It should be noted that low-speed centrifugation utilized in the approach from which solids were discarded prior to DNA extraction could result in loss of some spore-forming bacteria that have a high affinity for denatured milk protein solids.
This targeted, culture-independent approach was compared with corresponding culture-dependent approaches, using the August sample as a representative test case. The thermophilic sporeformer-enriched August sample (AT) from which DNA was extracted was also subject to culturing on BHI agar at thermophilic temperatures (ATP) before pooling colonies and extracting DNA. DNA from both AT and ATP samples was sequenced (see Fig. S5 in the supplemental material) and compared.
The results were different in that the dominant species identified following enrichment (i.e., without culturing on BHI agar) was T. thermosaccharolyticum, while this species was not detected in the cultured sample, potentially showing bias due to the agar and conditions used. T. thermosaccharolyticum has previously been detected in dried vegetables (37), and more recently, Thermoanaerobacterium spp. have been detected in the core microbiome of raw milk (38). T. thermosaccharolyticum is a thermophilic, anaerobic sporeformer and was previously classified as a member of the Clostridium genus, although it was subsequently reclassified (39). It is a known canned-food spoilage organism (40) that produces hydrogen and causes swelling in canned foods. This ability to efficiently produce hydrogen makes it a potentially important organism for sustainable biohydrogen production (41, 42). The previous nondetection of the species in dairy powders may relate to an inability to grow on the agar substrates conventionally employed by the dairy industry.
A decision was made to further investigate the potential B. cereus sequences identified to eliminate the potential that the taxonomic classifier was misassigning other members of B. cereus sensu lato and B. cereus sensu stricto, as members of this group are notoriously difficult to differentiate (43). Toxin profiling and PanPhlAn analysis confirmed that the strains had the potential to be pathogenic, while also highlighting differences between them. All except one of the samples that contained B. cereus showed complete alignment with both Nhe toxin-encoding genes. The exception was the February sample, from which one-half of the nhe L1 gene was absent. This may be due to misassembly or poor coverage of the genome, or it may reflect a natural mutation in this strain (Fig. S4). It should be noted that alignment to toxin genes infers potential ability to produce toxins but does not conclusively indicate functional toxin presence. Also, a total of 105 to 109 CFU of toxin-producing B. cereus is needed to cause food poisoning (15). The highly sensitive, qualitative nature of shotgun metagenomic analysis suggests that further steps will need to be taken to determine the risk associated with the detection of these and other toxin genes in food samples. PanPhlAn analysis showed that the 7 samples containing B. cereus contained reads that clustered with at least 5 strains of B. cereus by Euclidean distance, suggesting that the strains identified in the samples are different and not as a result of persistent contamination in the powder production facility or elsewhere.
The negative correlation between B. cereus and B. licheniformis/B. paralicheniformis raises the possibility of competition to determine the dominant species and the possible inhibition of the other species. The potential bacteriocin genes detected in all samples may impact the relationship dynamics observed and may have an antagonistic effect on some species currently or in future food products. Lichenicidin is a two-peptide lantibiotic previously shown to be active against pathogenic Gram-positive bacteria, including B. cereus (44), and we noted alignment to two peptides associated with lichenicidin in all four samples containing B. licheniformis. However, potential identification of two peptide genes does not suggest the presence of the whole bacteriocin gene cluster or infer correct posttranslational modifications to produce a functional bacteriocin. Positive correlations among virulence, disease, defense, dormancy, and sporulation highlight the need to be cautious of sporeformers in food products. A greater understanding of their relationships will aid the prevention of spoilage and pathogenic species that cause concern in food processing.
Disadvantages that need to be overcome in order to allow for the routine use of the sequencing technologies employed in this study primarily relate to the cost of analysis, which is currently too expensive for large-scale routine use. Additionally, there are challenges relating to assembly of genomes from shotgun metagenomic sequencing (22) and difficulties arising from insufficient accuracies associated, to different extents, with taxonomic classifiers (27). There are some solutions emerging, whereby new lower-cost, rapid sequencers are arriving on the market, with MinION (45) leading the way toward quick portable detection systems for microorganisms. Through the generation of more reference genome sequences and good-quality shotgun metagenomic sequencing, reference databases and the accuracy of results will improve.
This study highlights monthly diverging contamination patterns in whey powder production, which converged into 3 distinct mesophilic sporeformer population groups from 12 powder samples produced in the same production facility, and shows that the way in which the powders are treated postproduction (namely, incubation temperature postreconstitution) influences which bacteria germinate and become dominant. Shotgun metagenomics is a useful tool to delve deeper into the understanding of sporeformers and their relationships in food processing, although it brings with it its own set of caveats and need for guidelines for use and interpretation of results.
MATERIALS AND METHODS
Sample preparation and enrichment.Twelve (i.e., one per month for a calendar year) 5-g whey powder samples from a single supplier were suspended aseptically at a concentration of 10% (wt/vol) in sterile 1/4-strength Ringer's solution (Sigma-Aldrich). Each reconstituted sample was incubated at 80°C for 12 min to select for spore-forming bacteria, as previously described (3, 17). An aliquot was then plated onto brain heart infusion (BHI) agar before incubation at mesophilic (30°C) and thermophilic (55°C) sporeformer-suitable temperatures, and CFU · g−1 values were calculated.
Following spore selection, in which 10% (wt/vol) reconstituted samples were incubated at 80°C for 12 min, reconstituted samples were enriched in a manner consistent with that previously employed to select for low-abundance mesophilic sporeformers, but with incubation at 30°C for 48 h as opposed to at 32°C, as previously documented (3, 17). This is shown in Fig. S5 in the supplemental material.
To facilitate an investigation into the extent to which plating altered the identity of the populations of sporeformers detected, powder sourced in August (A) was treated as described above but with incubation at 55°C to select for thermophilic sporeformers (AT). Both of the enriched August samples were plated onto BHI agar before incubation at the corresponding mesophilic (30°C; AMP) and thermophilic (55°C; ATP) temperatures for 48 h (Fig. S5). DNA was extracted from the surface of these agar plates as described below.
DNA extraction and library preparation.A total of 50 ml of samples that were reconstituted, heat-treated, and enriched by incubation for 48 h at 30°C were centrifuged at 900 × g for 20 min. The resultant pellet was discarded, and the supernatant was centrifuged at 4,500 × g for 20 min. Pellets were washed in 1 ml of sterile 1/4-strength Ringer's solution and centrifuged at 13,000 × g for 2 min. Washing was repeated twice more, centrifuging for 1 min each time at 13,000 × g. DNA was extracted from the resultant pellet as described below. For DNA extraction from a combination of all colonies on the surface of agar plates, 5 ml of 1/4-strength Ringer's solution was added and colonies were washed off using a sterile Lazy-L spreader (Sigma-Aldrich). An aliquot of 4 ml of the resultant fluid was removed and centrifuged at 4,500 × g for 20 min. The resultant pellet was suspended in 1 ml of 1/4-strength Ringer's solution and centrifuged at 13,000 × g for 2 min. The pellet was washed twice more in 1 ml of Ringer's solution and centrifuged at 13,000 × g. Before the final centrifugation, 200 μl was removed from the mix to a sterile tube. This was centrifuged at 13,000 × g for 1 min. This pellet corresponded to 0.8 ml of culture from the original 4 ml.
All pellets were resuspended in 150 μl of 45 mg · ml−1 lysozyme and incubated at 37°C for 30 min. Samples were centrifuged at 13,000 × g for 1 min, and the supernatant was removed. From this point, the MoBio PowerFood DNA isolation kit protocol was followed, including the use of the alternative lysis step for difficult-to-lyse cells. DNA was eluted in 60 μl of 10 mM Tris-HCl and stored at −20°C. DNA was quantified and quality checked using the Qubit double-stranded DNA (dsDNA) high sensitivity assay kit (Bio-Sciences, Dublin, Ireland) and by visualizing on a 1% (wt/vol) agarose gel.
Samples were prepared for metagenomic sequencing according to Illumina Nextera XT library preparation kit guidelines and sequenced on an Illumina MiSeq sequencing platform at the Teagasc DNA Sequencing Facility with a 2 × 250 V2 kit, using standard Illumina sequencing protocols.
Bioinformatics pipeline.Raw metagenomic shotgun reads were checked for the presence of human and bovine reads filtered on the presence of quality and quantity and then trimmed to 170 bp with a combination of Picard tools (http://broadinstitute.github.io/picard/) and SAMtools (46). Kraken with a filter threshold of 0.2 (47) and SUPER-FOCUS (48) were used to determine microbial composition to species level and biological functions, respectively. MetaPhlAn2 (28) and Kaiju (49) were also utilized to determine microbial composition for comparison to Kraken results. Eleven B. cereus toxin-associated genes (Table 1) were downloaded from the NCBI database and aligned to sequence reads using Bowtie2 (50). Metagenomes were assembled into contigs using IDBA-UD (51), and toxin genes were aligned to these using Bowtie2 (50). Toxin genes were also aligned to contigs using Mauve version 20150226 (52)—more specifically, progressiveMauve (53)—using default parameters (default seed weight, determine locally collinear blocks [LCBs], full alignment, iterative refinement, min LCB weight = default, sum-of-pairs LCB scoring). PanPhlAn (54) was utilized to determine B. cereus strain profiles in each sample. A total of 32 complete B. cereus genome sequences were downloaded from the NCBI database and a pangenome was generated and compared to all samples. PanPhlAn outputs were visualized using GraPhlAn (55). Spearman correlations with Benjamini-Hochberg corrections for multiple comparisons were made in R using the Hmisc package and visualized using corrplot. Diversity was analyzed using the R vegan package, and data visualization was performed using ggplot2. Sequence reads were aligned to the Bagel 3 bacteriocin database (56) using DIAMOND (57) to determine the bacteriocin potential of the bacteria identified.
Accession number(s).Sequence data have been deposited in the European Nucleotide Archive (ENA) under the study accession number PRJEB24853.
ACKNOWLEDGMENTS
This research was funded by the Department of Agriculture, Food and the Marine (DAFM) under the FIRM project SACCP, reference number 14/F/883. Research in the Cotter lab is also funded by the SFI award “Obesiobiotics” (11/P1/1137).
Thanks to Orla O'Sullivan, Fiona Crispie, and Laura Finnegan of the Food Bioscience Department in Teagasc Moorepark.
FOOTNOTES
- Received 29 May 2018.
- Accepted 31 July 2018.
- Accepted manuscript posted online 3 August 2018.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.01305-18.
REFERENCES
- Copyright © 2018 American Society for Microbiology.
