Diversity of root-associated fungi of Rhododendron simsii in subtropical forests: fungal communities with high resistance to anthropogenic disturbances

Species of the Ericaceae or heath family are widely distribute in continental ecosystems and their special ericoid mycorrhizas (ERM) are considered beneficial to their survival and persistence in variable habitats. Currently, increasing anthropogenic disturbances and improper forest management are affecting subtropical forests of China where these native species located. These activities not only affect plant communities above-ground, but also impose pressures on microbial communities below- ground. In this study, root-associated fungal communities of Rhododendron simsii in four forest types under different anthropogenic disturbances were identified using an Illumina Miseq platform, i.e., old growth forests, secondary forests with one cutting (SEC I), secondary forests with two cuttings (SEC II), and Chinese-fir plantations (PLF). Intra- and inter-annual variations were analyzed by comparing samples taken in different seasons and years. The results show that: (1) over 1000 OTUs were found in hair roots with most from the division Ascomycota and Basidiomycota belonging to different functional groups; (2) while there were a few indicator OTUs specific to different forest types, seasons and years, the proportion of shared taxa was quite large, accounting for 44.9–79.4% of the total OTUs; (3) significantly positive correlations were found between disturbance sensitivity and temporal variations in common fungal orders, and both in major fungal orders were significantly different among fungal functional groups in which putative and possible ERM fungi were highly resistant to disturbances and low temporal variations. The high disturbance resistance and temporal persistence of putative ERM fungi may be essential for the successful adaptation of R. simsii in disturbed subtropical forests of China.


Introduction
Rhododendron simsii Planch. is popular as a landscape plant because of its beautiful flowers. It occurs naturally in China, Japan, Laos, Myanmar and Thailand. In most southern provinces of China, R. simsii populations occur in subtropical forests, along forest margins and in open upland thickets at 500-2700 m a.s.l. (Editorial Committee of Flora of China, Chinese Academy of Sciences 2006). However, their natural habitats are frequently disturbed by anthropogenic activities. Most Rhododendron species have typical ericoid mycorrhizas (ERM) (Fumiaki et al. 2003;Zhang et al. 2009;Tian et al. 2011), which are considered to be helpful for both the host plant and their fungal partners to be distributed widely in the world (Cázares et al. 2005;Bougoure et al. 2007;Fujimura and Egger 2012;Gorzelak et al. 2012).
Abstract Species of the Ericaceae or heath family are widely distribute in continental ecosystems and their special ericoid mycorrhizas (ERM) are considered beneficial to their survival and persistence in variable habitats. Currently, increasing anthropogenic disturbances and improper forest management are affecting subtropical forests of China where these native species located. These activities not only affect plant communities above-ground, but also impose pressures on microbial communities below-ground. In this study, rootassociated fungal communities of Rhododendron simsii in four forest types under different anthropogenic disturbances were identified using an Illumina Miseq platform, i.e., old growth forests, secondary forests with one cutting (SEC I), secondary forests with two cuttings (SEC II), and Chinesefir plantations (PLF). Intra-and inter-annual variations were analyzed by comparing samples taken in different seasons and years. The results show that: (1) over 1000 OTUs were found in hair roots with most from the division Ascomycota and Basidiomycota belonging to different functional groups;

3
Anthropogenic disturbances and particular types of forest management such as clear-cutting, selective cutting, planting of high-value species, not only influence plant communities above-ground but also underground fungal communities associated with plant roots. For example, ERM fungal communities under different vegetation types such as bogs, rough grazing and forest plantations (Hazard et al. 2014), and arbuscular mycorrhizal fungi (AMF) along an anthropogenic disturbance gradient (Gavito et al. 2008), show significant changes. The results of 378 case studies summarized by Shade et al. (2012) indicated that most soil microbial communities were sensitive to disturbance. In addition, only 23% of microbial communities recover based on their composition and/or function from 148 cases involved in resilience, (the rate of recovery after disturbance), suggesting that disturbance plays a major role in mycorrhizal fungal and soil microbial communities, but recovery depends on specific conditions. In this study, intra-and inter-annual variations of fungal communities associated with hair roots of R. simsii growing in Chinese subtropical forests were examined across different types of anthropogenic disturbances. Three hypotheses are proposed: (1) fungal communities associated with hair roots of R. simsii would be different along the disturbance gradient, across seasons and years; (2) different fungal taxa would show different responses to disturbance and temporal variations; and, (3) disturbance resistance and temporal persistence of fungal functional groups might be correlated, and ERM fungi would be more consistently associated with their R. simsii hosts than any other fungal groups.

Study site
Gutianshan National Nature Reserve in Zhejiang Province (GNNR) was selected as the sampling site for this study. Annual average temperatures are above 15 °C and Annual precipitation are more than 1750 mm. The vegetation of the site is subtropical evergreen broad-leaf forest (Yu et al. 2001), and the main soil type is subtropical red soil  Due to different disturbance histories, there are four distinct types of stands in the GNNR: Chinese-fir (Cunninghamia lanceolata (Lamb.) Hook) plantations (PLF), old growth forests (OGF), secondary forests with one cutting (SEC I), and secondary forests with two cuttings (SEC II). The PLF were planted 20 years ago after clear-cutting of secondary growth forests. The OGF are undisturbed forests located at the heart of the GNNR and have not been harvested or thinned in the past 100 years. Both SEC I and SEC II were clear-cut about 50 years ago, and SEC II was selectively cut 20 years ago. Since the previous disturbances, all stands have naturally regenerated. Twelve 1-ha (100 m × 100 m) plots were established with three replicates for each forest type.

Sampling procedure
Hairy root samples of four R. simsii plants were collected from each sample plot over six periods, March, June, September and December of 2012, and March 2013, and March 2014. There were 288 individual samples (4 samples × 12 plots × 6 sampling times). The fine roots were collected from soils at four directions around the taproots. After soaking and cleaning in sterile water, 20 (about 1-cm long) hair root segments were randomly selected and put separately into a 1.5 mL centrifugal tube for DNA extraction. A 200-g soil sample was collected for elemental analyses from the top 10-cm layer of soil around each sample plant.

DNA extraction and Illumina MiSeq sequencing
DNA was extracted according to Zhang et al. (2016). The ITS1 region special for fungi was selected to be amplified by PCR with 1723F and 2043R primers; the PCR reaction and amplicons purification procedures are based on Zhang et al. (2017). The Illumina MiSeq sequencing were completed by Shanghai Majorbio Co. Ltd. in China. Raw fastq files were demultiplexed and quality-filtered using QIIME (ver 1.7) following Zhang et al. (2016).
Open reference operational taxonomic unit (OTU) picking was done with pick_open_reference_otus.py using the default uclust method and ITS 12-11 dataset (97% similarity), and singletons were removed during OTU picking. Sin-gle_rarefaction.py in Qiime was used to generate an OTU table with even reads of 10,000 in each root sample. Rare fungal OTUs occurring in only one or two samples out of 72 were not considered in further analyses.
Mean relative abundances (MRA d ) of fungal orders under the four forest types were calculated and the disturbance sensitivity of each fungal order defined as the coefficient of variance (CV d ) of MRA d . Mean relative abundances (MRA t ) of fungal orders in six sampling times were also calculated and the temporal variation of each fungal order was defined as the coefficient of variance (CV t ) of MRA t of sampling times.

Statistical analysis
The differences population densities and mean DBH (diameter at breast height) of R. simsii between the four forest types was analyzed using one-way ANOVA with SPSS 22.0. The otu_category_significance.py in Qiime analyzed and determined the indicator fungal species for the different forest types, seasons, and years. The main influencing factors of root-associated fungal communities were identified with "Envfit" function. The "rda" function in the R package "vegan" performed principal components analysis (PCA), and the first three components (PC1, PC2, and PC3) were used as plant parameters in "envfit". The impacts of forest type, season and year on the community composition of root-associated fungi were evaluated with the "Adonis" function in the R package "vegan".

Populations of Rhododendron simsii in anthropogenic disturbance sites
The plant compositions in OGF and two secondary forests (SEC I and SEC II) are similar, but there are more large-diameter trees in the OGF. The dominant species (Chinese-fir) and plant composition in PLF are obviously different from the secondary and old growth forests. Mean DBH and population density of R. simsii in the four forest types are shown in Appendix S1. There were no significant differences in these two parameters among the forest types (P > 0.05), although the density was the lowest and mean DBH was the highest in the OGF.

Fungal OTUs in hair roots of Rhododendron simsii
Rarefaction curves of the observed fungal OTUs in hair roots of R. simsii are shown in Fig. 1. Mean observed fungal OTU richness at the sample depth of 10,000 reads were 331, 314, 329 and 268 in OGF, SEC I, SEC II and PLF, respectively. There were no significant differences between any two forest types, although the fungal richness in the PLF samples was the lowest.
The number of specific and shared OTUs of the stands under different anthropogenic disturbances, seasons and years is shown in the Venn diagram ( Fig. 2). Secondary forests shared more OTUs with OGF compared to the PLF (Fig. 2a). The shared OTUs in all stands were 408, accounting for 36.6% of the total number in all forest types and 85.2% of the total reads in each sample. The specific number of fungi were 18, 0, 5 and 24 for the OGF, SEC I, SEC II, and PLF, respectively.
The total number of fungal OTUs were 922, 796, 929 and 883 for the spring, summer, autumn and winter seasons, respectively (Fig. 2b). There were 431 shared OTUs in 2012 over all the seasons, while any two seasons shared over 610 OTUs. More seasonal specific OTUs (35) were found in the spring, compared with the summer, autumn, winter seasons (8, 17, and 12, respectively). From 2012 to 2014, the number of OTUs in the spring seasons was 610, 563 and 571, respectively (

Fungal community structures
The phyla and main classes in fungal community compositions under different forest disturbances are shown in Fig. 3. Ascomycota and Basidiomycota were the two dominant phyla in all stand types, accounting for 95.6-98.8% of the total reads. Other phyla, such as Zygomycota, Glomeromycota and Chytridiomycota were less than 1.2-4.4% of the total. Common classes of fungi were Leotiomycetes, Sordariomycetes, Dothideomycetes, Eurotiomycetes, and Pezizomycetes of the Ascomycota, as well as Agaricomycetes and Tremellomycetes of the Basidiomycota.
Principal component analysis (PCA) shows the community structure of root-associated fungi of R. simsii at different sampling times with different anthropogenic disturbances (Fig. 4). Community structure in spring and winter were similar in the OGF and SEC I, while they were similar in autumn and winter in SEC II and the PLF. The fungal community structure in the summer was always different than in any other seasons. Over the three sampling years, compositions of root-associated fungal communities were The results of the Adonis analysis shows that the effects of season, year, and forest type on the community composition of root-associated fungi of R. simsii were insignificant (P > 0.05) for seasons, (although the number of OTUs in summer was less), but significantly different for years and anthropogenic disturbances (P < 0.01) ( Table 1).

Factors affecting fungal community associated with Rhododendron simsii hair roots
Our results show that plant community, rather than soil parameters and climatic factors, had significant effects on root-associated fungal communities (Table 2). For plant communities, based on principal component 2(PC2), the proportions of AM (arbuscular mycorrhiza), ECM (ectomycorrhizal) and ERM plants were significant different. The plant species contributing significantly to PC2 included two species associated with arbuscular mycorrhizas, Cunninghamia lanceolata and Syzygium buxifolium Hook. et Arn, two dominant species in subtropical forests associated with ectomycorrhizal fungi, Castanopsis eyrei and Lithocarpus glaber (Thunb.) Nakai, as well as a common species associated with ericoid mycorrhizas, Vaccinium carlesii Dunn. These results indicate that both dominant trees and common neighbors associated with ERM may affect the fungal communities of R. simsii hosts. The AM: ECM: ERM ratios also show significant effects on the root-associated fungal communities. Soils or climatic parameters were insignificant on root-associated fungal communities of R. simsii.

Indicator fungal species of stands with different anthropogenic disturbances
There were five fungal OTUs showing significant preference to a specific forest type (Table 3). The number of  indicator species to OGF, SECI, SECII, and PLF were 3, 0, 0, and 2, respectively. There was no significant indicator for seasons. Eleven indicator species emerged in 2013 and 2014. Fungal species in families such as Herpotrichiellaceae 1, Myxotrichaceae 1, and Sebacinaceae 1 and orders such as Helotiales 1, Sebacinales 1 and Sebacinales 2 might be potential ericoid mycorrhizal fungi. Phialocephala fortinii was a putative DSE (dark septate endophyte), Tomentella sp. was a kind of ECMF, Penicillium sp. was a saprophytic species or pathogen, and Lecanicillium sp. was identified as a plant pathogen according to Tedersoo et al. (2014).

Disturbance sensitivity and temporal variations of functional fungal groups associated with hair roots of Rhododendron simsii
Fungal orders associated with hair roots of R. simsii were divided into five groups, e.g., putative ERMF, possible ERMF, putative ECMF, AMF, and fungal orders without identified functions (Table 4). These orders ranged from very sensitive to relatively tolerant of anthropogenic  disturbances and temporal variations. The order Helotiales had the lowest CV d (coefficient of variation to disturbance), indicating that it had the highest tolerance to logging and commercial planting. It also showed that AM fungi were more sensitive to anthropogenic disturbances, which may be due to the appearance of a large number of AMF's hosts introduced by tree planting in PLF. For temporal variations of fungal orders, CV t (coefficient of variation to time) values were (mean ± SE) 54. 0 ± 11.8, 76.5 ± 11.7, 63.5 ± 25.0, and 192.4 ± 33.1 for putative ERM fungi, possible ERM fungi, putative ECM fungi and AM fungi, respectively (Table 4). Significantly positive correlations were found between disturbance sensitivity and temporal variations (R 2 = 0.647, P < 0.001; Fig. 6).

Diversity of fungi associated with hair roots of Rhododendron simsii
Based on high through-put sequencing, Ascomycota and Basidiomycota were the dominant phyla in root-associated fungal communities of R. simsii, while Zygomycota, Glomeromycota and Chytridiomycota were present in relative lower abundance, similar to those observed in Vaccinium spp.; however, the 1115 OTUs in the roots of R. simsii were much less than those on roots of Vaccinium carlesii (over 5000 OTUs) in the same forests (Zhang et al. 2016). Some putative ERMF orders such as Helotiales and Sebacinales were also found on roots of other ericaceous plants (Xiao 1994;Read 1996;Piercey et al. 2002;Midgley et al. 2004;Mark 2006;Selosse et al. 2007;Walker et al. 2011;Weiß et al. 2011), and were also observed in this study. Besides putative ERM fungal taxa, other fungi, for example, DSE (Phialocephala fortinii), ECM fungi (species of Tomentella and Russula), AM fungi (species of Glomeromycota), saprophytic or pathogenic fungi (e.g. species of Penicillium and Lecanicillium), were found in these fungal communities and they were considered to play either positive or negative roles for their hosts in different environments. These fungal taxa were also observed in other research on ericaceous plants (Phuwiwat and Soytong 2001;Zare and Gams 2001;Allen et al. 2003;Walker et al. 2011;Wurzburger et al. 2011;Sun et al. 2012), and their functions on roots of ericaceous plants should be evaluated in further studies.

Determinants of root-associated fungal communities of Rhododendron simsii
PC2 and the AM: ECM: ERM proportions in the forest types, in which the mycorrhizal hosts belonged to different families, showed significant effects on root-associated fungal communities of R. simsii (Table 3). The most significant impact on the proportional changes of plants was the introduction of Cunninghamia lanceolata, the single dominant tree species that could stimulate the plant community to host more AM fungi (Bruelheide et al. 2011). Another impact was from clear-cutting and select-cutting (focused on local dominant tree species), and both would alter the fungal proportions maintained by native dominant tree species together with ECM fungi (i.e., Pinus massoniana Lamb. of the Pinaceae and Castanopsis eyrei and Castanopsis carlesii (Hemsl.) Hay. of the Fagaceae) and AM fungi (i.e., Schima superba of the Theaceae, Daphniphyllum oldhami (Hemsl.) Rosenth. of the Daphniphyllaceae), as well as ERM fungi (i.e., Vaccinium carlesii).
Compared with PC2, more ericaceous species such as Rhododendron ovatum and R. latoucheae, and other species such as Corylopsis glandulifera Hemsl. and Loropetalum chinense (R. Br.) Oliv.of the Hamamelidaceae family significantly contributed to PC3, however, the fungal communities on R. simsii roots have not shown obvious effects. Therefore, root-associated fungi of R. simsii were influenced by the combination of AM, ECM and ERM fungi on plant roots instead of more ERM hosts. It has also been welldocumented that fungal diversity and fungal community composition in soil and on plant roots are driven by plant community composition (Wu et al. 2012), and that high plant diversity and productivity support high fungal diversity (Schadt et al. 2003;Carney and Matson 2006;Dumbrell et al. 2011). From this point of view, it is easy to understand that the total number of OTUs in PLF was less due to its simple plant community composition and its low productivity compared with OGF and the two secondary forests. Fig. 6 Relationship between disturbance sensitivity and temporal variations of major functional fungal orders associated with hair roots of Rhododendron simsii (Yellow is putative ERM order, Green is possible ERM order, Blue is putative ECM order, Red is AM order, White is nonfunctional fungal order) Typical ECM fungi found on the roots of R. simsii and the AM: ECM: ERM proportions of plant community had significant effects on root-associated fungal communities of R. simsii (Tables 2, 3 and 4). Mycorrhizal connections formed by ERM fungi have been found between ECM and ERM plant species such as Quercus ilex L. (Fagaceae, with ECM) and Erica arborea L.(Ericaceae, with ERM) (Bergero et al. 2000), and Pinus sylvestris L. (Pinaceae, with ECM) and Vaccinium vitis-idaea L. (Ericaceae, with ERM) (Grelet et al. 2009(Grelet et al. , 2010. The same fungal association such as on Rhizoscyphus ericae (putative ERM) and ECM have also been observed on roots of ericaceous plants (Wurzburger et al. 2011;Wu et al. 2012;Zhang et al. 2016). This suggests that in a natural ecosystem, mycorrhizal fungi could form common mycelia networks that connect plants (Johnson and Gilbert 2015). Mycorrhizal networks formed between ECM and ERM hosts may be essential for species coexistence and ecosystem functioning in subtropical forests in China, and perhaps for forests in general. However, other tree species, such as exotics, introduced into local forests, would carry their own specific fungi and other microbes with them, resulting in fungal communities becoming even more complex.
From this investigation, soil parameters and climatic factors did not have significant effects on fungal community composition on roots of R. simsii (Table 2). In some studies, on fungal communities along nitrogen or pH gradients, nitrogen levels and pH had significant effects (Zijlstra 2006;Hofland-Zijlstra and Berendse 2010;Fujimura and Egger 2012). In our study of root-associated fungi of V. carlesii, soil organic carbon, total nitrogen content, total phosphorus content and ammonium nitrogen content significantly affected fungal community composition (Zhang et al. 2016). However, in this study, fungal communities on roots of R. simsii were insensitive to edaphic or climatic factors and more responsive to plant community compositions, indicating that different ericaceous species depend on different determinants for their root-associated fungi. This is further evidence for the popularity of R. simsii compared to other ericaceous plants.

Disturbance resistance and temporal variations of root-associated fungi
A large number of OTUs was shared by all stand types (Fig. 2), and the number of their specific OTUs was limited ( Fig. 2 and Table 3), indicating that root-associated fungal community composition and structure on R. simsii hair roots was adaptive and relatively stable under anthropogenic disturbances. According to all organisms that directly contribute to a particular functional process in an ecosystem (Allison and Martiny 2008), these shared OTUs were divided into different functional groups. High disturbance tolerance (or low disturbance sensitivity) and community resilience may result from taxa either surviving the disturbance or being able to rapidly re-colonize the disturbed area (Lekberg et al. 2012). The high density of R. simsii plants in secondary forests and in plantations after at least 20 years (Appendix S1), and putative ERM fungal taxa showed relatively high disturbance resistance, also explains why both host populations and fungal communities were strongly competitive in secondary succession. The proportion of functional fungus groups changed along the restoration period as SEC II, SEC I and OGF. Hutton et al. (1997) also observed that ericoid mycorrhizal colonization on disturbed sites returned to levels comparable to undisturbed sites after 12 years regeneration from disturbance, and the recovery was much faster than soil microbial communities in oak forests 100 years after anthropogenic disturbances (Fichtner et al. 2014). The gap between microbial taxa to resistance (insensitivity to disturbance), and resilience (the rate of recovery after disturbance), was larger (Lekberg et al. 2012). Furthermore, compared with resistance, resilience investigation was less. Shade et al. (2012) noted that only 23% of 148 resilient-examined cases on composition and/or function of soil microbes could return to a pre-disturbance state. Under a greater disturbance, soil microbes could never fully recover even if the recovery time was longer. ERM fungi and the recovery of the hosts was short. It is not clear whether functional mycorrhizal fungi of R. simsii in disturbed forests (PLF, SEC I and SEC II) would be recruited continuously into fungi communities as in OGF. But the relative stability and resilience of root-associated fungal communities, especially these functional fungi, must be essential for survival and growth of ericaceous plants in disturbed forests.
Significant temporal variations of root-associated fungal communities of R. simsii were observed and inter-annual (between 2 years) variations were higher than intra-annual (within a year) variations (Table 1). Significant correlations were found between disturbance sensitivity and temporal variation of fungal taxa (Fig. 6), and functional fungal groups showed different responses to CV d and CV t , implying that fungal taxa with greater resistance to disturbance would be observed more constantly in different seasons and years. Putative and possible ERM fungal orders were relatively higher resistant to disturbances and had lower temporal variations; however, non-essential taxa were occasionally observed with lower disturbance resistance and higher temporal variations. The study site experienced an extreme drought and high temperatures during the summer of 2013, and the fungal communities of the spring of 2014 were similar to those of the 2012 spring (Fig. 4). The distinction of fungal communities in spring of the 3 years and the flourish of indicator species for the year 2014 could be partially due to the inter-annual shifts of climatic factors and partially to the succession of fungi. Cázares et al. (2005) proposed that non-mycorrhizae prevailed in the earliest stages with some DSE, ECM and ERM plants, but the three would become equally dominant plants over time. Comparing early and mature successional stages with different vegetation in subtropical forests of Australia and the United States, Dickie et al. (2013) found that even though both stages had similar mycorrhizal and non-mycorrhizal hosts, in the mature stages only orchid mycorrhiza emerged. In OGF, populations of R. simsii tended to decline even though their mycorrhizal composition and structure were the best. In order to accompany dominant trees in subtropical forests with succession, R. simsii has significant support from their mycorrhizal associations, and sustainable and stable ERM fungi is necessary.

Conclusions
This study on root-associated fungal communities of R. simsii investigated over 1000 operational taxonomic units and determined that were mainly in the Ascomycota and Basidiomycota phyla. The results show the effects of time and different forest types under anthropogenic disturbances on fungal communities, although the proportion of shared operational taxonomic units remained large, accounting for 44.9-79.4% of the total number. The changes in plant communities with different mycorrhizal associations significantly affected root-associated fungal communities of R. simsii, while soil and climatic factors had negligible effects.
The disturbance sensitivity and temporal variations of major fungal orders were significantly different among fungal functional groups, and putative and possible ericoid mycorrhizal fungi were highly resistant to disturbances and low temporal variations. There were significantly positive correlations between disturbance sensitivity and temporal variations of common fungal orders. The high disturbance resistance and temporal persistence of putative ericoid mycorrhizal fungal taxa may be essential for successful adaptation of R. simsii in disturbed subtropical forests of China.