Soil Fungal Community Characteristics at Timberlines of Sejila Mountain in Southeast Tibet, China

Soil fungal community characteristics of alpine timberlines are unclear. In this study, soil fungal communities in five vegetation zones across timberlines on the south and north slopes of Sejila Mountain in Tibet, China were investigated. The results show that the alpha diversity of soil fungi was not different between the north- and south-facing timberlines or among the five vegetation zones. Archaeorhizomyces (Ascomycota) was a dominant genus at the south-facing timberline, whereas the ectomycorrhizal genus Russula (Basidiomycota) decreased with decreasing Abies georgei coverage and density at the north-facing timberline. Saprotrophic soil fungi were dominant, but their relative abundance changed little among the vegetation zones at the south timberline, whereas ectomycorrhizal fungi decreased with tree hosts at the north timberline. Soil fungal community characteristics were related to coverage and density, soil pH and ammonium nitrogen at the north timberline, whereas they had no associations with the vegetation and soil factors at the south timberline. In conclusion, timberline and A. georgei presence exerted apparent influences on the soil fungal community structure and function in this study. The findings may enhance our understanding of the distribution of soil fungal communities at the timberlines of Sejila Mountain.


Introduction
Alpine timberline, an ecological transition zone from forest to alpine krummholz, is one of the most important climate-driven ecological boundaries [1]. In ecotones, mountain forests gradually transition to scrubland or grassland habitats with an increase in altitude. Global warming is conducive to the regeneration and growth of high-altitude forests and has led to the upward climb of timberlines [2,3]. Concomitantly, the tree size, density and coverage at timberlines also change [2,3].
At timberlines, tree distribution results in various habitat patches with different climate [4,5], hydrology [6], soil and biological characteristics [7], providing a diverse microenvironment available for soil microorganisms. Soil-vegetation interactions of timberlines are strongly closed to soil fungal taxa with diverse ecological functions [8]. Tree species may also mediate the soil fungal community [9]; for example, dominant tree species at timberlines are cold-resistant coniferous trees, most of which are ectomycorrhizal. Ectomycorrhizal fungi promote nutrient absorption [10,11], the stress resistance [12,13] of trees and the establishment of seedlings [14,15] at timberlines. Ectomycorrhizal fungi rely on their hosts to spread, and a decrease in their richness with a rise in elevation can be observed at timberlines [8], but the decreased species richness is not necessarily linked to the decreased ectomycorrhizal symbiosis; for example, there is no host effect on ectomycorrhizal fungal richness, and even an increase in species richness with a host decrease can also be observed [16]. Therefore, the ectomycorrhizal species richness is not contrary to the presented important role of ectomycorrhizal symbiosis for certain tree species. While the ectomycorrhizal abundance increases with the host density, the dominant role of certain ectomycorrhizal taxa may result in a decrease in ectomycorrhizal species richness. Certainly, other potential factors also act on the richness of ectomycorrhizal fungi, such as the soil C/N [16] or pH of the timberlines [17]. Suillus, Tomentella and Cortinarius of Basidiomycota and Cenococcum of Ascomycota are dominant genera of ectomycorrhizal fungi at several timberlines all over the world [16,[18][19][20].
Saprophytic fungi are devoted to the decomposition of soil organic matter for promoting the soil nutrient cycle [21]. Soil saprotrophic fungi are usually more abundant under forests than shrubs at timberlines; for instance, Mucoromycota [22]. Saprotrophic fungi from Basidiomycota, Ascomycota and Zygomycota also occur frequently in soils from timberline forests, which are rich in wooden substrates [20]. These saprophytic fungi are believed to promote the turnover of soil organic matter in these harsh ecosystems together with ectomycorrhizal fungi [22].
At present, reports on the soil fungi of timberlines are very limited. We know that the changes in the soil fungal community occur along the elevation gradient due to the tree height [8], but tree species and obvious transitional differences in the number and spatial distribution of trees at different timberline zones with different vegetation characteristics also result in soil environmental heterogeneity [8,23]. For instance, in the Shenlongjia timberline ecotone, the shrubbery has a considerable amount of total and available phosphorus in the soil compared to the coniferous forest, but the soil pH and humidity levels are reversed [20]. Furthermore, in Changbai Mountains, the soil in the tree islands within the birch timberline ecotone has a significantly higher water content, total carbon and total nitrogen than open areas at equal elevation [23]. These soil heterogeneity effects related to timberline vegetation on the soil fungi community of timberlines need to be further explored.
Sejila Mountain is located in southeast Tibet, China, and is characterized by natural alpine timberlines (an altitude of over 4000 m). As the received solar radiation at the timberline on the south slope of Sejila Mountain is six times as much as that on the north slope, the heliophile Sabina saltuaria and shade-tolerant Abies georgei are the dominant tree species at the timberlines on the south and north slopes, respectively. Studies have confirmed the ectomycorrhizae associations between A. georgei and Basidiomycetes (e.g., Russulaceae, Cortinariaceae and Amanitaceae) [24,25], whereas no evidence for a symbiotic relationship between S. saltuaria or genus Sabina and fungi has been found so far, and even the distribution of soil saprotrophic fungi remains unknown at the timberlines of Sejila Mountain. It is urgent to note that the timberline in southeast Tibet has risen 69 m in the past 100 years [3]. Therefore, this research investigated the soil fungal community characteristics at the north-and south-facing timberlines of Sejila Mountain to describe soil fungal community characteristics at the timberlines and clarify the relationship between soil fungal communities and vegetation and soil factors. The fundings will provide important reference value for accurately understanding the role of soil fungi in the changing timberline ecosystems under climate change. Because dominant plants have a great impact on soil fungal community assembly [8,9], we assumed that the mycorrhizal fungi are mainly affected by host distribution, and therefore the diversity, composition and function of mycorrhizal fungi is expected to change with the decrease in mycorrhizal host species.

Study Site
This study was conducted at the Sejila Mountain (29 • 10 -30 • 15 N, 93 • 12 -95 • 35 E, 2200 m-5400 m a.s.l.), located in Nyingchi, Tibet, China [26]. The mountain comprises part of the dark coniferous forest in southeast Tibet. The Sejila Mountain is located in the transition zone between semi-humid and humid areas in southeast Tibet. The mountain extends from northwest to southeast, forming a large range of east-west slope. This region is characterized by sub-alpine temperate semi-humid climate, warm winter, cool summer and distinctive dry and wet seasons. The region has an annual average temperature of 2.0 • C-4.5 • C. The highest and lowest monthly average temperatures of the region are 11.1 • C (July) and −14.0 • C (January), respectively [27]. The region has an average annual relative humidity of 78.8% and average annual evaporation of 544 mm, accounting for 48.0% of the average annual precipitation (1134 mm). The precipitation in spring and summer accounts for 79.4% of the annual precipitation [27]. The soil in the region is mainly mountain brown soil and acid brown soil (pH = 4-6) [26].
The investigated timberlines of Sejila Mountain were divided into the timberline on the south slope (hereinafter referred to as south timberline) and the timberline on the north slope (hereinafter referred to as north timberline). In general, timberline ecotones consist of forest line, tree line and tree species line [28,29]. The forest line is the upper limit of the closed forest, with a canopy density ≤ 0.2. The tree line is the altitude at which the tree height is less than 5 m. The tree species line is the highest altitude of tree growth and comprises the upper limit for growth of isolated dwarf, curved trees. The closed forest is below the forest line and the scrubland is above the tree species line [30]. According to this, five vegetation zones along elevational gradients (<200 m) at the two timberlines were sampled: closed forest (F), forest line (FL), tree line (TL), tree species line (TSL) and shrub (S) [31]. In this study, the closed forests mainly consisted of Sabina saltuaria at the south timberline and Abies georgei at the north timberline. The alpine shrubs were dominated by Rhododendron tanastylum. This is a shrub widely spread under the tree canopies and the open areas of the timberlines, occurring together with Rhododendron aganniphum. The tree coverage and density of S. saltuaria and A. georgei significantly decreased from the closed forests to the shrubs, whereas the coverage of R. tanastylum increased. The geographical location of the study areas is shown in Figure 1. Information on the study sites is shown in Table 1.

Plot Setting and Vegetation Survey
In mid-August 2019, three 20 m × 20 m plots were set up at each of the five vegetation zones at each timberline. The replicate plots were placed at least 50 m apart perpendicular to the elevational gradient. A total of thirty plots were established in this study. Species and coverage of trees and shrubs in each plot were recorded. Tree individual numbers were also recorded.

Sample Collection
Soils were sampled in each plot in an S-shape manner to cover the entire plot [32]. After removing litter and organic layers, a soil drill (diameter, 50 mm) was used to collect 15 soil cores (0-20 cm) from each plot to increase sample representativeness [32]. The soil samples were sealed in sterilized bags, marked and transported to the laboratory. The soils were sieved through a 2 mm sieve to remove gravels and roots. The soils with same fresh weight were taken from the soil samples of the same plot and then evenly pooled to form 1 composite soil sample. A total of 30 composite soil samples (2 slope aspects × 5 vegetation zones × 3 repeat plots) were obtained. The composite samples were then divided into two parts: one part was stored at −80 • C for subsequent microbiological analysis and the other part was air-dried to analyze soil chemical characteristics. A volumetric soil sampler (100 cm 3 ) was used to collect an extra soil sample from each plot at a depth of 10 cm to analyze the soil physical characteristics.

Determination of Soil Physical and Chemical Properties
Soil porosity, soil bulk density and field holding water capacity were analyzed using the volumetric soil samplers following the methods described by State Forestry Administration [33]. The air-dried soil samples were sieved using a 0.15 mm sieve and then used for the analysis of soil pH, soil organic carbon, total nitrogen, total phosphorus, ammonium nitrogen and available phosphorus. The determination of these soil chemical indexes followed Lu's experimental manual [32]. An electrode method with a 1:5 soil:water ratio was used to determine soil pH using a pH meter (PHSJ-5, Leici, Shanghai, China). Soil organic carbon and total nitrogen contents were determined via the potassium dichromate oxidation and Kjeldahl methods. Total phosphorus and available phosphorus contents were determined via the molybdenum antimony ascorbic acid colorimetry method, whereas ammonium nitrogen content was determined via the indigo blue colorimetry method. PowerSoil ® DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, USA) was used to extract soil fungal genomic DNA from 0.5 g fresh soil sample following the manufacturer's protocol. The quality of genomic DNA was determined via electrophoresis using 1.2% agarose gel. The fungal ITS1 region was amplified using the universal primers ITS1F (5 -CTTGGTCATTTAGAGGAAGTAA-3 , coverage 89.8%) [34,35] and ITS2R (5 -GCTGCGTTCTTCATCGATGC-3 , coverage 93.7%) [35,36]. The PCR reaction system contained 4 µL of 5 × Fastpfu buffer, 2 µL of 2.5 mmol/L dNTPs, 0.8 µL of forward primer (5 mol/L), 0.8 µL of reverse primer (5 mol/L), 0.4 µL FastPfu Polymerase, 0.2 µL BSA and 10 ng template DNA, and was then supplemented with ddH 2 O to 20 µL. The PCR reaction conditions were: pre-denaturation at 95 • C for 3 min, 35 cycles × (denaturation at 95 • C for 30 s, renaturation at 55 • C for 30 s, extension at 72 • C for 45 s), extension at 72 • C for 10 min and incubation at 10 • C before sequencing. Each sample was run in triplicates. The PCR amplicons of different soil samples were homogenized to 10 nmol/L and mixed in equal volume for library construction.

Library Construction and Sequencing
The library was constructed using TruSeq ® DNA PCR-Free Sample Preparation Kit. The library was then quantified using qubit and Q-PCR methods. Library sequencing was conducted using Illumina NovaSeq 6000. Each sample was tagged based on the barcode sequence and PCR-amplified primer sequence. Reads of each sample were flashed (V1.2.11, http://ccb.jhu.edu/software/FLASH/ (accessed on 20 May 2023)) [37] and spliced after the interception of the barcode and primer sequence. The spliced sequence represented raw tags, which were filtered [38] to obtain high-quality tags (clean tags). Quality control for the tags was conducted using QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_ fastq.html (accessed on 20 May 2023)) [39]. The process was as follows: (1) tag interception: raw tags were cut from the first low-quality base site with continuous low-quality value (default quality threshold < 19) and the base number reached the set length (default length value, 3); (2) length filtering of tags: the data set of tags obtained by intercepting tags was further filtered to obtain a continuous high-quality base length less than 75% of the tags length. Tag sequences obtained after quality control were compared with the species annotation database to detect chimeric sequences (https://github.com/torognes/vsearch/ (accessed on 20 May 2023)) [40]. The final effective tags were obtained by removing chimeric sequences [41]. The details for specific data preprocessing and quality control are shown in Supplementary Table S1. The sequenced raw data were uploaded on the NCBI SRA (accession number; PRJNA785564) (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA785 564 (accessed on 20 May 2023)).

OTU Clustering and Taxonomic Annotation
High-quality sequences were clustered using UPARSE (V11, http://www.drive5 .com/uparse/ (accessed on 20 May 2023)) [42] into OTUs (operational taxonomic units) with 97% similarity threshold. The sequence with the highest frequency in an OTU was selected as the representative sequence of that OTU. The blast method in QIIME (V1.9.1) (http://qiime.org/scripts/assign_taxonomy.html (accessed on 20 May 2023)) [39] and UNITE (V8.3) databases (https://unite.ut.ee/ (accessed on 20 May 2023)) [43] were used to annotate OTUs taxonomically (kingdom, phylum, class, order, family and genus, species). Multiple sequence alignment was conducted using MUSCLE (V5, http://www.drive5.com/ muscle/ (accessed on 20 May 2023)) [44] to explore the phylogenetic relationship among all OTU representative sequences. The data for each sample were rarefied to the number of sequence reads obtained in the sample with the least amount of data (48,579 reads). After this, all the samples had saturating species accumulation curves (Supplementary Figure S1), indicating that the sequencing depth was sufficient and that fungal diversity in the soil samples was effectively covered by sequencing. Thus, subsequent fungal diversity and community structure analyses were conducted based on the rarefied data. The workflow for sequencing, quality control and OTU annotation is shown in Figure 2. FUNGuild is a functional annotation tool used to predict the functional composition of soil fungal communities of fungal amplicons obtained through high-throughput sequencing and other methods. However, it can only annotate information of fungal species on trophic-mode levels (pathotroph, symbiotroph and saprotroph). The modes can be further subdivided into multiple guides. Each guide comprises the species with similar absorption and utilization of environmental resources [45]. Herein, the ecological function of soil fungi was described based on the mode and guild levels.

Data Analysis
In order to determine differences between the two timberlines in soil fungi, the north and south timberlines were indicated by "1" and "2", respectively. Rarefaction curves were generated using R (V3.6.3) [46]. Venn diagrams were drawn using the R and VennDiagram (V1.7.1) package to show the quantitative distribution of soil fungal OTUs of different timberlines and vegetation zones [47]. The data were tested for homogeneity of variance before alpha diversity analysis. Levene's test was conducted using R and car (V3.0.12) to check whether the data variances of timberlines and vegetation zones were equal (p > 0.05) [48]. If necessary, the data were log-transformed to obtain equal variances. Two-way ANOVA was used to test differences in vegetation and soil factors, soil fungal diversity indices between the two timberlines and the five vegetation zones using R. The tested alpha diversity indices were: observed species, Shannon entropy of counts (Shannon), Simpson's index (Simpson), Chao1 richness estimator (Chao1), abundance-based coverage estimator (ACE), Faith's phylogenetic diversity metric (PD) and Good's coverage of counts (Goods_coverage). An alpha diversity package of the scikit-bio development team (http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.html#module-skbio.diversity .alpha (accessed on 20 May 2023)) was used to calculate the indices. The data are expressed as mean ± standard deviation (SD).
The R and ggplot2 (V3.3.5) package [49] were used to show the accumulative bar diagrams, representing the relative abundance of soil fungal phyla and predicted fungal modes of each vegetation zone. Heat maps were generated using the R and pheatmap (V1.0.12) package [50] to show the relative abundance of the top 10 fungal genera and fungal guilds of each vegetation zone. The sample differences in species complexity were evaluated using beta diversity analysis. Beta diversity of the weighted UniFrac distance was calculated using QIIME (V1.9.1). PERMANOVA tests (permutation = 999) based on Bray-Curtis distance were conducted using the R and vegan (V2.5.7) package to determine the soil fungal community structure differences [51]. The similarity in soil fungal

Data Analysis
In order to determine differences between the two timberlines in soil fungi, the north and south timberlines were indicated by "1" and "2", respectively. Rarefaction curves were generated using R (V3.6.3) [46]. Venn diagrams were drawn using the R and VennDiagram (V1.7.1) package to show the quantitative distribution of soil fungal OTUs of different timberlines and vegetation zones [47]. The data were tested for homogeneity of variance before alpha diversity analysis. Levene's test was conducted using R and car (V3.0.12) to check whether the data variances of timberlines and vegetation zones were equal (p > 0.05) [48]. If necessary, the data were log-transformed to obtain equal variances. Two-way ANOVA was used to test differences in vegetation and soil factors, soil fungal diversity indices between the two timberlines and the five vegetation zones using R. The tested alpha diversity indices were: observed species, Shannon entropy of counts (Shannon), Simpson's index (Simpson), Chao1 richness estimator (Chao1), abundance-based coverage estimator (ACE), Faith's phylogenetic diversity metric (PD) and Good's coverage of counts (Goods_coverage). An alpha diversity package of the scikit-bio development team (http://scikit-bio.org/docs/ latest/generated/skbio.diversity.alpha.html#module-skbio.diversity.alpha (accessed on 20 May 2023)) was used to calculate the indices. The data are expressed as mean ± standard deviation (SD).
The R and ggplot2 (V3.3.5) package [49] were used to show the accumulative bar diagrams, representing the relative abundance of soil fungal phyla and predicted fungal modes of each vegetation zone. Heat maps were generated using the R and pheatmap (V1.0.12) package [50] to show the relative abundance of the top 10 fungal genera and fungal guilds of each vegetation zone. The sample differences in species complexity were evaluated using beta diversity analysis. Beta diversity of the weighted UniFrac distance was calculated using QIIME (V1.9.1). PERMANOVA tests (permutation = 999) based on Bray-Curtis distance were conducted using the R and vegan (V2.5.7) package to determine the soil fungal community structure differences [51]. The similarity in soil fungal community structure of different timberlines and vegetation zones was analyzed using a principal coordinate analysis (PCoA) in R and ggplot2 package.
The relationships among fungal, vegetative and edaphic factors were assessed using redundancy analysis (RDA) and Pearson correlation. The data were standardized using zscore for data uniformity before RDA analysis. Multicollinearity factors (variance inflation factor > 10) for vegetation and soil were excluded from analysis. RDA (permutation = 999) was conducted based on the standardized data. Canoco (V5.0) was used for RDA analysis [52]. Pearson correlation was performed in R (V3.6.3) [46].

Vegetation and Soil Characteristics
The two-way ANOVA ( Table 2) showed that tree density (F = 5.21, p = 0.03) at the north timberline was significantly higher than that at the south timberline. Shrub coverage (F = 11.72, p < 0.01) in the shrublands surpassed that in the forests and forest lines. Moreover, tree coverage (F = 111.73, p < 0.01) and density (F = 25.14, p < 0.01) were considerably greater in the forests and lower in the shrublands compared to the other four vegetation zones, respectively. The interaction between the timberline site and vegetation zone had a significant effect on shrub coverage (F = 3.18, p = 0.04). Specifically, at the south timberline, shrub coverage in the forest was lower than that in the other vegetation zones, and increased gradually from forest to shrubland. Tree coverage and density, however, were highest in the forest compared to other zones, exhibiting a decline from forest to shrubland. The trend on the north timberline was similar, but shrub coverage in the shrubland was significantly higher than that in the forest and forest line. Nonetheless, tree coverage and density were substantially greater in the forest than in the other vegetation zones.

Soil
The two-way ANOVA (Table 3) indicated that the soil field water-holding capacity (F = 6.14, p = 0.02) at the north timberline was significantly higher than that at the south timberline, while the soil bulk density (F = 3.16, p = 0.04) in the forests and tree lines was significantly higher than that in the tree species lines. There was also a significant interaction between the timberline site and vegetation zone regarding soil bulk density (F = 4.00, p = 0.02). Specifically, the soil bulk density in the forest at the south timberline displayed a notably higher value than that in other vegetation zones, while the soil field water-holding capacity in the tree species line surpassed both the forest and shrubland by a significant margin. At the north timberline, only the tree line's soil bulk density displayed a significant elevation from other vegetation zones.
The soil pH (F = 45.21, p < 0.01) and available phosphorus content (F = 27.80, p < 0.01) were significantly higher at the south timberline when compared to the north timberline, whereas the soil ammonium nitrogen content (F = 13.84, p < 0.01) was significantly higher at the north timberline ( Table 3). The available phosphorus content (F = 3.56, p = 0.02) was significantly higher in the forests than in the other three vegetation zones, except for tree lines. Moreover, the interaction between the timberline site and vegetation zone had a significant effect on soil pH (F = 16.30, p < 0.01), total phosphorus (F = 5.04, p = 0.01) and ammonium nitrogen (F = 3.83, p = 0.02). The pH levels in the south timberline's forest and forest line were notably distinct from other vegetation zones, while the soil organic carbon concentration in the shrubland was considerably higher than that in the forest. Furthermore, the total phosphorus content of the tree line was significantly greater than shrub. On the other hand, the north timberline exhibited lower soil pH and available phosphorus levels in the forest compared to other vegetation zones. Additionally, the forest and tree line had a significantly higher soil ammonium nitrogen content than the tree species line and shrub.
The soil C/N (F = 87.27, p < 0.01) was significantly higher at the north timberline than at the south timberline, while shrublands had a significantly higher C/N (F = 36.02, p < 0.01) than the other vegetation zones (Table 3). Additionally, soil C/N (F = 17.68, p < 0.01) and C/P (F = 3.63, p = 0.02) were significantly affected by the interaction between the timberline site and vegetation zone. Notably, the C/N and C/P of shrub at the south timberline were significantly higher than those of other vegetation zones, while the N/P of shrub was significantly higher than that of the forest and tree line. Finally, the C/N of shrub at the north timberline was significantly higher than that of other vegetation zones, except for the forest line.

Soil Fungal Diversity
Alpha diversity indices were not significantly different between the timberlines and vegetation zones, and there was no interaction between these factors (F = 1.333-3.568, p = 0.073-0.262, Table 4, Figure S1).

Relative Abundance of Soil Fungal Phyla and Genera
Soil fungi were classified into 12 phyla, 43 classes, 104 orders, 201 families, 343 genera and 311 species. Notably, soil fungal communities of the two timberlines were dominated by Basidiomycota (12-80%) and Ascomycota (10-57%) (Figure 3). An obvious change in relative abundance from Basidiomycota to Ascomycota was observed with higher elevation at the north timberline. The relative abundance of Basidiomycota significantly decreased from forest (71%) to shrub (19%). On the contrary, the relative abundance of Ascomycota increased from forest (16%) to shrub (53%) (Figure 3). The relative abundance of Ascomycota was higher compared with that of Basidiomycota at all vegetation zones of the south timberline. Across all vegetation zones, the relative abundance of Mortierellomycota was also rather high (7-25%), especially at the south timberline.

Relative Abundance of Soil Fungal Phyla and Genera
Soil fungi were classified into 12 phyla, 43 classes, 104 orders, 201 families, 343 genera and 311 species. Notably, soil fungal communities of the two timberlines were dominated by Basidiomycota (12-80%) and Ascomycota (10-57%) (Figure 3). An obvious change in relative abundance from Basidiomycota to Ascomycota was observed with higher elevation at the north timberline. The relative abundance of Basidiomycota significantly decreased from forest (71%) to shrub (19%). On the contrary, the relative abundance of Ascomycota increased from forest (16%) to shrub (53%) (Figure 3). The relative abundance of Ascomycota was higher compared with that of Basidiomycota at all vegetation zones of the south timberline. Across all vegetation zones, the relative abundance of Mortierellomycota was also rather high (7-25%), especially at the south timberline.  The soil fungal communities at the two timberlines comprised different dominant fungal genera (Figure 4). Archaeorhizomyces (8-30%) of Ascomycota was the most dominant genus across all vegetation zones and was more abundant in the south timberline ecotone. The relative abundance of Russula of Basidiomycota ranged from 8% to 30% at the north timberline (except in the shrub zone (0.43%)), but only from 0.18% to 0.34% at the south timberline. The soil fungal communities at the two timberlines comprised different dominant fungal genera (Figure 4). Archaeorhizomyces (8-30%) of Ascomycota was the most dominant genus across all vegetation zones and was more abundant in the south timberline ecotone. The relative abundance of Russula of Basidiomycota ranged from 8% to 30% at the north timberline (except in the shrub zone (0.43%)), but only from 0.18% to 0.34% at the south timberline.

Soil Fungal Community Similarity
The PERMANOVA showed significant differences in soil fungal community composition (based on OTUs) between the two timberlines (R 2 = 0.109, F = 3.441, p = 0.001) and among the vegetation zones (R 2 = 0.082, F = 2.488, p = 0.002). These differences were verified through the PCoA explaining 41.22% of the total variation ( Figure 5). The soil fungal communities were distinguished by the PC1 axis (27.23%) based on timberline ( Figure 5). The fungal community compositions in the closed forest, forest line and tree line at the north timberline (all with a dominance of A. georgei) were similar, clustering at the left of the axis, whereas fungal communities at the south timberline were more similar to vegetation types with few or no A. georgei (tree species line and shrub) at the north timberline, clustering at the right of the axis. PC2 (13.99%) indicated the difference in fungal communities across the vegetation zones at the southern timberline. The fungal communities of forest and shrub at the south timberline were similar and were mainly distributed in the lower half of PC2, the fungal communities at the tree species line clustered in the middle and similar communities from other vegetation zones at the south timberline were scattered across PC2. However, the soil fungal communities of the forest, forest line and tree line at the north timberline mainly clustered in the middle of PC2.

Ecological Function of Soil Fungi
The Saprotroph mode (25-36%) and Soil_Saprotroph guild (16-30%) dominated at the south timberline, with little variation across the vegetation zones, whereas the Symbiotroph mode and Ectomycorrhizal guild were characteristic of the north timberline together with the Saprotroph mode and Soil_Saprotroph guild (Figures 6 and 7). The relative abundance of the Symbiotroph mode (56-3%) and Ectomycorrhizal guild (41-1%) decreased across the vegetation zones at the north timberline, whereas the relative abundance of the Saprotroph mode (11-35%) and Soil_Saprotroph guild (8-30%) increased across the vegetation zones (Figures 6 and 7).

Ecological Function of Soil Fungi
The Saprotroph mode (25-36%) and Soil_Saprotroph guild (16-30%) dominated at the south timberline, with little variation across the vegetation zones, whereas the Symbiotroph mode and Ectomycorrhizal guild were characteristic of the north timberline together with the Saprotroph mode and Soil_Saprotroph guild (Figures 6 and 7). The relative abundance of the Symbiotroph mode (56-3%) and Ectomycorrhizal guild (41-1%) decreased across the vegetation zones at the north timberline, whereas the relative abundance of the Saprotroph mode (11-35%) and Soil_Saprotroph guild (8-30%) increased across the vegetation zones (Figures 6 and 7). the south timberline, with little variation across the vegetation zones, whereas the Sym-biotroph mode and Ectomycorrhizal guild were characteristic of the north timberline together with the Saprotroph mode and Soil_Saprotroph guild (Figures 6 and 7). The relative abundance of the Symbiotroph mode (56-3%) and Ectomycorrhizal guild (41-1%) decreased across the vegetation zones at the north timberline, whereas the relative abundance of the Saprotroph mode (11-35%) and Soil_Saprotroph guild (8-30%) increased across the vegetation zones (Figures 6 and 7).

Relationship between Soil Fungal Community Characteristics and Vegetation and Soil Factors
The first two axes of the RDA explain 25.95% of the total variation (16.80% for the first axis and 9.15% for the second axis) (Figure 8 and Table 5). The observed species richness, Chao1, ACE and PD were negatively correlated with the vegetation zone (altitude), shrub coverage, soil field holding-water capacity, ammonium nitrogen and C/N, but were positively associated with the soil pH and bulk density. Negative correlations of the timberline site with the Shannon and Simpson indexes were observed. Ascomycetes, saprophytic and soil saprophytic fungi were negatively correlated with the timberline site, tree factors and soil ammonium nitrogen, and positively correlated with the soil pH, whereas the opposite was the case for Basidiomycetes, symbiotic and ectomycorrhizal fungi.

Relationship between Soil Fungal Community Characteristics and Vegetation and Soil Factors
The first two axes of the RDA explain 25.95% of the total variation (16.80% for the first axis and 9.15% for the second axis) (Figure 8 and Table 5). The observed species richness, Chao1, ACE and PD were negatively correlated with the vegetation zone (altitude), shrub coverage, soil field holding-water capacity, ammonium nitrogen and C/N, but were positively associated with the soil pH and bulk density. Negative correlations of the timberline site with the Shannon and Simpson indexes were observed. Ascomycetes, saprophytic and soil saprophytic fungi were negatively correlated with the timberline site, tree factors and soil ammonium nitrogen, and positively correlated with the soil pH, whereas the opposite was the case for Basidiomycetes, symbiotic and ectomycorrhizal fungi.
but were positively associated with the soil pH and bulk density. Negative correlations of the timberline site with the Shannon and Simpson indexes were observed. Ascomycetes, saprophytic and soil saprophytic fungi were negatively correlated with the timberline site, tree factors and soil ammonium nitrogen, and positively correlated with the soil pH, whereas the opposite was the case for Basidiomycetes, symbiotic and ectomycorrhizal fungi. However, Pearson correlation analysis further showed that the observed species richness, PD and Shannon and Simpson indexes were significantly positively associated with soil pH, but the Shannon and Simpson indexes were obviously negatively associated with soil ammonium nitrogen ( Figure 9C,F). At the north timberline, the coverage and density of trees were significantly negatively related to Ascomycetes, saprophytic and soil saprophytic fungi, and significantly positively related to Basidiomycetes, symbiotic and ectomycorrhizal fungi ( Figure 9A). The shrub coverage had contrary linkages compared with the tree factors. Therefore, the composition and function of soil fungal communities at the north timberline changed with the tree and shrub coverage and tree density. At the north timberline, the soil C/N and pH were negatively correlated with Basidiomycetes, symbiotic and ectomycorrhizal fungi, whereas ammonium nitrogen was positively correlated; the opposite was the case for Ascomycetes, saprophytic and soil saprophytic fungi, especially regarding soil C/N ( Figure 9D). Unlike the north timberline, there was no significant correlation between soil fungi, vegetation and soil factors at the south timberline ( Figure 9B,E). The Pearson correlation also indicated that the Chytridiomycota relative abundance significantly negatively responded to shrub coverage, soil field water-holding capacity, C/N and ammonium nitrogen, but significantly positively responded to soil bulk density and pH (Figures 8 and 9F). south timberline ( Figure 9B,E). The Pearson correlation also indicated that the Chytridiomycota relative abundance significantly negatively responded to shrub coverage, soil field water-holding capacity, C/N and ammonium nitrogen, but significantly positively responded to soil bulk density and pH (Figures 8 and 9F).

Discussion
The findings of the present study show that the composition and function of soil fungal communities between the north and south timberlines of Sejila Mountain had a strong site effect. The soil fungal community characteristics at the north timberline regularly changed with the tree coverage and density, but the tree species line and shrub were more similar to the vegetation zones of the south timberline regarding their traits. This indicates that the A. georgei presence also had an impact on the soil fungal community composition and function.

Distribution of Mycorrhizal Fungi and Its Relationship with Vegetation
The distribution of ectomycorrhizal fungi is strongly influenced by the tree species and their distribution [53]. In this study, A. georgei is a prominent species at the northern timberline ecosystem due to its preference for shade and moisture [54]. The ectomycorrhizal fungi predominantly occurred at the northern timberline and their abundance had a positive correlation with the density and coverage of A. georgei. The family Russulaceae is the dominant group of ectomycorrhizal fungi. Russula is the type genus of Russulaceae under Basidiomycetes. Its members can form ectomycorrhizae with trees of Larix, Picea, Pinus, Fagaceae, Salicaceae and Tiliaceae [55]. A study on Russulaceae fungi along the altitude gradient of 2600-4500 m in Sejila Mountain also showed that 70% of Russulaceae species are distributed under A. georgei forest, but they rarely appear in the high-altitude Rhododendron shrubs [24,25]. Several research studies have confirmed that Russula can form an ectomycorrhizal symbiosis with tree species of the Abies genus, such as Abies lasiocarpa [56], Abies balsamea [57], Abies koreana [58] and Abies firma [59]. This supports the notion that A. georgei is their potential host too. The number of A. georgei individuals has a great effect on the root density in soil, further affecting the relative abundance of ectomycorrhizal fungi [60]. This is the reason for the positive synergistic change in the relative abundance of Basidiomycetes, Russula, symbiotic and ectomycorrhizal groups with the tree coverage and density of A. georgei. This is basically consistent with our hypothesis of the host effect. Four species in Russula genus were identified in this study, but their relative abundances were extremely low. This implies that the high relative abundance of Russula was due to the presence of other unidentified species. The significantly negative correlation between the Shannon index and tree density at the north timberline indicates that the ectomycorrhizal fungi may inhibit other groups of fungi at the timberline. Owing to their remarkable cold tolerance, ectomycorrhizal fungi are considered to be crucial for seedling establishment, nutrient absorption and the distribution expansion of A. georgei [61]. In summary, the presence of A. georgei led to systematic changes in the soil fungal community characteristics.

Distribution of Saprophytic Fungi and Its Relationship with Vegetation
Soil saprotrophic fungal communities at the timberlines of Sejila Mountains are largely populated by the Archaeorhizomyces fungi. In this study, regarding the symbiosis relationship between S. saltuaria and mycorrhizal fungi, which has not been reported previously, here, it was also found that the low relative abundance of mycorrhizal fungi at the south timberline was dominated by S. saltuaria. However, there were plenty of Archaeorhizomyces belonging to the class Archaeorhizomycetes of Ascomycetes at the south timberline [62].
It is not clear enough whether the taxonomic status and ecological function of Archaeorhizomyces are a result of difficulties in isolation and culture. Only two species have been systematically described: Archaeorhizomyces finlayi [63] and Archaeorhizomyces borealis [62]. Archaeorhizomyces are non-mycorrhizal rhizosphere-associated fungi with saprotrophic activity. For example, in culture, A. finlayi grows slowly on both glucose and cellulose as a sole carbon source, indicating that it may be involved in decomposition [63]. At alpine timberlines and in other stressful ecosystems [22,[64][65][66], Archaeorhizomyces is a common fungal group, similar to our study. Archaeorhizomyces also often dominate soils of rainforest [67], temperature forests [68], boreal forests [69] and arctic tundra [70]. This explains the strongly positive associations among Ascomycetes, Archaeorhizomyces, saprophytic and soil saprophytic groups in this study ( Figure S3C).
At the north timberline, Archaeorhizomyces were negatively correlated with Basidiomycetes, Russula, symbiotic and ectomycorrhizal groups and A. georgei coverage and density, but positively correlated with shrub coverage ( Figure S3A). This may be related to the alleviative competition with Russula with the decrease in A. georgei hosts at the north timberline. Although the precise ecological niches of Archaeorhizomyces remain unknown [63], the higher relative abundance of Archaeorhizomyces within the ecotone implies that fungi of Archaeorhizomyces are more suitable for this transitional habitat lacking tree cover. One species of Archaeorhizomyces was identified in the current study, and the relative abundance was low. This also indicates that the high relative abundance of Archaeorhizomyces is attributed to the presence of other unidentified species. Overall, it is highly probable that Archaeorhizomyces collaborate with ectomycorrhizal fungi to facilitate the circulation of materials and energy at the timberlines of Sejila Mountain.

Relationship between Fungal Distribution and Soil Factors
In this study, timberlines on different slopes and diverse vegetation zones at varying altitudes show considerable variations in both physical and chemical properties of the soil [71,72]. These variations play a significant role in shaping the traits of soil fungal communities, together with timberline vegetation [73]. At the north timberline, due to a deficiency in base cations in the litters, A. georgei litters may lead to soil acidification [74], resulting in a significant negative correlation between soil pH and tree density, coverage and ectomycorrhizal fungi. Previous studies have confirmed that ammonium nitrogen promotes the symbiosis between ectomycorrhizal fungi and trees [75][76][77]. The soil total nitrogen pool at the vegetation zones with more trees at the north timberline ensured a continuous supply of ammonium nitrogen. A decrease in the relative abundance of ectomycorrhizae at the tree species line and shrub may be related to a significant reduction in N pool supply, in addition to a reduction in host density. The soil stoichiometric ratios indicate the mineralization capacity for soil carbon, nitrogen and phosphorus [78,79]. The soil C/N at the north timberline increased with a decrease in tree coverage and density, whereas N/P decreased. This indicates that soil nitrogen tended to be mineralized at the vegetation zones with more trees, which is conducive for the accumulation of soil ammonium nitrogen [80]. On the contrary, soil nitrogen mineralization weakened at the tree species line and shrub, and ammonium nitrogen accumulation was low. It is noteworthy that, at the south timberline, there is no obvious correlation between soil fungal community characteristics and vegetation and soil factors. We speculate that there may be other potential factors that play an important role in regulating soil fungal communities, which needs further investigation.

Conclusions
In the present study, the soil fungal diversity, community structure and ecological function of the north and south timberlines of the Sejila Mountain in China's southeastern Tibet along a 200 m elevation gradient were explored. The findings show that the site effect of timberlines and A. georgei presence exerted apparent influences on the soil fungal community structure and function at the two timberlines. These soil fungal community traits were related to coverage and density, soil pH and ammonium nitrogen at the north timberline, whereas they had no associations with the vegetation and soil factors at the south timberline.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/jof9050596/s1. Figure S1: Rarefaction curves for soil fungal communities at the south and north timberlines of Sejila Mountain; Figure S2: Venn diagrams for soil fungal OTUs at the south and north timberline of Sejila Mountain; Figure S3: Pearson correlation between vegetation factors and main fungal genera, phyla, modes and guilds at the timberlines of Sejila Mountain. Table  S1: Results for data preprocessing and quality control.

Data Availability Statement:
The sequenced raw data were uploaded on the NCBI SRA with an accession number PRJNA785564.