Microbial sequencing
DNA was extracted from 0.75 g of soil using the Power Soil DNA extraction kit (Qiagen) according to the manufacturer’s instructions. The primers ITS4ngs and ITS3mix targeting the ITS2 region of fungal genes (Tedersoo et al. 2015) and the primers 515FB and 806RB (Caporaso et al. 2012; Apprill et al. 2015; Paradaet al. 2016) targeting the V4 region of the 16S rRNA gene in bacteria were used. Presence of PCR product was verified using agarose gel electrophoresis. The PCR products were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter). Adapters and barcodes were added to samples using the Nextera XT DNA library preparation kit set A, B, and C (Illumina, San Diego, CA, USA). The final PCR product was purified again with AMPure beads, checked using agarose gel electrophoresis and quantified using a Nanodrop spectrophotometer before equimolar pooling. We pooled all fungal samples (192) from each time point in one Illumina Miseq PE250 run and divided the bacterial samples over two separate runs (96 samples each). With two time points analysed, this resulted in two MiSeq runs for fungi and four runs for bacteria. Libraries were sequenced at McGill University and Genome Quebec Innovation Center. Extraction negatives and a mock community were used and further sequenced in each sequencing run. This mock community consisted of 10 fungal species and associated bacteria, and was included to control for the potential variation between sequencing runs and to increase the accuracy of the bioinformatics analysis.
Bacterial sequences and fungal sequences were analysed using the PIPITS pipeline and the Hydra pipeline, respectively (Gweon et al. 2015; De Hollander 2017). In short, fungal sequences were paired using VSEARCH and quality was filtered using standard parameters of the pipeline. The ITS2 region was extracted using ITSx (Bengtsson-Palme et al.2013). Short reads were removed, and sequences were clustered based on a 97% similarity threshold using VSEARCH and fungal chimeric sequences were removed by comparing with the UNITE uchime database. The representative sequences were identified using the RDP classifier against the UNITE database (Abarenkov et al. 2010) and clustered further into phylotypes. For bacterial sequences, VSEARCH was used to pair and cluster sequences. For classification, the SINA classification was used with the SILVA database. Fungi were assigned to potential functions using FunGuild (Nguyen et al. 2016) and curated using in-house databases (Hannula et al. 2017; Tedersoo et al.2015; Mommer et al. 2018). For potential-plant pathogenic fungi, the target plant species was checked based on available literature (Watanabe 2018). The sequences created in this study are deposited to ENA with accession number PRJEB31856.
The data from the conditioning phase (May 2017) and two months after establishment of the responding phase (September 2017) were analyzed together. We first filtered out microbial taxa that were present in less than ten samples and had an abundance of less than 0.01%. For ITS data, the sequences derived from other organisms than fungi were removed and for 16S data, sequences originating from mitochondria and chloroplasts were removed. For both bacteria and fungi, samples with less than 1000 reads remaining were removed from the dataset and read numbers were normalized using total sum scaling (TSS).
Mock communities were used to inspect the filtering done for fungi and to compare the runs with each other. After filtering, we detected the same 13 fungal OTUs in all mock communities sequenced and the abundances of these OTUs in the mock communities in different sequencing runs were highly correlated (Pearson correlation, R2=0.97, p<0.001).