Ectomycorrhizal fungal communities associated with Larix gemelinii Rupr. in the Great Khingan Mountains, China

Study site and sampling

Field experiments were approved by the National Field Scientific Observation and Research Station of Forest Ecosystem in Greater Khingan Mountains, Inner Mongolia. This study was conducted in three coniferous forests typical for L. gemelinii located on the Greater Khingan Mountains in Inner Mongolia in northeast China. One site was located to the north at the Genhe Ecological Positioning Station (GH) and two other sites were to the south at the Hangguanliang Natural Reserves (HGL) and the Saihanwula Natural Reserves (SHWL) (Fig. S1). The three sites were all second forest ecosystems and approximately 60 years old. The forests were protected without any anthropogenic disturbances. L. gemelinii was the dominant EM plant species at each site with very few Pinus tabuliformis individuals. The study areas were mainly exposed to a cold-temperate continental monsoon climate. The mean annual temperature (MAT) ranged from −4.7 to 0.1 °C. The mean annual precipitation (MAP) ranged from 374 to 506 mm, according to the climate data from the WorldClim dataset, with 30-second spatial resolution (Hijmans et al., 2005).

The field work was carried out in July, 2017. Seven individuals were selected at each site. The plants were >10 m away from each other to ensure sample independence (Lilleskov et al., 2004). The DHB of the individuals collected in the sites were similar to avoid the influence of basal area on EM fungal community. The fine roots were excavated from each plant by tracing the roots from the base of the trunk at three points and merging them as one sample. The root samples were transferred to the laboratory in an icebox and stored at −20 °C until processing. In parallel with root sampling, a rhizosphere soil sample (∼200 g) from each individual was collected and pooled to form a composite sample at each site. In total, 21 root samples and three soil samples were collected across three sites. The soil samples were air-dried, sieved through a two mm mesh, then used for the analysis of soil properties. The latitude and longitude were recorded using a high-sensitivity GPS instrument (M-241, Holux Technology Inc., Taiwan, China). Climate conditions, including the MAT and MAP of each site, were derived from the WorldClim global climate dataset as mentioned above (Hijmans et al., 2005). Details on the geographical location and climate of each site are given in Table S1.

Analysis of soil physicochemical properties

A detailed description of the analysis protocol can be found in Wang et al. (2019b). Soil pH was measured using a FiveEasy pH meter (Mettler Toledo, Zurich, Switzerland) after mixing dried soil with distilled water at a 1:2.5 ratio (w/v). The total nitrogen (N) content was quantified by adopting the semi-micro-Kjeldahl method. The total phosphorus (P) and total potassium (K) were measured using an iCAP 6300 inductively coupled plasma spectrometer (Thermo Scientific, Wilmington, USA). Soil organic matter (OM) was determined by weight loss of 5 g dry soil after combustion of organic matter at 700 °C. The soil properties are presented in Table S1.

Molecular analysis

A detailed description of the molecular analysis can be found in Wang et al. (2019b). Root samples were carefully washed using tap water. The fine roots (<2 mm diameter) were cut into 1-2-cm-long pieces. The EM root tips were identified according to their specific morphological characteristics (e.g., emanating hyphae, shape, and color) under a stereomicroscope. Two hundred healthy EM root tips per sample were randomly selected, resulting in 4,200 EM root tips from 21 samples. The EM root tips were cleaned with sterilized distilled water and stored at −20 °C prior to DNA extraction.

The total genomic DNA was extracted using the modified cetyltrimethylammonium bromide (CTAB) method (Gao et al., 2013). A semi-nested PCR was adopted to amplify the ITS2 region of the rDNA in a Veriti 96-well Thermal Cycler (Applied Biosystems, Foster City, USA). First, the primers ITS1F (Gardes & Bruns, 1993) and ITS4 (White et al., 1990) were used to amplify the entire ITS region in 25 µl reaction mixtures including 1 U KOD-Plus-Neo DNA polymerase (Toyobo, Osaka, Japan), 2.5 µl of 10 × buffer, 2 mM of each dNTP, 25 mM MgSO₄, 10 µM of each primer, and c. 5 ng DNA template. The PCR temperature profile included an initial denaturation at 95 °C for 5 min, 20 cycles of 94 °C for 1 min, 58 °C for 50 s, and 68 °C for 30 s, then a final extension at 68 °C for 10 min. The PCR products were diluted 30 times and 1.5 µl of the resulting solution was used as the template for the second PCR amplification targeting the ITS2 region. The conditions for the second PCR amplification were identical to the first amplification, with the exception of our use of fITS7 (Ihrmark et al., 2012) and ITS4 primers equipped with unique barcode tags (12-base barcode sequences). For each sample, three independent replicates were conducted and pooled to represent the sample. The pooled PCR products of each sample were purified using the Wizard SV Gel and PCR Clean-Up System (Promega, Madison, USA). The DNA concentration of each purified sample was measured in a ND-2000 NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, USA) and then mixed with equimolar amounts (100 ng from each sample). The pooled products were sequenced using the paired-end (2 × 250 bp) option on an Illumina MiSeq PE250 platform at the Environmental Genome facility of the Chengdu Institute of Biology, Chinese Academy of Sciences, China.

Bioinformatics analysis

Raw sequence data were filtered using the QIIME platform (Caporaso et al., 2010). We discarded reads with ambiguous bases >6, mismatched primer or barcode sequences, length <250 bp, or an average quality score <20. The ITS2 region of the filtered sequences was extracted using the fungal ITSx software package (Bengtsson-Palme et al., 2013). Potential chimeras were identified and detected using the chimera.uchime command in Mothur 1.31.2 (Schloss et al., 2009) by referring the entries in the unified system for the DNA-based fungal species linked to the classification (UNITE) database (Koljalg et al., 2013). The high-quality and non-chimeric ITS2 sequences were clustered into operational taxonomic units (OTUs) at the cutoff of 97% using the UPARSE pipeline (Edgar, 2013), during which singletons were also removed. The representative sequence (most abundant) of each OTU was searched against the UNITE database (v. 8.2, release date: 02.04. 2020) using the basic local alignment search tool (BLAST) (Altschul et al., 1990). Fungal OTUs were identified and assigned to a taxonomic identity based on the criteria proposed by Tedersoo et al. (2014). The EM fungal OTUs were identified according to Tedersoo & Smith (2013) and Tedersoo et al. (2014) if they best matched the known EM fungal taxa and lineages. A random resampling procedure was conducted to construct a subset to a depth of 1,038 using the sub.sample command in Mothur in order to eliminate the effect of uneven sequence depths across samples on downstream community analysis. The representative sequences for each EM fungal OTUs were submitted to the European Nucleotide Archive (ENA) under accession nos. LR824644-LR824765. Information about the identified EM fungi is shown in Table S2.

Statistical analysis

All statistical analyses were performed in R v. 3.6.3 (R Development Core Team, 2019). The pairwise Euclidean distance among the sites was calculated based on the latitude and longitude of each site and then converted to the principal coordinates of neighbor matrices (PCNM) vector using the pcnm command in the PCNM package (Dray, Legendre & Peres-Neto, 2006). The positive PCNM eigenvector was retained for subsequent analysis. Accumulation curves of observed EM fungal OTUs in each site were drawn using the specaccum command in the vegan package (Oksanen et al., 2013). The EM fungal OTU richness, and Shannon and Simpson diversity indices were calculated using the diversity command in the vegan package (Oksanen et al., 2013). One-way analysis of variance (ANOVA) was used to examine the site effect on the three diversity indices, and Tukey’s honestly significant difference (HSD) test was adopted to perform the pairwise comparisons to detect significant differences among sites. All data were tested for normality and homogeneity of variance before ANOVA, and only fungal OTUs’ richness was log transformed. A Venn diagram was generated to analyze the shared and unique EM fungal OTUs among the three sites using the online software Venny (v.2.1.0, http://bioinformatics.psb.ugent.be/webtools/Venn/).

Taxonomic composition of the whole EM fungal community was demonstrated in a Krona chart using an online Krona tool (v.2.6, Ondov, Bergman & Phillippy, 2011), where circles from inside to outside represent different taxonomic levels. A distance matrix for the EM fungal communities (Hellinger-transformed OTU read data) was constructed using the Bray-Curtis method (Clarke, Somerfield & Chapman, 2006). The EM fungal community was visualized using non-metric multidimensional scaling analysis (NMDS) with the metaMDS command in the vegan package (Oksanen et al., 2013), and the ordiellipse function was used to fit the 95% confidence intervals of sites onto the NMDS ordination. We fit spatial, climatic, and soil variables to the ordination plot based on the environmental fitting test using the envfit command in the vegan package (Oksanen et al., 2013) to identify the main factors affecting the EM fungal community. We determined the significance of difference in the fungal communities among sites by permutational multivariate analysis of variance (PerMANOVA) using the adonis command in the vegan package (Oksanen et al., 2013). The mantel test was carried out to determine the relationship between the dissimilarity of the EM fungal community and the spatial distance using the mantel command in the vegan package (Oksanen et al., 2013). We used the non-parametric Kruskal-Wallis test for the data of all relative abundances that did not satisfy the normality and homogeneity of variance before and after square root or logarithmic transformation to further evaluate the variation of the relative abundances of dominant EM fungal lineages (>1% of total reads) across different sites. The Dunn’s test with Bonferroni correction was employed to detect the significant differences among sites using the posthoc.kruskal.dunn.test function in the PMCMR package (Pohlert, 2014). Additionally, random forest modelling (Breiman, 2001) was used to identify the main drivers of EM fungal OTUs richness and community composition by using the randomForest command in the randomForest package. Random forest analysis relies on the regression tree method using bootstrap samples of training data and random feature selection in tree construction Liaw & Wiener (2002). An increase in mean squared error (%InsMSE) was used to evaluate the importance of each predictor (Breiman, 2001). The importance and significance of each predictor on response variables were assessed using rfPermute and rp.importance commands in the rfPermute package (Archer, 2016), and the rf.significance command in the rfUtilities package was used to test for model significance (Evans & Murphy, 2018). Additionally, statistical analyses of fungal composition (NMDS, envfit analysis, PerMANOVA, and random forest analysis) were performed based on the presence/absence data (Sorensen distances) to ensure the robustness of our results.

In order to assess the distribution difference of the EM fungal OTUs across three sites, site/EM fungus occurrence analysis was conducted according to the study conducted by Toju, Tanabe & Ishii (2016). Briefly, a sample-level matrix (presence/absence data) was obtained by performing binarization on the original OTUs table (reads data with rows representing samples and columns representing OTUs). The sample-level matrix was transformed into a site-level matrix, in which rows represented sampling sites, columns represented fungal OTUs, and the cell entries indicated the numbers of root samples of specific site-fungus combinations. Randomized site-level matrices were generated using a shuffle-sample null model based on the sample-level matrix with 1,000 permutations. The occurrence of EM fungal OTUs for each site was evaluated based on a d’ interaction specialization index (Blüthgen et al., 2007) using the dfun command in the Bipartite package (Dormann et al., 2009). The standardized d’ values for each fungal OTU were generated using the formulas: d’ = [d’_observed − Mean (d’_randomized)]/SD (d’_randomized), in which d’_observed estimated the d’ of the observed site-level matrix while Mean (d’_randomized) and SD (d’_randomized) were the mean and standard deviation of the d’ values generated by 1,000 randomized site-level matrices. Similarly, the standardized d’ value for each site was estimated using the same procedure. We provided the standardized d’ value of the abundant EM fungal OTUs (>0.1% of total reads) as the estimation of rare OTUs is difficult. Additionally, we identified the site-fungus pairs exhibiting a strong presence by calculating the two-dimensional occurrence (2DO) based on the original and randomized site-level matrices mentioned above, following the equation: 2DO (i, j) = [N_observed (i, j) Mean (N_randomized (i, j))]/SD (N_randomized (i, j)), where N_observed (i, j) is the number of root samples in which a given combination of site and EM fungal OTU was detected in the original site-level dataset, and Mean (N_randomized (i, j)) and SD (N_randomized (i, j)) are the mean and standard deviation of the number of samples for the site–fungus pairs across 1,000 randomized matrices. The P values were corrected based on the false discovery rate (FDR) (Benjamini & Hochberg, 1995).

Fungal database summary

A total of 229,667 non-chimeric ITS2 sequences were filtered from 230,289 raw sequences after quality filtering. These high-quality sequences were clustered into 474 non-singleton OTUs (226,675), of which 141 (186,372 sequences) were identified as EM fungal. After rarefying all samples to the same number of sequences (1,038), 122 EM fungal OTUs were retained for downstream analysis. The 30 most abundant OTUs accounted for 89.9% of the EM fungal sequences (Fig. S2A) and 105 of 122 OTUs (86%) occurred in less than three samples (Fig. S2B).

EM fungal diversity

The accumulation curves of EM fungal OTUs in each site did not reach a plateau, indicating that additional sample collection would result in more undiscovered OTUs (Fig. 1A). One-way ANOVA showed that EM fungal OTUs’ richness (log transformed) and Shannon and Simpson indices significantly differed across three sites, and HSD turkey tests further indicated that the indices of GH were significantly higher than those in HGL and SHWL (Figs. 1B–1D). Venn diagrams showed the number of shared and unique OTUs of L. gemelinii across different sites and the results indicated that only six OTUs were shared by the three sites, meanwhile GH accounted the highest proportion of site-specific OTUs (65 OTUs, 53.3% of total OTUs), followed by HGL (20, 16.4%) and SHWL (8, 6.6%) (Fig. 2). Random forest analysis demonstrated that MAT, spatial PCNM, soil P, N, OM, and MAP were significantly responsible for the EM fungal OTUs’ richness. MAT was the best predictor variable, meanwhile 22.49% of variation was explained through modeling (R² = 0.222, P = 0.007, Fig. 3A). The linear model was adopted to estimate the relationship between EM fungal OTUs’ richness and MAT, and the results indicated that the OTUs’ richness of EM fungi significantly decreased with the increasing MAT (R² = 0.40, P = 0.001, Fig. S3).

EM fungal community composition

The Krona diagram indicated that Tricholoma, Tomentella, Piloderma, and Suillus were the dominant fungal groups at the genus level (Fig. 4). All EM fungal OTUs were assigned to 21 fungal lineages in our study, of which the most abundant were /tricholoma, /tomentella-thelephora, /suillus-rhizopogon, and /piloderma (accounting for 36.7%, 14.7%, 14.4%, and 12.1% of the total sequences) (Table S3 and Fig. S4). NMDS ordination revealed that EM fungal communities of the three sites were clearly separated, particularly the community in GH (Fig. 5). In addition, PerMANOVA further indicated significant differentiation of EM fungal communities across different sites (R² = 0.317, P = 0.01). An environmental fitting test demonstrated all spatial PCNM, climate conditions, and soil properties were significantly correlated with the EM fungal community (Fig. 5 and Table 1). Likewise, Random forest analysis revealed that geographic distance, MAT, soil N, P, and MAP were the most important determinants for dissimilarity in the EM fungal community and the variation explained by the model was 26.52% (R² = 0.264, P = 0.001, Fig. 3B). Mantel test showed that the community of EM fungi was significantly correlated with the geographic distance (R = 0.54, P = 0.001, Fig. S5). The results based on the presence/absence data (Sorensen distances) were similar to those mentioned above (Figs. S6–S8). Kruskal-Wallis tests indicated significant site differences in the relative abundance of dominant fungal lineages in /tricholoma (χ² = 7.59, P = 0.02), /tomentella-thelephora (χ² = 13.26, P = 0.001), /suillus-rhizopogon (χ² = 6.72, P = 0.02), /piloderma (χ² = 7.32, P = 0.02), /wilcoxina (χ² = 12.62, P = 0.002), /russula-lactarius (χ² = 13.75, P = 0.001) and /cortinarius (χ² = 13.15, P = 0.001), but not in /sebacina (χ² = 4.76, P = 0.09) (Fig. S9). For example, the relative abundance of /tricholoma was significantly lower in GH than in HGL and SHWL (Fig. S9A), however, the relative abundances of /tomentella-thelephora, /russula-lactarius, and /cortinarius were significantly higher in GH than in HGL and SHWL (Figs. S9B, S9F and Figs. S9H).

Occurrence analysis on site/fungus association indicated that eight of 42 (19%) abundant EM fungal OTUs significantly occurred in specific site, and all sites (100%) harbored EM fungal OTUs that showed different distributions. Twenty-two of 125 (17.5%) pairs of fungal OTU and sites exhibited remarkably strong occurrences (Fig. 6).