2.5. Statistical Analysis

The “soil total fungi” dataset contained the total OTUs detected by high-throughput sequencing. The soil fungal community dataset was divided into three subsets based on trophic modes: the “symbiotic fungi”, “saprophytic fungi”, and “pathogenic fungi” datasets (i.e., containing the OTUs of the trophic modes symbiotroph, saprotroph, and pathotroph, respectively.). The four fungal datasets (total, symbiotic, saprophytic, and pathogenic fungi) were analyzed separately. A Mantel test was used to check whether soil samples were independent or spatially autocorrelated [11]. The OTUs table was used for downstream community composition and diversity analysis. The relative abundance of soil fungal phyla levels, trophic modes, and functional guilds in each sample were calculated and ranked. Non-metric multi-dimensional scaling (NMDS) analysis based on the Bray–Curtis distance and analysis of similarities (ANOSIM) were carried out to examine soil fungal community composition dissimilarities among different vegetation types [58]. When exploring the relationship between altitude variation and soil fungal community composition dissimilarities, a generalized additive model was used within the ordisurf function in R package ‘vegan’ (version 2.6-2) [59]. The α-diversity was estimated using OTUs richness, the Simpson index, Pielou index, and Shannon–Wiener index [60]. To calculate β-diversity, the Jaccard index was used for the pairwise dissimilarity of species composition and the Bray–Curtis index was used for abundance-weighted dissimilarity. The β-diversity and its components (turnover and nestedness) were computed using the function beta.pair and beta.pair.abund in R package ‘betapart’ (version 1.5.4) [61]. Altitude distance (based on Euclidean distances) was computed using the function vegdist in R package ‘vegan’ (version 2.6-2).

To disentangle the relationship between soil fungal community composition and environmental variables, a distance-based redundancy analysis (db-RDA) and Monte Carlo permutation test (999 permutations) were performed. For the db-RDA analysis, the significance of a full model including all the explanatory variables was first tested and then the model was simplified by forward-model selection, using the function ordiR2step in R package ‘vegan’ (version 2.6-2). To evaluate the relationship between α-diversity and environmental variables, Pearson correlation analysis was used. The pairwise geographic distance (GEO) among sample plots was calculated using the R package ‘geosphere’ (version 1.5-14) [62], according to the geographic coordinates of each sample plot. Multiple regression on matrices (MRM) methods were used to examine the relative effects of geographic distance and dissimilarity in environmental variables (based on Euclidean distances) on β-diversity and its two components (based on the Bray–Curtis index). To exclude strong collinearity among environmental variables, environmental variables with high correlations were removed before the MRM test. To prevent the influence of data overfitting, two MRM tests were performed in R package ‘ecodist’ (version 2.0.9) [63]. The first MRM test was run to remove insignificant variables, then a second test was run with significant variables, and the model results for the second test were reported. All statistical analyses and plotting were conducted in R (version 4.1.1; http://www.r-project.org/ (accessed on 26 February 2022)).