Rhizosphere Microbiome is Associated with Yield and Medical Value of the Perennial Panax Notoginseng

Background: Rhizosphere microbiome play important roles in promoting plant growth. However, it is not well understood how rhizosphere microbiome were driven by medical plants during growth stages and whether they contribute the accumulation of medical values. Panax notoginseng is a perennial medicinal plant, which belowground biomass and saponin contents are the important indicators of its value. Here, we use high-throughput sequencing method to study the temporal dynamics of perennial P. notoginseng rhizosphere microbiomes and the relationship between the indicators and core rhizosphere microbiomes. Results: The results show that the diversity, composition and network structures of the bacterial and fungal communities are mainly driven by the developmental stages. And succession characteristics of bacterial and fungal diversity show similar parabolic patterns during the developmental stages. Enrichment and depletion of the bacterial and fungal communities are active at the 3-year root growth (3YR) stage. From samples collected at a large-spatial P. notoginseng production area at the 3YR stage, we obtained 639 bacterial and 310 fungal core operational taxonomic units (OTUs). Analysis of the data indicate that the microbiome diversity is related to the belowground biomass and total saponin contents. Some genera, such as Pseudomonas, Massilia, Sphingobium, and Phoma are positively correlated to the belowground biomass, and genera likely Staphylotrichum, Chaetosphaeria, and Podospora are positively correlated with total saponin contents. Additionally, we identied 36 microbial functions involving in plant-microbe and microbe-microbe interactions, nutrition acquisition, and disease resistance. They are related to belowground biomass and saponin contents. Conclusions: In short, this study provides a detailed and systematic survey of rhizosphere microbiome in P. notoginseng and reveals that P. notoginseng rhizosphere microbiomes are largely driven by the developmental stages, while the core microbiomes are related to the belowground biomass and saponins contents of the plant. The nding may enhance our understanding of the interaction between microbes and perennial plants and improve our ability to manage root microbiota for sustainable production of the herb medicine. beta-alanine metabolism are positively correlated with belowground biomass. Other functions such as fructose and mannose metabolism, cysteine and methionine metabolism are positively correlated with total saponin contents, which also shows positive correlation with the abundances of naphthalene degradation and aminobenzoate degradation. These data suggest that the P. notoginseng plant select functionally benecial microbial communities in the rhizosphere according to its own growth need. Taken together, we determine that core rhizosphere microbiomes associate with belowground biomass and saponin contents of P. notoginseng, revealing the adaptability of rhizosphere microbes regulating the important indicators of value in medicinal plants. Our data also demonstrate that climatic factors, edaphic factors, and rhizosphere microbiomes play synergetic roles in promoting belowground biomass and saponin contents of P. notoginseng.

understanding the succession in rhizosphere during plant growth stages is of great signi cance for improving agricultural practices [9].
Medicinal plants are essential for improving human health, approximately 75% of the population in developing countries mainly rely on herb-based medicines for health care around the world [10]. Many medicinal plants are perennial, including transitions from seedling, owering, to aging annually. In case of Panax ginseng, recent study indicates that the rhizosphere bacterial diversity decreases and fungal diversity increases at the root growth stage compared to vegetative, owering and fruiting stages [11]. Major active components of medicinal plants are the secondary metabolites, whose contents are in uenced by many environmental factors, such as climatic and edaphic factors [12,13]. Nevertheless, rhizosphere microbiomes directly impact the metabolome of host plants, thus in uencing the e cacy of herbal medicines [14]. When Matricaria chamomilla grow in different locations, the contents of secondary metabolites are different and depend on microbe composition [15]. Therefore, the succession characteristics of rhizosphere microbiome during the developmental stages of perennial medicinal plants could mediate the accumulation of secondary metabolism. However, it is not well understood to what extent a perennial medicinal plant during development stages can select a constant and bene cial rhizosphere community. It is of importance for starting a rigorous investigating the interaction and association between rhizosphere microbes and secondary metabolism which traits is crucial for medicinal plant breeding.
Panax notoginseng is an important medicinal plant known as "Nanguo Shencao (Miracle Plant from South China)", and used as the main raw material in Yunnan Baiyao and Xuesaitong due to its bloodinvigorating effects [10,16]. The therapeutic effects of P. notoginseng are mostly attributed to the plant's bioactive saponin constituents, including notoginsenoside R1 and ginsenosides [17,18]. P. notoginseng is a typical perennial plant mainly cultivated in the southwest of China [19]. As a perennial plant, it has vegetative, reproductive and root growth stages annually. It has been found that continuously cropping has reduced the number of rhizobacteria and fungal diversity [20,21]. However, it is unclear how rhizosphere microbial community assembles during annual growth stages. After application of biofertilizer (Burkholderia and Rhizobium), Panax ginseng yield was enhanced by 17.0-19.1% and ginsenosides (Rg1 and Rb1) contents were improved [22]. These results indicate that core rhizosphere microbes may interact with Panax plants to regulate their secondary metabolism. Typically, P. notoginseng has a 3-year-long crop cycle before harvest [23]. Therefore, it is valuable to study the relationship between the core rhizosphere microbiomes and the saponin contens at the harvest stage. This information relates to the content of P. notoginseng saponins and may help to use bene cial microbes to improve the quality and yield of the herb in the future.
We hypothesize that rhizosphere microbiomes are driven by developmental stages and core rhizosphere microbiomes are associated with the quality of P. notoginseng at the harvest stage. Here we utilize amplicon and metagenomic sequencing to analyze the pro les of diversity, composition, and potential functions of rhizosphere microbiomes during different development stages of P. notoginseng. Speci cally, we: (i) characterize the succession characteristics of rhizosphere microbiome in response to the development stages; (ii) identify the core microbiome and function related to the quality of P. notoginseng at the harvest stage.  (Fig. 1). Brie y, one-year-old seedlings were transplanted into each eld and then cultivated in strict accordance with the Good Agricultural Practices [24,25]. Every eld had three separate, 1.4×8.0 m 2 plots as replicates. Before transplanting, ten soil samples were collected from random locations in each plot and combined to generate one bulk soils (BL). Then, rhizosphere soil samples were collected at one of the three developmental stages of P. notoginseng in each year: (1) vegetative (V, May); (2) owering (F, Jul), (3) root growth (R, Oct), namely 2-year V (2YV), 2-year F (2YF), 2-year R (2YR), 3-year V (3YV), 3-year F (3YF), and 3-year R (3YR). At each sampling time, ten randomly selected healthy plants were removed from each plot, and mixed to generate one rhizosphere sample [26]. Thus, a total of 105 samples comprising 90 rhizosphere soil samples at six developmental stages and 15 bulk soil samples in ve elds were obtained, sieved (2 mm) and stored at -80 o C for DNA extraction.

Materials And Methods
Spatial sampling experiment at the harvest stage The belowground biomass and contents of saponins, the most important medicinal components of P. notoginseng, accumulate rapidly at the harvest stage (3YR), an active period of interaction between rhizosphere microbes and plants. To explore the potential contribution of core bacterial and fungal communities to the accumulation of belowground biomass and saponins contents of P. notoginseng, we also conducted a larger spatial scale sampling at the harvest stage in October 2017 across Yunnan and Guangxi provinces, the main production bases of P. notoginseng. Twenty-six sites extending from 22°2′N to 25°42′N and 100°6′E to 106°30′E were selected ( Fig. 1 and Additional le 2). Three 1.4×8.0 m 2 plots were sampled as replicates at each site. Ten healthy plants were randomly collected from each plot, with rhizosphere soil being combined and processed as described above. Besides, corresponding root samples were carefully washed to root biomass and medical value measurements as described below. Totally, 78 rhizosphere soil samples and P. notoginseng roots were obtained. All soil samples were homogenized, sieved (2 mm) and stored at -80 o C for DNA extraction.

Analysis of edaphic and climatic factors
Edaphic factors, namely soil pH, organic matter (OM), available phosphate (AP), available potassium (AK), and total nitrogen (TN) of soil samples from developmental stage experiment were measured using standard test methods2 [27]. For the spatial sampling experiment at the harvest stage, edaphic factors were also tested 2Climatic factors (MMT: Month mean temperature; MAP: Mean annual precipitation; MAT: Mean annual temperature) were acquired from the Worldclim database (www.worldclim.org/) using a geographic information system according to the longitude and latitude values of the collected samples (ArcGIS version 2.0). DNA extraction and amplicon sequencing DNA was extracted from soil samples using the FastDNA SPIN Kit for Soil (MoBio Laboratories, Inc., USA). Ampli cation and puri cation were performed as previously described (Additional le 1: Table S1) [28]. Sequencing was performed on the Illumina MiSeq PE 250 platform (Shanghai Biozeron Co., Ltd., China). Obtained paired-end sequences were demultiplexed and ltered using QIIME2 [29]. Then, chimeric sequences were removed using the USEARCH tool on the basis of UCHIME [30]. Sequences were split into operational taxonomic units (OTUs) with a 97% similarity level using the UPARSE pipeline afterwards [30]. Representative sequences of OTUs were assigned to taxonomic lineages using the RDP classi er against the SILVA database (release 132) [31] for bacteria and UNITE database (release 7.1) for fungi [32]. OTUs that have less than two reads or fail to be aligned to the database (i.e. Unclassi ed) were removed from datasets before further analysis. Finally, 16033 bacterial OTUs and 5147 fungal OTUs were obtained in the developmental stages experiment, while 42988 bacterial and 12796 fungal OTUs were identi ed in the spatial sampling experiment.

Metagenomics sequencing
To explore the response of functions to the quality of P. notoginseng, soil samples from nine sites (i.e., containing a wide variety of total saponin contents (ranging from 4.48% to 11.05%) were selected for shotgun metagenomic sequencing analysis. Metagenomic libraries were constructed using a TruSeq TM DNA PCR-free Sample Prep Kit (Illumina, USA) according to the manufacturer's instructions [33]. The metagenomic libraries were sequenced on a HiSeq 2500 sequencer (Illumina, USA), and 150 bp pairedend reads were generated.
Determination of belowground biomass and saponin contents in P. notoginseng After washing, root samples collected in the second experiment were weighed to obtain belowground biomass. The saponin contents of the P. notoginseng roots were analyzed using HPLC (Agilent 1260) analysis according to previous study [34].

Data analysis
Changes in microbiome composition and diversity during growth stages OTU richness and Shannon index were calculated to estimate the α-diversity of bacterial and fungal communities based on rare ed sequences number using the rrarefy and diversity functions in "vegan" package in R [35,36]. Linear least-square regression analyses with second-order term were performed to t α-diversities to developmental stages using "lm" function. Non-parametric Kruskal-Wallis method was used to further test the differences of α-diversities between different growth stages, followed by Nemenyi test for multiple comparisons using "PMCMR" package [37]. Redundancy analysis (RDA) was performed to assess the in uence of developmental stages and environmental factors on bacterial and fungal communities using rda function in "vegan" package [36]. Then constrained analysis of principal coordinate (CAP) was conducted to visualize the succession of bacterial and fungal communities based on Bray-Curtis dissimilarities using capscale and cmdscale functions in "vegan" package [38,39]. Permutational multivariate analysis of variance (PERMANOVA) was carried out to determine effect size and signi cance of developmental stages on microbial compositions using adonis function in "vegan" package [36]. Differentially abundant taxa among microbial communities obtained from different stages were estimated by tting a negative binomial generalized linear model in "edgeR" package [40]. Brie y, the calcNormFactors function was applied to normalize library size based on trimmed mean of M values method. Then the common and tagwise dispersions were obtained using the estimateDisp function with a design matrix. The glmFit function was performed to t a generalized linear model with a negative binomial distribution, and differential abundant OTUs were tested using likelihood ratio test (glmLRT function); P values were further corrected using "BH" method [41].

Network analysis of microbiomes at different developmental stages
To evaluate the potential variation in interactions among microbial taxa as plant growth, core microbiomes in the samples of each developmental stage were identi ed to construct association networks following two criteria as previous described: the top 10% abundant OTUs, and the OTUs present in more than 85% and 60% samples for bacteria and fungi, respectively [42]. The compositional robust method, SPIEC-EASI (SParse InversE Covariance Estimation for Ecological ASsociation Inference), was adapted to construct association networks using the selected bacterial and fungal core OTUs [43]. All networks were constructed using the SpiecEasi package (v 1.0.7). The neighborhood selection procedure (MB) was used to infer the optimal sparseness based on the Stability Approach to Regularization Selection (StARs) with the variability threshold of 0.05 [44]. Networks were generated and analyzed using the "igraph" package [45]. A series of metrics were calculated to characterize the topological attributes of nodes and networks, including average path length, modularity, degree, betweenness centrality, and closeness centrality. The fast greedy method was applied to detect subgroups in the networks and the modularity was calculated as criterion of the division using the cluster_fast_greey and modularity functions in "igraph" package [46]. The distance between different networks was calculated on the basis of the graphlet degree distribution (GDD) [47]. Multidimensional scaling (MDS) was performed to visualize the distances between networks using cmdscale function in "vegan" package [36]. The robustness of different networks to random and targeted "attack" (node removals) were further assessed using natural connectivity as the measure of graph stability [48]. The random removal of nodes (n=30 repetitions) evaluated the baseline robustness of the network. And targeted "attack" was performed following two strategies, in which nodes were removed in decreasing order of their betweenness or degree. The calculation of natural connectivity, "attack" robustness was performed using customized R script.
Core microbiomes in the spatial sampling dataset For the spatial sampling dataset, variable clustering was conducted to lter environmental factors to reduce collinearity effects among environmental variables using the "Hmisc" package []. Core bacterial and fungal OTUs were selected based on two criteria as described above. To explore the potential correlation between productive characters of P. notoginseng and core microbiomes, random forest (RF) analysis was performed using "randomForest" and "rfPermute" packages [50]. Bacterial and fungal compositions were represented by two axes of NMDS analysis on the core microbiomes [51,52]. Multivariate regression model and variance decomposition analysis were conducted to validate the results of RF analysis. Linear least-square regression analyses were conducted to assess the relationship between OTU richness of core microbiomes and the productive characters of P. notoginseng. Mantel test was conducted to estimate the correlation between β-diversity of core microbiomes and important indicators of medical value (i.e. belowground biomass and saponin contents) based on Bray-Curtis distance of microbiomes and Euclidean distance of indicators using mantel function in "vegan" package [36]. RF analysis was further performed to identify microbial taxa at genus level which potentially contribute to belowground biomass and total saponin contents. The numbers of marker genera were selected based on 20-fold cross-validation using rfcv function in "randomForest" package [50].

Metagenomics analysis
The raw reads from metagenome sequencing were used to generate clean reads by removing adaptor sequences, trimming, and removing low-quality reads (reads with an N base threshold of 10 and a minimum quality threshold of 20) at Biozeron, Shanghai, China. The clean reads were further trimmed using Sickle software, and trimmed reads shorter than 80 bp were discarded [53]. To identify and remove host-originated reads, the trimmed reads were mapped to P. notoginseng genomes using Bowtie 2 [54]. The Q20%, Q30%, and GC% variation ranges were 96.82-97.89, 91.10-93.42, and 58.03-60.94, respectively (Additional le 1: Table S2). The N50 and N90 values ranged from 778-1,297 bp to 543-593 bp, respectively, and the max contigs ranged from 67,400 to 1,151,948 bp (Additional le 1: Table S3). The metagenes were predicted using MetaProdigal (http://prodigal.ornl.gov/). Non-redundant metagenes (unigenes) were obtained using CD-HIT with an identity cut-off of 95% and a coverage cut-off of 90% (http://www.bioinformatics.org/cd-hit/) [55]. To obtain functional information for the unigenes, the protein sequences against the KO database were blasted using DIAMOND, and corresponding KEGG pathways as well as abundance were calculated against KEGG databases (release 93.0) using MinPath [56].

Structural equation model analysis
Structural equation model (SEM) was performed to further validate the casual relations among biotic and abiotic factors in the two datasets using "lavaan" package in R, respectively [57]. The axes of principal component analysis based on scaled variables represented edaphic factors, and axes of principal coordinate analysis based on Bray-Curtis dissimilarities were the representative of variations in microbial communities. For each dataset, group SEM were conducted to x path coe cients among environmental factors and host factors, but keep coe cients free for paths associated with bacteria and fungi. Goodness of t of two models were assessed using chi-square statistics and root-mean-square error of approximation (RMSEA).

Results
Diversity and composition of rhizosphere microbiomes driven by P. notoginseng developmental stages Succession characteristics of α-bacterial and fungal diversity were similar during the developmental stages of P. notoginseng based on the values of richness and Shannon indices ( Fig. 2a and 2b). Linear least-square regression with second-order term shows that the α-diversities of bacteria and fungi are parabolas with plant growth (Richness: R 2 =0.75, P<0.001 for bacteria; R 2 =0.52, P<0.001 for fungi. Shannon indices: R 2 =0.71, P<0.001 for bacteria; R 2 =0.40, P<0.001 for fungi). Compared to bulk soils (BL), values of P. notoginseng for both rhizosphere bacterial and fungal communities increase signi cantly. During the growth of P. notoginseng, bacterial α-diversities increase gradually in soils from 2-yearvegetative (2YV) stage to 2-year-root growth (2YR) stage. Nevertheless, the α-diversities of fungal communities peak in the 2YV stage, then uctuate mildly until 3YF stage. Interestingly, compared to 3year owering (3YF) stage, the richness and Shannon indices in soil bacterial and fungal communities at the 3-year-root growth (3YR) stage are remarkably reduced. These results indicate that the α-diversity in microbial communities is shaped by the developmental stage, and dynamic succession characteristics in bacterial and fungal communities are similar.
We used redundancy analysis (RDA) to preliminarily characterize the most in uential factors to the variance of microbial communities at the developmental stages ( Fig. 2c and Additional le 1: Table S4).
Permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis dissimilarities revealed that developmental stages signi cantly drove the variation within the bacterial (R 2 =37.19%, P<0.001) and fungal (R 2 =26.14%, P<0.001) communities. Constrained and unconstrained analysis of principal coordinates showed distinct separation among samples of different stages, and samples of 3YR stage could be separated from others ( Fig. 2d and Additional le 1: Figure S1). Bray-Curtis distances between rhizosphere soils and bulk soils samples also support that the plant developmental stages shape the bacterial and fungal communities, and Bray-Curtis dissimilarities increase with incremental growth (Fig. 2e). These results demonstrate that developmental stages drive succession characteristics in β-diversity of rhizosphere microbiomes.
Response of network structure to P. notoginseng developmental stages All networks obtained follow the power-law distribution ( Fig. 3 and Additional le 1: Table S5). Bacterial and fungal networks from rhizosphere soils are much larger than those from bulk soils, as indicated by more nodes and edges in Fig. 3a-3c. The MDS based on graphlet degree distribution distances show that bacterial and fungal communities signi cantly determine the disparate networks (R=0.97, P<0.001), and the distances among fungal networks constructed from different stages are much larger than that of bacterial networks (Fig. 3a-3b and Additional le 1: Figure S4). The number of nodes and edges in both bacterial and fungal networks constructed from different developmental stages show an increasing and subsequent decreasing trend (Fig. 3c). Modularity and average path length of rhizosphere bacterial and fungal networks decrease more than that of bulk soils (Fig. 3d-3e). As plants grow, modularity and average path length of the bacterial networks decrease slightly, while that of fungal networks increase. OTUs potentially served as keystone taxa in bacterial and fungal networks also changed with the plant growth (Additional le1: Figure S5). These results indicate that the network structure of rhizosphere microbiomes positively reacts to the developmental stages of P. notoginseng.
We further explored the robustness of different networks through using the natural connectivity (NC) as the index of graph stability (Fig. 3f). In bacterial graphs, networks obtained from the bulk soils are the most robust, while the network robustness from other stages decrease somewhat. In fungal graphs, networks representing the stages of 2YV and 3YV are the two most robust. Interestingly, compared to the network from bulk soils, the natural connectivity of fungal networks constructed from the rhizosphere soils decreases at a faster rate when sequentially removing nodes in the decreasing sequences of degree or betweenness centralities from the network. These results reveal that the transplantation of P. notoginseng reduces the stability of the network.
Core microbiomes related to belowground biomass and saponin content at the harvest stage over spatial scale To further determine the core rhizosphere microbiomes related to quality of P. notoginseng at the 3YR stage (harvest stage), 78 P. notoginseng samples were collected from a large-spatial production area. The belowground biomass, saponincontents and rhizospheric microbiomes of P. notoginseng show a regional speci city at the harvest stage (Additional le 1: Figure S6-S7 and Additional le 7). Additionally, we identi ed 639 bacterial OTUs and 310 fungal OTUs as core microbiomes (Fig. 4a-4b, Additional le 1: Figure S7 and Additional le 8). On average, the core bacterial microbiome accounted for 33.27%±9.07% of the sequences, while the core fungal community accounted for 54.07%±13.74%. The representative phyla include Proteobacteria (43.35%), Acidobacteria (13.81%) and Cyanobacteria (9.82%) in the core bacterial community. At fungal class level, Sordariomycetes (55.04%) and Mortierellomycetes (28.41%) dominate the core community. When a distance-based linear model and forward selection procedure were performed to identify the major environmental variables that shape the core microbial communities, soil pH and OM stand out in both bacterial and fungal communities (Additional le 1: Methods S1, Figure S8 and Table S6).
To relate microbiota to the belowground biomass and saponin contents of P. notoginseng, we identi ed 67 core bacterial genera as the important variables with signi cance level of P<0.05 ( Fig. 4g and Additional le 9). The relative abundances of some genera are signi cantly and positively correlated with belowground biomass (P<0.05), such as Pseudomonas, Massilia, Sphingobium, and Variovorax. While the abundances of Nitrospira and Comamonas present a negative correlation (P<0.05). Additionally, the total saponin contents exhibits a signi cant negative correlation with the abundance of some genera, such as Pseudomonas, Sphingobium, Massilia, and Variovorax.
Based on RF analysis, we also identi ed 53 fungal genera (as the important variables with signi cance level of P<0.05) (Fig. 4h and Additional le 9). The relative abundances of some genera are signi cantly and positively correlated with the belowground biomass (P<0.05), such as Solicoccozyma, Phoma, Fusidium, and Metapochonia, while the abundance of Staphylotrichum, Phomatospora, and Acremonium exhibit a negative correlation. Three genera are signi cantly and positively correlated with the total saponin contents (P<0.05), namely, Staphylotrichum, Chaetosphaeria, and Ganoderma. But three genera, including Devriesia, Phoma, and Metapochonia, show a signi cant negative correlation. These results indicate that the core bacterial and fungal genera are related to the belowground biomass and total saponin contents of P. notoginseng.
Functions of rhizosphere microbiome related to the belowground biomass and saponin contents in P. notoginseng We obtained a total of 4616 KEGG orthology (KO) functional categories and found that the highly abundant KOs mainly involve several KEGG functional categories of level 2, such as "carbohydrate metabolism", "amino acid metabolism", "lipid metabolism", "metabolism of terpenoids and polyketides", "xenobiotics biodegradation and metabolism", "membrane transport" and "signal transduction" (Fig. 5a and Additional le 10-11).
Some functions related to belowground biomass and saponin contents of P. notoginseng were obtained in 257 KEGG functional categories of level 3 (Fig. 5a, Additional le 1: Table S8 and Additional le 12). Among them, 36 functions involve in plant-microbe and microbe-microbe interactions, such as bacterial secretion system (ko03070), agellar assembly (ko02040), bacterial chemotaxis (ko02030), and twocomponent system (ko02020). Twenty-ve microbial functions mediate nutrition acquisition and plant growth, among which galactose metabolism (ko00052), phenylalanine metabolism (ko00360) and betaalanine metabolism (ko00410) show a positive correlation with belowground biomass (Fig. 5b and Additional le 1: Table S9). Meanwhile positive correlations are also observed between the abundance of fructose and mannose metabolism (ko00051), cysteine and methionine metabolism (ko00270) and total saponin contents. The abundances of functions mediating disease resistance, such as toluene degradation (ko00623), nitrotoluene degradation (ko00633) show a positive correlation with belowground biomass. Meanwhile, the abundances of naphthalene degradation (ko00626) and aminobenzoate degradation (ko00627) are positively correlated with the total saponin contents ( Fig. 5b and Additional le 1: Table S9).
Rhizosphere microbiomes driven by P. notoginseng developmental stages and core microbiome related to its root growth and saponin contents Structural equation model (SEM) was constructed to evaluate the effects of environmental factors and developmental stages on rhizosphere microbiomes (Fig. 6a). Month mean temperature (MMT), edaphic factors, and developmental stages were identi ed as the signi cant drivers that in uence the rhizosphere microbiomes (Fig. 6b). Developmental stages were found to have a signi cantly direct effect on the PCoA1 (bacteria, 0.76, P<0.001; fungi, 0.34, P<0.001), PCoA2 (bacteria, 0.19, P<0.05; fungi, 0.60, P<0.001) and bacterial richness (0.42, P<0.001). These data further reveal that P. notoginseng developmental stages drive its rhizosphere microbiomes.
Although edaphic factors, climatic factors and rhizosphere microbiomes collectively determine the variability in belowground biomass and total saponin contents of P. notoginseng (Fig. 6c), the core bacterial and fungal communities show driving effects. The core bacterial community has a remarkably direct effect on both biomass (PCoA1, 0.29, P<0.001) and content (richness, 0.28, P<0.05; PCoA4, 0.30, P<0.05). The fungal bacterial community also has a signi cantly and positively direct effect on the contents (richness, 0.15, P<0.05; PCoA2, 16, P<0.05) but a negative effect on the biomass (PCoA1, -0.27, P<0.01; PCoA2, -0.18, P<0.05). These data further suggest that core rhizosphere microbiomes are related to belowground biomass and total saponin contents.

Discussion
Here, we determined the spatiotemporal dynamics of diversity, composition, and network structures of fungal and bacterial microbiomes during the overall growth stages of P. notoginseng and evaluated its rhizosphere microbiomes at a spatial level to determine their core microbiomes and functions related to the belowground biomass and saponin contents. Our results shows that rhizosphere microbiomes are driven by P. notoginseng developmental stages, and variations of rhizosphere bacterial communities are also reported during the development stages of Arabidopsis thaliana [6]. Compared to bulk soils, the bacterial and fungal diversity exhibits increasing trends in rhizosphere soils of P. notoginseng. The αdiversity of rhizobacteria in rhizosphere soils of cotton (Gossypium hirsutum L.) were also increased compared with bulk soils [58]. Meanwhile, we noticed that bacterial and fungal diversity increased and then decreased during P. notoginseng growth, which is similar to the observations in P. ginseng [59] and Pseudostellaria heterophylla [60]. These data reveal that succession characteristics of bacterial and fungal diversity are similar as the P. notoginseng grew. Our constrained analysis of principal coordinates shows that there is a clear boundary between the bulk and rhizosphere samples. This indicates that compared with the bulk soil community, rhizosphere bacterial and fungal community structures have a host-selective effect, which is consistent with other reports [61,62]. We observed a sharp increase in Cyanobacteria in the soils of 3YR stage. This may be because Cyanobacteria have a strong ability to colonize the roots of plant and promote plant growth to meet the needs of the rapid root growth at this stage [63]. Putting together, these results indicate that P. notoginseng can select a subset of microbes and build up succession characteristics according to its growth needs, and the strength of its rhizosphere effect varies with its growth.
P. notoginseng has a complex network structure in the rhizosphere microbiomes. Compared to the bulk microbiomes, the network is both larger and more complex, a nding similar to the observation in the wild oat (Avena fatua) [64]. A large complex network may facilitate more interactions and niche-sharing. However, after transplantation, the modularity and average path length of the bacterial and fungal networks decreased, re ecting the homogeneity of the rhizosphere environment driven by the root exudates [65,66]. In other words, P. notoginseng transplantation reduces the stability of the network. In a previous study, antibiotic treatment in early life also reduced the networks stability of intestinal microbiota [67]. Thus network complexity, a previously undescribed property of the rhizosphere microbiome, appears to be a characteristic of this habitat. The release of exudates relies on the demand of plants, and varies in terms of plant developmental stages, thus leading to the dynamic changes of microbiomes [6,68]. Plant exudation, such as amino acids, sugars, phenolics, and sugar alcohols, changed with plant developmental stages [69][70][71]. Sugar exudation could accumulate during Acidobacteria growth [72]. High phenolic acid (-phenol 2,4-di-tert-butylphenol and vanillic acid) concentrations inhibited microbial biomass [73], and phenolics showed a negative correlation with phyla Actinobacteria [74]. Therefore, the rhizosphere microbiome of the 3-year root growth stage is distinguished from that of the other developmental stages of P. notoginseng possibly due to the in uence of the root exudates. P. notoginseng root in ate and grow rapidly at the 3-year root growth stages (harvest stage). Thus, it is particularly valuable to analyze the rhizosphere microbiomes at this stage to understand the microbiome functions of the perennial plant. Our data show that the core bacterial and fungal communities are related to the belowground biomass and total saponin contents of P. notoginseng. In foxtail millet, Bacillales, Actinomycetales, Rhizobiales, and Sphigobacteriales were reported to potentially correlate with the high plant productivity . In our study, abundances of some bacterial and fungal genera are signi cantly and positively correlated with the belowground biomass and total saponin contents of P. notoginseng. Many microbes, such as Bacillus altitudinis and Paenibacillus polymyxa, are bene cial for P. ginseng growth, ginsenoside accumulation, and morbidity reduction [76,77]. The application of Burkholderia and Rhizobium increased the root growth and ginsenoside contents in P. ginseng [22]. These microbes are not only associated with plant growth and secondary metabolites but also can maintain hormonal balance, control root development, promote nutrient acquisition, and prevent diseases of plant, so that they were known as plant bene cial microbes [78,79]. Therefore, it is important to understand rhizosphere microbes in order to strike a balance between high quality and sustainable production of the herb medicines.
Plant-microbe and microbe-microbe interactions have important effects on rhizosphere microbiomes [80]. This can be evaluated by the correlation between the abundance of functions related to environmental information processing and the belowground biomass. This type of functions commonly includes ABC transporters for exchange of complex molecules, two-component system and bacterial secretion system for intercellular signaling, as well as N metabolism pathway [81][82][83][84][85][86]. In our study, certain functions such as galactose metabolism, phenylalanine metabolism and beta-alanine metabolism are positively correlated with belowground biomass. Other functions such as fructose and mannose metabolism, cysteine and methionine metabolism are positively correlated with total saponin contents, which also shows positive correlation with the abundances of naphthalene degradation and aminobenzoate degradation. These data suggest that the P. notoginseng plant select functionally bene cial microbial communities in the rhizosphere according to its own growth need. Taken together, we determine that core rhizosphere microbiomes associate with belowground biomass and saponin contents of P. notoginseng, revealing the adaptability of rhizosphere microbes regulating the important indicators of value in medicinal plants. Our data also demonstrate that climatic factors, edaphic factors, and rhizosphere microbiomes play synergetic roles in promoting belowground biomass and saponin contents of P. notoginseng.

Conclusion
Our data show that developmental stage is a signi cant driving force affecting the rhizosphere microbiomes, and succession characteristics of bacterial and fungal diversity were similar and parabolic patterns during the developmental stages, re ecting the adaptability of the plant microbial communities to the changes in the plant development. Association analysis at the large-scale reveals that the core rhizosphere microbiome is related to the important medicinal value indicators of the perennial P. notoginseng. This work provides a comprehensive understanding in predicting the response of microbial communities and their potential functions in plant growth and accumulation of secondary metabolites. In addition to harness the power of microbiomes to evaluate herbal medicine quality and quantity, our data also provide valuable information in the guidance of microbial breeding for medical plants. The model of experimental samples collection. Blue and red boxes represent models of developmental stages and spatial sampling experiments, respectively. Bulk soils were collected from farmlands before P. notoginseng transplanting. Rhizosphere soil samples were collected from P. notoginseng plants.