High-Throughput 16S rRNA Gene Amplicon Analysis of Raw Milk Microbiota

Milk samples 3 (farm A) and 4 (farm B) subjected to culture-based analysis of 1,500 colonies were also used for amplicon sequencing.

As microbial loads in fresh raw milk are comparatively low and did not exceed log 5 cfu/mL in samples 3 and 4, bacterial cells were concentrated prior to DNA extraction. A volume of 150 mL was centrifuged at 8,000 × g for 20 min at 4°C to pellet bacterial cells. The supernatant consisting of a fat layer and skim milk was carefully removed, except for 10 mL skim milk, which were left to re-suspend the pellet. The suspension was transferred to a 50 mL tube and additionally centrifuged at 5,000 × g for 10 min at 4°C. Again, the fat layer was carefully removed.

Sedimented casein was removed according to Murphy et al. (2002). The cell suspension was divided into 10 aliquots of 1 mL and 300 μL 0.5 M EDTA (pH 8.0) as well as 200 μL TE buffer (pH 7.6) were added to each aliquot. Within approximately 1 min, the casein micelles disintegrated due to chelating of calcium ions indicated by clarification of the solution, which was centrifuged at 16,000 × g for 1 min at room temperature. The supernatant was discarded and the pellet re-suspended in 100 μL Ringer’s solution. All ten subsamples were then pooled and centrifuged again. The supernatant was discarded and the pellet re-suspended in 350 μl quarter strength Ringer’s solution.

The PathoProof DNA Extraction Kit (Thermo Fisher Scientific), including enzymatic lysis of microbial cells, was applied according to the manufacturer’s instructions to extract DNA, producing a final volume of 100 μL DNA extract per sample. As bulk tank milk contains high amounts of eukaryotic DNA originating from somatic cells of the cow, enrichment of bacterial DNA using the Looxster Enrichment Kit (Analytik Jena, Germany) was performed according to the manufacturer’s instructions resulting in 30 μL of DNA for each sample.

The concentration of total DNA was determined using the Qubit^® 2.0 fluorometer (Life Technologies Co.) using the dsDNA HS Assay Kit. In a further step, the concentration of bacterial DNA in both extracts was determined using quantitative real-time PCR. PCR mixtures contained 4 μL Phusion^® Buffer HF (Thermo Fisher Scientific Inc.), 1 μL dNTPs (20 nM), 2 μL each of primer 515F (5′-GTGCCAGCMGCGCGGTAA) and 806R (5′-GGACTACHVGGGTWTCTAAT) (10 pmol/μL), 0.1 μL Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific Inc.), 1 μL SYBR-Green (diluted 1:50,000), and 2 μL of DNA extract (diluted 1:10) in a final volume of 20 μL. After initial denaturation at 98°C for 30 s 40 cycles of 5 s denaturation at 98°C, 10 s primer annealing at 52.5°C, and 10 s elongation at 72°C were run. A DNA standard composed of somatic DNA extracted from bulk tank milk and different fractions (10, 1, and 0.1%) of a bacterial DNA extract of a Pseudomonas and Enterococcus raw milk isolate was used for quantification of bacterial DNA.

DNA quantification resulted in a total DNA concentration of 15 ng/μL for sample 3 of farm A and 10 ng/μL for sample 4 of farm B. Sample 3 contained 4.8% (0.7 ng/μL) bacterial DNA, sample 4 2.5% (0.25 ng/μL). DNA extracts were then standardized to contain 0.2 ng/μL bacterial DNA. Library preparation followed a two-step protocol (Berry et al., 2011) based on amplification of the V3–V4 region of the 16S rRNA gene using primers 341F (5′-CCTACGGGNGGCWGCAG) and 785R (5′-GGATTAGATACCCBDGTAGTC) (Klindworth et al., 2013). Due to the low concentrations of template microbial DNA, the PCR protocol was adjusted and the number of replicates as well as the number of cycles augmented. For the first PCR step, eight parallel PCRs (four barcodes in duplicate) were run for each sample, each reaction using 3.5 μL DNA extract containing 0.7 ng bacterial DNA. After initial denaturation at 98°C for 30 s, 30 cycles of 5 s denaturation at 98°C, 10 s primer annealing at 55°C, and 10 s elongation at 72°C were run. The second PCR step added barcodes (dual combinatorial indexing) and Illumina adaptors to the amplified fragments and was done using the same protocol but 2 μl of PCR product of step 1 as template and was run only for 10 additional cycles. A negative control (PCR blank) was included using PCR-grade water as template.

Duplicate PCRs for any given samples (i.e., having the same barcodes) were transferred into a 1.5 mL tube and purified using Agencourt AMPure XP Beads (Beckman Coulter Inc.). DNA concentrations were determined via Qubit^TM 2.0 fluorometer and adjusted to 2 nM. All libraries were pooled and sequencing was performed in paired-end mode (2 × 275 cycles) using a MiSeq platform (Illumina).

Raw reads were processed using the IMNGS pipeline (Lagkouvardos et al., 2016) based on UPARSE (Edgar, 2013). After demultiplexing, forward and reverse reads were merged and trimmed by five nucleotides on each end. Chimera filtering was done using UCHIME (Edgar et al., 2011). Operational Taxonomic Units (OTU) clustering (97% identity) was performed by USEARCH 8.0 (Edgar, 2010) and OTUs occurring at a relative abundance <0.05% in all samples were discarded. Taxonomical identification was assigned using the RDP classifier (Wang et al., 2007) and diversity analyses were done in Rhea (Lagkouvardos et al., 2017). Relative abundances of OTUs were normalized to account for differences in sequence depth and α-diversity was assessed on the basis of species richness as well as Shannon and Simpson diversity indices. As automated identification may be imprecise and often terminates at higher ranks, 487 OTU sequences not assigned to a genus were manually identified at the genus level (95% similarity) using EzBiocloud (Yoon et al., 2017). OTUs for manual identification were selected according to the following criteria: all OTUs that were identified only to the phylum level or above, all OTUs belonging to families detected by culturing, and the 100 most abundant OTUs in each sample. ß-diversity was computed for both the culture-dependent and -independent results based on generalized UniFrac distances (Chen et al., 2012) and visualized in a phylogram based on the Ward’s minimum variance method (Murtagh and Legendre, 2014).