Forest Age and Plant Species Composition Determine the Soil Fungal Community Composition in a Chinese Subtropical Forest

Fungal diversity and community composition are mainly related to soil and vegetation factors. However, the relative contribution of the different drivers remains largely unexplored, especially in subtropical forest ecosystems. We studied the fungal diversity and community composition of soils sampled from 12 comparative study plots representing three forest age classes (Young: 10–40 yrs; Medium: 40–80 yrs; Old: ≥80 yrs) in Gutianshan National Nature Reserve in South-eastern China. Soil fungal communities were assessed employing ITS rDNA pyrotag sequencing. Members of Basidiomycota and Ascomycota dominated the fungal community, with 22 putative ectomycorrhizal fungal families, where Russulaceae and Thelephoraceae were the most abundant taxa. Analysis of similarity showed that the fungal community composition significantly differed among the three forest age classes. Forest age class, elevation of the study plots, and soil organic carbon (SOC) were the most important factors shaping the fungal community composition. We found a significant correlation between plant and fungal communities at different taxonomic and functional group levels, including a strong relationship between ectomycorrhizal fungal and non-ectomycorrhizal plant communities. Our results suggest that in subtropical forests, plant species community composition is the main driver of the soil fungal diversity and community composition.


Introduction
Fungi are a highly diverse component of soil microbial communities. They play essential roles in many aspects of ecosystem development, functioning and stability [1]. As pathogens and symbionts, soil fungi influence the species composition and dynamics of plant communities [2,3]; as saprotrophs, they are an essential component of the soil food web influencing nutrient cycling and carbon sequestration [4]. About 100,000 fungal species have been described, but it has been estimated that there may be from 1.5 to 5.1 million extant fungal species [5]. Most investigations on fungal diversity and distribution have been accomplished in temperate, tropical and boreal regions (e.g. [6][7][8]). Consequently, there is a limited knowledge in the subtropics, although subtropical forest ecosystems are among the most prominent plant diversity hotspots in the holarctic realm [9].
It has been suggested that high plant diversity supports high fungal diversity [6,10,11], implying that subtropical ecosystems might harbor a high fungal diversity [12]. Recent studies on the fungal diversity in Chinese subtropical forests focused on fungal biomass [13], ECM sporocarp and mycorrhizal root analyses [14][15][16], and on the cultivable active microbial community [17], but there is still a gap of knowledge on the soil fungal diversity and the factors influencing the community composition.
In forest ecosystems, fungal diversity and community composition have been reported to be closely linked to numerous biotic and abiotic factors, such as stand age [17,18], elevation [19,20], plant diversity [6,21], plant productivity [2,22], and soil environment [23,24]. The impact of these factors on the soil fungal community has rarely been studied in combination, particularly in the very species-rich subtropical natural forests. In forest communities dominated by only one or a few tree species, such as temperate or boreal forests, it is expected that the soil fungal community is strongly dependent on the identity of a particular tree species. For example, it has been shown that different combinations of plant species have a paramount influence on soil fungal community composition [25]. Consequently, species identity effects of a particular tree species cannot be separated from diversity effects (compare [26]). In tropical forests with 6 to 9 different species per 10 sampled individuals [27] and in subtropical forests with 4.4 to 8.3 different species per 10 sampled individuals, as is the case in our study area (Helge Bruelheide, personal communication), plant species identity effects are expected to be an important factor influencing the fungal community.
In general, a high diversity in plant species is expected to result in high diversity levels at other trophic levels, as diverse plant communities offer a larger amount and larger heterogeneity of resources to consumers, parasites and symbionts [28,29]. In contrast to simply focusing on richness relationships, comparisons of communities at different trophic levels have the advantage that the data also accounts for structural factors and idiosyncratic links between species. The latter seems to be particularly important for the analysis of specificity in plant species -microbial community relationships [30,31], and suggests concordance in the composition of tree species and fungal communities.
We assessed the diversity and community composition of soil fungi using a natural species richness gradient along a chronosequence representing three forest age classes in a Chinese subtropical forest located in the Gutianshan National Nature Reserve (Zhejiang, China). The study plots were established in the framework of the Biodiversity and Ecosystem Functioning project (BEF-China) [32]. As soil parent material and other topographic variables along this chronosequence did only vary to a minor degree, we expected that most of the variation in the fungal community composition should be caused by the biotic components of the ecosystem. To obtain sufficient depth in our analyses, we employed massively parallel and targeted ITS rDNA pyrotag sequencing, an approach currently used in a number of studies investigating soil bacterial [33,34] and fungal communities [35,36]. The objectives of this study were to characterize fungal diversity in this subtropical forest ecosystem, to assess the fungal community composition among the forest age classes, and analyze the contribution of the soil environment, plant species richness and composition on the fungal community. In particular, we hypothesized that 1) fungal community composition differs among the forest age classes, 2) plant cover and variables related to plant biomass explain most of the variation in the soil fungal community composition, and 3) that there is a positive relationship between plant species and soil fungal community composition.

Ethics Statement
All necessary permits for the described field studies were issued by the Administration Bureau of the Gutianshan National Nature Reserve, Zhejiang, China.

Study Site and Soil Sampling
The study was conducted at the Gutianshan National Nature Reserve (NNR) in the Zhejiang Province in south-eastern China (29u8'18" -29u17'29" N, 118u2'14" -118u11'12" E). The NNR is approximately 81 km 2 in size and is located in a mountainous region with a typical subtropical climate. The mean annual temperature at the NNR is 15.3uC with a maximum of 38.1uC in July and a minimum of 26.8uC in January [37]. The annual mean precipitation is 1964 mm (calculated based on data from 1958 to 1986), occurring mostly between March and September [38]. Approximately 57% of the reserve is natural forest [38]. The NNR vegetation type is representative of typical subtropical forest ecosystems consisting of mixed evergreen broad-leaved species. A total of 111 woody species, including 24 ectomycorrhizal (ECM) and 87 non ECM tree species, from 41 families have been recorded at the study plots, with Castanopsis eyrei (Fagaceae), Daphniphyllum oldhamii (Daphniphyllaceae), and Schima superba (Theaceae) as the dominant woody species [32].
The study was carried out in 12 Comparative Study Plots (CSPs), which were randomly selected and stratified by successional age [32] (Fig. S1). Each CSP had a size of 30 m630 m, divided into 9 subplots with 10 m610 m size representing three age classes since natural regeneration (Young: 10-40 yrs; Medium: 40-80 yrs; Old: $ 80 yrs, with four replicate plots each). Abundance of all tree (.1 m in height) species was assessed in the entire plot, while herb species abundance (,1 m in height) was only surveyed in the central subplot. Forest floor litter biomass was determined with a PVC ring (19 cm in diameter) driven into the undisturbed litter layer in spring 2009. In each CSP, 4 litter cores were taken at randomly selected locations (total area of the 4 cores = 0.11 m 2 ) and pooled for dry weight measurement. Soil sampling was carried out between March and April 2009, where within each CSP soil samples were collected to 10 cm depth in each of the 9 subplots using an auger with a 10-cm diameter. The 9 subsamples were then bulked into one composite sample per study plot. Homogenized samples were immediately hand-sieved (#2 mm) in the field to remove stones, roots, macrofauna, and litter materials. Subsamples were frozen in liquid nitrogen, transported in an ice box and were stored at 220uC until molecular analysis. The analysis methods of soil characteristics and the vegetation data (presented in Table 1) are described in detail in [13].

DNA Extraction, Amplicon Library Generation and Pyrosequencing
Soil microbial genomic DNA was extracted from approximately 1 g (wet weight) of soil using a MoBio PowerSoil DNA Isolation Kit (MoBio Laboratories Inc. Carsbad, CA, USA) following the manufacturer's instructions. Fungal ITS rDNA amplicon libraries were produced using fusion primers [39] designed with pyrosequencing primer B, a barcode and the fungal specific primer ITS1F [40] as a forward primer and pyrosequencing primer A and the universal eukaryotic primer ITS4 [41] as a reverse primer. We used a set of 10 bp MID-barcodes provided by Roche (Roche Applied Science). ITS rDNA amplicon libraries were produced from a pool of two dilution levels (106 and 1006) from each soil DNA extract. The PCRs were performed in three replicate reactions per sample and per dilution to account for potentially heterogeneous amplification.
PCR reaction was carried out in a total volume of 50 ml containing 1 ml diluted DNA template, 1 ml 25 pmol of each of the two custom fusion primers, 25 ml Go Tag H Green Master mix (Promega), and nuclease free water (Promega). We used a touchdown PCR program with a denaturation at 95uC for 5 min followed by 10 cycles of denaturation at 94uC for 30 sec, annealing at 60-50uC for 45 sec (21uC per cycle), and extension at 72uC for 2 min; and then 30 cycles at 94uC for 30 sec, 50uC for 45 sec and 72uC for 2 min, then finalized by an extension step at 72uC for 10 min. PCR products were analyzed using 1.5% agarose gel and equimolar volumes of the amplified products of the expected size (ca. 600 bp) from the three positive replicate amplicons per sample were homogenized. The pooled products were gel purified using a Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany). The amount of DNA in the purified amplicons was quantified using a fluorescence spectrophotometer (Cary Eclipse, Agilant Technologies, Waldbronn, Germany). An equimolar mix of the 12 amplicon libraries was subjected to unidirectional pyrosequencing from the ITS1F end of the amplicons using a 454 titanium amplicon sequencing kit and a Table 1. Elevation of study sites, plant parameters and soil characteristics of the comparative study plots (CSPs) and their relationships with forest age class based on Pearson correlation analysis.

Bioinformatics Analysis
Multiple levels of sequence processing and quality filtering were performed as described in [42]. Briefly the 454 fungal ITS sequences were initially extracted based on 100% barcode similarity. Simultaneously, sequence reads with an average quality score of less than 25, read length of ,200 bp, ambiguous bases and homo-polymers of .8 bases were removed, barcodes and primers were trimmed using the split libraries script available in the Quantitative Insights In Microbial Ecology pipeline (QIIME) [43]. Based on preliminary sequence analysis, reads .450 bp were trimmed to a maximum read length of 450 bp using MOTHUR [44]. Sequences were then clustered and assigned to operational taxonomic units (OTU) using the QIIME implementation of CDHIT with a threshold of 97% pairwise identity after a pre-filtering step to remove identical sequences using default parameters. The most abundant representative sequences were selected and assigned to respective taxa according to NCBI taxonomy. The searches were based on blastn against the NCBI nucleotide database for fungi excluding uncultured and environmental sequences using the blast based taxonomic assignment script of the software for Cleaning and Analyzing Next Generation Sequences -CANGS [45].
Representative sequences of OTUs of the fungal ITS pyrotags assigned under the fungal kingdom were checked for chimeras using the chimera uchime algorithm using the same dataset as a reference, as implemented in MOTHUR. Subsequently a total of 2398 sequences identified as potentially chimeric and had less than 90% alignment length to sequences in the NCBI fungal reference database, and 1445 reads not assigned to the fungal kingdom were removed from the sequence dataset.
Finally, we found 14,136 sequences that were grouped into 1260 OTUs with variable number of reads per sample (Table S1). Consequently, the number of sequences per sample was normalized to the smallest sample size, 872 reads per sample, using the normalized shared command as implemented in MOTHUR.

Ectomycorrhizal Fungi Designation
The Ascomycotan and Basidiomycotan fungal OTUs were further identified as putative ectomycorrhizal (ECM) fungi at the genus level based on literature [46][47][48][49][50] as described in [50]. All ECM genera which were presented as non mycorrhizal (NM) by Tedersoo et al. [50] in contrast to Rinaldi et al. [48], were treated as NM. Furthermore, we did a manual NCBI blast search for the representative sequences of those genera reported to be composed of both ectomycorrhizal and saprotrophic species. Accordingly only those OTUs with blast similarities of .97% with sequences derived from mycorrhizal roots were maintained in our final putative ECM fungal dataset.

Statistical Analysis
The fungal OTUs were parsed by sample in order to calculate the abundance of fungal OTUs using the sequence count of each of the non-singleton OTUs as abundance value [42,51]. Based on a preliminary rank index analysis, we calculated dissimilarities between all pairs of samples using the log (x+1) transformed abundance data and Bray-Curtis dissimilarity coefficient in order to obtain an abundance based dissimilarity matrix. To assess the effect of singletons on the fungal community distribution, we calculated the non-metric multidimensional scaling (NMDS) ordinations with 20 random starts from the datasets with and without singletons. The correlation between the ordinations was tested using the Procrustes correlation analysis using the protest function [52] of the vegan package [53], where the significance of the congruence between any two ordinations was tested by a Monte Carlo procedure with 999 permutations. We found that the fungal community composition was not affected by the presence or absence of singletons (Procrustes correlation coefficient = 0.966, P,0.001, suggesting nearly identical ordinations). Similar analysis using the presence or absence dataset also showed that the NMDS ordinations were significantly correlated (Procrustes correlation coefficient = 0.97, P,0.001). Thus we performed the subsequent analyses using the dominant fungal community excluding singletons.
We used the function ANOSIM of the vegan package to explore the similarity of fungal community composition among the forest age classes. Alpha and beta diversity based fungal community compositions across the three age classes were compared using the abundance-based pair-wise Bray-Curtis dissimilarity and the Sorenson pair-wise dissimilarity matrix accounting for beta diversity using the betapart package [54]. The relative abundance of ECM fungi at the family level and the distribution of one of the most abundant ECM family at the OTU level across the three forest age classes was visualized by a heatmap using the function heatmap2 of the gplots package [55].
The vegetation and soil characteristics were checked for collinearity with the function varclus in the Hmisc package [56], using Spearman's rank correlation. The resulting set of noncollinear parameters was maintained for further analyses (Table 1 and Fig. S2). Based on preliminary rank index analysis of these exploratory variables, non-metric multidimensional scaling (NMDS) analysis was performed using the nmds function of the labdsv package [57] based on Bray-Curtis distances to visualize the distribution of the fungal communities. The function envfit of vegan was used to calculate a post-hoc regression of the individual environmental variables on the ordination scores. R 2 or the goodness-of-fit values and their significances were calculated using 999 random permutations. Pearson correlation coefficients were calculated to elucidate patterns of correlations of environmental variables using the function corr.test of the Psych package [58].
Model of multivariate analysis of variance was constructed using distance based redundancy analysis (dbRDA) using Bray-Curtis distances with the function capscale of vegan to determine the most influential environmental variables on the fungal community composition. Permutational multivariate analysis of variance (Permanova) was used to determine the interaction of the selected variables and assess their influence on the fungal community composition using the adonis function of vegan with the Bray-Curtis distances and 999 permutations.
The concordance of plant and fungal community ordinations was assessed using Procrustes correlation analysis as described above using plant and fungal NMDS ordinations. The contribution of individual soil characteristics to the plant and fungal community congruence was tested using ANOVA of Procrustes residuals. Based on the preliminary normality tests, in all the above analysis both plant species abundances and environmental variables were log (x+1) transformed. Statistical analysis was performed using R version 2.15.2 [59].

Pyrosequencing Data Analysis and Taxonomic Assignment
A total of 16,534 fungal ITS pyrotag reads with an average length of 450 bp were obtained from the 12 soil samples collected across the three forest age classes. After multiple levels of sequence processing, quality filtering, and sequence number normalization followed by a 3% dissimilarity clustering we found 1027 fungal OTUs including 457 (44%) singletons from 10,464 fungal ITS reads (Table S1). Taxonomic distribution of the 570 abundant fungal OTUs, excluding the singletons, showed distribution of the fungal community mainly across four fungal phyla with 50.2% belonging to Basidiomycota, 37.3% to Ascomycota, 2.3% to Zygomycota, and 0.4% to Glomeromycota. We were able to assign 73.7% and 68.8% of these OTUs to the family and genus levels, respectively. The remaining 9.8% OTUs were grouped as unidentified fungi indicating the presence of unknown fungal taxa in this particular subtropical forest ecosystem.
Within the two major fungal phyla, the Ascomycota comprised of 213 OTUs, which were mainly members of the subphylum Pezizomycotina (180/213). A total of 35 ascomycetous families were identified among which the Trichocomaceae, Myxotrichaceae, Elaphomycetaceae, Hypocreaceae, and Herpotrichiellaceae were the most represented. From the total of 286 OTUs identified as members of the phylum Basidiomycota, the majority (282 OTUs) was assigned to the subphylum Agaricomycotina. A total of 36 basidiomycetous families were found, with Russulaceae, Sebacinaceae, Cortinariaceae, Thelephoraceae, Amanitaceae, Tricholomataceae, and Clavulinaceae being the most abundant ones.
It is also noteworthy that, within the two dominant fungal phyla, 10% of the 213 ascomycetes and 75% of the 286 basidiomycetes fungal OTUs were identified as potential ectomycorrhizal (ECM) fungi. The ECM fungal OTUs (234 OTUs) represented 22 families (5 belonging to Ascomycota and 17 to Basidiomycota) and 34 genera (5 Ascomycota and 29 Basidiomycota). The relative abundance of ECM fungal communities at the family level showed that Russulaceae, Thelephoraceae, Helotiaceae, Amanitaceae, Sebacinaceae, Entolomataceae, Cortinariaceae, Elaphomycetaceae and Ceratobasidiaceae were the most abundant families across the three forest age classes. At the genus level, Russula, Tomentella, Amanita, Thelephora, Entoloma, Sebacina, and Cortinarius were the most abundant (Fig. S3 a, b).

Fungal Community Assemblage across the Three Forest Age Classes
The presence or absence of singletons showed no significant effect on the fungal community ordination at the fungal Kingdom level (see Material and Methods). These ordinations were also consistent when confined to the phylum Ascomycota (Procrustes correlation coefficient = 0.742, P,0.01), phylum Basidiomycota (Procrustes correlation coefficient = 0.929, P,0.001), and ECM fungi (Procrustes correlation coefficient = 0.985, P,0.001) levels tested. Analysis of similarity of the dominant fungal community assemblages revealed significant differences of the fungal community across the three forest age classes (ANOSIM, R = 0.345, P,0.01), which was consistent in the phylum Basidiomycota (ANOSIM, R = 0.24, P,0.05) and ECM (ANOSIM, R = 0.36, P,0.05) fungal communites.
Both abundance and incidence based pair-wise fungal community dissimilarity analysis indicated that fungal alpha diversity (adonis R 2 = 0.235, P, ,0.01) and beta diversity (adonis R 2 = 0.231, P,0.05) were significantly correlated with forest age. Hierarchical cluster analysis of these dissimilarity matrices also revealed a clear cluster of the young forest plots (Fig. 1 a, b). Consistently our results showed the presence of age class related fungal community distribution in this forest ecosystem. We found a total of 20, 15 and 27 ascomycetous and 60, 43 and 49 basidiomycetous unique fungal OTUs in the young, medium and old forest age classes, respectively. The ten most abundant and unique ascomyceteous and basidiomycetous fungal OTUs in each forest age class are presented in Table S2.
Complete hierarchical cluster analysis based on the relative abundance of the ectomycorrhizal fungal families and visualization using a heatmap also indicated age class related distribution of some ECM families (Fig. 2a). For example, members of the Lyophyllaceae and Gomphidiaceae were detected only in the young age class forest while Russulaceae (142 OTUs) and Thelephoraceae (80 OTUs) were the most abundant groups in all the three forest age classes. Similar visualization of ECM fungal OTUs in both families indicated the presence of age class related ECM fungal OTUs (the heatmap for the family Russulaceae is depicted in Fig. 2b). All detected ascomycetous ECM fungi appeared generally across the forest age classes, while some basidiomycetous ECM fungi were found only in the young or shared between the medium and the old forest age classes (Table S3).

Relationships between Fungal Communities and Environmental Variables
NMDS analysis followed by environmental variable fitting to assess the relationship of individual variables to the ordination plot Figure 1. Complete hierarchical clusters based on pair-wise dissimilarity matrices derived from (a) abundance based Bray-Curtis dissimilarity matrix (b bray ) and (b) incidence based Sorenson dissimilarity matrix accounting for beta diversity (b sor ). Both cluster diagrams showed that the observed significant effect of forest age to the fungal alpha (adonis R 2 = 0.235, P, ,0.01) and beta diversity (adonis R 2 = 0.231, P,0.05) differences is mainly attributed to the fungal communities of the young age class forests. doi:10.1371/journal.pone.0066829.g001  indicated that the fungal community composition was significantly related to forest age and to plant and soil parameters ( Table 2, Fig. 3a). In contrast to the ascomycetous communities the basidiomycetous and ECM fungal community ordinations were influenced by forest age and plant parameters (Table 2, Fig. 3b).
The dbRDA based model selection, however, indicated SOC and elevation to be the most important factors shaping the fungal community composition (F = 1.63, P,0.01) ( Table 3). Marginal tests also showed that SOC and elevation were significantly related to the fungal community composition (SOC -F = 1.84, P,0.01; elevation -F = 1.42, P,0.05). We also found a significant interaction of forest age with SOC and elevation in the fungal (adonis F = 1.37, P,0.05) and ascomycetous (adonis F = 2.01, P = 0.002) communities (Table 3). These results indicated the role of a particular or specific group of environmental drivers that shape the diversity and composition of fungal communities.

Relationship between Fungal and Tree Communities
The correspondence of the fungal and tree community as revealed by the Procrustes correlation analysis showed a significant relationship across the study sites, at the level of the fungal  kingdom (Procrustes correlation coefficient = 0.664, P,0.01).
Consistently we found a significant correlation of the plant community to the Basidiomycotan and ECM fungal communities. Similar analysis based on ECM and non-ECM tree communities demonstrated a significant agreement with the fungal community. Contrary to our expectations the concordance between ECM fungi and ECM tree communities was not significant (Table 4, Fig.  S4). ANOVA based tests to assess the influence of soil characteristics (Table 1) on the relationships between fungal and tree community assemblages showed no significant effect of the soil parameters on the observed concordance between tree and fungal communities (Bonferroni corrected P.0.05).

Fungal Diversity and Taxonomic Assignment
In this study we found evidence for a highly diverse soil fungal community in sub-tropical forests of China. Despite our stringent sequence quality filtering and normalization of sequence reads per sample we found 1027 fungal OTUs (including 457 singletons). The number of fungal OTUs we found in this forest ecosystem, without sequence read normalization as presented in many of the published works, is relatively higher than previously reported fungal diversity records of 1077 fungal OTUs from Quercus spp. ectomycorrhizas [60], 47 arbuscular mycorrhizal fungal OTUs from root samples of 10 plant species [61], and up to 1000 fungal OTUs per 30,000 reads from European spruce and beech forest soil [35].
The fungal communities were predominantly members of the phylum Basidiomycota, followed by Ascomycota and Zygomycota. This observation is in accordance with results from European forests [35]. In contrast, in Australian subtropical forests Curlevski et al. [62] showed a dominance of Ascomycota over Basidiomycota and Zygomycota, while He et al. [63] reported Zygomycota as the dominant phylum. Such diverse results may indicate lack of a general global pattern in distribution and abundance of the main fungal phyla in forest soils. Although this study has a limitation due to the relatively low number of final sequences per sample, the observed differences with other studies could also be attributed to differences in DNA extraction protocols, differential amplification of different PCR primers and PCR conditions used in the different studies. This suggests a need for larger scale investigations with standardized approaches [64].
In the investigated Chinese subtropical forest ecosystem dominated by non-ECM trees, which accounts for 78.4% of the plant species, we found 75% of the basidiomycetous and 10% of the ascomycetous OTUs to be putative ECM fungi. The result corroborates the observed prevalence of basidiomycetous ECM fungi based on ECM root analysis in a subtropical forest in Dujiangyan [16] and supports the observation that ECM fungi tend to have higher diversity in mixed forests with a higher proportion of non-ECM tree species as compared to forests composed of exclusively or dominantly ECM trees (e.g. [65]). Here, we observed that Russulaceae and Thelephoraceae were the most abundant ECM fungal families. This observation has also been documented in tropical [7,8,11], subtropical [16,66,67], Mediterranean [68], and temperate [16,30] forest ecosystems, suggesting a global pattern of distribution for these two ECM fungal families.

Relationships between Fungal Communities and Forest Age Classes
Consistent with previous studies we found a strong influence of forest age on the fungal community composition [16,18,69]. Our results showed significant differences in the fungal community of the kingdom fungi, phylum Basidiomycota and ECM fungi between the three age classes. The observed clear cluster of the young forest plots (Fig. 1) is mainly related to the decline in the mean Shannon and Simpson diversity indices of the fungal and ECM communities from young to old age class forests. This indicates that some species might become dominant in the fungal community with time. This observation is consistent with that of Keizer and Arnolds [70] and Wallander et al. [18], who suggested that the ECM fungal diversity increases fairly rapidly in the first 30-40 years of forest development and gradually decreases to an intermediate or rather constant level. However, confirmation of this pattern of change would require a longer successional sequence [17,71,72].
Detailed analysis of the ECM fungal communities at the family and genus levels also revealed distinct fungal communities in the three age classes. The observed pattern of ECM fungal distribution indicates the importance of forest age in subtropical forest ecosystems and suggests succession of fungal communities during forest development. This pattern forms a strong contrast to the plant communities, where among the 148 tree species found in our study forest, only one shrub (Photinia glabra) and two woody species (Glochidion puberum, Platycarya strobilacea) were significantly concentrated in the oldest and youngest forests respectively [32].

Relationships between Fungal Communities and Environmental Variables
We found a strong relationship between individual environmental variables such as tree layer cover, woody plant biomass, litter layer, dead wood cover, sand and SOC with the fungal community composition. These observations corroborate with previously reported positive correlations between the index of plant primary production with fungal diversity and community composition [73]. In accordance with the known role of soil fungi in decomposition of leaf litter and dead wood in forest ecosystems [22,74,75], we also found a significant contribution of the litter and dead wood cover on the basidiomycetous and ECM fungal community composition. Analysis of the role of these environmental variables on the fungal community composition accounting the effects of other parameters indicated that forest age, elevation, and SOC are the most important variables. This result is in line with our observation of the strong correlation of forest age with tree layer cover, woody plant biomass, and dead wood cover, while SOC is correlated with sand and litter layer (see Table S3). Our finding of elevation as one of the important variables shaping fungal community composition is in agreement with Bahram et al. [20]. Considering the strong correlation between forest age, herb layer cover, and herb species richness with the elevation of the study plots and the suggestion of Bahram et al. [20] to take elevation as a proxy for environmental variables such as precipitation and temperature, our result indicates, the need for further large-scale study considering more climatic variables.

Relationship between Fungal and Tree Communities
In accordance with our hypothesis, the Procrustes correlation analyses revealed a remarkable congruence between tree and fungal communities, except for the ECM tree and ECM fungal communities. The absence of concordance between ECM fungal and ECM tree communities in this study contradict the reports from boreal forests where close linkages between these groups were documented [25,76]. This observed lack of correlation could be attributed to our sampling strategy. First although we sampled 9 subsamples per plot, bulking of the samples would have homogenized the fungal diversity and reduced the plant species effect on the fungal community. Second since ECM fungal distribution shows relationship with soil horizons [77], our sampling depth might have restricted the number of plant species specific ECM fungal species recovered in this study. Nevertheless, our results strongly support previous observations in tropical [8] and subtropical [14] forests, and indicated that the spatial pattern of plant communities, proportion of ECM and non-ECM plants as drivers of the fungal community composition.

Conclusion
This study provides evidence of high soil fungal diversity in this subtropical forest ecosystem. In general our findings suggested the existence of plant and soil related effects on soil fungal community composition. Forest age, elevation of the study plots and SOC were found to be the most important factors influencing the fungal community composition. The observed significant correspondence between tree and fungal communities and the relative contribution of the ECM and non-ECM tree species communities on the fungal community composition suggests the existence of unidentified specific links of resources provision and utilization between plants and rhizosphere fungal communities in this subtropical forest. Elucidating these links will require experimental testing in manipulated systems that vary in plant species diversity, which is the main focus in the ongoing BEF China project.