Influence of the slope aspect on the ectomycorrhizal fungal community of Quercus variabilis Blume in the middle part of the Taihang Mountains, North China

Songpo Wei · Yijing Song · Liming Jia

Abstract The slope aspect is one of the most critical topographic factors in mountainous areas. Little is known, however, about the effect of the aspect on the ectomycorrhizal(ECM) fungal community. Additionally, we know very little about the composition of ECM fungal communities associated with Quercus variabilis, which is widely distributed in East Asia. In this study, we compared the richness, community composition, and exploration types of ECM fungi associated with Q. variabilis between predominantly south- and north-facing slopes in the Taihang Mountain, North China for the first time. DNA was extracted from the root tips of Q. variabilis, and Illumina MiSeq sequencing was used to identify ECM fungi. In total, 168 operational taxonomic units belonging to 28 genera were detected, and the ECM community was found to be dominated by Russula, Inocybe, Tomentella, Scleroderma, and Cortinarius. Compared with the north-facing slopes, the ECM communities on the south-facing slopes had higher diversity. The community composition and exploration types were directly affected by the slope aspect. Also, the aspect-induced edaphic variables, such as total phosphorus,total nitrogen, total potassium, pH, and soil water content,were important sources of variation in ECM fungal richness and distributions of exploration types. Different genera tended to be distributed in various slope aspects. Cenococcum, Genea, and Clavulina were significantly enriched in north-facing slopes, while Geopora, Helvelosebacina, Scleroderma, Gyroporus, Astraeus, Boletus, Tricholoma, Hebeloma, Cortinarius and unclassified Thelephoraceae were more abundant in south-facing slopes. Hydrophobic ECM fungi were obviously enriched in the south-facing slope, but there was no statistical difference between hydrophilic among the south- and north-facing slopes. Our study deepened our knowledge of the aspect-driven variation in ECM fungal communities associated with Q. variabilis.

Keywords Ectomycorrhizal community · Illumina sequencing · Quercus variabilis · Slope aspect · Edaphic factors

Introduction

Forests are complex ecosystems that provide a multitude of resources to people (Lindberg et al. 1997). And ectomycorrhizal (ECM) symbiosis is omnipresent in global forests (Franklin et al. 2014). Additionally, an estimated 20,000-25,000 ECM fungal species exist in the forest ranges around the globe(Brundrett 2009), promoting plant growth by improving water and nutrient absorption, providing tolerance to environmental stresses, and playing an essential role in forest nutrient cycling(Smith and Read 2008; Kumar and Atri 2018). Because the performance benef its of ECM for plants vary among fungal species or strains (Hoeksema 2010; Johnson et al. 2012), it is necessary to understand ECM fungal community composition.

ECM fungal communities are affected by various biotic and abiotic factors (Essene et al. 2017). Due to the specificity of many ECM fungi to tree genera (Bruns et al. 2002), host plants almost determined the types of mycorrhizae (e.g., Ectomycorrhiza, Orchid Mycorrhiza) (Kennedy et al. 2018). And many studies have shown that the host plant distribution always affects the ECM richness and community composition significantly (Burke et al. 2009; Bahnmann et al. 2018; Saitta et al.2018; Veach et al. 2018). Moreover, apart from host species identity, the ECM fungal community structure and richness could be influenced by the host genotypes (Wang et al. 2018),host tree age (Zhang et al. 2014), host tree diameter at breast height (Burke et al. 2009), and even health status of plants(Ishaq et al. 2018). Some studies, especially those conducted under large-scale conditions, have suggested that host plants may have the strongest effect on the community composition of ECM fungi (Tedersoo et al. 2012; Polme et al. 2013), even other factors such as elevation, climatic, and edaphic parameters were almost negligible compared with the role of the host types(Saitta et al. 2018). However, at the same time, many studies on ECM fungal communities have suggested that the effects of abiotic environmental factors cannot be ignored (Essene et al. 2017; Glassman et al. 2017; Suz et al. 2017). Actually, a range of edaphic and terrain factors, including the pH (Tedersoo et al. 2014; Veach et al. 2018), moisture content (Jarvis et al.2015; Nickel et al. 2018; Reis et al. 2018), phosphorus (Twieg et al. 2009; Burke 2015; Fan et al. 2018; Zavisic et al. 2018),nitrogen (Kjoller et al. 2012; Van Strien et al. 2018), potassium(Bauman et al. 2016), soil carbon (Jarvis et al. 2015; Addison et al. 2018; Veach et al. 2018), C:N ratio (Coince et al. 2013;Maghnia et al. 2017; Bahnmann et al. 2018), elevation (Jarvis et al. 2015; Matsuoka et al. 2016; Bowman and Arnold 2018),slope steepness (Matsuoka et al. 2016), and slope aspect (Scattolin et al. 2008), have been shown to be directly or indirectly associated with ECM fungal communities. And, there is also an intergeneric variability in the fungal responses to the above factors (Tedersoo et al. 2003; Morgado et al. 2015).

For mountain areas, the slope aspect is undoubtedly an important environmental factor equal to the elevation and slope steepness. The nature of the differences in the slope aspect is the differences in the net solar radiation received(McCune and Keon 2002). Additionally, existing research has shown that the slope aspect characteristics directly determine the hydrothermal (Tedersoo et al. 2012; Jarvis et al. 2013,2015; Morgado et al. 2015) and light conditions (Kummel and Lostroh 2011), forming a different microenvironment(Sternberg and Shoshany 2001; Fang et al. 2004) that, in turn,affects the biotic community and edaphic properties (Chen and Huang 1997), as well as the tree species composition(Sternberg and Shoshany 2001). Some studies have demonstrated that the slope aspect has a significant impact on the community composition and OTU richness of arbuscular mycorrhizal (AM) fungi (Chu et al. 2016; Liu et al. 2017; Chai et al. 2018) and total soil microbe (Frey et al. 2016; Ai et al.2018; Xue et al. 2018). It could be inferred that the difference in the slope aspect may shape different ECM community structures. However, unfortunately, compared with the elevation and slope steepness, the relationship between the aspect and ECM community is not well studied (Scattolin et al. 2008;Geml 2019). At the same time, some studies considered the differences between microbial communities from different aspects to be relatively minor (Zhang et al. 2010; Schlatter et al. 2018; Xu et al. 2018). Those contradictory results indicated that the effect of the aspect on microbial communities is more complicated than expected, and the relationship between the aspect and ECM fungi needs further research.

The Taihang Mountains, which form the boundary between the Shanxi plateau and North China Plain, was marked by dry,stony, and poorly fertile soils. To restore vegetation damaged by long-term human activities, the Taihang Mountains were selected as a major afforestation region in China presently(Delang and Yuan 2015). Considering that mycorrhizal seedlings have been widely used in afforestation worldwide, and its effectiveness was conf irmed many years ago (Querejeta et al. 1998), mycorrhization should be able to play an essential role in vegetation restoration. However, owing to the lack of research, little information is available about the ECM fungus associated with plants in the Taihang Mountains (Wang et al.2012; Wei et al. 2018). In this area, the slope aspect variance always leads to significant differences in site conditions, such as, most notably, north-facing slopes typically retaining more water, accumulating more organic matter and being more productive than south-facing slopes (Yang et al. 1993; Yang et al.2006). Therefore, it is likely that the slope aspect substantially impacts ECM fungal communities.

Oaks, the most diverse genus of ECM host plants, are an essential component of temperate forests in the Northern Hemisphere (Garcia-Guzman et al. 2017). However, their ECM communities have been scarcely studied compared with those of conifers (Smith et al. 2007). Approximately 21 species belonging to oaks have been studied worldwide (Garcia-Guzman et al.2017), but the ECM fungal community of Q. variabilis, which is widely distributed in East Asia, has not been well investigated(Garcia-Guzman et al. 2017; Wei et al. 2018). As a widespread oak species, Q. variabilis is also a major constituent of deciduous broad-leaved warm temperate forests in northern China(Chen and Huang 1997), and it is widely used for soil and water conservation, especially for the afforestation of difficult site conditions. We have had a preliminary understanding of the morphological features of ECM root tips of Q. variabilis (Wei et al.2018). Additionally, it is necessary to conduct in-depth research on the ECM fungal community associated with Q. variabilis.

Here, we collected ECM root tips and soil samples from eight different slopes during the growing season and explored the richness and composition of the ECM fungal community of Q. variabilis Bl. in different slope aspects based on the rDNA sequencing of samples from ECM root tips using the Illumina MiSeq platform. We investigated three issues: (1) the distribution of ECM fungi at the different slopes; (2) whether the slope aspect has a significant direct effect on ECM fungal diversity and community composition; and (3) whether the ECM fungal richness and community composition were related to edaphic factors that were associated with the slope aspects.

Materials and methods

Site description

The study site was located in the central part of Taihang Mountains (37° 7′ N, 113° 58′ E), near Xingtai, in Hebei Province (Fig. 1), China. The sites have a warm, temperate and semi-humid climate characterized by the synchronization of high temperatures and ample precipitation, with terrain comprising hills that are 591-652 m above sea level.The average annual precipitation and air temperature from 1996 to 2016 were 558.7 mm and 12.8 °C, respectively. The rainfall distribution is uneven from season to season, with over 80% of the rainfall concentrated in the summer months(July and August). The growing season for typical deciduous species in this region is from April to October. As is typical of stony hillside stands, the soils are derived from gneiss with high gravel contents. The site is Q. variabilis-pure secondary forests (approximately 30-40 years old).

Roots and soil sampling

Eight 20 × 20 m plots (Table 1) were established in pure Q. variabilis stands with different slope aspects to investigate the ECM community structure, and the plots were separated at least 100meters apart, the aspect (ASP) and slope (SLO) were measured using a compass. The sampling dates were staggered to ensure that all of the samples could be treated in a timely fashion.

Table 1 General characteristics of sampling points

All the root system samples were collected from June to August 2016. Fifteen trees were selected randomly in each plot. For each tree, one approximately 20 × 20 × 20 cm soil block containing the roots of Q. variabilis was collected. To avoid collecting the roots from other plants, all roots were traced from the trunk of the selected trees. Each sampled root system was placed into a plastic bag and was taken to the laboratory immediately and stored at 4 °C until analysis. A total of 120 root samples was collected. The roots were washed carefully with tap water, and then all the living ectomycorrhizal root tips were separated from the dead roots by examination at 0.7-11.5 × magnification (OLYMPUS SZX16, Japan)with attention to their visual appearance and elasticity (Aerts et al. 1989). All the root systems samples were processed within ten days. Next, every five random samples from the same plots were mixed as one root sample and were storedat − 20 °C before DNA extraction. Finally, 24 root system samples from all eight plots were analyzed in this study.

The soil sample collection was accompanied by root system samples collection. For each tree, more than 200 g of soil was collected from around the area of the root system samples and was dried in the open air. In total, 120 soil samples were collected. At each plot, every five soil samples were mixed equally together as one soil sample, resulting in 24 mixed soil samples for follow-up chemical analysis.

Soil chemical analysis

The soil pH was measured using a pH meter (INESA INSTRUMENT, China) after mixing 10 g of the soil sample with 25 mL of 1 M potassium chloride (KCl). The soil organic matter (OM)content was determined using the Walkley-Black method of quantifying the amount of oxidizable soil carbon (C) content via reaction with acidic dichromate (Cr2O72− ), with the OM(g kg −1) equaling the C (g kg -1 × 1.724) (Nelson and Sommers 1996). To determine the total soil nitrogen (TN) content, the soil samples were melted by sulfuric acid (H2SO 4) and hydrogen peroxide (H2O 2 ), and then the TN was determined by a continuous flow analytical system (AA3; SEAL Analytical,UK). To determine the total soil potassium (TK) and total soil phosphorus (TP) contents, the soil samples were melted using sodium hydroxide (NaOH) at 720 °C, and then the TK and TP were determined by atomic emission photometry (SpectrAA 220 Atomic Absorption Spectrometer; VARIAN, USA)and the molybdate method using a UV spectrometer (Agilent 8453, USA), respectively. The available soil phosphorus (AP)was extracted with 0.5 M sodium bicarbonate (NaHCO 3 ), and the P concentration was determined colorimetrically using the molybdate method with a UV spectrometer (Agilent 8453,USA) (Olsen et al. 1954). The available potassium (AK)was extracted with 1 M ammonium acetate (CH3COONH 4 ;pH = 7.0), and the K concentration was quantified by atomic emission photometry (SpectrAA 220 Atomic Absorption Spectrometer; VARIAN, USA). The mass water content of the soil(SWC) was determined using the oven-drying method.

DNA extraction and polymerase chain reaction (PCR)

For each plot, all the living ectomycorrhizal root samples from the same sampling date were thoroughly mixed; the 24 mixed root samples were then ground in liquid nitrogen. The total genomic DNA in the roots was extracted using a modified CTAB method (Doyle and Doyle 1990), and the extracted DNA was purified using an E.Z.N.A soil DNA kit (Omega Biotek, USA) according to the manufacturer’s protocol (Murray and Thompson 1980; Mundra et al. 2015). The ITS1 region of the fungal rDNA gene was amplified by PCR using primers ITS1F (CTT GGT CAT TTA GAG GAA GTAA) and ITS2 (GCT GCG TTC TTC ATC GAT GC) (White et al. 1990; Gardes and Bruns 1993; De Beeck et al. 2014). The PCR amplification was conducted in a total volume of 20 μL containing 4 μL of 5 × FastPfu Buffer (TransGen Biotech, China), 2 μL of dNTPs(2.5 mM), 1.6 μL each of forward/reverse primer (5 μM),0.4 μL of FastPfu Polymerase (TransGen Biotech, China),0.2 μL of BSA, and 10 ng of template DNA, with sufficient ddH2O to bring the final volume to 20 μL. The DNA amplification was performed in an ABI GeneAmp 9700 thermocycler(PE Applied Biosystems, USA) using the following thermalcycling parameters: an initial denaturation at a temperature of 95 °C for 3 min, followed by 32 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s, and then a final extension at 72 C for 10 min. The amplicons were gel purified using the Axy-Prep DNA gel extraction kit (Axygen Biosciences, USA) and were quantified using QuantiFluor TM -ST (TBS-380, Promega,USA). After the purified amplicons were pooled in equimolar amounts, they were PE sequenced (2 × 250) on an Illumina MiSeq platform according to standard protocols obtained from Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China).The raw sequences obtained from this study were deposited in the NCBI Sequence Read Archive (SRP149624).

Sequence processing

The sequencing data were processed using QIIME (Caporaso et al. 2010) with the following criteria: Only sequences > 200 bp in length with an average quality score > 20 and without ambiguous calls were included in the subsequent analysis. Reads which could not be assembled were discarded. The sequences were clustered into OTUs at a 97% similarity threshold using UPARSE (Edgar 2013), and chimeric sequences were identified and removed using UCHIME. They were then assigned taxonomically using the UNITE database 7.0 (Nilsson et al. 2018).

We classified each OTU into an ecological guild by FUNGuild(Nguyen et al. 2016), all the OTUs that were categorized as ECM fungi with a conf idence level designated as “highly probable” or“probable” were included in this analysis (Nguyen et al. 2016).After filtering for ECM fungal taxa, the OTU table was subsampled to obtain the minimum sequences to normalize for different sequencing depths using micca 1.2 (Albanese et al. 2015). All ECM taxa that could be identified to genus level were also classified into different exploration types and hydrophilic or hydrophobic properties according to previous studies (Agerer 2001, 2006).

The alpha diversity of the samples (richness, Shannonwiener index, and Shannon evenness) was calculated using EstimateS 9.1.0 (Colwell 2013). To simplify the estimation of the relationship between the aspect and ECM fungi, the azimuth of the aspect was transformed to northness using the cosine function (Roberts 1986) as follows:

Data analyses

Except for the part described above, all the other statistical analyses were performed in R 3.4.4 (R Core Team 2018).The relative abundance of each ECM genera (families) was calculated as the number of genera (families) sequences divided by the total number of ECM sequences. According to the data of the slope aspect, all eight plots were divided into two groups, the north- and south-facing slope, by Ward’s method cluster analysis using the “hclust” command.The correlations between ECM fungal diversity indices and edaphic factors were evaluated using Spearman’s rank correlation by “rcorr” command in the “Hmisc” package (Frank and Harrell 2018; Schlatter et al. 2018).

The differences in the total OTU richness, richness and abundance of hydrophobic/hydrophilic ECM fungi between slope aspects were evaluated using generalized linear models with the “glm” command, and they were modeled following Poisson distribution. The relationship between the diversity indices (Shannon-Wiener index and Shannon evenness) and the aspect was analyzed by linear models. The prediction of the data on the north-facing (assuming a slope aspect of 0°)and south-facing (assuming a slope aspect of 180°) slopes utilized the “predict” command.

We used variation partitioning based on redundancy analysis (RDA, ‘rda’ command in the “vegan” package) to quantify the contribution of the aspects, slopes and edaphic variables to ECM richness and exploration types (Oksanen et al. 2017). The proportion of variation in ECM richness and exploration types explained by slope aspect and edaphic variables calculated independently by “permu.hp” command in the “rdacca.hp” package (Lai 2019). The Mantel tests were conducted in R using the “mantel” command from the“vegan” package (Oksanen et al. 2017) to determine the relationship between the Euclidean distance of the slope aspect and Hellinger distance of the ECM fungal composition.

The Bray-Curtis metric of the ECM fungal communities from different sites was tested using one-way analysis of similarities (ANOSIM) by calculating 999 random permutations using PAST 3 (Hammer et al. 2001). Using the same metrics,we visualized the differences in the ECM community structure using nonmetric multidimensional scaling (NMDS) ordination as implemented by the “vegan” and “ggplot2” packages(Wickham 2016; Oksanen et al. 2017). LefSe analysis was run using the online version at https://hutte nhowe r.sph.harva rd.edu/galax y/, with the default settings as follows: ɑ = 0.05 for the factorial Kruskal-Wallis test among different aspects(no subclasses), the threshold on the logarithmic LDA score for discriminative features was set at 2.0 (Segata et al. 2011).

In this study, we treated the aspect as a categorical factor (north-facing and south-facing) for NMDS analysis,ANOSIM analysis, and LefSe analysis while the aspect was treated as a continuous variable (transformed azimuth) for linear models, generalized linear models, and RDA analyses.

Results

Taxonomic representation

In total, 492 OTUs with 921,009 gene sequences were obtained. The number of OTUs accurately ascribed to taxonomic ranks decreased at each rank level, with 412 OTUs(83.70%), 386 OTUs (78.5%), 334 OTUs (67.9%), and 280 OTUs (56.9%) identified at the class, order, family, and genus levels, respectively. 21 OTUs were animal pathogens,51 OTUs were saprotrophs, 11 OTUs were endophytes, one OTU belonged to arbuscular mycorrhizal fungi, 6 OTUs were fungal parasites, 11 OTUs were plant pathogens, 168 OTUs belonged to ectomycorrhizal fungi, and 223 OTUs couldn’t be sure of functional classification.

OTUs that belonged to ECM fungi were selected, the samples were then subsampled to 28,614 sequences per sample to normalize for different sequencing depths. Finally, 168 OTUs with 686,736 sequences from 24 mixed root samples were included in further analysis. Among the 168 OTUs,121 OTUs belonged to 28 genera, and 47 OTUs, which were identified at the family level belonged to 6 families.According to the relative abundance, the 20 and 100 most abundant ECM OTUs accounted for 84.5% and 96.4% of the total ECM sequences, respectively. Basidiomycetes was the dominant phylum, and it accounted for 98.0% of the total sequences, and only 2.0% represented ascomycetes. For richness, 149 of the 168 OTUs (88.7%) belonged to basidiomycetes, and 19 OTUs (11.3%) belonged to ascomycetes.

At the genus level, Inocybe (32 OTUs), Tomentella (31 OTUs), Russula (7 OTUs), Cortinarius (6 OTUs), Helvella(6 OTUs) were the five most OTU rich genera. And Russula(364,937 sequences), Inocybe (76,157 sequences), Tomentella(49,550 sequences), Scleroderma (30,326 sequences), Cortinarius (14,177 sequences) were the five most abundant genera. ECM fungi in Thelephoraceae that were not identified at the genus level also had high OTU richness (38) and abundance (89,293 sequences). The composition of the ECM fungal communities from different plots was prof iled according to the relative abundance based on the sequence numbers (Supplementary Table S1).

ECM fungal diversity and edaphic factors across different aspects

Through cluster analysis by Ward’s method, eight plots were divided into two groups: the south-facing slope (180°-230°)and north-facing slope (290°-345°) (Table 1).The soil characteristics showed variability between the various aspects(The detailed data was shown in Supplementary Table S2).According to the linear models, TN, OM, and SWC correlated positively with ASP; TP, TK, pH correlated negatively with ASP. Thus, TN, OM, and SWC were significantly higher in north-facing plots, and TP, TK, pH were significantly higher in south-facing plots. However, the slope aspects had no significant effects on AP, AK, and C/N ratio (Table 2).

Based on the generalized linear models (for OTU richness) and linear models (for Shannon and Shannon evenness), the ASP had significant effects on OTU richness(P< 0.001), the Shannon-Wiener index (P = 0.002), and Shannon evenness (P = 0.011). The predicted results showed that the diversity indices of ECM communities on the south-facing slope were 1.6 times (OTU richness), 1.6times (Shannon-wiener) and 1.4 times (Shannon evenness)higher than those on the north-facing slope, respectively(Fig. 2). ECM fungal diversity indices were positively related to TP (R = 0.44, 0.45, 0.44;P< 0.05, for richness,Shannon, and Shannon evenness diversity, respectively),pH (R = 0.41, 0.50, 0.47,P< 0.05, for richness, Shannon and Shannon evenness diversity, respectively). Shannon and Shannon evenness diversity was negatively related to TN(R = − 0.45, − 0.45,P< 0.05, respectively). Moreover, there was a positive correlation between OTU richness and TK(R = 0.44,P< 0.05).

Table 2 Linear model results presenting the effects of slope aspects(cosine transformed) on edaphic factors

Variation in ectomycorrhizal fungal composition among aspects

At genus level, 16 genera were detected in both the northfacing and south-facing slopes (Fig. 3 a). Five genera presented in the south-facing slopes that were not detected in the north-facing slopes were Tricholoma, Hebeloma,Sphaerosporella, Butyriboletus, and Geopora. And Cenococcum, Timgrovea, Hysterangium, Clavulina, Otidea, Chloridium and Strobilomyces were only found on the northfacing slopes. At the OTU level, 50 OTUs were recovered only from south-facing slopes, 53 OTUs occurred only on north-facing slopes, and 65 OTUs were detected on both slopes (Fig. 3 b).

The NMDS analysis of all 24 mixed samples, which was based on the Bray-Curtis distance, demonstrated a clear difference in the ECM fungal composition from samples taken from different aspects (Fig. 4). And the ANOSIM analysis based on the Bray-Curtis similarity index also revealed that the aspects had a significant effect on the ECM community at the OTU level (R = 0.121, P = 0.001). Also, the Mantel test revealed that the Hellinger distances of the ECM fungal composition were significantly correlated with the ASP(R = 0.129, P = 0.009)

Fig. 2 Predicted (mean ± SE) ECM fungal OTU richness, Shannon-Wiener index, and Shannon evenness across north-facing (0°) and southfacing (180°) slopes

Fig. 3 Venn diagrams of detected OTUs (a) and genera(b) on the two different slope aspects

ASP was a significant contributor to variance in ECM fungal richness (18.2%, P = 0.001) when the slope aspect was calculated independently with RDA (Fig. 5 a). As showed by the biplot, Cenococcum had a higher richness on the north-facing slopes, while Gyroporus, Cortinarius had a higher richness on the south-facing slopes. When considering only edaphic factors, nine edaphic factors (TP, AP,TK, AK, OM, TN, C/N, pH, and SWC) explained 21.9%(P = 0.004) of the variance in ECM fungal richness. Among them, TP, TK, TN, and pH were significant source of variation in ECM fungal richness (Fig. 5 b, Table 3), but the remaining variables, AP, AK, OM, C/N, and SWC played less important roles (Table 3).

LEfSe analysis was used to test for significant taxonomic differences among the north- and south-facing slopes(Fig. 6). And the results indicated that the distribution of many groups at different taxonomic levels could be significantly affected by aspect. Ascomycota and Basidiomycota appeared as the biomarker of north-facing and south-facing slopes, respectively. At the genus level, Cenococcum,Genea, and Clavulina were significantly enriched in northfacing slopes, while Geopora, Helvelosebacina, Scleroderma, Gyroporus, Astraeus, Boletus, Tricholoma, Hebeloma, unclassified Thelephoraceae and Cortinarius were more abundant in south-facing slopes.

Exploration types of ECM fungi across different aspects

When the factor was calculated independently with RDA(Fig. 7 a), ASP was a significant contributor to variance in distributions of exploration types (19.6%, P = 0.002). ECM of medium-distance fringe (MDF) exploration type and long-distance (LD) exploration type showed tendencies to associate with south-facing slopes. However, slope aspects had less influence on the distributions of short-distance (SD)exploration type, contact (C) exploration type, and mediumdistance mat (MDM) exploration type. When considering only edaphic factors, all edaphic factors explained 33.9%(P = 0.010) of the variance in distributions of exploration types. TP, TN, and SWC were significant source of variation in exploration types (Fig. 7 b, Table 4), but the remaining variables, such as AP, TK, AK, OM, C/N, and pH, played less important roles (Table 4). As shown in Fig. 7 b, MDF and LD exploration types tend to be distributed in areas with a high level of TP, and SD exploration type tends to be distributed in areas with a high level of TN and SWC. And C exploration type tends to be distributed in areas with a high level of SWC.

Fig. 4 NMDS analysis of the ECM fungal communities associated with Quercus variabilis

Fig. 5 RDA ordination biplot indicating the influence of aspect (a)and edaphic (b) variables on the richness of genera. Only the 15 most OTU-rich genera (AMA Amanita, BOL Boletus, CEN Cenococcum, CLA Clavulina, COR Cortinarius, GYR Gyroporus, HEL Helvella, INO Inocybe, LAC Lactarius, RUS Russula, TOM Tomen-tella, TUB Tuber, BOC Boletaceae unclassified, SEB Sebacinaceae unclassified, THE Thelephoraceae) and significant variables (ASP slope aspect, TP total phosphorus, TN total nitrogen, TK total potassium, pH) were shown in the figure

Table 3 Proportion of variation in ECM fungal richness explained by edaphic factors

Generalize linear model analysis indicated that both OTU richness and abundance of hydrophobic ECM fungi were negatively correlated with cosine transformed aspect(P< 0.001 and P = 0.009, for OTU richness and abundance,respectively). It meant that hydrophobic ECM fungi had higher richness and abundance on the south-facing slopes(Table 5). However, the OTU richness and the abundance of hydrophilic ECM fungi were not significantly associated with the aspect (P = 0.156 and P = 0.083, for OTU richness and abundance, respectively).

Discussion

Fig. 6 The cladogram based on LEfSe analysis (P < 0.05, LDA score > 2) showing the significantly different abundant of ECM fungi among the north- and south-facing slopes. The color nodes represent ECM taxa that are significantly enriched in the corresponding slope aspect. Yellow nodes denote that the taxa are not significantly higher in any aspects. From the center to the periphery, the taxa represent the kingdom, phylum, class, order, family, and genus levels. Names of biomarkers are given for the phyla and class levels in the figure, and the names of other biomarkers are represented by letters

In this paper, we investigated the ECM fungal communities of Q. variabilis, which has not been well studied previously. Among the ECM fungi that were associated with Q. variabilis in this region, 13 genera such as Tomentella,Cenococcum , Cortinarius, Inocybe, Russula, Lactarius,Peziza, Tuber, Amanita, Boletus, Clavulina, Scleroderma,and Tricholoma were frequently detected (more than half of all 17 sites worldwide, Supplementary Table S3) in ECM fungal communities of other oaks grown in Asia (Wang et al.2012; Toju et al. 2013; Zhang et al. 2013; Huang et al. 2014;Zhang et al. 2014; Geng et al. 2016; He et al. 2016), Europe(Courty et al. 2008; Leski et al. 2010; Trocha et al. 2012;Lancellotti and Franceschini 2013; Scattolin et al. 2014; Suz et al. 2014) and North America (Smith et al. 2007; Morris et al. 2008, 2009; Dickie et al. 2009; Moser et al. 2009).In our study, the above 13 genera accounting for 82.5% of the total sequences and 57.7% of the total OTU richness.These findings indicated that the ECM fungi associated with Q. variabilis were similar to other oaks in Northern Hemisphere at the genus level. However, some genera, such as Chloridium, Genea, Geopora, Helvella, Otidea, Sphaerosporella, Astraeus, Gyroporus, Hebeloma, Hysterangium,Strobilomyces, Timgrovea, and Tylopilus were less common in ECM fungal communities of other oaks worldwide. Interestingly, as a biomarker of the south-facing slopes in this region, Gyroporus had never been reported in ECM communities of other oaks worldwide. Based on this result, we assumed that Gyroporus might be an endemic genus that was associated with the Q. variabilis in this region. Although the high-throughput sequencing method used in this study is powerful, its imitations (amplifying only one of the two ITS regions) make it difficult to identify most of the found taxa down to the species level. More endemic ECM fungal species may be detected due to changes in sequencing methods. In addition, Five genera (Russula, Inocybe, Tomentella,Scleroderma, and Cortinarius) with the highest abundance among all 34 taxa in this region accounted for over 77.9%of the total sequences and 45.8% of the total OTU richness.Thus, the composition of the ECM fungal communities in this region was also in line with the theories of previous researchers that the composition is very OTU rich with a few species being abundant and many species being rare(O’Hanlon 2012).

Fig. 7 RDA ordination biplot indicating the influence of aspect (a)and edaphic (b) variables on the exploration types. C contact exploration, SD short-distance exploration, MDM medium-distance mat exploration, MDF medium-distance fringe exploration, LD long-dis-tance exploration. Only significant variables (TP total phosphorus;TN total nitrogen; SWC soil water content, ASP slope aspect) were shown in the figure

Table 4 Proportion of variation in distributions of exploration types explained by edaphic factors

Table 5 Generalize linear model results presenting the effects of slope aspects (cosine transformed) on OTU richness and abundance of hydrophobic and hydrophilic ECM fungi

Undoubtedly, the slope aspect had a significant impact on the ECM fungal community composition and richness in this study. The predicted results based on generalized linear models showed that the south-facing slopes have a higher OTU richness than the north-facing slopes. This finding was in line with studies focusing on AM fungi in the arid ecosystems of Inner Mongolia and the Tibetan Plateau in China (Chu et al. 2016; Liu et al. 2017). However, the situation on the oak stands growing in Hungary is quite contrary to our results. In that study, although the OTU richness of saprotrophs and plant pathogens changed significantly with the slope aspect, the OTU richness of the ECM fungi was not affected by the slope aspect (Geml 2019). It seemed that changes in host plants in some studies may also lead to uncertainty in results (Liu et al. 2015; Chu et al. 2016; Chai et al. 2018; Geml 2019). Considering the strong influence of the aspect on the distribution of host plants (Sternberg and Shoshany 2001) and the impact of host plants on the ECM fungal communities (Tedersoo et al. 2012), investigating the effect of the slope aspect on ECM communities under conditions of unf ixed host plants may make it difficult to obtain “pure” effects of the slope aspect. In our study, a single host plant and the relatively even-aged forest were chosen to minimize the effect of host plants on the ECM community and highlight the impact of the slope aspect, and it also allowed us to study the effects of the aspect on ECM fungal communities with less interference from host plants changes (Jarvis et al. 2015).

Some previous studies have focused on AM or ECM fungi and suggested that the slope aspect had weak direct effects on AM fungus community composition or OTU richness but was affected by aspect-induced changes in edaphic and plant communities (Scattolin et al. 2008; Chu et al. 2016;Liu et al. 2017). However, in this study, as a continuous variable, the cosine-transformed slope aspect (value from− 1 to 1) clearly showed a significantly negative correlation with the diversity indices (i.e., the ECM communities on the south-facing slope had a higher OTU richness, Shannon-Wiener index, and Shannon evenness). And slope aspect was a significant contributor of variance in ECM fungal richness and distributions of exploration types (Figs. 5 a, 6 a). On the basis of our findings above, it can be concluded that aspect has a significant direct effect on ECM fungal communities.Our results were supported by the study of Geml (2019),which showed that there was a significant linear relationship between the differences in aspect (degrees) and distance of fungal community composition (Bray-Curtis dissimilarity). Thus, the slope aspect was more similar to continuous variables such as altitude and was different from the spatial variation caused by site effects.

RDA revealed that edaphic factors, such as TP, TN,TK, pH, and SWC, were essential sources of variation in ECM fungal richness and distributions of exploration types(Tables 3 and 4). And the five edaphic factors above were significantly affected by the slope aspect (Table 2). These findings suggest that some of the effects of slope aspect on ECM fungal richness and distributions of exploration types may be explained by aspect-induced differences in edaphic factors.

The diversity of ECM fungi decreased with the increase of TN. This finding was also reported by Lilleskov et al.(2002) and De Witte et al. (2017). However, some studies found that soil nitrogen was positively correlated with ECM richness (Schlatter et al. 2018), or had no significant effect on ECM fungal community (Coince et al. 2013; Addison et al. 2018). According to the results of RDA in this study,ECM exploration types showed particular tendencies to soil TN levels. These results may explain the different responses of ECM fungi to the N level.

TP played an essential role in shaping ECM fungal communities and exploration types. These results reflect those of Burke (2015) who also found that ECM community structure significantly correlated with edaphic parameters,especially P. Figure 7 b plots the enrichment of LD and MDF exploration types in the higher TP conditions. Many studies had conf irmed that ECM fungi could enhance P mobilization via many mechanisms (Teste et al. 2016), and the stronger ability of LD exploration type to forage P than SD and C types (Lilleskov et al. 2019) might be the reason for its higher abundance in high TP conditions.

Consistent with the previous studies, SWC and pH were critical edaphic factors that determine ECM communities and exploration types (Tedersoo et al. 2014; Nickel et al.2018; Van Geel et al. 2018; Veach et al. 2018). The above two factors were both related to the slope aspect. For the SWC, under a small scale, there will be no significant difference in precipitation between plots. Therefore, the slope aspect maybe the most important factor that influences SWC(Mendez-Toribio et al. 2016) in present study.

According to the LEfSe analysis (Fig. 6), different ECM fungal genera showed particular tendencies to associate with a specific aspect, and this finding was in line with previous studies (Reverchon et al. 2012; Morgado et al. 2015). LEfSe analysis also allowed us to obtain the extent to which the principal ECM taxa were affected by the slope aspect in this region. Functional differences may explain the inter-genera variability in the richness and abundance of the ECM fungi.Unlike Geml (2019) who did not find a significant difference between exploration types among south- and north-facing slopes, the analysis of the ECM exploration types (Fig. 7 a)and hydrophobicity of the identified ECM fungi (Table 5)indicated that aspect significantly affected the exploration types of ECM in this study.

The LD and MDF exploration types, which were often considered to be hydrophobic ECM (Agerer 2006), were more likely to be distributed on the south-facing slopes in our study (Fig. 7 a, Table 5). This group consists most notably of genera Cortinarius, Astraeus, Boletus and Scleroderma. Among these genera, Scleroderma had been reported to have a high ability to persist during droughts (Mrak et al.2017). Thus, it might be speculated that some genera with higher abundance on the south-facing slope than on the north-facing slope, especially the taxa with the biomarker of the south-facing slope, are likely to have the similar droughtresistant ability.

Contrary to hydrophobic ECM, hydrophilic ECM was almost unaffected by aspect. Most of the SD and C exploration types were considered to belong to hydrophobic ECM(Agerer 2006), and it represented by genera such as Russula and Cenococcum . In this study, Russula had the highest relative abundance among all genera, and the distribution of Russula was independent of the slope aspect or edaphic factors, which showed the close relationship between Russula and Q. variabilis, and this relationship was also reflected in other ECM communities associated with oaks (Reis et al.2018).

One unanticipated finding was that Cenococcum, which is widely distributed throughout the world and is considered to be well-suited to various site conditions (Dickie 2007),only occurred on north-facing slopes in the present study.The work of Kummel and Lostroh (2011) revealed that host plants living under lower light conditions have a higher abundance of Cenococcum. Based on this observation,the lower light levels in the north-facing slopes might, to a certain extent, explain the distribution of Cenococcum.Unfortunately, the sequence abundance in the sequencing method we used might be an inaccurate proxy for the actual abundance within a mixed community (Palmer et al. 2018).To accurately analyze the relationship between an individual taxon, quantitative fluorescence PCR (QF-PCR) with genusspecific primers may be used to conduct quantitative analyses of a particular ECM genus in future research.

Slope aspect is an independent factor, and also represents a set of environmental factors. But the correlation between aspect and other environmental factors may pose a certain challenge for analysis. To unbundle the edaphic factors from the slope aspect, the large sample size is needed in future research. We concentrated only on partial abiotic factors in this study. According to the RDA, slope aspect and edaphic factors could explain 18.2% and 21.9% of the variation in the ECM richness, respectively, and could only explain 19.6% and 33.9% of the exploration types differentiation,respectively. These results indicate that factors other than the parameters measured in the present study might have a significant effect on the ECM community. Many studies have focused on abiotic factors, which are susceptible to slope aspects, such as soil temperature (Ebel 2012), light conditions (Kummel and Lostroh 2011), and microenvironment (Zumsteg et al. 2013). Besieds, the characteristics of symbiosis between ECM fungi and host plants should not be ignored. Future research should focus more on biotic factors that may be affected by the slope aspect, such as the state of the host’s root growth (Jumpponen et al. 2010), the host plant metabolic activities (Buee et al. 2005; Voriskova et al.2014), the interspecific ECM competition (Kennedy 2010),and the presence of a helper bacterial community (Mujic et al. 2016).

Conclusions

This is the first study that has investigated the effects of slope aspect on ECM fungal communities associated with Q. variabilis in the Taihang Mountain, China. The ECM community was found to be dominated by Russula , Inocybe,Tomentella, Scleroderma, and Cortinarius. The OTU richness, Shannon-Wiener, and Shannon evenness of ECM communities on south-facing slope were 1.6 times, 1.6 times and 1.4 times higher than those of the north-facing slope,respectively. Our study also emphasized the direct effects of the slope aspect on the ECM richness and exploration types. Simultaneously, the aspect-induced edaphic variables, such as TP, TN, TK, pH, and SWC, were important sources of variation in ECM fungal richness and distributions of exploration types. Slope aspect and edaphic factors could explain 18.2% and 21.9% of the variation in the ECM richness, respectively, and could explain 19.6% and 33.9%of the variation in the distributions of exploration types,respectively. At the genus level, Cenococcum, Genea, and Clavulina were significantly enriched in north-facing slopes,while Geopora, Helvelosebacina, Scleroderma, Gyroporus,Astraeus, Boletus, Tricholoma, Hebeloma, Cortinarius and unclassified Thelephoraceae were more abundant in southfacing slopes. Hydrophobic ECM fungi were obviously enriched in the south-facing slope, but there was no statistical difference between hydrophilic among the south- and north-facing slopes. This study has deepened our knowledge of the changes in ECM communities associated with Q. variabilis caused by the slope aspect.

AcknowledgementsWe are grateful to Shanghai Majorbio Biopharm Technology Co., Ltd for providing the free online platform of Majorbio I-Sanger Cloud Platform for us to perform some data processing. And we thank American Journal Experts (AJE) for English language editing.