Different patterns of belowground fungal diversity along altitudinal gradients with respect to microhabitat and guild types

Background Fungi are key components of belowground ecosystems with various ecological roles in forests. Although the changes in the richness and composition of belowground fungi across altitudinal gradients have been widely reported, only a few studies have focused on the microhabitat types along altitudinal gradients. Results Here, we analyzed the effect of altitudinal gradients on the fungal communities of belowground microhabitats (roots and soil) and their different ecological roles. Root and nearby soil samples were collected from 80 Pinus densiora trees in 16 localities at various altitudes across the country, and their community composition was analyzed using DNA metabarcoding of the ITS2 region on the Illumina MiSeq platform. We found that ectomycorrhizal (ECM) fungi and ericoid mycorrhizal fungi were more abundant in the soil, whereas the relative abundance of endophytic and saprotrophic fungi was higher in the roots. Altitude negatively affected the species richness of root-inhabiting and ECM fungi but did not inuence that of non-ECM fungi. However, the composition of ECM fungi was less inuenced by altitude than that of non-ECM fungi. Conclusion Our results demonstrate that microhabitat types and altitudinal gradients affect the diversity, composition, and function of fungal communities associated with a single plant host and contribute to a better understanding of plant-associated fungal communities.

Furthermore, in some studies on multiple vegetation with different hosts, the identity of vegetation had a greater effect on the community than abiotic factors [28,33]. Only a few studies have focused on a single host that minimize the effect of host identity on the altitudinal patterns of the fungal community [29,[34][35][36]. However, these studies mostly investigated the composition changes in the fungal community and not the species richness.
To our knowledge, there are only a few reports on how altitudinal gradients are related to fungi in different belowground microhabitats [37] or different guilds [38]. Although these studies provided valuable insights of the relationship between altitudinal gradients and fungal communities, the relationship between altitudinal gradients and the fungal community of different microhabitats or guilds remain unknown.
To this end, we examined how altitudinal gradients affect belowground fungal communities in different microhabitats (root and soil) and different guilds (ECM and other non-ECM guilds). We selected the widely distributed ectomycorrhizal tree Pinus densi ora (red pine) as a target species to minimize the host effect and focus on altitude and microhabitat [39]. In the present study, the composition and diversity of fungal communities were compared between microhabitats and altitudinal gradients, and the results were interpreted according to the functional characteristics of fungal guilds. We expected a signi cant change in community compositions across altitudinal gradients in both microhabitats and guilds, especially in the fungal communities that are not directly associated with P. densi ora. In species richness level, we hypothesized that plant-associated fungi (i.e. ECM or root-inhabiting fungi) would be more affected by altitudinal gradients due to change in host tness along the altitude. Our ndings would offer a more comprehensive understanding of the fungal interactions between root and soil microhabitats. A clear understanding of these changes can help to predict the reactions of microbial diversity to changing environments and climate [26,40].

Sample Collection
Sampling was conducted in 16 well-protected, P. densi ora-dominant South Korean forests from 2019-2020 ( Fig 1A, Table 1). The sampling sites were sorted into three groups (high, mid, and low) according to their altitude with respect to the vegetation distribution patterns of the nearby area [41,42]. We randomly selected ve mature P. densi ora trees with similar diameter at breast height (DBH) per site that were at least 20 m apart [43]. The trees were estimated to be 30-40 years old with a DBH of 32.36 ± 12.17 cm. To avoid including the roots from other plants, roots were manually traced from the trunk of P. densi ora after removing the litter layers around the trees. In total, root and soil samples were collected from 80 P. densi ora trees. For each tree, we sampled at least two lateral roots with mycorrhizal root tips from different directions and two soil samples from each root sample in different directions.
The root and soil samples were stored in plastic bags and transported in an icebox. In the laboratory, the soil samples were sieved, and those from the same tree were pooled into one sample and stored at −20°C until DNA extraction. Root samples were stored at 4 °C prior to surface sterilization and DNA extraction. Sampling was permitted and accompanied by the Korea National Arboretum Authority.
Hereafter, we refer to the two sample types (root and soil) as "microhabitats" or "belowground microhabitats," which will be treated as categorical variables.

DNA extraction
For the root samples, we removed the soil and debris from the surface by gently shaking them in distilled water to prevent damage to the root. After washing, all ne roots with fresh-looking root tips were separated using sterilized scissors and the surface of the separated roots were sterilized by submerging in 3% hydrogen peroxide for 1 min. Surface-sterilized root samples were rinsed three times with sterile distilled water to remove possible contaminants from the soil. Root sterilization was conducted within a week of sampling. After sterilization, we pooled the root samples from the same tree into one sample and air-dried them on sterilized lter paper placed on a clean bench overnight. Dried ne roots were nely ground in an autoclaved mortar using liquid nitrogen. products were separated on a 1% agarose gel and puri ed using an ExpinTM PCR SV kit (GeneAll Biotechnology, Seoul, South Korea). The second round of PCR was conducted using the same protocol, but cycles were reduced to 10, and the unique identi er sequences were attached following the Nextera XT index kit protocol (Illumina, San Diego, CA, USA). Products from the second PCR were checked and puri ed as described above. The concentration of each amplicon library was measured using NanoDrop2000 (Thermo Fisher Scienti c, Waltham, MA, USA) and then pooled in equimolar quantities. Sequencing was performed on an Illumina MiSeq platform at Macrogen (Seoul, South Korea).

Sequence data analysis
Raw sequencing data were processed using QIIME v. 1.8.0 [45]. After merging paired-end sequences and quality ltering (Q ≥ 20), sequences were clustered into operational taxonomic units (OTUs) with 97% similarity using Vsearch v. 2.6.2 [46]. Chimeras were identi ed and removed using the 'uchime_denovo' method in Vsearch [46] with the reference database from UNITE 8.2 [47] and Seoul National University Fungal Collection (SFC). Representative sequences of OTUs were assigned to the most abundant sequences. For taxonomic assignment, the BLAST algorithm and reference database from UNITE 8.2 [47] and SFC were used following the criteria of Tedersoo et al. [20]. Non-fungal sequences and rare OTUs (≤10 reads or ≤5 samples) were ltered out to remove potential sequencing errors [48]. The postclustering curation method was performed in the default setting to merge co-occurring 'daughter' OTUs with more abundant 'parent' OTUs [49] in addition to a manual BLAST search in NCBI Genbank for better taxonomic resolution. Among the 160 samples, seven were discarded due to low sequence reads (two root samples and ve soil samples; root: GY03P-4R, GW04P-3R; soil: CC01P-3S, GY01P-4S, GY03P-5S, GB02P-2S, and GW02P-3S).
For the trophic mode assignment, we used the FUNGuild database [50] and manually checked the guild classi cation of OTUs based on the literature and eld experience. We used OTUs with proper reference or con dence rankings of 'probable' or 'highly probable' for further analysis of the ecological guilds. Then, we sorted the OTUs into six guilds (ectomycorrhiza, ericoid mycorrhiza, endophytes, saprotrophs, saprotrophs + endophytes, and saprotrophs + plant pathogens) and rare ed the samples to a minimum number of sequences (total: 5,600; ECM: 1,100; non-ECM: 1,000). Alpha-diversity indices (the number of OTUs, chao1 richness, Shannon diversity, and good coverage) were calculated in QIIME for each functional group.
Statistical analysis R v.3.6.1 was used for further statistical analysis [51]. The relationship between alpha-diversity indices and altitude gradients was examined using the Pearson correlation method with the cor.test function and linear regression test with the lm function in R. We compared the relative abundances (sequence reads of ECM OTUs/sequence reads of all OTUs in the sample) of the root and soil samples from each altitudinal range using t-test. Linear discriminant analysis effect size (LEfSe) analysis was conducted to discover taxa that are found in speci c microhabitats and/or altitudinal range using Huttenhower lab galaxy site with default options (LDA > 2.0, p < 0.05; https://huttenhower.sph.harvard.edu/galaxy/) [52]. The differences between communities were analyzed using the Jaccard dissimilarity index and assessed using non-metric multidimensional scaling (NMDS) methods in the phyloseq R package [53].
Permutational multivariate analysis of variance (PERMANOVA) statistical tests were performed to determine the differences in community composition with sampling sites, altitude, and microhabitat (root and soil) using 'adonis' function in vegan R package with 999 permutations and Jaccard dissimilarity index [54]. UPGMA (average linkage clustering) dendrograms were constructed based on the Jaccard distance to identify the connection of root and soil samples. To determine the changes in the fungal community composition along altitudinal and geographic distances (pairwise distance between two sampling sites), Mantel correlations were calculated using Jaccard dissimilarity metrics with Spearman's rank correlation and 9,999 permutations in the vegan R package [54]. As our sampling was conducted at multiple sites across various altitudes, we added geographic distance in the analysis, along with altitude, to compensate for the effect of geographic distance on community composition. Geographic distances between sampling sites were calculated using the GPS coordinates. The R package ggplot2 was used to visualize the results [55].

Results
Overall characteristics of fungal communities The fungal community structure associated with P. densi ora was determined for 78 root samples and 75 soil samples from 16 locations in South Korea after removing samples with insu cient sequence reads ( Figure 1A; Table 1). In total, 11,150,431 rRNA sequence reads were obtained after quality ltering, and 10,712,967 reads after ltering singletons and non-fungal sequences. The average length of sequence reads was 386.01 ± 48.18 base pairs. We found 12,927 OTUs after post-clustering, and 2,206 OTUs remained after ltering rare OTUs (≤10 reads or ≤5 samples) and normalization. From each microhabitat, 1 831 OTUs and 1 928 OTUs were found in the roots and soil, respectively. On average, 195.86 ± 53.13 OTUs and 222.76 ± 71.42 OTUs were found in each root and soil sample, respectively. Good's coverage ranged from 0.978 to 0.998 with an average of 0.989 ± 0.003, indicating that our sequencing depth was su cient to represent most of the fungal OTUs (Additional le 1: Table S1).

Fungal community composition
The composition of belowground fungi varied between the root and soil microhabitats. At the genus level, Tomentella, Geminibasidium, Suillus, Oidiodendron, and Cortinarius were the most abundant taxa ( Figure  1B, Additional le 2: Table S2). Among the ectomycorrhizal genera, Suillus and Phellodon showed higher relative abundance in roots than in soil. In contrast, the relative abundance of Tomentella, Oidiodendron, and Cortinarius was higher in the soil than in root. Meanwhile, non-ectomycorrhizal fungi also showed different relative abundances according to the microhabitats. The genera Mycena, Trichoderma, and Basidiodendron were more abundant in roots than in soil, whereas Gemnibasidium, Umbelopsis, and Oidiodendron were more abundant in soil (Additional le 2: Table S2). Among the plant pathogenic fungi, Venturia was mostly found in root microhabitats but in relatively low abundance compared to other guilds. The relative abundance of Geminibasidium, Umbelopsis, and Mycena showed signi cant differences across microhabitats according to LefSe analysis ( Figure 1B, Additional le 3: Figure S1).
The composition of belowground fungal communities also changed across altitudinal gradients. Suillus was more abundant in root samples collected at high altitudes, whereas Geminibasidium was more abundant in soil samples from low altitudes. Several genera, such as Mortierella and Tricholoma, showed opposite altitudinal patterns between root and soil samples (Fig. 1B, Additional le 2: Table S2). According to the results of LefSe analysis, in both microhabitats, Trichoderma sp. (OTU 31) was found at high altitudes, whereas Geminibasidium and Mycena were signi cantly abundant at low altitudes. Among ectomycorrhizal fungi, Astraeus and Cortinarius were the discriminating taxa of root and soil samples respectively from low altitudes (Additional le 4-5: Figure S2-3).
In terms of ecological function, saprotrophs and ECM were the major guilds in both microhabitats.
Saprotroph (p < 0.01), endophytes (p < 0.001), and saprotroph + plant pathogens (p < 0.001) were more abundant in roots, whereas the ECM (p < 0.001) and ericoid mycorrhiza (p < 0.001) were more abundant in the soil (Fig 2A). The relative abundance of saprotroph + endophytes was not signi cantly different between the microhabitats. Relative abundance of Saprotroph + endophyte guild was higher in the soil than in the root at low altitudes, and the pattern was opposite at high altitudes (p < 0.05). Similarly, the relative abundance of ECM was signi cantly higher in the soil than in the roots only at high altitudes (Fig.  2B, p < 0.001).

Change in fungal diversity across microhabitats and altitudinal gradients
We found that microhabitat type had a signi cant effect on the community composition of belowground fungal communities (Fig. 3, Table 2, Additional le 6: Table S3). However, their in uence differed according to the ecological groups. According to the results of the PERMANOVA test, the effects of microhabitat were higher in non-ectomycorrhizal fungi (R 2 = 0.083) than in ectomycorrhizal fungi (R 2 = 0.018) or the entire fungal communities (R 2 = 0.048). UPGMA analysis showed that the composition of fungal communities between the two microhabitats was much more similar in ECM fungal communities than in non-ECM fungal communities (Additional le 7: Figure S4).
Altitude had a signi cant effect on community composition, but the effect was different between ECM and non-ECM fungal communities (Fig. 3, Additional le 6: Table S3). The community of non-ECM fungi was more strongly in uenced by altitude than that of ECM fungi, regardless of the microhabitat type (Additional le 6: Table S3). Similar to the results of the PERMANOVA test, those of the Mantel test revealed that altitude exhibited signi cant distance-decay patterns with fungal communities. The negative correlation between altitudinal difference and fungal community dissimilarity was stronger in the soil than in the roots in all groups (Additional le 8: Figure S5). In contrast, the similarity of fungal communities was more negatively correlated with geographic distances in roots than in soil (Additional le 9: Figure S6). Overall, both microhabitat and altitude were signi cant factors that affected the fungal community composition. In most microhabitats, non-ectomycorrhizal fungal communities were more affected by altitudinal differences than the total fungal or ectomycorrhizal fungal communities.
The relationship of the alpha-diversity indices of fungal communities associated with P. densi ora and altitude varied according to microhabitats and ecological groups (Fig. 4). In roots, signi cant negative relationships were found between alpha-diversity indices and altitude in both groups. However, in soil, different patterns were observed between the ECM and non-ECM communities. Although a signi cant negative correlation between altitude and alpha-diversity of ectomycorrhizal and overall fungal communities was detected, non-ECM communities exhibited no signi cant relationships with altitude. *calculated based on Jaccard distance and relative abundance at the OTU-level with 999 permutations (*p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001).

Discussion
We found that altitudinal gradients and microhabitats were both signi cant factors affecting the richness, diversity, and composition of fungal communities. However, their in uence and pattern differed between both of microhabitats and ECM / non-ECM communities. To our knowledge, this is one of the few studies that covered the altitudinal changes in fungal communities across different microhabitats and guilds. It is important to note that in contrast to previous studies, ECM fungal richness showed signi cant change along altitudinal gradients.
Similar to the ndings of previous studies [37,56], we observed that microhabitat signi cantly affects belowground fungal community composition. However, in our study, the difference between microhabitats was much smaller in ECM communities than in non-ectomycorrhizal communities. The similar composition of ECM communities is likely due to the ECM mycelia in the soil as most ectomycorrhizal fungi form extraradical hyphae in the soil environment around colonized roots [2]. Although ectomycorrhizal fungi in both microhabitats mostly overlapped, species richness and diversity were consistently higher in the soil than in the root. Several factors may explain the differences between microhabitats. First, inactive spore banks or relic DNA from dead mycelia may in uence these differences in alpha diversity. Spore banks are formed by spores from nearby fruiting bodies that are spread by wind. Dormant spores stored in the soil can initiate associations if the right host plant is present [57] and can survive in soil for several years even in the absence of the right host [58]. Similarly, DNA fragments from dead fungi have been detected for several years [59]. Because we could not differentiate between active and inactive fungi using NGS, DNA from inactive fungi might have contributed to higher ectomycorrhizal fungal richness in soil in our data. Second, ECM fungi colonizing non-ectomycorrhizal hosts [60][61][62] or herbaceous plants [7] have been reported in previous studies. Mycelia from other plant hosts could have also in uenced ECM richness in soil.
On the other hand, the composition of non-ECM fungi exhibited higher variation than ECM fungi between microhabitats. The clear functional separation may explain this difference. Most non-ECM genera that were abundant in roots were endophytes or plant pathogens, such as Phialocephala, Mycnea, and Venturia. Phialocephala is a well-known endophyte found in the roots of various host trees, including P. densi ora [63]. On the other hand, several in vitro studies suggested that Mycena can colonize ECM plant roots and promote growth of host plant [64,65]. Venturia is a plant pathogenic fungus that is widely distributed in the Northern temperate area [66]. However, the presence of Venturia in the roots of P. densi ora has not been reported previously; therefore, further studies are needed to understand their relationship. In the soil, ericoid mycorrhizal fungi from Oidiodendron were more abundant than in the root. The presence of Oidiodendron in soil around the roots of ectomycorrhizal trees may be attributed to nearby understory vegetation or their saprobic activity in soil [67].
In addition, compositional changes were observed along the altitudinal gradients (Table S1 and Figure  S5). Distance decay of ECM fungal communities following altitudinal gradients has been observed in most studies, irrespective of host identity control [30,[68][69][70]. Interestingly, the effect of altitudinal gradients was weaker in the ECM than in the non-ECM groups of both microhabitats in our study. As our study focused on a single host species, the ECM community may have been less affected at the regional scale. Bahram et al. [28] have reported that host tree species was the most important factors that affected the ECM community composition. In contrast, free-living fungi, such as saprotrophs, in soil are expected to be more affected by geographic distance than host-dependent fungi because the former are relatively more functionally overlapped [71] and sensitive to species pools in the area [69,72,73].
Similarly, belowground endophytes in roots have been reported to be relatively host-independent [74] and dispersal-limited, which may account for the differences in altitudinal and geographic patterns. However, other environmental factors such as annual temperature or soil chemical properties around individual trees could have also resulted in these differences [28,35].
Species richness and diversity decreased along altitudinal gradients in the root fungal communities of P. densi ora, whereas altitudinal patterns of alpha diversity differed between ecological groups in the soil fungal communities. In previous studies, a decline in ECM species richness has been observed along altitudinal gradients in both microhabitats [28,75]. However, studies with controlled host species have reported no relationship between ECM diversity and the altitude in the roots of Fagus sylvatica [29] and in the soil of P. sylvestris [35,36]. This could be due to the characteristics of the host tree (P. densi ora) in our study; however, further studies are needed to clarify the effect of altitude on fungal richness associated with P. densi ora. Although the alpha diversity of the non-ECM group decreased at high altitudes in the root microhabitat, there was no signi cant change in soil. Decreased growth and photosynthetic productivity of P. densi ora has been reported at high altitudes [76]. The lowered productivity of host trees might have affected the richness and diversity of associated fungi, i.e., ECM and root-associated non-ECM groups [77]. The non-signi cant relationship between altitude and alpha diversity of the non-ECM community in soil may suggest that the decomposition of organic matter and/or productivity of understory vegetation is less affected by altitude [78].

Conclusion
Overall, our study provides a basis for understanding belowground fungal community compositions associated with P. densi ora in temperate forests of South Korea along altitudinal gradients. Our results indicate strong correlations between richness, diversity, and altitude in plant-associated fungal communities, i.e., non-ECM fungi in roots and ECM fungi, but not in non-ECM fungal communities in soil.
Microhabitat types and guilds of fungal taxa are important factors in determining fungal diversity patterns along altitudinal gradients, which highlights the need to study various ranges of microhabitat types to understand the factors that in uence root-inhabiting fungi and ECM fungi. By extending the targets to diverse types of plant hosts and microhabitats, these ndings can provide a solid foundation for plant microbiome management.