Processing of microbial community datasets

Generated sequencing data were initially quality filtered with the Trimmomatic tool version 0.36 (Bolger et al., 2014). Low quality reads were truncated if the quality dropped below 15 in a sliding window of 4bp. Subsequently, all reads shorter than 100bp and orphan reads were removed. Remaining sequences were merged, quality-filtered and further processed with USEARCH version 10.0.240 (Edgar, 2010). Filtering included the removal of reads shorter than 250 bp or longer than 490 (fungi) or shorter than 400 or longer than 470 bp (bacteria) as well as the removal low quality reads (expected error > 1) and reads with more than one ambitious base.

Processed sequences of all samples were dereplicated, concatenated, and obtained unique sequences were denoised and clustered into zero-radius operational taxonomic units (zOTUs) with the unoise3 algorithm implemented in USEARCH version 10.0.240 (Edgar, 2010). All OTUs consisting of one single sequence (singletons) were removed. Subsequently, remaining chimeric sequences were removed using UCHIME (Edgar et al. 2011) in reference mode with the SILVA SSU Ref NR 99 132 database (Quast et al., 2013) as reference data set for bacteria and the QIIME release of the UNITE database version 7.2 (Kõljalg et al., 2013)

161

for fungi. Filtered sequences were mapped on remaining unique sequences to determine the occurrence and abundance of each unique sequence in every sample. To assign taxonomy of bacteria and fungi, chimera-free sequences were classified by BLAST alignment against the most recent SILVA database (Quast et al., 2013) and the most recent UNITE database (Kõljalg et al., 2013), repectively, with an e-value threshold of 1e-20. Concatenated sequences of all sequences were mapped on the final set of unique sequences to calculate the evenness and abundance of each unique sequence in all samples. All non-bacterial or non-fungal zOTUs were removed based on their taxonomic classification in the respective database. Final zOTU tables for bacteria and fungi are provided in Tables S4 and S5, respectively. Only zOTUs occurring in more than one sample were considered for further statistical analysis.

Samples with less than 22 (bacteria) and 16 (fungi) sequences per sample were removed prior statistical analysis, resulting in 229 samples for bacteria, and 231 samples for fungi.

2.12 Statistical Analysis

All statistical analyses were performed using R version 3.4.0 (R Core Team, 2016) and the packages therein. Differences were considered as statistically significant with p ≤ 0.05.

Differences in alpha or beta diversity as well as sequencing depth with regard to cropping system, treatment and inoculation (yes/no) were tested by a Kruskal-Wallis test. There were no significant differences of the mean sequencing depths between the intercropping or monocropping systems, treatments or inoculation. In consequence, zOTU tables were not rarefied as recommended by McMurdie and Holmes (2014).

A variety of alpha diversity indices (Richness, Shannon index of diversity and Michaelis Menten Fit) were calculated using the R-packages picante version 1.6-2 (Kembel et al., 2010) and drc version 3.0-1 (Ritz and Streigbig, 2016). Sample coverage was estimated using the Michaelis-Menten Fit calculated in R. For this purpose, richness and rarefaction curves were calculated utilizing the picante package (Kembel et al., 2010). OTU tables were rarefied using the rrarefy function in vegan version 2.4.4 and samples with less than 17,177 (soil bacteria), 4,172 (rhizosphere soil bacteria), 508 (root bacteria), 22 (leaves bacteria), 6,105 (soil, fungi), 596 (rhizosphere soil, fungi), 61 (root, fungi) and 16 (leaves, fungi) sequences were removed prior alpha diversity analysis. Richness and diversity were calculated using the specnumber and diversity functions, respectively. The Michaelis-Menten Fit was subsequently calculated from generated rarefaction curves using the MM2 model within the drc package (Ritz and Streibig, 2016). All alpha diversity indices were calculated 10 times. The average from each iteration was used for further statistical analysis. Final tables

162

containing bacterial and fungal richness and diversity are provided in Tables S6 and S7, respectively.

Differences in community composition were investigated by permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001) based on Bray-Curtis distance matrices using 999 permutations. A significant p-value in PERMANOVA for beta-diversity can be driven by true biological differences, differences within treatment (variance) or both (Anderson, 2001). In case of significant p-values in PERMANOVA, we tested for differences in homogeneity using permutational analysis of multivariate dispersions (PERMDISP, Anderson, 2006) with 999 permutations. NMDS, PERMANOVA and PERMDISP were run using functions; metaMDS, adonis and betadisper, respectively, in the R package vegan (Oksanen et al., 2016). To investigate in differences between cropping cropping regimes, pairwise Adonis with p-value adjustment “BH” based on Bray-Curtis distances were used (Martinez Arbizu, 2017). Bacterial and fungal communities were tested separately. Differences with regard to crop species were tested after exclusion of unplanted soil samples. The effect of cropping systems and inoculation on diversity and richness of fungi and bacteria in all investigated compartments were analyzed separately to avoid spatial pseudoreplication. Global effects (calculated for both sampling times together) of plant compartment and crop species on fungal and bacterial communities were tested with strata = pot, as we had pseudoreplicated data. The two sampling dates were also analyzed separately as to assess whether the observed effects at harvest 1 would be maintained at harvest 2.

Data (including plant and soil parameter as well as alpha-diversity) were tested for normal distribution with shapiro.test Test and homogeneity of variance with leveneTest function with the package car version 2.1-5 (Fox and Weisberg, 2011). Differences between measured environmental and plant parameters were calculated with Kruskal-Wallis test followed by multiple comparing using dunnTest with Benjamini-Hochberg p-value adjustment or Tukey’s post hoc test using the HSD.test function in the R package FSA version 0.8.17 (Ogle, 2016) and agricolae version 1.2-8 (De Mendiburu, 2014), respectively.

To identify zOTUs highly associated to application with respect to plant compartment and harvest, multipattern analyses were applied. For that purpose, bacteria and fungi were investigated using the multipatt function from the IndicSpecies package version 1.7.6 (DeCáceres and Legendre, 2009). Only fungal and bacterial zOTUs found in at least three samples were used. The biserial coefficients (R) with a particular plant species and treatment were corrected for unequal sample size using the function r.g (Tichy and Chytry, 2006).

163 3. Results and Discussion