Saprotrophic fungal diversity predicts ectomycorrhizal fungal diversity along the timberline in the framework of island biogeography theory

Host and site

B. ermanii Chamiss, a deciduous broad-leaved tree species, occurs naturally in open spaces within subalpine and boreal forests in Northeast China, Japan and the Russian Far East.^40,41,42 B. ermanii is also a typical tree species of the timberline, because it can resist to severe frost, strong wind and low temperature by ecophysiological and ecomorphological flexiblity.^43,44,45 In addition, a rich mycobiome of B. ermanii has been reported,^16,46,47 which may facilitate the survival and spread of the host plant.

On the northern slope of Changbai Mountain, Northeast China, B. ermanii grows over a broad elevation range from ca. 1700 to 2100 m a.s.l., and forms pure stands along >300 m vertical belt below the timberline.^42,48 The establishment of Changbai Nature Reserve (CNR) in 1960 gave strict protection to the Erman’s birch (B. ermanii) forests.⁴⁹ All our sampling was located within the core region of the CNR. The region has a typical continental temperate monsoon climate, with higher elevations experiencing lower temperature and greater precipitation.⁵⁰ Soils of the Erman’s birch forests are Permi-Gelic Cambosols, whereas the soils beyond the upper and lower limits of B. ermanii zone are Permafrost cold Cambosols and Umbri-Gelic Cambosols, respectively. Climate change, particularly rising temperature, threatens the populations of B. ermanii and the whole Erman’s birch forest ecosystem in Changbai Mountain.^50,51

Field sample collection

We sampled fine roots and neighboring soils for B. ermanii individuals at six elevation-related habitats of B. ermanii on 2–8 September, 2018 (Fig. 1a). These six elevation habitats include the upper limit (2069–2116 m), tree islands (1997–2042 m), treeline (1949–1992 m), pure stands (including two sub-sites isolated by over 2 km; 1900–1926 m), ecotone of dark coniferous forests and Erman’s birch forests (1742–1765 m) and lower limit (sparse individuals in coniferous forests; 1688–1706 m), respectively. In each habitat, 14 trees were randomly selected, and all trees were located more than 20 m apart from other sampled trees to ensure independence of each sample. Neighboring soils and fine roots were collected according to the protocols of Yang et al.²⁸ and Lankau and Keymer,⁵² respectively. Briefly, with the trunk as center and the DBH as distance from the stem, we collected four soil cores (diameter = 3.5 cm, depth = 10 cm) after removal of litter and mixed them as a single composite soil sample (Fig. 1b). We excavated surface soils near the base of trees and traced roots from the base to terminal fine roots in three directions. The fine roots of three directions were combined as a composite fine root sample, and each raw sub-fine-root section was nearly 6 cm wide and 8 cm long (Fig. 1c, d). All the samples were brought back to the laboratory with ice bags within 8 h. Soil was sieved through a 2-mm mesh and divided into two subsamples: one was stored at 4 °C to determine the soil properties, whereas the other was stored at −40 °C for subsequent DNA extraction. Fine roots were rinsed with sterile water and cut into 1.5-cm segments: one subsample (ca. 80%) was stored at 4 °C to determine root biochemical traits, and the other subsample (ca. 20%) was stored at −40 °C for subsequent DNA extraction. In the field, tree height, canopy diameter, DBH, elevation, slope, latitude, and longitude of each sampled tree were recorded. In total, 84 fine roots and 84 neighboring soils were collected.

Measurement of soil properties and root traits

We measured 28 soil properties, including soil pH, moisture, conductivity, dissolved organic carbon, dissolved organic nitrogen (DON), ammonium nitrogen, nitrate nitrogen, total carbon, total nitrogen, total phosphate, total potassium, total calcium, total magnesium, total manganese, total iron, total aluminum, available phosphate, available potassium, available calcium, available magnesium, available manganese, available iron, available aluminum, C/N ratio, C/P ratio and the proportions of clay, silt and sand. The measurement methods of soil pH, moisture, ammonium nitrogen, nitrate nitrogen, total carbon, total nitrogen and total content of other elements followed our recent study.²⁸ In addition, soil conductivity was determined with a soil to water ratio of 1:5 by conductivity meter (Mettler Toledo FE30, Shanghai, China). Mehlich 3⁵³ and three-acid-system (nitric acid, perchloric acid, and hydrofluoric acid) were used to extract the available and total content of elements, respectively. Total and available content of phosphate, potassium, calcium, magnesium, manganese, iron, and aluminum were measured using an ICP Optima 8000 (Perkin-Elmer, Waltham, MA, USA). The proportions of clay, silt, and sand were measured by Laser Particle Sizer LS13320 (Beckman, Brea, CA, USA).

Eighteen root traits, including root total carbon (RTC), root total nitrogen (RTN), root phosphate, root potassium, root calcium, root magnesium, root manganese, root iron, root aluminum, root C/N ratio, root N/P ratio, lignin, cellulose, hemicellulose, soluble sugar, soluble protein, free amino acid (FAA) and free fatty acids (FFA), were also measured. Specifically, RTC and RTN were determined with a carbon–hydrogen–nitrogen (CHN) elemental analyzer (2400 II CHN elemental analyzer; PerkinElmer, Boston, MA, USA). Root phosphate, root potassium, root calcium, root magnesium, root manganese, root iron, and root aluminum were measured in ICP Optima 8000 (Perkin-Elmer, Waltham, MA, USA). Soluble protein was measured by a dye-binding assay.⁵⁴ FAA was analyzed by the amino acid analyzer L-8800 (Hitachi, Tokyo, Japan) with leucine as the standard sample. FFA was determined by NEFA FS kits (Diasys, Holzheim, Garman) and the automatic biochemical analyzer AU680 (Olympus, Tokyo, Japan). The measurement methods of lignin, cellulose, hemicellulose, and soluble sugar followed that of our previous study.¹⁶

Calculation of distance to forest edge

The location of each tree individual was determined by latitude and longitude. A high-resolution map (treecover2000) on global forest cover at a spatial resolution of 30 m was used as a base map.⁵⁵ In the map, the areas where forest cover was more than 30% were defined as the “mainland” in the IBT framework and shown as the green grids in ArcGIS (Fig. S1). This standard referred to the proposal of Convention on Climate Change Kyoto.⁵⁶ Then, we calculated the minimum distance of each tree to the neighboring forest edge (i.e., green grids) by using the function Near of the Proximity tool box in ArcGIS 10.0 (ESRI, Redlands, CA, USA).

Sequencing and bioinformatics

Soil total DNA was extracted from 0.5 g of soil by using FastDNA^® Spin kit for Soil (MP Biomedicals, Solon, Ohio, USA). Total DNA of fine roots was extracted from 0.3 g of plant tissue by using Qiagen Plant DNeasy kits (Qiagen, Hilden, Germany). PCR procedures, including primers (ITS1-F: CTTGGTCATTTAGAGGAAGTAA, ITS2: GCTGCGTTCTTCATCGATGC) and conditions were described in our previous studies.^16,28 The PCR products of all samples were normalized to equimolar amounts and sequenced on the Illumina MiSeq PE300 platform of the Majorbio Company, Shanghai, China.

We first merged the paired-end reads using FLASH.⁵⁷ QIIME 1.9.0⁵⁸ and Cutadapt 1.9.1⁵⁹ were applied for quality filtering, trimming, and chimera removal. Altogether 8,238,146 sequences passed quality filtering (parameters: minlength = 240; maxambigs = 0; phred quality threshold = 30). ITSx 1.0.11 was used to remove the flanking small ribosomal subunit (SSU) and 5.8 S genes,⁶⁰ leaving the ITS1 region for further analyses. The putative chimeric sequences were removed using a combination of de novo and reference-based chimera checking, with the parameter –non_chimeras_rentention = union in QIMME.⁶¹ The remaining sequences were then clustered into operational taxonomic units (OTUs) at 97% similarity threshold by using USEARCH.⁶² Singletons were also removed during the USERCH clustering process. Fungal taxonomy was assigned to each OTU by using the Ribosomal Database Project Classifier with minimum confidence of 0.8.⁶³ The UNITE v.8.0 (http://unite.ut.ee) release for QIIME served as a reference database for fungal taxonomy.⁶⁴ The OTU table was then curated with LULU, a post-clustering OTU table curation method, to improve diversity estimates.⁶⁵

After removing non-fungal sequences, the final data set included 7,849,126 fungal sequences covering 6663 OTUs in 168 samples (minimum 4662; maximum 70,779; mean 46,721 sequences per sample). The rarefaction curves of the average observed OTU number are shown in Fig. S2. FUNGuild was used to assign each OTU to a putative functional guild, and the assignments with confidence ranking “possible” were assigned as “unknown” as recommended by the authors.⁶⁶ We further modified the assignment of EcM fungi (the subject in the present study) and their lineages according to.³¹ For some OTUs that were simultaneously assigned to endophytic, saprotrophic, or pathogenic fungi, we considered these as endophytes in roots and saprotrophs in soil samples.

Statistics

All statistical analyses were conducted in R 3.5.2⁶⁷ and AMOS 21.0 (AMOS IBM, New York, USA). In order to analyze the alpha diversities of soil fungi and the three most dominant guilds (viz., EcM, endophytic and saprotrophic fungi) at the same sequencing depth, the data set was subsampled to 4662 reads with 30 iterations. The mean number of observed OTUs was used to represent the diversities of total fungi, EcM fungi, endophytic fungi, and saprotrophic fungi, as previously implemented in.^15,68 Numbers of EcM fungal lineages and saprotrophic genera, families, orders, and classes of each sample were also calculated based on the same subsampling.

First, linear and quadratic regression models were used to determine the effect of elevation on diversities of total fungi, EcM fungi, endophytic fungi, and saprotrophic fungi. The model with lowest Akaike’s information criterion (AIC) value was selected. In order to account for spatial effects, linear mixed-effects models (LMMs) were fitted using the lme4 package⁶⁹ to analyze the variation in diversities of total fungi, EcM fungi, endophytic fungi, and saprotrophic fungi along the elevation gradient with latitude and longitude as random factors. Corrected Akaike Information Criterion (AICc) for small data sets was used to identify the best mixed-effects model from linear and quadratic polynomial models. The significance of each LMM was tested by the function Anova in the car package.⁷⁰ Marginal (m) and conditional (c) R² were calculated by the function r.squaredGLMM in the MuMIn package.⁷¹ Marginal R² (R²_m) represents the variance explained by fixed effects, whereas conditional R² (R²_c) represents the variance explained by both fixed and random effects.

Second, to test the application of IBT on EcM fungal diversity, DBH and RTC of B. ermanii were chosen as proxies of island area and energy, respectively, whereas DFE was chosen as the proxy of island isolation (i.e., island distance to mainland). Linear regression models were used to assess the species-area, species-energy, and species-isolation relationships. In order to account for spatial effects, LMMs were also used for these independent relationships with latitude and longitude as random factors as described above. Classical power-law function models were used to identify the species-area relationship for EcM fungal diversities in roots and soils using z-values in the formula S = CA^{z 72} to compare EcM fungi with macroorganisms in previous studies. Furthermore, ordinary least squares (OLS) multiple regression models were performed to identify the relative contributions of DBH, RTC, and DFE on pattern of EcM fungal diversity when considering other predictive variables. Here, five spatial vectors (PCNM1-5) with significant positive spatial autocorrelation (Fig. S3) were obtained by the principal coordinates of neighbor matrices (PCNM) method,⁷³ and added into OLS multiple regression models to consider the possible geographic effect. EcM fungal diversities in roots and soils, 28 soil properties, 18 root traits, elevation, slopes, DBH, tree height, canopy diameter, and DFE were standardized (average = 0 and SD = 1) before the OLS multiple regression analysis. AIC was used to identify the best OLS multiple regression model, as implemented in the MASS package.⁷⁴ Variance inflation factor (VIF) was calculated for each model by the function vif in the car package. We used the criterion VIF < 3 to adjust to multicollinearity of predictive variables. The function forward.sel in the packfor package⁷⁵ was implemented to estimate the relative contributions of each predictive variable on the variation in EcM fungal diversity. The best OLS multiple regression models in this step are called as OLS multiple regression models #1 to distinguish these from the following models.

Third, to test the possible effect of biotic interactions on EcM fungal diversity, diversities of endophytic and saprotrophic fungi in soils and roots were added into the OLS multiple regression models. These resulted in best OLS multiple regression models #2. The significant differences between OLS multiple regression models #1 and #2 were identified by the function anova in the stats package.⁶⁷ In addition, partial least squares regression (PLSR) was performed to identify the effect of biotic interactions (i.e., saprotrophic fungal diversity) on EcM fungal diversity, as implemented in the pls package.⁷⁶ In order to test the evolutionary source effect within the fungal kingdom, PLSR was also used to identify the relationships between EcM fungal lineages and saprotrophic fungal taxa at the genus, family, order, and class levels, respectively.

Finally, we used the integrated model (i.e., structural equation modeling (SEM)) to combine different ecological frameworks, including IBT, biotic interactions, the climate-driven hypothesis and effects of soil properties and root traits, to predict EcM fungal diversity along the timberline. Specifically, (1) we built a SEM theoretical model on EcM fungal diversity (Fig. S4), in which elevation, island isolation, island area, island energy, biotic interactions, and soil properties acted as ‘predictive aspects’. (2) Based on the best fitted models among OLS multiple regression models #1 and #2, we added the corresponding predictive variables to SEM: Our first step was to assign the predictive variables from OLS multiple regression models to each ‘predictive aspects’ in SEM by 1:1 (e.g., DFE as island isolation and DBH as island area). Because there were more than one of soil properties and root traits screened in the best OLS multiple regression models, we then continued to add other predictive variables into SEM stepwise. (3) AIC was used to screen the best SEM model among alternative models. The total, direct and indirect effects of each ‘predictive aspects’ were interpreted by summing the standardized path coefficient (SPC). Alternatively, variation partitioning analysis (VPA) was performed to identify the shared and independent contributions of IBT, biotic interactions, and other predictors (i.e., elevation and soil properties) on EcM fungal diversity by the function varpart in the vegan package.⁷⁷

Source: Ecology - nature.com