The characters of root-derived fungi from Gentiana scabra Bunge and the relations with their habitats

Gentiana scabra Bunge (gentian) is one of the best-known traditional Chinese medicinal plants and its dry roots and rhizomes are usually used as medicine. To our knowledge, the main active metabolites of gentian, such as gentian polysaccharide and gentiopicroside, determine the quality and pharmaceutical effect of gentian. The root-derived fungi including endophytic fungi and rhizosphere fungi are key factors which affect the secondary metabolism and the utilization of soil nutrients of gentian. Therefore, it is necessary to investigate the characteristics and dynamics of root-derived fungi related to gentian, as well as their relationships with the quality of gentian. The population, diversity, and the dominant type of endophytic and rhizosphere fungi in gentian were determined by using ITS rRNA gene amplicon and sequencing methods. The potential influences of different habitats on fungi communities and the correlation between fungi communities, the main active metabolites of gentian and soil physicochemical properties were analyzed by statistical methods. The population and diversity of endophytic and rhizosphere fungi varied with both habitats and growth stages, showing a significant difference. Among them, the predominant genera of endophytic fungi were Lecidella, Leptoshaeria and unclassified_p_Ascomycota with the relative abundance of 35.1%, 9.0% and 7.3%, respectively, while the predominant genera of rhizosphere fungi were complicated. Compared to endophytic fungi, the diversity of rhizosphere fungi was obviously affected by soil physicochemical properties, including pH, water content, alkali hydrolysis nitrogen (AN), available phosphorus (AP), etc. Additionally, we found that the accumulation of gentiopicroside was positively correlated with endophytic fungi of Epicoccum and rhizosphere fungi of Mortierella, Solicoccozyma, Talaromyces and Trichoderma. The accumulation of gentian polysaccharide was negatively correlated with endophytic fungi of Lenzites, Mucor, Myrothecium and Saccharomycopsis and rhizosphere fungi, such as Botrytis, Cadophora, Cladophora, Didymela, Fusarium, etc. In this work, the differences on population, diversity, and the dominant type of endophytic fungi and rhizosphere fungi in gentian were revealed. The relationship between the main active metabolites of gentian and root-derived fungi, and the relationship between soil physicochemical properties and root-derived fungi were clarified, respectively. It is believed that our work will guide to explore strategies on improving the quality of gentian by regulating soil factors and root-derived microbial community structure with microorganisms.


Introduction
Gentiana scabra Bunge (gentian), as a traditional Chinese medicinal herb, exerts anti-inflammatory, antioxidant, anti-viral, and antimicrobial properties by main chemical components of gentian polysaccharide and gentiopicroside (Xu et al. 2020;Joksic et al. 2021;Zhao et al. 2019). As the demand increases, the source of wild gentian has been seriously destroyed owing to over-harvesting (Liu et al. 2020b) and the cultivated ones have now become the principal commodity. However, there are differences on the quality of gentian cultivated in different habitats, what we're interested in is the role that root-derived microbes play on the differences.
As a holobiont (Taule et al. 2021;Hassani et al. 2018), there were accompanied with various microbial communities and formed plant-microbial complexes during different developmental stages of gentian (Araujo et al. 2019). The root-derived fungi, including endophytic and rhizosphere fungi have a close symbiotic relationship with their hosts. The endophytic fungi produced antibiotics, enzymes, and other bioactive compounds to enable their host plants to survive in competitive habitats (Bogas et al. 2022). The rhizosphere fungi are often related to root exudates, and many organic compounds secreted by plant roots provide energy for rhizosphere fungi. Meanwhile, the rhizosphere fungi can provide nitrogen and phosphorus in return for their hosts, significantly promote the utilization efficiency of soil nutrients by plants, thus promote plant growth and even improving plant quality . Research on the endophytic fungi isolated from medicinal plants suggested that the fungal mycelia and fermentation broth contained the same metabolites as the host plants (Uzma et al. 2018;Rajamanikyam et al. 2017;Masi et al. 2018), which meant that there were relationships between plant metabolites and root-derived fungi. Hence, it is necessary to investigate the characteristics and dynamics of root-derived fungi, as well as their relationships with the quality of gentian.
To the best of our knowledge, many previous studies have shown that external factors, such as geographic location (Coller et al. 2019), soil physicochemical properties (Arévalo-Gardini et al. 2020), and climate (Cavicchioli et al. 2019), have a substantial impact on the diversity and community composition of endophytic and rhizosphere fungi Suryanarayanan and Uma Shaanker 2021). Soil physicochemical properties, including pH, the contents of available phosphorus, alkaline hydrolysis nitrogen and organic matter led to a direct difference in the composition of rhizosphere fungi community (Jing et al. 2020;Xie et al. 2020). Since rhizosphere fungi could shift and colonize in plant roots, they have served as the important sources of endophytic fungi for plants (Kerdraon et al. 2019;Zheng and Gong 2019), resulting in differences on the community composition of endophytic fungi (Zhang et al. 2018). Additionally, the community of endophytic fungi was regulated by the host plant, which was called as host selectivity (Schmidt et al. 2019), varied in different developmental stages of the plant (Wassermann et al. 2019;Jia et al. 2020). Compared with the seedling stage, plants in the mature stage clearly had a relatively stronger ability to regulate the community of endophytic fungal. For instance, Zhang et al. (2018) demonstrated that the community structure of rice endophytic fungi varied with the time that there was residence in the field, and finally, the community structure tended to be stable as the rice matured. Thus, we wondered that whether there were correlations between soil physicochemical properties and root-derived fungi during the cultivation of gentian, and if these factors could affect the accumulation of main active metabolites (gentian polysaccharide and gentiopicroside).
Based on the former facts, samples of roots and rhizosphere soils of gentian were collected from three Geo-Authentic Product habitats, including Liaoning 1 3 Vol.: (0123456789) (LN), Jilin (JL), and Yunnan (YN) province, China in May, July and September of 2019, respectively. Some specific objectives of this work were as follows: (1) to determine the population, diversity, and the dominant type of endophytic and rhizosphere fungi in gentian; (2) to analyze the potential influence of different habitats on the structure of endophytic and rhizosphere fungi; (3) to investigate the link between soil physiochemical properties and the structure, composition and diversity of endophytic and rhizosphere fungi communities; (4) to study the correlation between both fungi communities and the main active metabolites of polysaccharide and gentiopicroside in gentian. The results of our study will provide valuable information for better understanding of the relationship between root-derived fungi and the quality of gentian and will promote growth and improve the quality during planting by fungal resources.

Sample collection and processing
The samples of gentian roots and rhizosphere soil were collected from three Geo-Authentic Product habitats including Fushun in Liaoning Province (LN, 125°2′ E, 42°11′ N, 290 m above the sea level), Dunhua in Jilin Province (JL, 128°26′ E, 43°46′ N, 388 m above the sea level), Baoshan in Yunnan Province (YN, 99°4′ E, 25°22′ N, 1845 m above the sea level), respectively. The sampling time points were in May (the period that gentian roots boomed), July (the period that gentian vigorously grew) and September (the period that gentian matured) of 2019, respectively. Before sampling, four replicate plots (10 m × 10 m each) were established, and a total of five sampling areas of 1 m × 1 m at the four corners and center of each plot were set. Next, plant roots with attached soil were randomly sampled at above five plots from a depth of 0-20 cm to make one mixture sample (Schmidt et al. 2019;Zhang et al. 2021). Fresh samples were placed on ice and immediately transported to the laboratory where rhizosphere soil samples were collected from the soil adhering to the root crowns, and then passed through a 2 mm sieve and then homogenized . Moreover, the rhizosphere soil that was used for DNA extraction was freeze-dried and then stored at −80 °C, while the soil that was used to determine the physicochemical properties was naturally air-dried. The roots of gentian were divided into three parts: one part was stored at −20 °C for DNA extraction; one part was stored at −4 °C to determine root vitality, while the other part was naturally dried to determine the gentiopicroside and polysaccharides.
The roots of gentian used for DNA extraction were washed with running tap water and then washed with sterile distilled water to fully clear the mud and dried on sterile filter paper. The root tissues were immersed in 75% ethanol for 2 minutes followed by 5% (v/v) sodium hypochlorite for 5 min. After each treatment, the samples were washed five times using sterile water (Ling et al. 2020). Finally, sterile water from the last treatment was placed on Luria-Bertani agar (LB) medium and cultured at 28 °C for 72 hours. There was no observable colony growth, indicating that the root surface had been sterilized.

DNA extraction, PCR amplification and Illumina sequencing
The pretreated roots were ground to powder with liquid nitrogen, and the genomic DNA of the endophytic fungi was extracted using an M5 Plant Genomic DNA Kit (Mei5 Bio) following the manufacturer's instructions, Genomic DNA of the rhizosphere fungi was extracted from rhizosphere soil with a FastDNA Spin Kit for Soil (MP Biomedicals, Solon, OH, USA) according to the manufacturer's instructions. The quality and quantity of DNA were checked on a Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). The ITS region of fungal ribosomal DNA was amplified using the primers ITS1F and ITS2R (Cullen et al. 2018). Three technical replicates were established per sample, and each sample was amplified by PCR using a 20 μL reaction with 5 M of each primer. The PCR system included 10 × buffer (2 μL), 2.5 mM dNTPs (2 μL), 5 μM forward and reverse Primer (0.8 μL), 2.5 units/μL Taq Polymerase (0.2 μL), 20 mg/mL BSA (0.2 μL) and Template DNA (10 ng), adding ddH 2 O to 20 μL. The PCR reactions were performed by a GeneAmp PCR System 9700 (Thermo Fisher Scientific) that used the following cycling conditions: initial denaturation step at 95 °C for 3 min (1 cycle); 30 cycles at 95 °C for 30 s, 50 °C for 30 s, 72 °C for 45 s, a final extension step at 72 °C for 10 min (1 cycle). The amplification products were checked on a 2.0% agarose gel and purified using the AxyPrep-DNA Gel Recovery Kit (AXYGEN, Union City, CA, USA), following the manufacturer's instructions. Afterwards, a sequencing library was generated by the addition of an Illumina sequencing adaptor to the product using an Illumina TruSeq DNA Library Preparation Kit (San Diego, CA, USA), following the manufacturer's instructions. Finally, the library was sequenced on an Illumina MiSeq (PE300) platform at Shanghai Majorbio Bio-pharm Technology Co., Ltd.

Bioinformatics analysis
Raw sequences were analyzed using Quantitative Insights into Microbial Ecology (QIIME) version 1.9.1 to eliminate low-quality sequences (Caporaso et al. 2010), including the sequences that had an average quality score < 20, without valid primer sequences or barcode sequences, or that contained ambiguous bases or non-fungal DNA. Forward and reverse reads were subsequently joined using the FLASH version 1.2.11 with a minimum overlap of 10 bp (Magoc and Salzberg 2011). Of the filtered sequences, the potential chimeras were subsequently checked using the UCHIME commanded in MOTHUR version 1.30.2 (Edgar et al. 2011). The non-chimeric sequences were clustered into the operational taxonomic units (OTUs) at 97% similarity based on the UPARSE pipeline using USEARCH version 7.0 after discarding all the singletons (Edgar 2013). The OTUs were taxonomically classified using the RDP Classifier version 2.11 against the NCBI Nucleotide Sequence Database with a confidence threshold of 70% (Wang et al. 2007;Schoch et al. 2020). To compensate for different sequencing depths, samples were rarefied to an even depth according to the Minimum sample sequence, leaving a total of 5419 fungal OTUs for additional analysis.

Plant physiological and biochemical indicators and soil physicochemical properties
The content of gentiopicroside in gentian was determined using reversed phase high-performance liquid chromatography (RP-HPLC) with an external standard calibration curve method based on the Pharmacopoeia of the People's Republic of China. Gentian polysaccharides were analyzed by the anthrone-sulfuric acid method, using glucose as a standard (Yu et al. 2019). The root vitality of gentian was measured using triphenyltetrazolium chloride (TTC) (Yan et al. 2020). The plant height of gentian was measured using a ruler. The air-dried and sieved (<2 mm) soil samples were used to determine the soil physicochemical properties, including pH, available phosphorus (AP) and alkaline hydrolysis of nitrogen (AN). In particular, the soil pH was measured using a pH meter in a soil water suspension (1: 2.5 w/v) after shaking for 30 min (Lei et al. 2020). The soil moisture and AN were determined using gravimetric analysis and alkali diffusion methods, respectively Xu et al. 2021a). The AP was determined using the sodium hydrogen carbonate solution-Mo-Sb anti spectrophotometric method .

Microbial community analysis
The Shannon and Simpson indices were calculated for endophytic and rhizosphere fungi at the OTUs level as a measure of alpha diversity using MOTHUR version 1.30.2. The differences on the population and diversity of both fungi between the different regions and growth stages were tested based on a one-way analysis of variance (ANOVA) using GraphPad Prism version 7.0 (San Diego, CA, USA) Nguyen et al. 2020). Non-metric multidimensional scaling (NMDS) of Bray-Curtis dissimilarity matrices was used to identify differences between microbial communities sampled from LN, JL and YN in May, July and September 2019, respectively (Yao et al. 2019). The distance matrices were calculated using QIIME version 1.9.1, and the analysis and visualization were accomplished using the vegan package (Callahan et al. 2016). An analysis of similarity (ANOSIM) was conducted on the pairwise differences among groups based on 999 permutations to statistically support the results (Veach et al. 2019).
The composition of fungi community was analyzed at the phylum and genus level, respectively. Moreover, to observe the differences on the community structure of endophytic and rhizosphere fungi in different regions and growth stages more intuitively, the fungal relative abundance was analyzed at the Class level using the ggplot package . Only the top 10 Classes were analyzed, and the relative abundance <0.01 or unclassified fungi were clustered into 'other' group. Additionally, the Spearman's 1 3 Vol.: (0123456789) rank correlation analysis was conducted to ascertain the correlations among fungal communities and environmental factors, including plant physiological and biochemical indicators and soil physicochemical properties, using the pheatmap package at the Genus level. The top 30 genera were analyzed. Furthermore, Linear Discriminant Effect Size (LEfSe) analysis was used to detect biomarkers that existed in different regions and growth stages using LEfSe (http:// hutte nhower. sph. harva rd. edu/ galaxy). The non-parametric Kruskal-Wallis (KW) sum-rank and the Wilcoxon rank-sum tests were used to determine the difference in species abundance between different groups, then LDA linear discriminant analysis was applied to estimate the influence of these different genera on the difference between the test groups (Moreno-Arrones et al. 2019). The LDA threshold was set as 2 and 4 for the endophytic and rhizosphere fungi, respectively.

Co-occurrence network and statistical analysis
The co-occurrence networks for the LN, JL and YN samples were constructed to provide an insight into the structure and ecological interactions of microbial communities using the routine CoNet in Cytoscape 3.8.0 (Marasco et al. 2018). Spearman correlation coefficients were used to generate the co-occurrence networks. The edge-specific permutation and bootstrap score distributions were conducted to calculate the statistical significance of the co-occurrence and mutual exclusions. The data was then re-normalized for each permutation and provided a null distribution that captured the similarity introduced by the composition. The unstable edges were then filtered. In each of these networks, the nodes that represent genus and edges represent significant co-occurrence relationships. Other topological properties for the network include the average degree, clustering coefficient, average path length, and betweenness centrality. Software of Gephi version 1.9.2 was used to explore network properties and visualize networks (Xun et al. 2021).
An ANOVA was performed to assess the differences on plant physiological and biochemical indicators (gentiopicroside, polysaccharide, plant height and root vitality), soil physicochemical properties (pH, AP and AN) using SPSS 25.0 software (IBM, Inc., Armonk, NY, USA) among different groups.

Difference and dynamic analysis of physiological and biochemical indicators of gentian and soil physicochemical properties
Gentiopicroside and gentian polysaccharides, as the main active metabolites of gentian, varied significantly in different regions and growth stages. In detail, the content of gentiopicroside ranged from 74.0 to 97.1 mg‧g −1 in LN, 53.6 to 80.2 mg‧g −1 in JL, and 74.4 to 108.6 mg‧g −1 in YN (Table 1, P < 0.05). The average content of gentiopicrin in JL was relatively low, and the contents of gentiopicroside in LN and YN were significantly higher in July (P < 0.05). Gentian polysaccharides showed a decreasing trend between the habitats of LN and JL during the period of May to Sept, and the content ranged from 15.3 to 21.2 mg‧g −1 and 13.0 to 18.2 mg‧g −1 , respectively (Table 1). The content of gentian polysaccharide in YN was relatively stable and kept at a higher level, reached 29.3 mg‧g −1 (Table 1). Furthermore, the root vitality (RV) of gentian in the three habitats decreased from May to July and then increased from July to Sept (Table 1). The height of gentian showed an increasing tendency during the period of May to Sept in habitats of LN and JL, and the height of gentian in YN was significantly lower than that in LN and JL (Table 1, P < 0.05). The pH, available phosphorus (AP) and alkaline hydrolysis of nitrogen (AN) of rhizosphere soil of gentian varied with habitats and growth stages. The pH of rhizosphere soil was weakly acidic and relatively stable during May to Sept in all three habitats, and the soil pH in LN and JL was slightly higher than that in YN ( Table 2, P < 0.05). The content of AP in LN significantly increased with the growth of gentian (P < 0.05), while the variation trend of the content of AP in YN was opposite, and the content of AP in JL was stable during May to Sept (Table 2, P < 0.05). The content of AN in LN and JL showed the increasing tendency and reached the maximum of 141.46 mg‧g −1 and 162.24 mg‧g −1 in July, respectively (Table 2, P < 0.05). However, the AN content in YN showed a decreasing trend from May to September (Table 2, P < 0.05).
The community composition of the endophytic and rhizosphere fungi differed in different habitats and growth stages 5419 fungal isolates were obtained in our work, which belonged to 1 domain, 1 kingdom, 8 phyla, 32 classes, 104 orders, 280 families, 630 genera, 1238 species. The Non-metric multidimensional scaling (NMDS) ordination revealed that the community composition of endophytic and rhizosphere fungi was significantly different and varied with growth habitats and growth stages of gentian (Fig. S1). For endophytic fungi, the operational taxonomic units (OTUs) could not be separated completely into diverse clusters according to the regions, whereas the OTUs in May were separated from those in July and September (Fig. 1A). This result indicated that the composition of endophytic fungi community in gentian was more affected by the variation of growth stages. Unlike the community of endophytic fungi, the OTUs were separated into different clusters in terms of different habitats for rhizosphere fungi, showing that the composition of rhizosphere fungi community was more affected by growth habitats of gentian (Fig. 1B).
The predominant Phyla of endophytic fungi were Ascomycota and Basidiomycota across all the regions, the relative abundances were 75.2% and 3.5%, respectively (Fig. S2). At the genus level, Leptosphaeria and Trichoderma were the predominant endophytic fungal genera in May for all habitats, and the relative abundances were 25.5% and 25.9%, respectively. Moreover, Botrytis (20.3%) and Leptodontidium (4.4%) were specifically enriched in LN, while Lecidella (10.5%) and Cladosporium (8.6%) were enriched in YN. However, in July and September, the compositions of endophytic fungi communities were similar, and Lecidella and unclassified_p_Ascomycota were the predominant genera among all the habitats, the relative abundance of Lecidella were 58.5% in LN, 57.1% in JL and 63.2% in YN, respectively. The relative abundance of unclassified_p_Ascomycota were 10.2% in LN, 11.7% in JL and 10.8% in YN, respectively (Fig. 1C).
The predominant Phyla of rhizosphere fungi were Ascomycota, Basidiomycota and Mucoromycota in all three habitats, and the relative abundances were 47.5%, 17.4% and 6.2%, respectively (Fig. S2). During the period of May, Lecidella and unclassified_p_Ascomycota were the predominant genera with the relative abundances of 52.5% and 9.5% in LN, and Trichoderma (4.3%) and Solicoccozyma (4.3%) were the predominant genera in JL and YN. Moreover, Fusarium, unclassified_f_Mortierellaceae and Epicoccum were enriched in JL, and the relative abundances were 5.5%, 8.4% and 16.5%, respectively. Saitozyma was abundant in YN with the relative abundance of 5.5%. In July, Fusarium was enriched in LN and JL, and the relative abundances were 9.1% and 50.4%, respectively. Unclassified_f_Mortierellaceae (14.0%) was abundant in LN with the relative abundance of 14.0%. In addition, Trichoderma (8.0%), Humicola (6.2%) and unclassified_f_Pyronemataceae (11.8%) were concentrated in LN, Botrytis (5.4%) was abundant in JL, while Saitozyma (12.0%) was enriched in YN. The rhizosphere fungal composition of LN and JL in September had few changes compared with that in July. However, Russula became the dominant genus of YN, the relative abundance was 87.6% (Fig. 1D).
Diversity of endophytic and rhizosphere fungi varied in different habitats and dynamically changed with the variation of gentian growth stages Based on the Shannon and Simpson diversity indices, the alpha diversity of endophytic and rhizosphere fungi differed in growth stages of gentian. For endophytic fungi, the diversity was significantly lower in May, and then increased rapidly in July and Sept in LN and JL, while in YN, the alpha diversity of endophytic fungi remained at a higher level during the growth stages of gentian ( Fig. 2A, C). For rhizosphere fungi, the diversity in May was significantly higher than those in July and Sept in the three habitats, and the lowest diversity was observed in YN in Sept (Fig. 2B, D). Moreover, the alpha diversity of endophytic and rhizosphere fungi also varied with growth surroundings. The diversity of endophytic fungi in YN was significantly higher than those in LN and JL in May (P < 0.05), and significantly lower than those in LN and JL in July (P < 0.05) but the diversity was turned to be consistence with each other in Sept. The diversity of rhizosphere fungi in LN was significantly higher than those in YN and JL in May (P < 0.05), and significantly higher than that of JL in July (P < 0.05). However, there was no significant difference between LN and JL in Sept, but the diversity of rhizosphere fungi in LN and JL were significantly higher than that of YN (P < 0.05). Additionally, the differences between endophytic and rhizosphere fungi on the number of observed OTUs based the sobs diversity were compared, and the results indicated that the number of endophytic fungi significantly increased (P < 0.05), while the number of rhizosphere fungi decreased from May to Sept. The numbers of rhizosphere fungi were significantly higher than endophytic fungi in May (P < 0.05) and were significantly lower than endophytic fungi both in all three habitats in July and Sept in all three habitats (Table S1, P < 0.05).
Distinctive abundance of endophytic and rhizosphere fungi caused the differences between habitats and growth stages Linear Discriminant Effect Size (LEfSe) analysis was used to identify distinctive taxa at the genus level. As shown in Fig. 3, distinctive specific taxa were enriched in different regions owing to the influence of geographical factors. The endophytic fungi Gorgomyces and Candida were specifically enriched in LN, and Ilyonectria and Hymenoscyphus were abundant in JL, while g_unclassified_c_Tremellom ycetes was enriched in YN (Fig. 3A). For the rhizosphere fungi, Humicola and Saitozyma were specifically enriched in LN and YN, respectively, while Fusarium, Epicoccum, Botrytis and Cadophora were abundant in JL (Fig. 3B).
Furthermore, some taxa of endophytic and rhizosphere fungi also specifically enriched with the growth of gentian in different habitats. In LN, the endophytic fungi like Humicola and Leptodontidium were specifically enriched in May, Acremonium and Mucor were specifically enriched in July, and Verticillium was specifically enriched in Sept (Fig. 4A). The rhizosphere fungi like Humicola, Trichoderma, Metarhizium and Talraomyces were specifically enriched in Sept (Fig. 4B). In JL, the endophytic fungi like Leptosphaeria and Trichoderma were specifically enriched in May, and Nigrospora and Talaromyces were specifically enriched in Sept (Fig. 4A). The rhizosphere fungi like Solicoccozyma and Epicoccum were specifically enriched in May, Phoma and Fusarium were specifically enriched in July, and Botrytis and Leptodontidium were specifically enriched in Sept (Fig. 4B). In YN, the endophytic fungi like Russula, Cladosporium and Epicoccum were specifically enriched in July, and Lecidella was specifically enriched in Sept (Fig. 4A). The rhizosphere fungi Hirsutella was specifically enriched in May, and Saitozyma and Penicillium were specifically enriched in July (Fig. 4B).
Interactive effects on the community structures between endophytic and rhizosphere fungi To further investigate the changes in specific taxa during different growth stages, the relative abundance of endophytic and rhizosphere fungi in the same region were compared at the Class level. There was some trade-off relationship between endophytic and rhizosphere fungi (P < 0.001). From May to Sept in LN, the relative abundances of Tremellomycetes, Sordariomycetes, Leotiomycetes and Dothideomycetes decreased in the roots of gentian, and increased in the rhizosphere soil. The relative abundances of Lecanoromycetes and Ascomycota increased in roots, and the opposite trend was found in the rhizosphere (Fig. 5A, D). In JL, the relative abundances of endophytic fungi of Sordariomycetes, Leotiomycetes and Eurotiomycetes decreased, and Lecanoromycetes increased from May to Sept, whereas these fungi as rhizosphere fungi changed in opposite manners (Fig. 5B, E). Moreover, in YN, the relative abundance of endophytic fungi of Agaricomycetes decreased, but Lecanoromycetes increased. However, as the rhizosphere fungi, the dynamic changes of these fungi were in contrary to endophytic fungi (Fig. 5C, F).
We also found that some fungi existed in roots and rhizosphere soil of gentian showed same change tendency. For example, the relative abundances of Dothideomycetes in JL (Fig. 5B, E), Tremellomycetes, Sordariomycetes, Mucoromycetes, Leotiomycetes, Eurotiomycetes and Dothideomycetes in YN (Fig. 5C, F), synchronously decreased from May to September.

Characters of interspecific relationship based on the co-occurrence networks
To investigate the differences in community structures of endophytic and rhizosphere fungi and reveal the relationships between genera, the co-occurrence networks were constructed based on the Spearman correlation coefficient for each region (P < 0.01), respectively. The co-occurrence networks of rhizosphere fungi were more complex but had a lower degree of modularity than those of the endophytic fungi. Detailed analysis indicated that the network topological coefficients including modularity and graph density, as well as the nodes and links of the network were obviously lower in endophytic fungi than those in rhizosphere fungi. However, the average clustering coefficient and average path length followed the opposite trend (Table 3).
The co-occurrence networks were also different in the growth regions of gentian. The co-occurrence network of endophytic fungi in JL was more complex than those in LN and YN. Furthermore, the degree, eccentricity, closeness centrality and betweenness centrality of the top 10 nodes in co-occurrence networks of three regions apparently differed. It was notable that the degree of unclassified_k_Fungi, Lecidella and unclassified_p_Ascomycota ranked as the top three in the co-occurrence networks of all regions, which indicated that they were the keystone taxa for the endophytic fungi of Gentian (Tables S2, S3, S4, Fig. S3A, B, C). In detail, the three key endophytic fungi taxa in LN were positively correlated to Acremonium, Verticillium, Lenzites, Monosporascus, Epicoccum and Penicillium, and negatively related to Leptospharria, Wallemia, Leptodontidium, Trichoderma and Humicola (Fig. 6A). The key taxa in JL were positively correlated with Sarocladium, Tylopilus, Lenzites, and Verticillium and were negatively related to Penicillium, Trichoderma, and Botrytis (Fig. 6B). Furthermore, the main taxa in YN were positively related to Mucor, Nigrospora, and Monosporascus, and were negatively correlated with Cladosporium (Fig. 6C).
The geographical location showed a more significant influence on the co-occurrence networks of rhizosphere fungi. The key taxa for rhizosphere fungi between the genera of each region differed. For example, the degree of Penicillium, Chaetomium and Metarhizium ranked as the top three in the rhizosphere fungal network of LN, and Fusarium, unclassified_f_Herpotric hiellaceae and unclassified_o_Helotiales were the top three genera in JL, while unclassified_o_Helotiales, unclassified_p_Ascomycota and Penicillium ranked as the top three in YN (Tables S5, S6, S7).
As for the relationship between genera in rhizosphere fungal networks, Penicillium was positively correlated with unclassified_f_Mortierellaceae, unc lassified_f_Lasiosphaeriaceae and Chaetomium in the LN, while Cryptococcus was negatively correlated with the latter three. In addition, Chaetomium was negatively correlated with Mortierella and unclassified_o_Hypocreales, and Metarhizium was positively related to Didymella, Humicola, and Solicoccozyma. Moreover, Trichoderma and Epicoccum were a negative influence for each other (Fig. S3D).  In the co-occurrence network of JL, Fusarium was positively related to Pseudeurotium, Mortierella and unclassified_f_Mortierellaceae, and negatively correlated with Humicola, Solicoccozyma and Talaromyces. In addition, both unclassified_f_Herpotrichiella ceae and unclassified_f_Helotiales were found to be negatively correlated with Epicoccum (Fig. S3E). In the YN network, unclassified_o_Helotiales positively correlated with Metapochonia and Trichoderma, and negatively related to Saitozyma, and Aspergillus. In addition, Solicoccozyma also had a positive relationship with Trichoderma (Fig. S3F).
Some specific fungal taxa correlated with physiology and biochemistry indicators of gentian and soil physicochemical properties Spearman correlation heatmaps were used to investigate the correlation between gentian quality factors, soil physicochemical properties and the endophytic and rhizosphere fungi of gentian (P < 0.05).
In detail, endophytes of Epicoccum, unclassified_f _Hyaloscyphaceae and unclassified_k_Fungi had a positive correlation with the content of gentiopicroside in gentian, while Botrytis was negatively related to gentiopicroside (P < 0.05). The content of gentian polysaccharide had a significantly positive correlation with Hyaloscyphaceae (P < 0.001), while significantly negatively related to Mucor (P < 0.01), Lenzites, Myrothecium and Saccharomycopsis (P < 0.05). In addition, most of the endophytic fungi, such as Epicoccum, Lenzites, Mucor, and Saccharomycopsis were significantly positively correlated with the height of gentian (P < 0.001). However, the correlation between endophytes and root vitality of gentian was not significant. As for soil physicochemical properties, the significant correlation between endophytes and soil pH was negative, and there were positive correlations between most of endophytes and available phosphorus (AP) of soil (Fig. 7A). Compared with the endophytic fungi, the rhizosphere fungi were more closely correlated with gentian physiological and biochemical indicators and soil physicochemical properties. Particularly, Mortierella, Solicoccozyma (P < 0.05), Talaromyces and Trichoderma (P < 0.01) were found to be positively related to the contents of gentiopicroside, while most rhizosphere fungi like Botrytis, Cadophora, Cladophora, Didymella, and Fusarium were the most significantly negatively correlated with the contents of gentian polysaccharide (P < 0.001). Additionally, there were significantly positive correlations between most of rhizosphere fungi, such as Botrytis, Cadophora, Fusarium, Humicola, and Leptodontidium, and the height of gentian (P < 0.05), as well as the pH and AP of the soil. The most of significant correlation between rhizosphere fungi and soil alkaline hydrolysis of nitrogen (AN) was negative (Fig. 7B).

Discussion
The community composition of endophytic fungi had spatial and temporal difference, and mainly affected by the selectivity of the host plant. Our research showed that the abundances of Agaricomycetes, Ascomycota and Lecanoromycota significantly increased in the roots of gentian, but decreased in the rhizosphere from May to September, which indicated that these rhizosphere fungi could shift from the rhizosphere soil to roots of gentian, resulting in the differences on the composition of endophytic fungi. This study also showed that the community structures of endophytic fungi of gentian in the three habitats were obviously different in May and July, but tended to be consistent in September, which indicated that the community of endophytic fungi was more affected by soil physicochemical properties and climate change, which was also a result that endophytic fungi could adjust themselves to help the host plants for adapting different surroundings (Suryanarayanan and Uma Shaanker 2021). Root vitality referred to the absorption and synthesis ability of roots, which directly affected the growth and nutritional status of plants. From May to September, the height and root vitality of gentian were increased, showing that the host selectivity increased with the growth of gentian, which was the main factor affecting the composition of endophytic fungi.
However, the community structure of rhizosphere fungi was significantly influenced by soil physicochemical properties and slightly regulated by the plant. In our work, most of the rhizosphere fungi with high relatively abundances were inversely related to alkaline hydrolysis of nitrogen (AN) to a large extend, which was the reason why the abundance of rhizosphere fungi in YN was lower than that in LN and JL. As shown in the result, most of the rhizosphere fungi with high relatively abundances were positively correlated with available phosphorus (AP), that means the plant root recruited the rhizosphere fungi which were AP-depending during the growth of gentian Sasse et al. 2018;Phour et al. 2020). Furthermore, the composition and abundance of rhizosphere fungi community were significantly different in May and July in LN and JL and the rhizosphere fungi in YN remained difference with LN and JL, these two results illustrated that the complex and changeable soil physicochemical properties were the dominant regulator for the rhizosphere fungi (Veach et al. 2019;Zhang et al. 2018). While the composition and abundance of rhizosphere fungi community slightly regulated by the plant for its tended to be consistent from July to September in LN and JL.
There were great differences on the co-occurrence networks of endophytic and rhizosphere fungi. For the Nodes, Links and Average Clustering Coefficient, the networks of rhizosphere fungi were relatively more than that of endophytic fungi, showing that the networks of rhizosphere fungi were more complex. However, the average degree and graph density of the networks of endophytic fungi were relatively higher than that of rhizosphere fungi, and the modularity was relatively lower than that of rhizosphere fungi, indicating that the scale of the networks of endophytic fungi were relatively small, but the connections between species in the networks of endophytic fungi were more closely. Additionally, the compositions of co-occurrence networks of endophytic and rhizosphere fungi were different from all three habitats, showing that the environment factors play a leading role on the construction of co-occurrence network for both endophytic and rhizosphere fungi. On one hand, the environment factors such as climate and soil physicochemical properties, significantly affected the species composition of rhizosphere fungi (Arévalo-Gardini et al. 2020;Xie et al. 2020;Park et al. 2020;Deltedesco et al. 2020;Dickie et al. 2013), and then affected the co-occurrence network. On the other hand, we knew that the new roots of gentian grew in May, and gentian plant entered a vigorous growth stage in July, then matured in September. During the early growth stage of gentian in May, the community of endophytic fungi were more affected by the environment factors (Zhang et al. 2018;van der Heijden and Hartmann, 2016;Nelson 2018), and the cooccurrence network differed with the variation of fungal species. With the maturity of gentian (from July to September), the host selective effects on endophytic fungi were enhanced and played a leading role in recruiting fungi (Grudzinska-Sterno et al. 2016;Miguel et al. 2019;Cavaglieri et al. 2009). Although the community of endophytic fungi in different habitats tended to be consistent, the co-occurrence network was different. In short, the environment factors regulated the interrelationship among fungi and were more important to the co-occurrence networks than the host selectivity of gentian.
A substantial amount of previous research showed that Epicoccum was widely found in the internal environment of the plant and played an important role in promoting the secretion of some compounds (Silva-Valderrama et al. 2021;Bian et al. 2021;Taguiam et al. 2021;Lee et al. 2020;Ogórek et al. 2020) and inhibiting the plant pathogens. Thus, Epicoccum nigrum was newly discovered properties as a biocontrol fungus. Moreover, Botrytis, as a common plant pathogen, was highly devastating for plant health (Petrasch et al. 2019;Veloso and van Kan 2018;Wang et al. 2021). In this work, the co-occurrence network suggested that the relationship between these two genera was found to be mutual inhibitory, and a correlation analysis showed that gentiopicroside was positively correlated with the endophyte Epicoccum, and negatively correlated with Botrytis. The cooccurrence network also suggested that Epicoccum was negatively related to Leptodontidium, while Botrytis was positively related to it. Therefore, we thought that the effect of gentian was finally shown as the result of antagonism of Epicoccum, Botrytis and Leptodontidium, and Epicoccum was likely to play a role through the network interaction of the these three or more nodes. Furthermore, Mortierella, as the plant growth-promoting fungi from soil, was widely applied in the agricultural inoculation (Dourou et al. 2017;Zhao et al. 2020). Solicoccozyma, Talaromyces, and Trichoderma are regarded as beneficial soil fungi (Zhao et al. 2021a, b;Yu et al. 2021;Singh and Balodi 2021), which can inhibit Fusarium (del Carmen Orozco-Mosqueda and Santoyo 2021), reduce disease, and then promote the Vol.: (0123456789) growth of gentian, thus, increasing the accumulation of gentiopicroside. Mucor and Myrothecium were the common plant pathogens and decreased the yield of commercial crops (Cossus et al. 2021;Choinska et al. 2020;Lakornwong et al. 2019), and Lenzites is regarded as a wood-rotting fungus and participate in the biodegradation of plant biomass (Knezevic et al. 2013;Wu et al. 2012;Xu et al. 2021b). As we know, the primary metabolites were connected to the plant stress response, and the plant could resist the biotic stress by re-shaping the primary metabolism (Tarkowski et al. 2020;Bezrutczyk et al. 2018). Our results suggested that the accumulation of gentian polysaccharide was negatively correlated with Lenzites, Mucor, Myrothecium and Botrytis and positively correlated with Coniochaeta. We inferred that Mucor, Myronthecium and Botrytis might influence the growth of plant to reduce the content of polysaccharides, while Lenzites might degrade the polysaccharides, thereby leading to a low content of gentian polysaccharide. It was readily apparent from the above analysis that there was almost no direct correlation between the differential fungi, marker fungi or network core fungi from different habitats and quality of gentian. Therefore, it can be hypothesized that a few species may play a direct role to the accumulation of plant metabolism products, while more key fungi rely on a microbial interaction network to exert function.

Conclusions
In this study, we characterized the population and diversity of root-derived fungi community in gentian and investigated their correlation with plant physiological and biochemical indicators and soil physicochemical properties, respectively, for the first time. Our results showed that the fungal diversity associated with gentian was considerable, showing a high diversity. The composition and structure of endophytic and rhizosphere fungi communities were affected by different habitats and growth stages. The dominant endophytic and rhizosphere fungi orders differed dramatically. Particularly, the compositions and structure of endophytic fungi community were mostly affected by the selectivity of the host plant. However, rhizosphere fungi were significantly influenced by soil physicochemical properties and slightly regulated by the plant. Additionally, we found that the accumulation of gentian metabolites, such as gentiopicroside and polysaccharide were affected by some specific taxa but more affected by the interaction of microbial communities. High-throughput sequencing methods and the correlation analysis developed in this work can be used for more systematically and comprehensively estimating the fungal communities of traditional Chinese medicines in future studies.