A detailed explanation of the ExStream system can be found in Mack et al. (2022). In summary, stream mesocosms were connected to the adjacent stream Bieber, which provided them with a constant water flow. The stream Bieber is part of the Rhine-Main-Observatory,¹ a Long-Term Ecological Research site in Germany (Haase et al., 2016; Mirtl et al., 2018). Each mesocosm was set up using substrate and organisms from the stream. A random subset of the mesocosms was exposed to either the insecticide chlorantraniliprole (Coragen, DuPont), increased fine sediment concentration, or both. Both insecticides and fine sediment are known key stressors of aquatic environments introduced into streams by agricultural runoff. The stressors were induced using a 4×2 factorial design by adding 0.2 μg/L, 2 μg/L, and 20 μg/L (acute stressor phase, 4 days) or 0.02 μg/L, 0.2 μg/L, and 2 μg/L (reduced stressor phase, 17 days) of the insecticide and 450 mL of fine sediment (<2 mm) to the mesocosms. Each possible combination of stressor levels was replicated eight times in addition to eight control mesocosms that did not receive any stressor, resulting in 64 mesocosms.

The goal of the ExStream experiment was to evaluate the individual and combined effects of the applied stressors on biodiversity and organic matter decomposition in streams. To investigate organic matter decomposition, cotton strips were added to all mesocosms. Cotton strips are mainly made of cellulose, which is a major source of carbon in stream ecosystems. Therefore, analyzing the biofilm on the cotton strips allowed the analysis of the diversity of microbial communities degrading organic matter.

The experiment was divided into a colonization phase (days −21 to −1) and a stressor phase (days 0 to 21). Two cotton strips were added to each of the 64 mesocosms on day −17 (128 in total) and recovered after 28 or 35 days, respectively for more information on the phases and cotton strip addition and recovery see Mack et al. (2022). Four cotton strips were washed away during the experiment, so 124 cotton strips were recovered in total. A 2-cm-long piece of each cotton strip was cut off and transferred into a ZR BashingBead Lysis Tube (0.1 & 0.5 mm) pre-filled with 1 mL of DNA/RNAShield (Zymo Research, Freiburg, Germany) using sterile laboratory gloves, forceps, and scissors. The samples were transferred to a laboratory, stored at −20°C, and then homogenized using a bead mill homogenizer (MM 400, Retsch, Haan, Germany) at 1,800 rpm for 30 min. 300 μL of each lysate were processed for amplicon sequencing at the University of Duisburg-Essen, Germany, and the remainder of each lysate was shipped to the University of Guelph, Canada, on dry ice and processed for metagenomics and total RNA-Seq.

Amplicon sequencing was carried out following the workflow described by Buchner et al. (2021). All subsequent processing steps were completed on a Biomek FX^P liquid handling workstation (Beckman Coulter, Brea, CA, United States). Briefly, replication of the samples was carried out before DNA extraction by transferring 60 μL from the bead-beating tubes to deep-well plates pre-filled with 133 μL of TNES buffer (50 mM Tris, 400 mM NaCl, 100 mM EDTA, 0.5% SDS, pH 7.5) and 6 μL of Proteinase K (10 mg/mL) following incubation for 3 h at 55°C for complete lysis of the samples. DNA was extracted using a modified version of the NucleoMag Tissue kit Macherey Nagel, Düren, Germany; for modifications see Buchner et al. (2021). Extraction success was verified using a 1% agarose gel.

The PCR for the amplicon library was performed using a two-step PCR protocol following Zizka et al. (2019). Samples were amplified in a first-step PCR using the Qiagen Multiplex Plus Kit (Qiagen, Hilden, Germany) with a final concentration of 1x Multiplex Mastermix, 200 mM of each primer [515F & 806R for 16S (Caporaso et al., 2011) and ITS3-CS1 & ITS4-CS2 for ITS-2 (Frey et al., 2016)], and 1 μL of DNA, and filled up to a total volume of 10 μL with PCR-grade water. The amplification protocol was: 5 min of initial denaturation, 25 cycles of 30 s denaturation at 95°C, 90 s of annealing at 50°C for 16S and 55°C for ITS-2, and 30 s of extension at 72°C, finished by a final elongation step of 10 min at 68°C. For subsequent demultiplexing, each of the PCR plates was tagged with a unique combination of inline tags (Supplementary File S1).

The first-step PCR results were cleaned up with magnetic beads. The PCR product was mixed with clean-up buffer (2.5 M NaCl, 10 mM Tris, 1 mM EDTA, 20% PEG 8000, 0.05% Tween 20, 2% carboxylated Sera-Mag SpeedsBeads (Cytiva Life Sciences, Marlborough, MA, United States), pH 8) at a 0.8x ratio and incubated for 5 min, washed two times with wash buffer (10 mM Tris, 80% EtOH, pH 7.5) for 30 s, dried for 5 min at RT and finally eluted in 40 μL of elution buffer (10 mM Tris, pH 8.5).

During the second-step PCR, samples were amplified with a final concentration of 1x Multiplex Mastermix, 1x Coralload Loading Dye, 100 mM of each primer, and 2 μL of the first-step product. Cycling conditions were the same as in the first-step PCR except for 61°C as annealing temperature and a decreased cycle number of 20. In the second-step PCR, each of the 96 wells was individually tagged so that the combination of the in-line tag from the first-step PCR and the index-read of the second-step PCR yielded a unique combination per sample. PCR success was verified using a 1% agarose gel.

PCR products were normalized to equal concentrations with normalization buffer (same as clean-up buffer, but with only 0.1% beads) following the same protocol as the clean-up after the first step but with a ratio of 0.7x and an elution volume of 50 μL. All normalized products were pooled in the final libraries in equal parts. The libraries were concentrated using a silica-membrane spin column (Epoch Life Science, Missouri City, TX, United States) by mixing 1 volume of the library with 2 volumes of binding buffer (3 M Guanidine Hydrochloride, 90 EtOH, 10 mM Bis-Tris, pH 6) for the binding step (1 min centrifugation, 11,000 x g), 2 washing steps (30 s centrifugation, 11,000 x g) with wash buffer and a final elution (3 min incubation at RT, followed by 1 min centrifugation at 11,000 x g) with 100 μL elution buffer. Library concentrations were quantified on a Fragment Analyzer (High Sensitivity NGS Fragment Analysis Kit; Advanced Analytical, Ankeny, United States). The libraries were then sequenced using the Illumina MiSeq platform with 2 lanes for each library with a paired-end kit (V2, 2×250 bp for 16S and V3, 2×300 bp for ITS) at CeGat (Tübingen, Germany).

DNA and total RNA were separately extracted from samples in 96-well plates using the NucleoMag DNA/RNA Water kit (D-MARK Biosciences, Toronto, Canada) that includes magnetic beads. Instead of using a magnetic plate to separate magnetic beads from buffers, we used the Magnetic Bead Extraction Replicator (V&P Scientific, San Diego, United States), which allows for the transfer of all magnetic beads from one lysate/buffer/elution plate to another without the need to remove the supernatant from individual wells.² The RNA extraction protocol involved a 25-min-long rDNase incubation step to digest DNA. Since the 96-well plates were open during the entire extraction, which posed a contamination risk, we added one negative extraction control to each row of each plate by replacing lysate with pure water. All extractions were performed under a sterile hood. DNA/RNA concentrations of all extracts and all negative extraction controls were measured using a Qubit fluorometer with the dsDNA HS Assay Kit and the RNA HS ASSAY Kit, respectively (Thermo Fisher Scientific, Burlington, Canada).

DNA and RNA libraries of all samples and negative extraction controls were prepared for metagenomics and total RNA-Seq using the NEBNext Ultra II DNA Library Prep Kit for Illumina and the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina, respectively (New England Biolabs, Whitby, Canada). For RNA library preps, we did not perform mRNA enrichment or rRNA removal and instead processed the entire RNA. The RNA library prep kit has a default insert size of 200 bp, and we chose an insert size of 150–350 bp for the DNA library preps to keep insert sizes approximately consistent. After library prep, we randomly selected 8 DNA sample libraries, 3 negative DNA extraction control libraries, 7 RNA sample libraries, and 4 negative RNA extraction control libraries and sent 2.5 μL of each to the AAC Genomics Facility at the University of Guelph, Canada for analysis on an Agilent Bioanalyzer 2,100 system (Agilent Technologies, United States) to confirm successful library preps and check for contaminations in negative extraction control libraries. After consultation with the sequencing facility (Center for Applied Genomics, Hospital for Sick Children, Toronto, Canada), we cleaned up all DNA and RNA libraries following the DNA/RNA library prep kit manual to remove primer dimers and unincorporated primers.

We pooled 5 μL of each DNA and RNA library for sequencing, respectively, including negative extraction controls. We pooled equal volumes instead of equal concentrations because this pooling strategy allows for an equal relative sequencing depth per sample as opposed to an equal total sequencing depth. That way, the relative number of reads per sample mirrored the relative amount of DNA/RNA, avoiding an over- or underrepresentation of samples with higher or lower DNA/RNA amounts. Size distributions of the DNA and RNA library pools were assessed with a bioanalyzer by the sequencing facility, and the average fragment size was 386 bp for the DNA library pool and 436 bp for the RNA library pool. Both pools were paired-end (2×100 bp) sequenced in a 50:50 ratio on a single lane of a NovaSeq 6,000 SP flowcell.

Raw data of the sequencing runs were delivered demultiplexed by index reads. Further demultiplexing by inline tags was done with the Python script “demultiplexer”.³ Sequences were subsequently processed with APSCALE v1.4 (Buchner et al., 2022) using default parameters. Paired-end reads were merged using vsearch v2.21.1 (Rognes et al., 2016). Primer sequences were trimmed with cutadapt v3.5 (Martin, 2011). For 16S sequencing, only sequences with a length of 252 ± 10 bp were retained, and for ITS-2 sequencing, only sequences with a length ranging from 240 to 460 bp were retained. Only sequences with an expected error of 1 passed quality filtering. Reads were dereplicated and singletons were removed. For OTU generation, sequence clustering was performed with a similarity threshold of 97%, and for ESV generation, denoising was carried out with an alpha value of 2 and a minimum size of 8 as implemented in vsearch. Before taxonomic assignment, the resulting OTU and ESV tables were filtered for potentially biased sequences using the LULU algorithm (Frøslev et al., 2017) implemented in APSCALE.

Subsequently, only OTUs and ESVs found in both replicates of the same sample were summed up for all samples. After this initial data filtering, reads still left in the negative controls were subtracted from OTUs or ESVs, respectively, to generate final OTU and ESV tables. Taxonomic assignment was performed using DADA2 with default parameters in combination with the database SILVA 138.1 designed for DADA2 (McLaren and Callahan, 2021) for 16S sequences and the database UNITE (Abarenkov et al., 2021) for ITS-2 sequences, respectively.

In an earlier study, we investigated 672 combinations of bioinformatic tools to identify the best-performing combination to process and taxonomically annotate microbial mock community datasets (Hempel et al., 2022). Based on these results, we processed both metagenomics and total RNA-Seq data as follows: we used Trimmomatic v0.39 (Bolger et al., 2014) to trim the leading and trailing low-quality nucleotides of each read by cutting reads if the average quality of nucleotides in a sliding window of size 4 was below a PHRED score of 20. After trimming, we excluded reads shorter than 25 nucleotides and error-corrected reads using the error-correction module of the assembler SPAdes v3.14.1 (Bankevich et al., 2012). Then we assembled the reads into scaffolds using MEGAHIT v1.2.9 (Li et al., 2015) with the parameter ‘presets’ set to ‘meta-large’ to adjust k-mer sizes for the assembly of large and complex metagenomes. All other parameters were set to default. Subsequently, we mapped reads to assembled scaffolds to determine the abundance of each scaffold using BWA v0.7.17 (Li and Durbin, 2009) with default parameters. We processed mapped reads using the function coverage of samtools v1.10 (Li et al., 2009) to obtain the mean per-base coverage for each scaffold. For taxonomic annotation, we used the SILVA132_NR99 SSU and LSU reference databases (Quast et al., 2013) in combination with kraken2 v2.1.1 (Wood et al., 2019) using default parameters. The setup of the kraken2 database for SILVA required manual adaptations, which are described in the Supplementary material. All code utilized is available on GitHub.⁴

Taxon abundances/P–A represented independent features, and we defined the dependent feature as the combinations of applied insecticide level (none, low, medium, high) and fine sediment addition (normal fine sediment concentration, increased fine sediment concentration) for each sample, resulting in eight classes that were predicted by the machine learning algorithms. Since correlated independent features add noise, we removed them by applying the SULOV (Searching for Uncorrelated List of Variables) algorithm using the function FE_remove_variables_using_SULOV_method of the Python module featurewiz v0.1.55,⁷ which identifies all pairs of highly correlated independent features (features with a Pearson correlation coefficient of >0.7 or < −0.7 by default), determines their Mutual Information Score (MIS) to the dependent feature, and keeps the independent feature with the highest MIS for each highly correlated feature pair.

Each ExStream mesocosm was sampled at two time points as part of the cotton strip assay, which meant that samples consisted of highly related paired samples, i.e., two samples of the same mesocosm. When splitting the data sets into train and test sets, we ensured that paired samples were assigned to the same training and test sets to avoid data leakage between the sets.

Initially, we applied a 90:10 train-test split to the datasets (109 train samples, 12 test samples) and performed training and testing without repetition, but due to large discrepancies between train and test scores, we changed the train-test split ratio to 80:20 (97 train samples, 24 test samples) and repeated both training and testing splits three times in total. During each repetition, we randomly selected 12 pairs (24 samples) of highly related samples for the test dataset and trained and tested all models across all datasets with the same randomly selected 12 sample pairs per repetition.

For feature selection, we used Recursive Feature Elimination to select the 20 most important features using the function RFE from scikit-learn with a DecisionTreeClassifier as the estimator.

It is generally recommended to test multiple machine learning algorithms (Greener et al., 2022), which is why we selected eight machine learning algorithms to predict stressor classes: k-Nearest Neighbors (KNN), Linear Support Vector Classification (LSVC), Logistic Ridge Regression (Ridge), Logistic Lasso Regression (Lasso), Multilayer Perceptron (MLP), Random Forest (RF), Support Vector Classification (SVC), and XGBoost (XGB). For thorough descriptions of these algorithms in a biological context see Greener et al. (2022) and Ghannam and Techtmann (2021).

All algorithms, except XGBoost, are available in scikit-learn. To run the XGBoost algorithm, we used the Python module xgboost v1.6.1 (Chen and Guestrin, 2016), which is compatible with scikit-learn. To optimize hyperparameters while avoiding overfitting, we performed Bayesian hyperparameter optimization with 10-fold cross-validation using the function BayesSearchCV of the Python module scikit-optimize v0.9.0.⁸ The function is compatible with scikit-learn and builds a performance probability model for given hyperparameters, which is used to select the most promising hyperparameters through iterative performance evaluations. While not every possible hyperparameter combination is tested that way, this approach provides a good trade-off between optimization results and runtime. Model prediction performance was evaluated using the Matthews Correlation Coefficient (MCC), which ranges from −1 to 1, where 1 means perfect predictions/performance, 0 means prediction performance as good as random guessing, and − 1 means all predictions are wrong, and increments between −1 and 1 can be interpreted in the same way as increments of the Pearson correlation coefficient. All hyperparameters tested can be found in the publicly available code⁹ and Supplementary File S2. The optimized hyperparameters were then used to train models on the entire training dataset, and model performances to predict classes of the testing dataset were evaluated using the MCC. During training on the entire dataset, learning curves were generated using the learning_curve function from scikit-learn. This process was repeated three times, as described above, and the mean average and standard deviation (SD) of the training and test MCC scores across the three repetitions were determined.

We tested each possible combination of taxonomic datasets (ITS-2, 16S, 16S + ITS-2, metagenomics, and total RNA-Seq), clustering or denoising methods (OTU, ESV; only applicable to amplicon sequencing data), taxonomic levels (phylum, class, order, family, genus, and species), data types (abundance, P–A), feature selection (with feature selection, without feature selection), and classification algorithms (KNN, Lasso, LSVC, Ridge, MLP, RF, SVC, and XGB), resulting in a total of 1,536 evaluated combinations.

For amplicon sequencing data, OTU clustering significantly improved SPP while ESV denoising significantly decreased SPP. This observation is in contrast to the emerging recommendation to denoise amplicon sequences into ESVs (Callahan et al., 2017; Knight et al., 2018). Studies comparing OTU clustering and ESV denoising approaches did not yet reach a consensus, showing that either both approaches lead to similar results (Glassman and Martiny, 2018; Vera-Gargallo et al., 2019; Kang et al., 2021), ESV denoising outperforms OTU clustering (Caruso et al., 2019; Tapolczai et al., 2019; Joos et al., 2020), or vice versa (Roy et al., 2019; Tedersoo et al., 2022). Our results support the latter, although more similar studies are required to determine if clustering or denoising is more appropriate for machine-learning-based environmental predictions using microbial communities.

In general, a higher taxonomic resolution provides a better picture of microbial communities, but our results show that the species level correlated worse with SPP than genus, family, order, and even class levels. For ITS-2 sequencing and omics-based methods, the high number of detected taxa at the species level might have added more noise than value to the data. This is indicated by the significantly positive impact of feature selection on SPP, i.e., the limitation of the number of included taxa. However, for 16S sequencing, feature selection had no impact on SPP while the species level still negatively correlated with SPP. This result may be related to the number of sequences that could not be assigned to the species level and were consequently dropped. The lower the taxonomic level considered, the harder it is to annotate taxonomy due to the lack of reference sequences in databases, and the more sequences are dropped from the downstream analysis. In microbiome amplicon sequencing studies, the taxonomic resolution is usually limited to the genus level due to the difficulty in designing primers that resolve microbial communities at the species level (Knight et al., 2018). Metagenomics allows for taxonomic resolutions at the species level or even strain level, but this requires sufficient sequencing depth (Knight et al., 2018). Dropping sequences from the analysis is equivalent to a loss of information, which could have decreased SPP at the species level. It is also possible that correlations between taxa and environmental variables are higher at lower taxonomic levels because lower taxonomic groups can be overall ecologically coherent, i.e., share similar physiologies, while higher taxonomic groups can be ecologically incoherent and have very different physiologies (Philippot et al., 2010; Choe et al., 2021; Auladell et al., 2022). Once reference databases have been extensively expanded and most sequences can be taxonomically annotated, it will be possible to determine if the lack of reference sequences or ecological incoherency of species explains lower SPP at the species level.

We were surprised that the data types (abundance/P–A) did not have an impact on SPP, given that many studies focus on methods to improve abundance estimates from HTS data (Dillies et al., 2013; Gloor et al., 2017; Weiss et al., 2017; Pereira et al., 2018). The difference in abundance and P–A data lies in the weight of the taxa; in P–A data, abundant and rare taxa are weighted equally, making the data more sensitive to noise but also to subtle differences in community composition. Using simulated data, Koh et al. (2019) demonstrated that P–A data is more powerful when taxa associated with an environmental variable are rare while abundance data is more powerful when those taxa are abundant. However, a large-scale morphological study on benthic invertebrates showed that ecological status classifications based on abundance and P–A data showed only minor variations (Buchner et al., 2019). In a microbial context, multiple HTS studies showed similar correlations of both abundance and P–A data with environmental variables (Muletz Wolz et al., 2018; Knowles et al., 2019; Farinella et al., 2022), while some studies showed that correlations differed between data types (Kask et al., 2020; Tavalire et al., 2021). These results indicate that the impact of data types might depend on the studied environmental variables, but if further research shows that both data types have similar predictive power for environmental assessments, as our results suggest, then P–A data could be used exclusively in future environmental assessment studies. This would avoid the rather complex and partially disagreeing statistical methods required when working with compositional data, i.e., HTS abundance data (Dillies et al., 2013; Gloor et al., 2017; Weiss et al., 2017; Pereira et al., 2018). Furthermore, if abundance and P–A data generate similar results, then the often-stated advantage of metagenomics to generate abundance data free from target PCR bias (Knight et al., 2018; Khachatryan et al., 2020) would become irrelevant, which would decrease the value of omics-based approaches in comparison to amplicon sequencing.

Feature selection can be applied to microbial data to remove noninformative, noisy, or redundant features (Ghannam and Techtmann, 2021). This is generally recommended because the high number of observed features can increase the risk of overfitting, which is described as the “curse of dimensionality” (Oudah and Henschel, 2018). However, feature selection goes against the proposed idea that a more holistic picture of environmental microbial communities is beneficial for predicting environmental variables, as it reduces the number of taxa included in prediction models. Our results suggest that feature selection improves SPP overall and especially for metagenomics, while the SPP of 16S sequencing was not impacted by feature selection. This indicates that the increased biodiversity coverage of omics-based methods might in fact not be beneficial for machine learning predictions and that datasets covering a lower number of taxa, as generated by amplicon sequencing, might result in more accurate and precise predictions. It should be noted, though, that ITS-2 sequencing detected approximately as many species as total RNA-Seq, and feature selection did increase the SPP of ITS-2 sequencing, showing that amplicon sequencing can also be significantly impacted by feature selection. Furthermore, the sequencing depth of metagenomics and total RNA-Seq in our study was very low, which could have influenced the impact of feature selection. If similar studies with a sufficient sequencing depth come to the same conclusion that omics-based methods in fact detect too many taxa for accurate and precise environmental assessments and require feature selection, then this would strongly tip the balance in favor of amplicon sequencing.

Machine learning algorithms had a substantial impact on SPP, and even when applying two different algorithms to the same data set, the resulting MCC could range from 0.38 to −0.05. This illustrates the importance of testing multiple machine learning algorithms, which is recommended in general (Greener et al., 2022). One of the most commonly applied machine learning classification algorithms for HTS data is RF (Smith et al., 2015; Frühe et al., 2020; Hermans et al., 2020; Lanzén et al., 2020; Dully et al., 2021; Ghannam and Techtmann, 2021; Marcos-Zambrano et al., 2021), which reveals which feature contributed most to a prediction. Other popular algorithms are XGB, Support Vector Machines (which include SVC and LSVC), Logistic Regression, and KNN (Ghannam and Techtmann, 2021; Marcos-Zambrano et al., 2021; Greener et al., 2022). However, among those algorithms, RF and (L)SVC did not significantly correlate with SPP in our study, while XGB and KNN significantly negatively correlated with SPP and only logistic regression, specifically Lasso and Ridge, significantly positively correlated with SPP. Linear algorithms have the lowest flexibility among all popular machine learning algorithms, since they assume only linear relationships, and while other algorithms can assume non-linear relationships, which increases their flexibility and is often considered beneficial for the analysis of large and complex data, this was not the case for our study. In contrast, MLP, which represents a simple neural network (NN) with the highest flexibility among all algorithms tested in our study, performed overall the best after Lasso and Ridge and specifically the best for omics-based methods that generated the largest datasets. NNs are currently among the most powerful machine learning algorithms for the analysis of extremely large data, and their impact is so significant that an entirely new field of research emerged around NNs, called deep learning (Greener et al., 2022). To unfold their potential, NNs require large amounts of samples that usually go beyond the number of samples generated in a single biological study. However, thousands of sampling sites are monitored for routine environmental assessments, and once the broad application of omics-based methods becomes more affordable, it will be interesting to see if NNs are required for good SPP based on omics data or if less complex machine learning algorithms will be sufficient or even more appropriate.

Overall, our study shows that data-processing methods should be chosen carefully since they can have a high impact on SPP and that methods resulting in the single best SPP are not necessarily the most appropriate overall. Therefore, we conclude that it is advisable to explore multiple sequencing and, in particular, data-processing methods to maximize prediction performance.