Collection of 14 batches of seedlings of annual A. sinensis in Min County, Tanchang County, Zhang County, Lintan County and Zhuoni County in Gansu Province, China. The seedlings were planted in Min County Medicinal Plants Growing Technology Extension Centre (34°22′30” N, 104°53′20” E; black soil, pH ≈ 8.0) with same measures of cultivation and management during the planting period, and it was cultivated to grow for 180 days then harvest. Three samples for each batch were randomly selected, each sample was comprised of 5 healthy plants, and 42 samples were taken for each growth stage. The seedlings and mature roots were both authenticated by Dr. Hui Yan. Carefully, uproot fresh plants and shake loose bulk soil that clings to the roots. Remove the rhizosphere soil which tightly attached to roots with a sterile brush. Fresh plants and soil samples were stored at -80 °C refrigerator for analyses.

According to manufacturer’s approach, total DNA was extracted from per sample using the HiPure Soil DNA Kits (Magen, Guangzhou, China). The bacterial sequence was amplified with the primers 806R (5’-GGACTACHVGGGTATCTAAT-3’) and 341F (5’-CCTACGGGNGGCWGCAG-3’). PCR reactions were performed in triplicate in 50 μL mixture containing 3 μL of 25 mM MgSO₄, 5 μL of 2 mM dNTPs, 5 μL of 10 × KOD Buffer, 1.5 μL of each primer (10 μM), 100 ng of template DNA, and 1μL of KOD Polymerase. The bacterial V3-V4 hypervariable region of 16S rRNA was amplified by PCR (94°C for 2 min, followed by 30 cycles at 98°C for 10 s, 62-66°C for 30 s, and 68°C for 30 s, and a final extension at 68°C for 5 min) using primers. Related PCR reagents were from TOYOBO, Japan. According to the manufacturer, the PCR amplicons were extracted from 2% agarose gels and purified using the AMPure XP Beads (Beckman Agencourt, USA), and quantified using ABI StepOnePlus Real-Time System (Life Technologies, Foster City, USA). Purified amplicons were pooled in equimolar and paried-end sequenced (PE250) on an Illumina platform in accordance with the standard protocols.

FASTP (Chen et al., 2018) (version 0.18.0) mainly used for the quality control of the raw reads which containing less than 50% of bases with quality (Q-value)>20 and more than 10% of unknown nucleotides (N). The splicing was done by FLASH (Magoc and Salzberg, 2011) (version 1.2.11). The clean tags were grouped into operational taxonomic units (ASVs) by UPARSE (Edgar, 2013) (version 9.2.64) according to the similarity ≥ 97%. The representative OUT sequences was analyzed by RDP classifier (Wang et al., 2007) (version 2.2) based on SILVA (Pruesse et al., 2007) database (version 132), were classified into organisms according to the confidence threshold of 0.8. Diversity indices of data were performed in QIIME (version 1.9.1) and R packages (version 2.5.3). Statistic analysis of Anosim test, Tukey’s HSD test, Welch’s t-test, Kruskal-Wallis H test, Wilcoxon rank test and Adonis was calculated in R project Vegan package (version 2.5.3).

The ACQUITY™ UPLC system was used to perform the medicinal components of A. sinensis with an A Thermo Syncronis C₁₈ (2.1 mm × 100 mm, 1.7 μm) column. The 0.1% formic acid in chromatographic (A) and acetonitrile (B) were used as mobile phase. The gradient elution were as follow: 0 ~ 2 min, 5% ~ 10% B; 2 ~ 10 min, 10% ~ 40% B; 10 ~ 13 min, 40% ~ 40% B; 13 ~ 19 min, 40% ~ 50.8% B; 19 ~ 23 min, 50.8% ~ 90% B; 23 ~ 24 min, 90% ~ 90% B; 24 ~ 24.5 min, 90% ~ 5% B; 24.5 ~ 26 min, 5% ~ 5% B The detection wavelength are 260 nm, 280 nm and 320 nm. The medicinal components (ferulic acid, chlorogenic acid, Z-ligustilide, senkyunolide A, senkyunolide H, senkyunolide I, n-butylphthalide, coniferyl ferulate) were measured by comparing with the calibration curves. The 8 reference compounds bought from National Institutes for Food and Drug Control (China), Shanghai Macklin Biochemical Co., Ltd (China) and Nanjing Jin Yibai Biological Technology Co., Ltd (China) respectively. The Waters ACQUITY™ Synapt Q-TOF mass spectrometer equipped with an electrospray ionization (ESI) source. ESI-MS spectra was generated in both positive (ESI⁺) and negative (ESI^-) ion modes through scanning from 100-1000 Da. The flow rate was 0.4 mL·min^-1, column temperature was kept at 35°C, with 2 μL injection volume. The optimized conditions include: capillary voltage 3 kV for both positive ion mode and negative ion mode; source temperature 120°C; desolvation gas temperature 350°C; desolvation gas flow 600 L/h; sample cone voltage 30 V; cone gas flow 50 L/h; extraction cone voltage 2.0 V; secondary collision energy 25 ~ 45 V. Quantitative mass spectrometric analysis was performed using Xevo Triple Quadrupole MS (Waters Corp., Milford, MA USA) equipped with electrospray. ESI-MS spectra were obtained by using MRM mode.

Using IBM SPSS Statistics 19 and R (4.2.0) to perform all statistical analyses. The mean and standard deviation were calculated by multiple comparison analysis and analysis of variance (ANOVA) for statistical tests. The univariate approach depends on t-tests (or their nonparametric alternatives). Spearman correlation coefficients were calculated using R v 4.2.0. RStudio v 2022.02.1 and SIMICA v14.1 were employed to perform PCA and OPLS-DA. Differences between groups were considered significant when p < 0.05. TBtools v 1.09876 (Chen et al., 2020) was employed to perform heatmaps and hierarchical clustering. Statistical analysis of rhizosphere microorganisms taxonomic differences between groups using STAMP v 2.1.3. Cytoscape v3.9.1 and MetScape v3.1.3 were employed to visualize the co-occurrence networks of microbiota and metabolome data. Results were built and optimized using OriginPro 2021.

The depth of sequencing in this experiment has basically covered all the species in the sample, as shown in the Figure S1 . As shown in Table 1 , the indices representing alpha-diversity were estimated as community richness and evenness. The diversity indices were lower in G compared with M, however, no significant differences in bacterial richness or diversity between M and G, as indicated by Sobs, Shannon, Simpson, Chao1, Pielou, and Pd. There are significant differences in coverage of ASV with low abundance between M and G, as indicated by Goods coverage index. The ACE index was significantly lower in G compared with M, showing that the richness and evenness of community in M were significantly higher than in G. NMDS was utillzed to analyze the comparability of the assemblage of bacterial communities between group M and G at the ASV level ( Figure 1 ). Group M was located in the left part of the graph, and the distribution of different samples was discrete. The group G was gathered at the right, and the distribution of different samples was centralized. The results showed a distinctive difference in bacterial community structures between group M and G.

The species composition of rhizosphere microorganisms was analyzed at phylum and genus level. Figure 2A shows Proteobacteria, Acidobacteria, Actinobacteria, Bacteroidetes, Planctomycetes, Verrucomicrobia, Gemmatimonadetes, Chloroflexi, Firmicutes, Patescibacteria are most abundant at the phylum level, representing about 96.6% of the microorganisms detected in the 84 soil samples. Compared with group M, the abundance of Proteobacteria, Gemmatimonadetes, Patescibacteria, and Firmicutes have significantly (p <0.05) increased while the abundance of Actinobacteria, Bacteroidetes, Verrucomicrobia has significantly (p <0.05) decreased in group G. Proteobacteria, Acidobacteria, Actinobacteria were the most three abundant in the rhizosphere of A. sinensis. The relative abundance of Actinobacteria declined, showed a negative response to plant growth. The relative abundance of Proteobacteria increased, showed a positive response to plant growth. The relative abundance of Acidobacteria showed no significant difference during the plant growth. This observation is consistent with many other successional studies: including Arabidopsis thaliana (Roller et al., 2016), wheat (Bulgarelli et al., 2015), rice (Lu et al., 2006), switchgrass (Mao et al., 2014), maize (Peiffer et al., 2013) and Avena barbata (Zhalnina et al., 2018). It is indicating that the reconstruction of rhizosphere microbial communities had a general rule in different plant species and soil types. Research showed (Zhalnina et al., 2018), rhizosphere microbes that respond positively to plant growth are predicted to have longer generation times based on codon-usage bias, which means that their genomes are characterized by slower growth rates. Since slower-growing organisms can have higher substrate utilization efficiency (Chaparro et al., 2013), growth efficiency may be preferred over growth rate in the rhizosphere. Top 10 phyla with significant differences between group M and G are shown in Figure 2B .

At the genus level, the dominant genera were Pseudomonas, Sphingomonas, RB41, Flavobacterium, Pedobacter, Candidatus Udaeobacter, Stenotrophomonas, Acinetobacter, Gemmatimonas, Chthoniobacter. The hierarchical clustering of TOP 35 genera showed that there were significant differences in microbial community structure between group M and G as shown in Figure 3A . The dominant genera also showed significant differences. Compared with group M, the abundance of Acinetobacter, Gemmatimonas, and Sphingomonas has significantly (p <0.05) increased while the abundance of RB41, Pedobacter, Candidatus Udaeobacter, Chthoniobacter has significantly (p <0.05) decreased in group G as shown in Figure 3B .

Untargeted metabolomics facilitates botanical metabolomics studies through efficient high-throughput screening techniques. At present, LC-MS/MS has been widely used for metabolomics studying (Tsugawa et al., 2021). In this study, UPLC-QTOF-MS/MS method was used to characterize metabolites in A. sinensis samples between M and G. Principal components analysis (PCA) ( Figure 4A ) and the OPLS-DA model ( Figure S2 ) showed that the metabolic profiles varied significantly (p < 0.05) between groups, and the samples were clustered into two groups significantly. The hierarchical clustering heatmap ( Figure 4B ) showed the differences between group M and G, suggesting that the synthesis of secondary metabolites pattern at seedling stage and collection period of medicinal material are distinct. Further analysis ( Table S1 ) combined with a hierarchical clustering heatmap ( Figure 4B ) revealed that group G had a higher level of organic acids, such as dicaffeoylquinic acid, chlorogenic acid, genipic acid, and ferulic acid. In addition, the contents of comuside, phenylethyl primeveroside and QUIZALOFOP-ETHYL, 2-dicyclohexylphosphino-2’,6’-dimethoxybiphenyl and [12]-gingerdione in group G were higher than M, while the levels of phenylpropanoids, phthalide and other organoheterocyclic compounds, flavonoids, amines, and fatty acid were lower compared to group M. These results suggest that different growth stages of A. sinensis have different secondary metabolite synthesis patterns, which indicate that changes in the composition of the secondary metabolites over time may contribute to the observed successional patterns in the rhizosphere microbiota ( Figure 1 ). Metabolites with the Variable Importance for Projection (VIP) value greater than 1 and the p-value less than 0.05, combined with s-plot (covariance greater than 0.05) were screened as markers that contributing to grouping. A total of 29 markers were identified. ( Table 2 ) (Zhang et al., 2019).

To further explain the cause of differences in metabolites of A. sinensis between group M and G, phthalides and organic acids were comprehensively analyzed. The result showed that the contents of medicinal components differed remarkably (p < 0.05) between groups. Specially, the level of coniferyl ferulate, senkyunolide A, and ligustilide of A. sinensis in group M were significantly (p < 0.01) higher than in group G ( Figure 5 ). The level of chlorogenic acid in group G was significantly higher than in group M (p < 0.001).

The utilization of specific chemical components of exudates by microorganisms has potential metabolic differences, which may be crucial for bacterial success in the rhizosphere (Zhalnina et al., 2018). Thus, the chemical composition of plant secondary metabolites synthesis and plant exudates may represent key means of shaping the microbial composition of the rhizosphere. The hierarchical clustering heatmap showed that the differential metabolites were clustered into two groups ( Figure 4C ), and the distribution and type of differential metabolites are highly consistent with Figure 4B . The results showed that the rhizosphere microbial community have different assembly patterns between group M and G, and the rhizosphere microbial community had a key contribution to grouping and quality-related factors of A. sinensis. Combined Figures 4B, C , we found that Actinobacteria, Proteobacteria and Bacteroidetes were predominant in synthesizing secondary metabolites such as phenylpropanoids, phthalide, and other organoheterocyclic compounds, flavonoids, amines, and fatty acid in the seedling stage or early growth stage of A. sinensis. Proteobacteria, Acidobacteria and Firmicutes were predominant in the synthesis of organic acid in the drug-producing stage or later growth stage of A. sinensis. Aromatic organic acids, such as cinnamic acid, have been found to shape plant rhizosphere microbes and influence plant-microbe interactions (Sasse et al., 2018; Cotton et al., 2019). Therefore, the synthesis of medicinal components and differential secondary metabolites is closely related to the growth stage of A. sinensis. The synthesis and secretion of secondary metabolites also affect the assembly of the rhizosphere microbial community.

In the case pf medicinal plants, specific microorganisms may be directly related to the biosynthesis of medicinal components of the host plant. For example, Lysobacter was identified strongly associated with gene CYP72A154, which was required glycyrrhizic acid biosynthesis of Glycyrrhiza uralensis Fish (Zhong et al., 2022). Therefore, to further elucidate the relation between differential microbial community, differential secondary metabolites and medicinal components, two interactive networks were constructed ( Figures 6 , 7 ). The significance test and correlation of differential microbial community, differential secondary metabolites and medicinal components were showed in Figures S3 , S4 .

There showed a complex correlation between differential microbial community and differential secondary metabolites ( Figure 6 ). For differential secondary metabolites, the flavonoids mulberroside F showed the greatest correlation with the differential bacterial members ( Figure 6 ), and the phenol 4-Tert-Amylphenol, the amines dehydrophytosphingosine, the organic acid genipic acid and pantothenic acid, the organoheterocyclic loliolide are also important nodes that correlate with the differential bacterial members. 4-Tert-Amylphenol, loliolide and mulberroside F were negatively correlated with JGI 0001001.H03, Candidatus Solibacter, Aquicella and Bryobacter, while genipic acid and comuside was positively correlated with Candidatus Solibacter and Bryobacter, genipic acid was positively correlated with Candidatus Solibacter (r > 0.7). Various rhizodeposits may differentially influence the composition of the rhizosphere microbiome composition (Pascale et al., 2020). Recently studies showed that selected secondary including flavonoid, coumarin, benzoxazinoid, phytohormones, and triterpenes affect the succession of rhizosphere microorganisms of host plants. (Pang et al., 2021). In this paper, we classified the differential secondary metabolites of A. sinensis according to its chemical structure, and made a further Venn comparative analysis on the primary, and secondary metabolites of carbohydrates and glycosides, organic acid, flavonoids, phthalides, and other organoheterocyclic compounds, phenol ( Figures 8A, B ). The results showed that carbohydrates and glycosides, organic acid, and phenol were positively correlated with Candidatus Koribacter (r > 0.7), carbohydrates and glycosides, organic acid were positively correlated with Candidatus Solibacter (r > 0.7), flavonoids, phenol, Phthalide, and other Organoheterocyclic compounds were positively correlated with Pedobacter and Bosea (r > 0.7), flavonoids and organic acid were positively correlated with Iamia (r > 0.7), flavonoids and phenol were positively correlated with Thermomonas (r > 0.7).

Flavonoids, phenol, phthalide, and other organoheterocyclic compounds were negatively correlated with JGI 0001001.H03, Candidatus Solibacter, Aquicella (r > 0.7), carbohydrates and glycosides, organic acids were negatively correlated with Clostridium sensu stricto 12, flavonoids and organic acid were negatively correlated with Burkholderia.Caballeronia.Paraburkholderia (r > 0.7), flavonoids and phenol were negatively correlated with Ralstonia (r > 0.7). From the above results, we can see that flavonoids, phthalides, and other organoheterocyclic compounds have consistent correlation with specific bacterial communities; carbohydrates and Glycosides, and organic acids have consistent correlation with the other specific bacterial communities, which have opposite trend with the former. In short, flavonoids, organic acids, carbohydrates and glycosides, phthalide, and other organoheterocyclic compounds are strongly related to microbes, which may be the main driving factors of microbial community assembly in the rhizosphere of A. sinensis.

There showed a complex correlation between differential microbial community and medicinal components ( Figure 7 ). Chlorogenic acid, senkyunolide A, and n-butylidenephthalide showed the greatest correlation with rhizosphere bacterial communities among 8 medicinal components. Chlorogenic acid was positively correlated with Rhodobacter, JGI 0001001.H03 (r > 0.7), while negatively correlated with Aridibacter, Arthrobacter, Cellulomonas, Luteolibacter (r > 0.7). Senkyunolide A was positively Kineosporia, Arthrobacter, Devosia, Sphingorhabdus, and Brevundimonas (r > 0.7). n-butylidenephthalide was positively with Arthrobacter, while negatively correlated with Alkanindiges, JGI 0001001.H03, Chthonomonas (r > 0.7). Interestingly, the Venn comparative analysis showed that Arthrobacter was positively correlated with senkyunolide A and n-butylidenephthalide, while negatively correlated with chlorogenic acid (r > 0.7) ( Figure 8C ). And JGI 0001001.H03 was positively correlated with chlorogenic acid, while negatively correlated with n-butylidenephthalide (r > 0.7) ( Figure 8C ). In the previous analysis, the level of chlorogenic acid in group G was significantly higher than in group M (p < 0.001), while the level of senkyunolide A was significantly lower than in group M (p < 0.001). And the content of n-butylidenephthalide in group M was also higher than in group G. Therefore, Arthrobacter and JGI 0001001.H03 may be the key microorganism related to the quality of A. sinensis. Arthrobacter is known for its nutritionally versatile nature and wade prevalence in stressful environments, with efficient survivability under high attitude stress conditions (Mukhia et al., 2021). And the strains of Arthrobacter also have the ability to degrade Dibutyl phthalate (DBP), dimethyl phthalate (DMP), 4-chlorophenol efficiently (Kim et al., 2008; Wang et al., 2019; Liu et al., 2020). Therefore, we speculate that Arthrobacter may be able to provide energy to A. sinensis decomposing harmful substances in high altitude and cold environment, and promote the synthesis of medicinal components such as n-butylidenephthalide and senkyunolide A in the early growth of A. sinensis.