Strains employed to subsequently create a model community were obtained from the Teagasc Dairy Production Centre (DPC) culture collection. These strains were Bacillus velezensis DPC 3313, Escherichia coli DPC 6912, Lactococcus lactis DPC 4268, Pseudomonas proteolytica DPC 6056 and Staphylococcus haemolyticus DPC 5987 and represent bacteria found in the bovine milk microbiome (Quigley et al., 2013). Strains were grown on tryptic soy agar (TSA; Thermo Fisher, Ireland) overnight at 30°C (B. velezensis and P. proteolytica) or 37°C (E. coli, L. lactis and S. haemolyticus) to generate fresh cultures. Cells were diluted in 10 ml phosphate-buffered saline (PBS; Thermo Fisher) to obtain a bacterial suspension of a density equivalent to approximately 10⁷ cells/mL for each strain. Equal volumes of strains were combined to generate a 5-strain model community for spiking into milk samples. This mixture was used as a live model community and a dead model community was generated through heat inactivation at 95°C for 15 min (resulting in an absence of subsequent growth on TSA).

Store-bought UHT skim milk was used for milk spiking studies in order to ensure a low load of background microorganisms. Two different viability conditions were used for spiking; live (consisting of 1 ml of the live-cell model community) and dead (consisting of 1 ml of the dead-cell model community). Cells were spiked into 10 ml of milk in triplicate and left to incubate at room temperature overnight. Unspiked milk samples and both the live and dead bacterial model communities, which did not undergo any treatment and were used for spiking, were used as controls. Milk samples were then centrifuged at 4,500 × g for 15 min at 4°C, the supernatant was discarded, and the cells pellets were subjected to two washing steps whereby the pellets were resuspended in sterile PBS and centrifuged at 13,000 × g for 1 min, after which the supernatant was discarded. Samples used for Shotgun and PMA-Shotgun metagenomic analysis were resuspended in 1 ml sterile PBS while, for samples used for metatranscriptomic- and 16S rRNA-based analysis, pellets were resuspended in 500 μl RNALater and stored at -80°C before extraction.

Each 1 ml sample was divided into 2 × 500 μl samples, one of which was subjected to a PMA (PMAxx dye, Biotium, CA, United States) treatment and the other without PMA treatment. Preliminary investigations were performed to optimise the final PMA concentration, incubation and light exposure time with the quantification of total bacteria performed using the Femto Bacterial DNA Quantification kit (Zymo Research, CA, USA). The final parameters for PMA treatment involved the use of a final PMA concentration of 20 μM and incubation in the dark at room temperature and light exposure with the PMA-Lite™ LED photolysis device (Biotium, USA) for 30 min. Following PMA treatment, both PMA-treated and PMA-free samples were subjected to DNA extraction using the MolYsis complete5 kit (Molzym GmBH & Co. KG, Bremen, Germany), with 50 μl of DNA eluted for downstream sequencing. The MolYsis kit was used as it had been found to significantly improve microbial sequencing depth for milk samples (Yap et al., 2020). gDNA was quantified using the Qubit dsDNA HS assay kit (Invitrogen) and stored at at -20°C before library preparation. Total bacteria levels were again quantified using the Femto Bacterial DNA Quantification kit (Zymo Research).

A total of 18 samples (6 PMA-treated, 6 non-PMA-treated, 6 controls) were prepared for shotgun metagenomic sequencing according to Illumina Nextera XT library preparation kit guidelines, with the use of unique dual indexes for multiplexing with the Nextera XT index kit (Illumina). Following indexing and clean up, samples were pooled to equimolar concentration of 1 nM. Samples were sequenced on an Illumina NextSeq 500 sequencing platform with a V2 kit, at the Teagasc DNA Sequencing Facility, using standard Illumina sequencing protocols.

Whole metagenome amplification was performed on samples using multiple displacement amplification (MDA) with the REPLI-g UltraFast Mini kit (Qiagen, West Sussex, United Kingdom), according to the manufacturer’s instructions. After amplification, DNA was quantified using the Qubit dsDNA BR assay kit (Invitrogen) and normalised to 400 ng in 7.5 μl. Two libraries were prepared, with each library containing 9 samples (6 samples (PMA-treated or non-PMA-treated) and 3 controls) per flowcell. The Rapid Barcoding kit (SQK-RBK004, Oxford Nanopore Technologies, United Kingdom) was used to prepare the libraries. For this, DNA was tagmented and barcodes were attached to fragments before pooling, followed by clean-up using Ampure XP beads (Beckman Coulter) to concentrate the pooled library directly prior to sequencing. DNA was sequenced using a GridION (release 19.12.6) with single flow cells for each library (R 9.4.1) that were primed prior to loading libraries with MinKNOW (core 3.6.5) and integrated basecalling by Guppy (3.2.10).

Samples to be used for 16S rRNA and metatranscriptomic-based analysis were subjected to RNA extraction using a TRIzol chloroform protocol with on-column DNase purification with the PureLink RNA Mini kit (ThermoFisher Scientific). Briefly, samples suspended in RNALater were first centrifuged at maximum speed (17,000 × g) for 1 min at 4°C, after which the supernatant was discarded. 1 ml of TRIzol reagent was added to each sample and mixed through pipetting up and down to disperse the pellet, followed by incubation for 10 min at 60°C. To each sample, 200 μl of chloroform was added and the tubes were shaken vigorously by hand for 15 s before incubation at room temperature for 3 min. The samples were then centrifuged at 12,000 × g for 15 min at 4°C, during which the mixture separates into a lower, red phenol-chloroform phase, an interphase, and a colourless upper aqueous phase which contains the RNA. Approximately 400 μl of upper aqueous phase of RNA was transferred to a RNase-free tube where an equal volume of ice-cold 100% ethanol was added for precipitation and vortexed for even mixing. The samples were then purified through on-column DNase treatment following the manufacturer’s protocol with final elution at 30 μl in RNase-free water. RNA concentration was quantified using the Qubit RNA HS assay kit (Invitrogen) and quality checked using an Agilent Bioanalyzer with the RNA 6000 pico assay kit (Agilent). RNA was stored at -80°C for use in downstream library preparation steps.

For 16S methods, 20 μl of RNA for each of the 9 samples (6 samples and 3 controls) were used for cDNA synthesis with the Superscript IV First Strand Synthesis System (Invitrogen) according to the manufacturer’s instructions, using 2 μM of the reverse primer of the V3-V4 region of the 16S rRNA gene that are commonly used in Illumina 16S sequencing. The cDNA generated was stored at -20°C overnight before use in library preparation.

Five microliters of the cDNA generated was used in library preparation according to Illumina 16S Metagenomic Sequencing Library Preparation guidelines, with the use of dual indexes for multiplexing with the Nextera XT index kit (Illumina). Following indexing and clean up, samples were pooled to equimolar concentration of 20 nM. Samples were sequenced on an Illumina MiSeq sequencing platform with a V3 kit, at the Teagasc DNA Sequencing Facility, using Illumina sequencing protocols.

The cDNA generated was normalised to 10 ng in 10 μl for use in library preparation with the 16S Barcoding kit (SQK-RAB204, Oxford Nanopore Technologies) according to the manufacturer’s instructions. Briefly, amplification of the 16S gene was done using barcodes, followed by library clean up using Ampure XP beads (Beckman Coulter) and pooling to 10 ng/μL. Sequencing adapters were attached before the loading of libraries on a single primed flow cell (R 9.4.1) and sequenced on the GridION (release 19.12.6) with MinKNOW (core 3.6.5) and integrated basecalling by Guppy (3.2.10).

For metatranscriptomics, rRNA depletion was carried out using the QIAseq FastSelect–5S/16S/23S kit (Qiagen) according to the manufacturer’s instructions with 20 μl of each of the 9 samples (6 samples and 3 controls). The first step of RNA fragmentation was performed at 89°C for 7 min. After bead clean up, rRNA-depleted RNA was stored at -80°C before use in library preparation.

First strand cDNA synthesis was performed using the NEBNext Ultra II Directional RNA Library Prep Kit (Brennan & Co., Dublin, Ireland) with 10 μl of rRNA-depleted RNA. Library preparation was performed according to the manufacturer’s instructions. Seven cycles were used for PCR amplification during indexing of adapter ligated DNA. The quality and quantity were measured using a high sensitivity DNA chip on the bionanalyzer (Agilent) and Qubit HS dsDNA kit (Invitrogen), respectively. An additional clean up step was done due to the presence of primer and adapter dimers before the quality and quantity was measured. Samples were pooled to equimolar concentration of 4 nM before sequencing on the Illumina NextSeq 500 using a mid-output (2 × 75 bp kit) at the Teagasc DNA Sequencing Facility, using Illumina sequencing protocols.

rRNA-depleted RNA was normalised to 100 ng in 7.5 μl for use in library preparation with the Direct cDNA Native Barcoding kit (SQK-DCS109, Oxford Nanopore Technologies) with the native barcoding expansion kit (EXP-NBD104, Oxford Nanopore Technologies) according to the manufacturer’s instructions. Briefly, reverse transcription and strand switching was done to prepare full-length cDNAs followed by barcode ligation. Barcoded samples were pooled prior to adapter ligation before a final clean up with Ampure XP beads (Beckman Coulter) before loading the library onto a single primed flow cell (R 9.4.1) and sequenced on a GridION (release 19.12.6) with MinKNOW (core 3.6.5) and integrated basecalling by Guppy (3.2.10).

Quality checks and adapter trimming for Illumina methods (shotgun, PMA-shotgun and metatranscriptomics) were done with FastQC (v. 0.11.8; Andrews, 2010) and cutadapt (v. 2.6; Martin, 2011) and host reads were aligned to the bovine genome (Bos taurus) and removed with Bowtie2 (v. 2.4.4; (Langmead and Salzberg, 2012). Quality checks for ONT methods (shotgun, PMA-shotgun and metatranscriptomics) were done with MinIONQC (Lanfear et al., 2019) and Nanoplot (v. 1.28.2; De Coster et al., 2018), with adapter trimming done with porechop (v. 0.2.4). Host reads were aligned with Minimap2 (v. 2.17-r974) to the bovine genome (Bos taurus) (Li, 2018). Sortmerna (v. 2.1b) was used for rRNA removal prior to the analysis of metatranscriptomic data (Kopylova et al., 2012). Taxonomic classification for all metatranscriptomics, shotgun and PMA-shotgun methods was performed with Kraken2 (2.0.7) using the Genome Taxonomy Database (release 89) which contains Bacteria and Archaea (Parks et al., 2018, 2020, 2022; Wood et al., 2019). For the Illumina 16S rRNA method, quality control, adapter trimming and denoising was done with Qiime2 (2021.2) with cutadapt and DADA2 (Martin, 2011; Callahan et al., 2016; Bolyen et al., 2019). For the ONT 16S rRNA method, analysis was performed using the q2ONT pipeline,¹ which uses porechop (v. 0.2.4) to demultiplex reads and trimmomatic (0.38) to discard reads shorter than 1,400 bp and crop reads to that length. Qiime2 (2021.2) was used for the remaining steps, where sequences were dereplicated and chimeric sequences were filtered out with vsearch at 85% identity before reads were aligned with mafft and highly variable positions were masked and filtered out. For both 16S rRNA methods, taxonomy was assigned with the Qiime2 trained classifier using the Silva 138 database (Quast et al., 2012).

Performance metrics were calculated in R (4.1.2; R Core Team, 2015), with the relative abundances of a model community strain correctly (TP) and incorrectly (FP) quantified by the method, along with the relative abundances of other taxa present in the sample correctly (TN) and incorrectly (FN) quantified. The following metrics were used:

Accuracy is the correct quantification of model community and other taxa across all observations, A = TP + TN/TP + TN + FP + FN

Diversity analysis was done with the vegan package (Oksanen et al., 2015), with beta diversity calculated as Bray-Curtis metrics, visualised in a principal coordinate analysis plot. The “adonis” function from the vegan package was used to calculate the permutational analysis of variance (PERMANOVA) to determine differences in composition of the community between groups of samples (number of permutations = 999). Hierarchical clustering was done with dendextend (Galili, 2015) using Bray-Curtis distances. Data was cleaned, analysed and visualised in R with ggplot2, tidyverse and ggpubr packages (Wickham, 2011; Wickham et al., 2019; Kassambara, 2020).

In general, the model community taxonomy was determined with one classifier to ensure consistency across methods. This classifier, Kraken2, was previously found to efficiently characterise the milk microbiome when compared to other classifiers (Yap et al., 2020). In addition, a trained classifier in Qiime2 was used for the analysis of 16S rRNA data (Illumina-sequenced 16S rRNA: i-16 s and ONT-sequenced 16S rRNA: o-16 s). Comparing across methods, the taxonomic identity of the model community was consistent to genus level in all cases (Table 1). This was also the case at species level, where assignment at species level was possible, with the exception that reads generated from Illumina sequencing were assigned as Staphylococcus hominis, while the same model community strain was classified as a Staphylococcus haemolyticus when ONT data was analysed. However, this is not unexpected as S. hominis is classified under the S. haemolyticus group of bacteria (Becker et al., 2014). ONT-generated metatranscriptomic data and the short read length 16S rRNA data (regardless of sequencing platform) could not be resolved beyond genus level. Indeed, Staphylococcus was not detected from within the ONT-derived metatranscriptomic data, most likely due to insufficient sequencing depth combined with low levels of gene expression by the strain (Supplementary Table S1). For downstream analysis, genus level data was used.

Total bacteria qPCR was performed to quantify the bacteria in each sample and the log number of cells/mL of each microorganism per sample was estimated using a combination of qPCR results and relative abundance values. The abundances of cells recovered from Live samples were greater than those in the heat-treated samples (Supplementary Figure S1). Some reads were classified into genera other than those corresponding to the 5 spiked strains (Supplementary Table S2), and some of these were also detected in the unspiked milk samples (Supplementary Table S3). The reads classified as others were not prioritised in the further analysis due to the particular focus on differentiating between the members of the model community in the Live and Dead samples.

To evaluate the performance of each method, several performance metrics such as accuracy, precision, sensitivity and F-score were calculated. The comparisons were made between the abundances of the live and dead model community controls that had been sequenced and analysed with the various methods. Overall, most approaches showed a high level of accuracy, which is the fraction of correctly identified taxa (model community or other) relative to all observations, with a median of 0.939 ± 0.046 (Figure 2). However, the F-score, which gives a sense of how precise and sensitive the method is, revealed that only half of the methods scored at least 0.5, indicating issues in precision or sensitivity in these instances. The two PMA-shotgun sequencing and analysis methods yielded outputs that were most precise, while the 16S rRNA methods were the most sensitive, including when separated into both live and dead spiked samples (Figure 2, Supplementary Figure S2). ONT-sequenced 16S rRNA (o-16 s) was the most accurate and had the best overall F-score, while metatranscriptomics (i-metaT) and PMA-shotgun (i-pma) had the highest accuracy and F-score out of all the Illumina-based methods, respectively.

Beta diversity analysis was employed to determine the differences between the outputs generated from the different approaches taken, and several sample clusters were apparent (Figure 3A). Samples first clustered by library type, whereby the outputs from DNA-based analysis (Shotgun and PMA-shotgun) were significantly distinct from those generated using the RNA-based approaches (metatranscriptomics and 16S rRNA; PERMANOVA, R² = 0.457, p < 0.01), and also by sequencing methods, with Illumina-sequenced and ONT-sequenced outputs clustering apart (PERMANOVA, R² = 0.195, p < 0.01). Differences in nucleic acid input and sequencing methods used were also apparent in the hierarchical clustering of the Bray-Curtis distances (Figure 3B). The composition of taxa varied across library types, with Escherichia and Lactococcus being found at greater abundance when DNA-based methods were used, while Pseudomonas was the most active on the basis of RNA-based methods (Figure 3C). In particular, the Illumina-sequenced 16S rRNA data points clustered away from those generated using other methods. This was also obvious in the taxonomic composition where abundances of other taxa besides the model community was lowest in i-16S compared to other methods (Figure 3C).

Among the DNA-generated dataset, the effect of PMA was clear, with lower DNA and overall calculated cell counts observed for samples that underwent PMA treatment (Figure 4). Notably, although the Shotgun methods do not exclusively select for viable cells, the accuracy of those methods were comparable to the PMA-shotgun methods (Figure 2). Precision was higher in PMA-shotgun methods compared to the Shotgun methods, while sensitivity was higher for i-pma and lower for o-pma when compared to i-shotgun and o-shotgun, respectively. In terms of the composition, most of the abundances of model community taxa in PMA-shotgun samples were lower in total cells/mL than those in Shotgun samples, apart from Escherichia cells/mL, which was higher in PMA-shotgun methods than Shotgun for both live and dead spiked samples (Figure 3C). The cell counts of Bacillus and Pseudomonas were generally lower, with the exception of i-shotgun where both taxa were found in abundances greater than the other DNA-based methods. Lactococcus was the most abundant of the model community taxa quantified, which was consistent across the two methods, with exception of o-shotgun, which detected higher abundances of Staphylococcus. Staphylococcus was clearly detected in the Shotgun samples spiked with the dead cells whereas levels were negligible in the PMA-shotgun samples. Between sequencing methods, distinct separation was found between Illumina and ONT methods (Figures 3A,B) and, in terms of composition, ONT methods detected higher abundances of other non-model community taxa, compared to the Illumina methods (Figure 3C).

Although there were no obvious differences in the accuracy of the methods across the different RNA-based approaches, the 16S rRNA method outperformed the metatranscriptomics approaches with respect to precision and sensitivity (Figure 2). In general, methods were more sensitive than precise, except for o-metaT which had poor sensitivity. Staphylococcus reads were detected in low abundances across all RNA-based methods, while read numbers corresponding to Pseudomonas were consistently found at high abundance (Figure 3C). Similar to the DNA-based methods, clear separation between library types was seen between Illumina- and ONT-based methods (Figure 3). As mentioned above, Illumina-sequenced 16S rRNA was the most distinct, clustering separately from other methods and with the lowest proportion of other taxa identified compared to the other RNA-based methods (Figures 3A,B). It was the only method that had detectable levels of Bacillus in live samples and its samples spiked with the dead-cell model community had a taxonomic profile that differed from the other samples spiked with dead cell concentrations (Figure 3C).