Rhizosphere soil samples of P. notoginseng were collected at five fields in Yunnan Province of China during a 2-year growth from February 2016 to October 2017 (Figure 1; Supplementary Table S1). Briefly, 1-year-old seedlings were transplanted into each field and then cultivated in strict accordance with the Good Agricultural Practices (Heuberger et al., 2010; Zhang et al., 2021b). Every experimental field had three separated 1.4 m² × 8.0 m² plots as replicates. Before transplanting, 10 soil samples were collected from random locations in each plot and combined to generate one bulk soil (BL). The rhizosphere soil samples were then collected at the three developmental stages of P. notoginseng in each year: (1) vegetative (V, May); (2) flowering (F, Jul); and (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, 10 randomly selected healthy plants were removed from each plot and mixed to generate one rhizosphere sample (Dong et al., 2018a). A total of 105 samples comprising 90 rhizosphere soil samples at six developmental stages and 15 unplanted BL samples in five fields were obtained, sieved (2 mm), and stored at −80°C for DNA extraction.

Edaphic factors, namely, soil pH, organic matter (OM), available phosphate (AP), available potassium (AK), and total nitrogen (TN) were measured for the soil samples using standard test methods (Supplementary Table S1). The monthly mean temperature (MMT) of sampling time was 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) (Supplementary Table S1).

Total DNA from soil samples was extracted using the FastDNA SPIN Kit for Soil (MoBio Laboratories Incorporation, California, USA) to characterize the bacterial and fungal communities. 515F (5′-GTGCCAGCMGCCGCGG-3′)/907R (5′- CCGTCAATTCMTTTRAGTTT-3′) primer pair was used to amplify the V4–V5 region of bacterial 16S rRNA gene, and ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′)/ITS2R (5′- GCTGCGTTCTTCATCGATGC-3′) primer pair was used to amplify the fungal ITS region (Mueller et al., 2014; Jiao et al., 2019). Sequencing was performed on Illumina MiSeq PE 250 platform (Biozeron Corporation Ltd., Shanghai, China). The obtained paired-end sequences were demultiplexed, merged, and filtered using QIIME2 and USEARCH (Edgar, 2013; Bolyen et al., 2019). Chimeric sequences were then removed using the USEARCH (Edgar et al., 2011). The remaining sequences were split into operational taxonomic units (OTUs) with a 97% similarity level using the UPARSE pipeline (Edgar et al., 2011). Representative sequences of OTUs were assigned to taxonomic lineages using the Ribosomal Database Project (RDP) classifier against the SILVA database (release 132) for bacteria and UNITE database (release 7.1) for fungi (Quast et al., 2013; Nilsson et al., 2019). OTUs with less than two reads or fail to be aligned to the database (i.e., unclassified) were removed before further analysis. Finally, 16,033 bacterial operational taxonomic units (OTUs) and 5,147 fungal OTUs were obtained. The raw sequence data were uploaded to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the bioproject number PRJNA559079. The relative abundance of OTUs was used in subsequence analyses except for the estimation of α-diversities and edgeR.

Operational taxonomic unit richness and Shannon index were calculated according to the rarefied sequence number to estimate the α-diversities of bacterial and fungal communities (bacteria: 44,332, fungi: 38,847) using the “rrarefy” and “diversity” functions in “vegan” package in R (Oksanen et al., 2019; Team RC, 2019). Linear least-square regression analysis with the second-order term was performed to fit the α-diversities to developmental stages using “lm” function. Nonparametric Kruskal–Wallis method was employed to further test the differences of α-diversities between different growth stages, followed by the Nemenyi test for multiple comparisons. Permutational multivariate ANOVA (PERMANOVA) was carried out to determine the effect size and significance of developmental stages on microbial compositions using “adonis” function in “vegan” package (Oksanen et al., 2019). Differentially abundant genera among microbial communities obtained from different stages were estimated by fitting a negative binomial generalized linear model in “edgeR” package (Robinson et al., 2010). After genera with CPM less than 100 in three samples were removed, “calcNormFactors” function was applied to normalize the library size according to the trimmed mean of the M value method. Common and tagwise dispersions were obtained using the “estimateDisp” function with a design matrix. The “glmFit” function was employed to fit a generalized linear model with a negative binomial distribution, and differential abundant genera were tested using the likelihood ratio test (glmLRT function). P-values were further corrected using the “BH” method (Benjamini and Hochberg, 1995).

A cross-domain co-occurrence network was first constructed for each developmental stage based on correlations to reveal the variations in potential interaction patterns of rhizosphere microbiomes during P. notoginseng growth. Only bacterial and fungal OTUs with a relative abundance >0.01%, and spearman correlations > 0.7 or < −0.7 with the false discovery rate (FDR) corrected P < 0.01 were used to construct the network (Zhang et al., 2021b). A set of node-level and network-level topological properties were then calculated in the “igraph” package in R (Ma et al., 2017). The node-level feature set included degree centrality (the number of adjacent edges), betweenness centrality (the number of shortest paths going through a node), closeness centrality (the number of steps required to access all other nodes from a given node), and eigenvector centrality (the tendency for a certain node to share connections with other nodes that link to many other taxa). These node topologies potentially represent the role of nodes in maintaining co-occurrence patterns (Jiao et al., 2017; Zhang et al., 2021b). The network-level feature set consisted of complexity measures (i.e., graph density, average path length, and clustering coefficient) and stability measures (i.e., natural connectivity) (Jun et al., 2010; Ma et al., 2016; Fan et al., 2018; Zhang et al., 2021a). Wilcox rank-sum test was performed to compare the node-level properties between bacterial and fungal nodes. Bacterial and fungal subnetworks at each developmental stage were further extracted using the “igraph” package in R.

Redundancy analysis was performed to assess the influence of developmental stages and environmental factors on bacterial and fungal communities using the “rda” function in the “vegan” package (Oksanen et al., 2019). Constrained analysis of principal coordinates (CAP) was then conducted to visualize the succession of bacterial and fungal communities based on Bray–Curtis dissimilarities using “capscale” and “cmdscale” functions in the “vegan” package (Oksanen et al., 2019).

A structural equation model was also established to further explore the potential causal relationships among biotic and abiotic factors using the “lavaan” package in R (Rosseel, 2012). The representative of variations in microbial communities were the axes of principal component analysis based on scaled variables represented edaphic factors and the axes of principal coordinate analysis (PCoA) based on Bray–Curtis dissimilarities. The goodness of fit of the two models was assessed using chi-square statistics and root mean square error of approximation (RMSEA).

The succession characteristics of bacterial and fungal α-diversities during the developmental stages were analyzed using OTU richness and Shannon indices (Figures 2A,B). Linear least-square regression with second-order term showed that the α-diversities of bacteria and fungi were parabolic with plant growth (richness: R² = 0.75, P < 0.001 for bacteria; R² = 0.52, P < 0.001 for fungi. Shannon indices: R² = 0.71, P < 0.001 for bacteria; R² = 0.40, P < 0.001 for fungi). Compared with those of BLs, the values of P. notoginseng for rhizosphere bacterial and fungal communities increased significantly. During P. notoginseng growth, the bacterial α-diversities increased gradually in soil from the 2YV stage to the 2YR stage. Meanwhile, the fungal α-diversities peaked at the 2YV stage and then fluctuated mildly until the 3YF stage. Interestingly, compared to the 3YF stage, the richness and Shannon indices of soil bacterial and fungal communities at the 3-year-root growth (3YR) stage were remarkably reduced. These results revealed the α-diversities of rhizosphere communities first increased and then decreased with plant growth and the dynamic succession characteristics were similar for bacterial and fungal communities.

Permutational multivariate ANOVA (PERMANOVA) based on Bray–Curtis dissimilarities revealed that the developmental stages significantly drove the variations in the bacterial (R² = 37.19%, P < 0.001) and fungal (R² = 26.14%, P < 0.001) β-diversities (Figure 2C). Constrained analysis of principal coordinates showed a distinct separation among different stages, and the samples of the 3YR stage could be separated from the others. Bray–Curtis dissimilarities between rhizosphere soils and BLs increased with incremental growth (Figure 2D). These results indicated that the biodiversity and community structure of rhizosphere bacterial and fungal communities at different stages of host plant development might be significantly related to the succession of host plant characteristics.

The microbial compositions in the rhizosphere soil of P. notoginseng also changed with the developmental stages (Figures 2E,F; Supplementary Figures S1,S2). Compared with those in BLs, the relative abundance of bacterial genera, namely, Bacillus, Paenibacillus, and other genera from Firmicutes largely decreased (Figure 2E). At the 3YR harvest stage, the relative abundance of Cyanobacteria and certain bacterial genera, namely, Dyadobacter, Dysgonomonas, Burkholderia, Delftia, Flavobacterium, Lechevalieria, and Hyphomicrobium increased significantly. Some fungal genera such as Cladosporium, Epicoccum, Leptosphaerulina, Phoma, Cladophialophora, Exophiala, Fusarium, Plectosphaerella, and Ilyonectria were enriched at the 3YR stage (Figure 2F). These results revealed that the composition variations of bacterial and fungal communities were active at the 3YR stage.

According to the cross-domain co-occurrence network constructed from each developmental stage, the variation trend of node number was similar to that of α-diversities (Figures 3A, 4A). However, the number of edges of rhizosphere networks was lower than that of the network obtained from BL, and the 3YR network exhibited the lowest number of edges (Figure 3B). The network-level topological properties, namely, average degree, graph density, natural connectivity, and clustering coefficient, also decreased after plant cultivation (Figure 4A). The average path lengths of rhizosphere networks were higher than those of BL networks (Figure 4A). In addition, the topologies of the 3YR network showed dramatic changes compared with those in previous stages (Figures 3A,B).

For the proportions of different edge types, the proportion of negative edges decreased in the rhizosphere soil compared with that in the BL (Figure 3C). At the 3YR stage, the ratios of bacteria–bacteria negative (BBN) edges and bacteria–bacteria positive (BBP) edges declined and those of fungi-associated links, especially fungi–fungi positive (FFP) edges increased. In the BL network, the node-level properties of bacterial and fungal nodes only exhibited slight differences (Wilcoxon rank-sum test, P > 0.05), except for the eigenvector centrality (Wilcoxon rank-sum test, P < 0.05; Figure 3D). After plant cultivation, the topological centralities of bacterial nodes became significantly higher than those of fungal nodes until the 3YF stage. The fungal nodes showed remarkably higher centrality values than the bacterial nodes at the 3YR stage (Wilcoxon rank-sum test, P < 0.05), except for the closeness centrality (Wilcoxon rank-sum test, P = 0.87).

For the single-domain subnetworks, the topological characters of bacterial networks showed consistent trends with the cross-domain networks (Figure 4B; Supplementary Figure S3). However, the fungal network exhibited an opposite trend at the 3YR stage, that is, the graph density, natural connectivity, and clustering coefficient of the fungal network at the 3YR stage increased compared with those of the networks of previous developmental stages. These results showed that the network variations of bacterial and fungal communities were active at the 3YR stage.

Monthly mean temperature (MMT), edaphic factors, and developmental stages were identified as significant drivers that influenced the rhizosphere microbiomes. RDA was used to preliminarily characterize the effects of environmental factors and developmental stages on rhizosphere microbiomes during plant growth (Figures 5A,B; Supplementary Table S2). The results showed that the developmental stages were the most important variables in bacterial (6.76%, P < 0.001) and fungal communities (5.26%, P < 0.001).

A SEM was then constructed to evaluate the most influential factors of the variance in microbial communities during plant growth (Figure 5C). The developmental stages were found to have a significant direct effect on 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). An indirect positive effect of development stages on PCoA1 (R = 0.39, P < 0.001) was also found via their impacts on richness (R = 0.420, P < 0.001). These data further emphasized that the developmental stages of P. notoginseng could significantly drive the rhizosphere microbiomes.