Bioinformatics

Adapter removal, quality trimming and demultiplexing using molecular identifier (MID) sequences were performed by the sequencing company (Fasteris). The base quality of the sequence reads was checked using FastQC⁷⁸. A split-sample filtering protocol for Illumina amplicon sequencing was used for two technical replicates per DNA sample⁷⁹. The raw sequences were quality filtered⁸⁰ (PRINSEQ-lite v.0.20.4) to remove reads with an average Phred quality score below 25. The paired-end reads were assembled and quality filtered with PANDASeq⁸¹ (v2.10). All reads with uncalled bases, an assembly quality score below 0.9, a read overlap below 20, or a base with a recalculated Phread-score below 1 were removed. After dereplicating chimeras were identified and filtered using UCHIME⁸² (v7.0.1090) with default settings. Sequences that were not present in both sample branches were discarded⁷⁹. Operational taxonomic units (OTUs) were created swarm with default settings⁸³ and amplicon sequencing variants (ASVs) were created with DADA2⁸⁴ (1.18.0). Sequences can barely be tracked down to represent distinct species or even strains, thus, the use of taxonomic units (OTUs or ASVs) is mandatory to perform analyses. The tools were used within the modular bioinformatic pipeline Natrix⁸⁵. Sequences wereBLAST (Basic Local Alignment Search Tool) aligned against the National Center for Biotechnology Information (NCBI) nucleotide database (nt) (downloaded Feb. 2020). OTUs with less than 100 reads in total, occurring in less than 10 samples, having a BLAST e-value larger than 1 × 10^–5 or a percentage of identity less than 90% were discarded. Sequencing depth waschecked with rarefaction curves and samples that did not reach saturation were manually curated. Namely, the samples NP1_336h, NP5_168h, ZnNP5_168h, Zn5_336h, C2_336h were removed before analysis.

True diversities and functional group abundances were computed from rarefied reads (function rrarefy from the vegan R Package⁸⁶, version 2.5.7⁸⁷). Reads were rarefied to the mean read amount (419,158) from all samples. True diversities are based on the Simpson index using the R package RAM⁸⁸ as well as Pielou`s evenness. True diversity, also known as Hill numbers or effective number of species, were chosen as a measurement for diversity as it is not a non-linear diversity index (e.g. Shannon index, Simpson index) but is suitable for linear comparisons between samples. True diversity and OTU richness was statistically compared against the control using Wilcoxon tests.

We assigned the nutritional mode to OTUs by using their BLAST based percentage identity values to species/taxa where the nutritional mode is known. We accepted individual thresholds for correctly resolving the assigned taxonomy to the phylum, clade, class, order, genus, species and strain level, while referring to Table ^S6. For the phylum/clade/class/order level we accepted a percentage identity of 90%, for the family level 92%, for the genus level 94.5%, for the species level 98.7% and for the strain level 99%⁸⁹. These thresholds were derived from research on bacteria and archaea and were used here for eukaryotes, since similar, universally usable, thresholds are not available for the V9 18S region in eukaryotic microorganisms. Differential abundant OTUs were calculated by their log₂fold changes using the R Package DESeq2⁹⁰ on raw count data with an accepted adjusted p value below 0.01. The design formula accounts for effects caused by individual treatments and not by temporal effects (design = ~ time + treatment).

For the non-metric multidimensional scaling (NMDS) we replaced zeroes in our dataset based on a Bayesian-multiplicative replacement (cmultRepl, R zCompositions package⁹¹) and calculated the Aitchison distance, as we are dealing with compositional data⁹². Aitchison distance is used as a community dissimilarity indicator. This distance matrix was used to compare the treatments among each other using a Benjamini–Hochberg p-adjusted pairwise adonis (adonis padj). The same test wasapplied to functional community compositions of individual treatments. Figures were created with R and Adobe Illustrator^17,87.