Analysis of eld-grown Panax ginseng rhizosphere microbiome provides clues for preventing red skin root syndrome

Ginseng red skin root syndrome (GRS) is one of the most common ginseng diseases. It leads to a severe decline in ginseng quality and seriously affects the ginseng industry in China. However, as a root disease, the characteristics of GRS rhizosphere microbiome are still unclear. The amplicon sequencing technology, combined with bioinformatics analysis, was used to explore the relationship between soil ecological environment and GRS. There were signicant differences in the diversity and richness of soil microorganisms between the rhizosphere with different degrees of disease, especially between healthy ginseng (HG) and heavily diseased groups. We also found that bacterial communities underwent multiple changes between complex stability and simple instability in different ginseng rhizospheres through the established interaction networks. The GRS group also had more competition with each other and ecological niche separation than the HG group. The fungal community's stability decreased signicantly in the early stages of the disease, followed by the formation of a stable and complex fungal community. The GRS groups signicantly increased interspecies cooperation and ecological niche overlap in the fungal network than the HG group. Microbes closely related to potential pathogenic fungi were also identied according to the interaction network, which provided clues for looking for biological control agents. Finally, the Distance-based redundancy analysis (dbRDA) results indicated that total P (TP), available K (AK), available P (AP), catalase (CAT), invertase (INV) are the key factors that inuence the microbial communities. This study collectively analyzed the changing characteristics in ginseng rhizosphere and provided the basis for soil improvement and biological control of eld-grown ginseng.


Introduction
Ginseng (Panax ginseng Mayer), mainly distributed in northeast China and South Korea, is known as the "King of Herbs" because of its substantial medicinal value. Because of the continuous cropping effect of ginseng, the method of cutting and planting ginseng was often used in China before. For a long time, we have lost a lot of forest resources for the ginseng industry. And now, this type of cultivation is no longer allowed. Due to the limitations of planting areas and planting patterns, it is more important to improve the quality and yield of ginseng for the sustainable development of the ginseng industry. Many diseases that occur during the growth process of ginseng are the biggest obstacles to our obtaining high yield and quality ginseng.
GRS has been reported as a disease caused by Al, Fe, and microbial stress in ginseng cultivation. The condition is characterized by rust-stained reddish-brown areas on the epidermis of ginseng roots. Areas may grow more signi cant as the planting years increase. In severe cases, areas can occupy more than 80% of the epidermis of ginseng roots. Unlike root rot disease, severe symptoms of red skin roots are also rarely associated with root rot. In the stem and leaf part of ginseng, however, there is no prominent characteristic. Previous studies have found that the application of chitosan induces red skin roots, and a variety of phenolic compounds and elements accumulate in red skin tissues [1]. In particular, the accumulation of Al and Fe may promote the accumulation of phenolic compounds and the activation of enzymes related to their oxidation. The activity of various antioxidant substances and antioxidant enzymes in red skin tissues is signi cantly increased, which can prevent phenolic compounds from being oxidized [2]. The microorganisms produce pectinase, cellulase and ligninase that damage the cell walls of ginseng roots, causing red skin symptoms that are exacerbated by the application of Fe 3+ [3]. The presence of pathogenic fungi has been reported, and non-biological factors in the red skin formation process seem to be closely related [4][5][6][7][8][9].
There are close interactions between plants, soil and microorganisms [10]. Healthy plant growth is dependent on the ability of microorganisms to promote soil material circulation and nutrient conversion [11,12]. The plants' root diseases are closely related to soil microbial community and physical and chemical properties, and there are a variety of potential pathogens in soil [13]. Therefore, to study a plant root disease, it is necessary to understand the microbial community and physical and chemical properties of rhizosphere.
The results of previous studies showed that the rhizosphere of diseased ginseng was separated from that of healthy ginseng, and there are considerable differences in the dominant bacteria and fungi genera in the two rhizospheres. Some bene cial microorganisms were signi cantly reduced in the rhizosphere of diseased ginseng, and the potential pathogen Ilyonectria was also identi ed [6]. And the diversity of bacteria and fungi in the rhizosphere of diseased ginseng was decreased [14]. There have also been studies exploring rhizospheres of plants with varying degrees of red skin and have found a strong correlation between metal element accumulation and red skin symptoms. According to the experimental results, Ilyonectria may not be necessary for red skin root [15].
In this study, we classi ed the ginseng rhizosphere into ve degrees based on the severity of the disease.
The 16S rRNA gene and ITS region were applied to explore the ve degrees rhizosphere microbial community composition on the Illumina PE250 platform. Nutrient composition and enzyme activity in rhizosphere were also examined, and their relationship with microbial changes was explored. This study elucidated the evolution rule of ginseng rhizosphere with the aggravation of GRS and provides the basis for preventing and treating this disease in eld production.

Results richness and diversity of microbial communities
In the bacterial analysis based on 16S rRNA, the OTUs number in the GRS2 group was signi cantly higher than that in the HG group, GRS1 group, and GRS3 group (P<0.05). The Observed species of the GRS2 group is also substantially higher than HG and GRS3 (P<0.05). Besides, The Shannon index and Simpson index of the GRS2 group are both signi cantly higher than that of other groups, and the Chao1 index of GRS2 group is also signi cantly higher than that of the HG group and GRS3 group (P<0.05) ( Table 1).
Based on the ITS rRNA analysis of fungi, there was no signi cant difference in OTUs numbers among the ve groups. The Observed species and Chao1 index of the GRS2 group were signi cantly higher than that of the GRS3 group (P<0.05). Also, the Shannon indices of the HG, GRS1, and GRS2 groups were considerably higher than those of the GRS3 and GRS4 groups (P<0.05). The Simpson index was signi cantly higher in the HG group than in the GRS3 and GRS4 groups (P<0.05) ( Table 1).
Differences and composition of microbial communities NMDS analysis based on bacterial OTUs level showed no signi cant separation of GRS2, GRS3, and GRS4 groups in all ve groups of samples, while HG and GRS1 groups were separated from other groups (Fig. 1A). Cluster analysis showed that samples from the HG and GRS3 groups were signi cantly clustered to form independent clusters. The samples were then analysed for bacterial abundance, belonging to 46 phyla, 246 families, and 593 genera. In the analysis of phylum level, Proteobacteria is the most dominant phylum (average 50.2%), followed by Acidobacteria (average 17.6%), Bacteroidetes (average 7.1%) Gemmatimonadetes (average 6.0%), and Actinobacteria (average 5.9%) ( Fig. 2A).
In terms of the fungi community, by NMDS analysis, we found the signi cant separation of HG and GRS2 groups from GRS3 and GRS4 groups (Fig. 1B). UPGMA cluster analysis showed that HG, GRS3, and GRS4 samples were signi cantly clustered and formed independent clusters. Furthermore, based on the ITS gene sequences annotation, we identi ed 12 phyla, 197 families, and 323 genera in the fungal community. At the phylum level, Ascomycota is the most dominant (average 53.8%), followed by Mortierellomycota (average 16.0%), Basidiomycota (average 7.6%), and Glomeromycota (average 1.3%) (Fig. 2B).

Biomarkers and different levels of taxa change
Through LEfSe analysis, we found six bacterial taxa in the HG group with LDA scores greater than 4, including Acidobacteria and Actinobacteria. In the GRS1 group, LDA scores greater than 4 were Alphaproteobacteria and Rhizobiales. Phyla level with LDA scores greater than 4 in GRS3 was Proteobacteria and Bacteroidetes. Potential biomarkers of the GRS4 group were selected as Rizobiaceae and Achromobacter ( Fig. 3A and Fig. 3B). Then, we conducted MetaStat analysis (from the phylum level to the species level) to understand further the bacterial community changes with the disease severity of ginseng (Fig. S2). The results showed that the bacterial abundance in rhizosphere at different levels varied signi cantly with varying degrees of disease. Notably, at the phylum level, compared with the HG group, the abundance of Proteobacteria in the four red skin groups was signi cantly increased, and Bacteroidetes has a greater abundance in GRS3 and GRS4 groups. Besides, the abundance of Actinobacteria, Chloro exi, and Gemmatimonadetes were lower in the severe groups (GRS3 and GRS4) than in the healthy and mild groups (GRS1 and GRS2). Interestingly, the genera with a relatively high abundance were concentrated in HG, GRS3, and GRS4 group, and there were also signi cant differences between the two severe groups (GRS3 and GRS4).
LEfSe analysis results for fungi showed that, in the HG group, Fusariun, Chaetomium, and Glomeromycota had LDA scores over 4. In the GRS1 and GRS2 group, Mortierellomycota and Humicola were screened out, respectively. Besides, there were more taxa with LDA scores above 4 in the GRS3 group, and they were concentrated in Ascomycota, including Nectriaceae, Herpotrichiellaceae, and Ilyonectria ( Fig. 3C and Fig. 3D). For the MetaStat analysis at the phylum level, the results showed that the relative abundance of Mortierellomycota, Glomeromycota, and Mucoromycota decreased with increasing ginseng red skin area ( Fig S3). Besides, among several fungi with greater relative abundance at the genus level, the GRS3 and GRS4 groups had a much greater abundance of Ilyonectria than the other three groups. In contrast, the abundance of Trichosporon and Paraglomus decreased with the increase of red skin area. The Humicola abundance of GRS2 and GRS3 groups was higher than the other three groups (Fig. S3).

Network analysis of microbial communities
Based on spearman correlation, microbial networks in ve groups from HG to GRS4 were established. Through microbial networks, we can nd that the keystone phyla of bacteria and fungi are Proteobacteria and Ascomycota, respectively ( Fig. 4 and Fig. 5). In the bacterial network, the number of nodes in the ve groups did not change signi cantly. However, the number of total links in the ve groups was not stable, reaching the maximum in the GRS4 group (3600). All four red skin groups increased their proportion of negative links compared to the HG group (Table S3). Among other topological properties, the average path length of the GRS4 group was lower than that of other groups. The modularity of red skin groups was more down than the HG group, especially the GRS4 group. Also, the GRS4 group clustering coe cient and average degree were higher than the other groups (Table S3).
For fungi, the GRS4 group had signi cantly more links than the other groups, and the four red skin groups had more proportion of negative links than the HG group (Table S3). Also, from HG to the GRS4 group, the modularity value decreases gradually. The average path length in the GRS4 group also decreased signi cantly. In contrast, the clustering coe cient and average degree increased for the GRS4 group (Table S3).
In addition, through the interaction network, we obtained the microbe that may be strongly correlated with GRS-related fungus (Ilyonectria, Cylindrocarpon, and Fusarium) mentioned in the previous studies (Table  S4) [8] [5] [16,17]. The results showed that the network relationship of the three fungi was different in the rhizosphere of ginseng with different degree of disease, and the interaction relationship may be more in GRS4 group.

Soil properties and their relationship with microbial communities
After determining soil physical and chemical properties, we found that the pH of the two severe groups (GRS3 and GRS4) decreased signi cantly compared with the other three groups (Fig. S4). Signi cant changes in some nutrients were also observed among the sample groups (Fig. S5). The levels of AP and TP were lower in the two severe groups than in the HG group and mild groups (GRS1 and GRS2), and the AK content of GRS1, GRS3, and GRS4 groups was signi cantly lower than that of HG and GRS2 groups ( The dbRDA results for fungi and bacteria demonstrate the relationship between microorganisms and environmental factors. In Fig. 6A, axis 1, and axis 2 accounts for 75.81% and 11.51% of the total variance. TP, AK, CAT, AP, and pH were signi cantly correlated with bacterial phylum composition. In the fungi plot, axis 1 and axis 2 account for 71.55% and 14.71% respectively, and INV, CAT, and AP are most closely related to the community composition of fungi (Fig. 6B).

Discussion
GRS is one of the most concerned ginseng diseases at present. In recent years, GRS has appeared in large areas in some ginseng producing areas in China, seriously damaging the ginseng industry. As a root disease, a thorough understanding of its rhizosphere characteristics is essential [18]. In this study, to provide basis for soil improvement and disease control, we made a more detailed division of the disease (from HG to GRS4) to comprehensively reveal the pattern of rhizosphere changes in this ginseng disease's background.

Changes of microbial diversity and richness in rhizosphere
Through the statistics of the Shannon index and Simpson index of each group, we found that with the increase of ginseng disease degree, the change of bacterial community diversity in rhizosphere was rstly increased and then decreased. Observed species and Chao1 index indicated that the richness of the rhizosphere bacterial community uctuated in these ve degrees, and the highest richness was found in the GRS2 group (Table 1). This result seems to imply that GRS2 (red skin area 25 to 50%) is a grade of concern because of the turning point of changes in the ginseng rhizosphere's bacterial diversity.
For fungi, the fungal diversity in the severe groups (GRS3 and GRS4) was signi cantly lower than that in the HG group, while the richness did not change signi cantly (Table 1). Overall, the fungal diversity in the GRS's rhizosphere was affected, especially in severe groups.

Changes of microbial community composition in rhizosphere
To better re ect the nonlinear structure of the ecological data and the inter-and intra-group differences in the samples, we performed NMDS analysis based on OTUs level (Fig. 1). The bacteria in HG and GRS1 groups were separated from the other three groups, which means that with the development of the disease, the bacteria in rhizosphere had a continuous change process and tended to be stable at the later stages (from GRS2 to GRS4) of development (Fig. 1A). For the NMDS analysis of fungi, we still found that the HG group was signi cantly separated from the GRS3 group and the GRS4 group. Interestingly, the groups were not separated by the degree of disease development (Fig. 1B). The possible reason is that the fungal community develop uctuating rather than gradual as the disease progresses. Further, by clustering analysis, we also found that the bacterial HG group was individually clustered; meanwhile, the other GRS groups were not incredibly signi cantly clustered, con rming the bacterial community's instability in the rhizosphere of the diseased ginseng ( Fig. 2A). In terms of fungi, in addition to the HG group, the GRS3 group and GRS4 group also had obvious clustering, which seemed to indicate that the fungal microbial community in the severe groups again formed the stable state different from the HG group. (Fig. 2B). Based on the above, the bacterial community of ginseng rhizosphere appears to be more sensitive than the fungal community at the beginning of the disease. On the contrary, the fungal community has more stability compared to the bacterial community.

Characteristics of microbial community taxa variation in rhizosphere
Proteobacteria is the most dominant phylum among bacteria, and its abundance increases with the degree of disease ( Fig. 2A and Fig. S2). The increase of Bacteroidetes abundance in GRS3 and GRS4 groups (Fig. S2) may be due to its good degradation ability, which is more suitable for survival in the ercely competitive soil [19]. The abundance of Actinobacteria, Chloro exi, and Gemmatimonadetes was reduced in severe groups (GRS3 and GRS4) (Fig. S2). Actinobacteria contribute to OM's the decomposition [20], but the ecological functions of Chloro exi and Gemmatimonadetes were unclear. In the LEfSe analysis, 24 biomarkers were obtained and distributed in four groups (Fig. 3A). We also identi ed Acidobacteria as a potential biomarker for the HG group, which has been reported to degrade plant-derived OM speci cally and is more abundant in soils of enriched plants [21]. In the GRS1 group, Rhizobiales, a clade of Alphaproteobacteria, were screened out and reported to x nitrogen or as a pathogen [22]. Most biomarkers in the GRS3 group belongs to Proteobacteria, and the other belong to Bacteroidetes. Although described as phytopathogenic, the GRS4 group biomarker, Rhizobiaceae, is not known whether it is related to rust roots [23].
The fungal community's dominant phyla are Ascomycota, Mortierellomycota, and Basidiomycota (Fig.  2B). In severe groups, the dominance of Mortierellomycota decreased (Fig. S3). Besides, there are Glomeromycota and Mucoromycota that decrease in abundance as the disease progresses. In the present study, we still found changes in the abundance of Ilyonectria in different groups at the genus level. We also found the fungi genus (Trichosporon and Paraglomus), whose abundance decreases with disease development (Fig. S3). We also obtained speci c biomarkers for each group through LEfSe analysis, such as Glomeromycota, Fusarium, and Chaetomium of the HG group (Fig. 3C). In the GRS3 group, biomarkers all belong to Ascomycota, the soil ecosystem's main fungal decomposer (Fig. 3D) [24,25]. Among them, except Ilyonectria, the other biomarkers were not associated with GRS, although they have been reported to be associated with plant diseases [26,27].
In general, we used statistical methods to look for changes in the composition of rhizosphere microbial communities with different disease degrees. Although taxa were found to be correlated (positively or negatively) with disease severity, their ecological function was unknown.

Changes of the microbial interactions in the networks
Co-occurrence networks can be constructed to analyze the interaction and co-existence patterns among different microorganisms, which is crucial to our further understanding of the changes in rhizosphere with GRS development [28]. Fig. 4 shows that the phylum with the highest connectivity in the ve groups' bacterial community is Proteobacteria, which is consistent with previous study [29]. Therefore, with the disease's development, the keystone phylum has not changed in ginseng rhizosphere. However, it should be noted that Proteobacteria accounts for a larger proportion in GRS3 and GRS4 group. Since Proteobacteria could exploit labile carbon sources and produce extracellular polysaccharides to bond sand particles, this may mean that the rhizosphere of the severe groups is more solid, while has a stronger nutrient metabolism [30,31]. Further, from the obtained topology parameters, we found that links and average degree were uctuating (from HG to GRS4). The average path length of the GRS4 group was signi cantly lower than the other four groups (Table S3). Thus, with the disease's development, the bacterial community in the ginseng rhizosphere underwent multiple changes from complex and stable to simple and unstable. In addition, by counting the ratio of positive to negative links, we can also predict that there may be more interspeci c competition and ecological niche separation in the GRS groups than in the HG group (Table   S3) [32].
For the fungal networks, we can nd that the keystone phylum of rhizosphere was Ascomycota, unchanged at all disease development stages. It can be seen that in ginseng rhizosphere, Ascomycota was essential for resisting the harsh environment and maintaining system stability [33]. Through the statistics of topology parameters, we found that from HG to GRS4, the number of links, clustering coe cient, and the average degree were the lowest in the GRS1 group, and then gradually increased (Table S3). Therefore, in the early stage (GRS1) of ginseng disease, fungal communities' interaction and stability decreased. Interestingly, after this, the soil's fungi gradually formed a complex and stable community (GRS4). From the change in the ratio of positively and negatively links, it appears that the competitive relationships and the taxa in the same ecological niche were severely affected in the GRS1 group (Table S3). Furthermore, the negative correlation in GRS groups was lower than that in the HG group. In previous studies, GRS was closely related to fungi [8] [34]. Here, based on the interaction network, we identi ed microorganisms that interacted closely with the reported potentially pathogenic fungi. Of these, Simplicillium, which is extremely negatively correlated with Ilyonectria, is reported as a biological control agent [35]. In addition, although these microorganisms, such as Musicillium and Arachnotheca, have not been fully studied, the establishment of the networks also provides a basis for the future search for microorganisms that have antagonistic interactions with pathogens and, thus, for biocontrol development.

Changes of physicochemical properties and their relationship to microbial communities
Changes in the soil's physicochemical properties directly determine the structural composition of the microbial community [36,37]. In this study, pH, nutrient composition, and enzyme activity of rhizosphere were measured from HG to GRS4 group. From the GRS2 group, the pH of rhizosphere decreased signi cantly (Fig. S4). Compared with other groups, AP and TP in rhizosphere of the severe groups (GRS3 and GRS4) were reduced considerably. Also, AK's content was different among the groups, while there was no signi cant difference in other nutrients (Fig. S5). In addition, the activities of three enzymes (CAT, INV, and PHO) were signi cantly decreased in the rhizosphere of the severe groups (GRS3 and/or GRS4), while the activity of the URE was not signi cantly changed (Fig. S6). The results indicated that the development of GRS was correlated with the physicochemical properties of rhizosphere. The reasons for the above changes in physicochemical properties may be related to the occurrence of diseases affects the production of enzymes and the conversion of nutrients in the soil or increases their consumption [38,39]. It is also associated with the release and accumulation of secretions from the root of diseased ginseng [2].
To explore the microbial-environmental linkages and gain insight into the key factors in uencing microbial communities' changes, we used dbRDA to examine microbial communities and environmental factors [40,41]. In the dbRDA results for bacteria, TP, AK, CAT, and AP were the main factors determining the distribution of bacterial community composition and taxa changes (Fig. 6A). In the dbRDA analysis of fungi, INV, CAT, and AP were signi cant factors affecting the structure change of the fungal community (Fig. 6B). This study suggests that AP and CAT may focus on future research on soil improvement due to their close correlation with both bacterial and fungal communities.

Conclusions
we studied the changes of soil ecological environment of ginseng rhizosphere with different red skin degree in a same ginseng farm. With the development of the red skin degree, the bacterial community's diversity rst increased and then decreased, and the richness also uctuated. The fungal community's diversity decreased obviously during the severe disease, but the richness was not affected. The rhizosphere of ginseng at different stages showed obvious differences, especially between HG and severe groups (GRS3 and GRS4). Besides, we also explored the changes of microbial taxa in different classi cation levels from HG to GRS4, as well as potential biomarkers for different groups. We established the interaction networks of rhizosphere microorganisms under this disease's background, and the keystone phylum of different groups of bacteria and fungi remained unchanged, namely Proteobacteria and Ascomycota. As the disease progressed, the bacterial community underwent multiple changes from complex and stable to simple and unstable, and the GRS group also had more intercompetition and ecological niche segregation than the HG group. The fungal community's stability decreased signi cantly in the early stages of the disease, followed by forming a stable and complex fungal community. In contrast to bacteria, the GRS groups signi cantly increased interspecies cooperation and ecological niche overlap in the fungal network than the HG group. In the network, we also obtained several microorganisms that are closely related to the potential pathogenic fungi, which can provide the basis for biological control. Soil TP, AK, AP, CAT, INV are the key factors in uencing the microbial communities. Overall, our study comprehensively analyzed the variation characteristics of the ginseng rhizosphere, providing evidence for disease control and soil improvement.

Site description and sample collection
We found that the ginseng grown at a ginseng farm showed more red skin root in Hunchun city (42.86′N and 130.37′E) in northeast China after previous investigations [42]. This area has a temperate oceanic monsoon climate. The ginseng farm had an average annual rainfall of 757 mm and an average yearly temperature of 3.5°C during ginseng cultivation (2015 to 2019). And this farm is the rst to grow ginseng, the application of pesticides strictly comply with the"Ginseng safe production technical speci cation of pesticide application (DB22/T 1233-2019)".
The sampling date is July, 2020. All soil samples were collected from the above ginseng farm with ginseng (5-year-old). The rhizosphere was divided into ve grades according to the degree of disease of ginseng, namely healthy ginseng (HG), rust root area 0 to 25% (GRS1), rust root area 25 to 50% (GRS2), rust root area 50 to 75% (GRS3) and rust root area greater than 75% (GRS4). Six rhizosphere samples were collected at each grade, a total of 30 samples. The rhizosphere was collected according to Riley and Barber's standards [43,44]. In brief, dig out the complete ginseng roots from the soil, gently shake off the large blocks, and then obtain the ginseng root's rhizosphere. Each soil sample was sieved through a 2 mm plastic mesh. Then, one part of each sample was stored at -80℃ for DNA extraction, and the other part was air-dried at ambient temperature for determination of pH, enzyme activity, and soil nutrients.
DNA extraction, PCR ampli cation, and Illumina NovaSeq sequencing Total genome DNA from samples was extracted using the CTAB method [45,46]. DNA concentration and purity were monitored on 1% agarose gels. According to the concentration, DNA was diluted to 1ng/µL using sterile water.
Mix the same volume of 1X loading buffer (contained SYB green) with PCR products and operate electrophoresis on 2% agarose gel for detection. PCR products were mixed in equidensity ratios. Then, the mixture of PCR products was puri ed with the Qiagen Gel Extraction Kit (Qiagen, Germany).
Sequencing libraries were generated using TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, USA) following the manufacturer's recommendations, and index codes were added. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scienti c) and Agilent Bioanalyzer 2100 system. At last, the library was sequenced on an Illumina NovaSeq6000 platform, and 250 bp paired-end reads were generated at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China).

Data analysis
Paired-end reads were assigned to samples based on their unique barcode and truncated by cutting off the barcode and primer sequence. Paired-end reads were merged using FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/), a high-speed and accurate analysis tool, which was designed to merge paired-end reads when at least some of the reads overlap the read generated from the opposite end of the same DNA fragment, and the splicing sequences were called raw tags [48]. Quality ltering on the raw tags was performed under speci c ltering conditions to obtain the high-quality clean tags according to the QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html) quality-controlled process [49,50]. The tags were compared with the reference database (Silva database, using UCHIME algorithm (UCHIME Algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html) to detect chimera sequences, and then the chimera sequences were removed [51,52]. Then the effective tags nally obtained. Detailed quality control information is shown in Table S1 and Table S2.
Sequences analysis was performed by Uparse software (Uparse v7.0.1001, http://drive5.com/uparse/) [53]. Sequences with≤97% similarity were assigned to the same operational taxonomic units (OTUs). Furthermore, all samples' rarefaction curves also indicate the reliability of the sampling depth (Fig. S1). The representative sequence for each OTU was screened for further annotation. For each representative sequence, the Silva Database (http://www.arb-silva.de/) was used based on the Mothur algorithm to annotate taxonomic information [54]. To study the phylogenetic relationship of different OTUs, and the difference of the dominant species in different samples (groups), multiple sequence alignment was conducted using the MUSCLE software (Version 3.8.31, http://www.drive5.com/muscle/) [55]. OTUs abundance information was normalized using a standard of sequence number corresponding to the sample with the least sequences. Subsequent analysis of alpha diversity and beta diversity were all performed basing on this output normalized data.
Rhizosphere pH, nutrients, and enzyme activity The rhizosphere's pH value was measured with a pH meter/potentiometer under the soil: water ratio of 1:2.5. The soil organic matter (OM) content was measured by the potassium dichromate external heating method. The total N (TN) was determined by the Kjeldahl method; the TP was determined by the alkali fusion molybdenum antimony colorimetry (China HJ 632-2011); the total K (TK) was determined by the sodium hydroxide fusion ame photometric method. The AP was determined by NaHCO 3 extraction molybdenum-antimony colorimetry; the available N (AN) was determined by the alkali diffusion method; AK was determined with the ammonium acetate extraction fame photometric method [56].
The activities of CAT, invertase (INV), urease (URE), and phosphatase (PHO) in the soil of ginseng were measured. The activity of CAT was determined by KMnO 4 titration; soil INV activity was determined by 3,5-dinitrosalicylic acid colorimetry; URE activity was determined by indophenol blue colorimetry; PHO activity was determined by disodium phenyl phosphate colorimetry method.

Statistical analysis
Alpha diversity was used to analyze the microbial community diversity in the sample and assess the microbial community's species richness and diversity differences. Five indices were selected to identify microbial alpha diversity: the observed OTUs number, the observed species number, the Chao1 index, Simpson index, and the Shannon index [47]. Beta diversity was used for comparative analysis of microbial community composition in different samples. Cluster analysis was preceded by Non-Metric Multi-Dimensional Scaling (NMDS) analysis was performed to re ect inter-group and intra-group relationships. Unweighted Pair-group Method with Arithmetic Means (UPGMA) Clustering was performed as a hierarchical clustering method to interpret the distance matrix. The Wilcoxon rank-sum test was used to analyze differences in alpha diversity and beta diversity among multiple groups. All of the indices in our samples were calculated with QIIME and displayed with R software.
Metastat analysis was used to carry out permutation test among groups to study the taxa with signi cant differences among groups in R software, and False Discovery Rate (FDR)-corrected P value (Q value) <0.05 was considered statistically signi cant. Besides, to nd statistically different biomarkers, the linear discriminant analysis (LDA) effect size (LEfSe) analysis was used to identify the differentially abundant taxa. LDA score greater than 4 was used as a potential biomarker.
Spearman correlation analysis was used to study the mutual change relationship between environmental factors and microbial taxa, and P<0.05 was considered a signi cant difference [57]. Furthermore, dbRDA analysis was used to re ect the relationship between community structures and environmental factors [58]. Both analyses were performed using the vegan package in R software.
To explore the microbial taxa with close interaction, we calculated the Spearman correlation coe cient for all samples to obtain the taxa correlation coe cient matrix. Then, the co-occurrence network was obtained through ltering. The ltering conditions were as follows: (1) remove connections with the correlation coe cient less than 0.6, (2) lter out self-connection of nodes, and (3) remove connections with node abundance less than 0.005%.
Environmental factor data were analyzed using the SPSS software (IBM Corporation, Armonk, NY, USA), and the results were expressed as the arithmetic mean value ± standard deviation. The differences in the means were compared by the Tukey test at P<0.05.

Declarations
Ethics approval and consent to participate: Not applicable.
Consent for publication: Not applicable.
Availability of data and material: Not applicable.
Competing interests: The authors declare that they have no competing interests. Supplementary Information Table S1. Statistical results of bacterial sequencing data processing. Table S2. Statistical results of fungal sequencing data processing.        Tables   Table 1 Diversity of the 16S rRNA gene-based bacterial and ITS rRNA gene-based fungi communities.

Figure 2
Unweighted Pair-group Method with Arithmetic Mean (UPGMA) clustering analysis with Weighted Unifrac distance matrix and the relative abundance of each sample and group at the phylum level. A) bacteria; B) fungi.

Figure 3
The linear discriminant analysis (LDA) effect size (LEfSe) analysis. A) LDA scores of bacteria differential taxa (LDA score>4); B) diagram of bacterial differential taxa; C) LDA scores of fungi differential taxa (LDA score>4); D) diagram of fungi differential taxa.

Figure 5
Overview of fungi networks. Different nodes represent different genera; the size of nodes represents the average relative abundance of the genus; the nodes of the same gate have the same color; the color of the lines between nodes corresponds to the positive and negative correlation (red is positively correlated; blue is negatively correlated).