Changes of Key Soil Factors, Biochemistry and Bacterial Species Composition during Seasons in the Rhizosphere and Roots of Codonopsis pilosula (tangshen)

: Codonopsis pilosula is a medicinal and edible herb with a rich nutritional value. In Gansu Province, China, its production quality and yield differ during the four seasons. Here, we investigated the differences in the microﬂoral composition and metabolic functions in the rhizospheric soil and roots of C. pilosula during the four seasons, and we also analyzed their dynamic and synergistic effects on C. pilosula growth and carbohydrate content change. The C. pilosula samples were analyzed for plant physiology, microﬂoral composition and metabolic functions in the rhizospheric soil and roots using high-throughput sequencing technology. Environmental indices including soil physiochemistry and meteorological conditions were also determined by the coupled chromatography–spectroscopy technique. The results revealed that the C. pilosula growth was affected by temperature, precipitation and light intensity, with the bacterial structures and functions of the soil and root samples showing obvious seasonal changes. Due to the diversity of microbial composition and community metabolic function, and the synergistic effect of microbial and environmental factors, there are signiﬁcant differences in stress resistance, physiological status and metabolites of C. pilosula in different seasons. Furthermore, the change in seasons was signiﬁcantly correlated with the quality and yield of C. pilosula . This study provides a scientiﬁc basis for soil improvement and the reﬁnement of local Radix C. pilosula cultivation methods.


Introduction
Codonopsis pilosula, which belongs to Campanulaceae [1], is a perennial herb with rich medicinal value [2]. It contains proteins, al-kyne glycosides, polysaccharides, amino acids, minerals and other active ingredients, which are beneficial to the human body [3][4][5]. As a traditional Chinese herbal medicine, C. pilosula has been planted in many countries and regions, mainly in Central Asia, East Asia and Southeast Asia [6,7]. Gansu Province, China, is one of the main producing areas of C. pilosula [8]. In recent years, due to the deterioration of the ecological environment and the changes in local soil texture, temperature, rainfall and other environmental conditions, the quality and yield of C. pilosula have been affected to some extent [9]. Hence, how to improve the quality and yield of C. pilosula and how to further explore the effects of biotic or abiotic stress on plant rhizosphere microecology, plant growth, development and accumulation of metabolites has become a hot research topic at present. Many scholars have devoted themselves to research on the environmental conditions and microorganisms affecting growth and development of C. pilosula [10,11], and some main pathogenic bacteria for C. pilosula root were obtained and identified [11]; meanwhile, beneficial soil microbes promoting seed germination, plant growth and photosynthesis of C. pilosula have also been obtained [12]. Internal and external factors affecting the growth of C. pilosula have been independently explored by many different research groups. Zhou et al. [13] analyzed the changes in soil conditions of C. pilosula under continuous cropping and found that continuous cropping would reduce the main nutrients such as N, P and K in the soil environment where C. pilosula grows, thus hindering the growth and development of C. pilosula. Wojciech [14] and Xie et al. [15] evaluated and improved the planting technology of C. pilosula, which provided a scientific and theoretical basis for improving the quality and yield of C. pilosula. Additionally, it was proven that the change in the growth of C. pilosula could be affected by various internal and external factors. So far, few studies have focused on seasonally driven changes in soil physicochemical properties and changes in microbial community structure in specific C. pilosula cultivation environments. There are no relevant reports on the impact of seasonally driven key factors on the quality and yield of C. pilosula [16].
This study mainly focused on C. pilosula collected in Huichuan area as a specific object. The characteristic changes in rhizospheric and non-rhizospheric physicochemical factors during the four seasons were determined and evaluated based on high-throughput sequencing technology and bioinformatics analysis methods. Furthermore, we analyzed the changing characteristics of the community structures and metabolic functions of rhizosphere soil bacteria and root endophytic bacteria as well. The purpose of this work was: (i) to explore the dynamic relationship among the change in the rhizospheric soil bacterial community, the change in the endophytic bacterial community and the internal indicators of C. pilosula in different seasons. (ii) to analyze the effects of key environmental factors in different seasons on the growth and development of C. pilosula. Ultimately, this study provides a scientific and theoretical basis for improving the soil fertility of C. pilosula in different periods and optimizing the planting methods and management measures of C. pilosula.

Sampling Site Information and Sample Collection
The sample plot in this study is Huichuan Town, Gansu Province (103 • 57 E-104 • 1 E; 35 • 6 N-35 • 9 N), which is one of the main producing areas of C. pilosula in China. C. pilosula was planted in sandy soil and grew slowly. It emerged from March to April and entered a period of rapid growth from mid-June to mid-October [17]. The aboveground part of C. pilosula withered in late October and then entered a dormant period with a daily average temperature of 18-20 • C [18]. Spring samples were collected on 3 April 2017, summer samples on 30 June 2017, autumn samples on 3 September 2017, and winter samples on 10 December 2017. The 5-point sampling method was employed to collect the samples: first determine the midpoint of the diagonal as the central sampling point, and then select four points on the diagonal which are equal to the distance from the central sample point as the sample point. The samples that included rhizosphere soil, bulk soil and roots of C. pilosula were collected in spring (S), summer (X), autumn (A) and winter (W), and each sample included three replicates (a, b and c). The collected samples were packaged and placed in sterile, sealed bags and, finally, stored in a refrigerator at −4 • C for later use. The locations of the sampling sites are shown in Figure S1.

Determination of Physicochemical Indicators of Soil and Plants
The soil physicochemical indicators measured in this study include: (i) Chemical properties: soil moisture content (WC) was measured by conventional drying and weighing methods, and electrical conductivity (EC) and pH were measured by electrode method [19,20]. (ii) Nutrient content: total nitrogen (TN) was determined by using the Kjeldahl method [21]. Available nitrogen (AN) was determined by the alkaline solution diffusion method [22]. Total phosphorus (TP) was determined by the molybdenum antimony anti-colorimetric method [23]. Determination of available phosphorus (AP) was conducted by the NaHCO 3molybdenum antimony anti-colorimetric method [24]. Total potassium (TK) was determined by the sodium hydroxide fusion method [25]. Determination of available potassium (AK) was carried out by NH 4 Ac-flame spectrometry [26]. The organic matter (OM) content was determined using the potassium dichromate method [27]. (iii) Soil enzyme activity: 5 g of soil was added to a 50 mL triangle bottle with 1mL of methylbenzene solution. After 15 min, 10 mL of 10% urea and 20 mL of pH 6.7 citrate buffer were added into the above mixture, shaken and incubated at 37 • C for 24 h. Then, the reaction solution was collected, and the urease activity was determined by the indophenol blue colorimetric method [28]. A total of 2 g of soil was added to a 50 mL triangle bottle with 15 mL of 8% sucrose solution and 5 mL pH 5.5 phosphate buffer, shaken and incubated at 37 • C for 24 h. Then, the reaction solution was collected and sucrase (SUC) was determined using the 3,5-dinitrosalicylic acid colorimetric method [29]. A total of 2 g of soil was added to a 100 mL triangle bottle with 40 mL distilled water and 5 mL 0.3% H 2 O 2 solution, sealing reacted for 20 h. Then, 1 mL of saturated aluminum-potassium-vanadium solution was added into the above solution, and catalase (CAT) was determined by potassium permanganate titration [30]. A total of 1.0 g of soil was added to a 10 mL protein extract solution, extracted for 15 min at 4 • C. Then, 1 mL extraction solution was obatained, and protease (PRO) was determined using a modified ninhydrin colorimetric method [31].

Identification of Rhizosphere Soil Bacteria and Root Endophytic Bacteria
Genomic DNA was extracted from the soil samples stored at −80 • C according to the MOBIO PowerSoil ® DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA USA), and the V3-V4 region of 16S rDNA was amplified with specific primers. The primer sequences were (5 -3 ): 341F: CCTACGGGNGGCWGCAG and 806R: 2GGACTACHVGGGTATC-TAAT. PCR was conducted with a designed amplification system (Table S1) under the amplification conditions as shown in Table S2. Then the PCR amplification products were recovered, and the products were quantified with Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The purified products were mixed in equal amounts, the sequencing adapters were ligated and the sequin library was constructed. Finally, the library construction and sequencing steps were performed by Gene Denovo Biotechnology Co., Ltd., Guangzhou, China. The Miseq PE250 platform of the Illumina company was used for sequencing.
In this study, the raw data were quality control filtered and then spliced and assembled to obtain valid sequences by FLASH (Fast Length Adjustment of Short, v 1.2.11), QIIME (Quantitative Insights Into Microbial Ecology, v 1.8.0) and the UCHIME algorithm (http://www.drive5.com/usearch/manual/uchime_algo.html, accessed on 17 July 2018), respectively. The microbial ecology pipeline software UPARSE (v 9.2.64_i86linux32) was used to cluster the 16S rRNA sequences at 97% similarity to different operational taxonomic units (OTUs). In this study, the above contents were analyzed by diversity analysis, community structure analysis, community function prediction and environmental factor correlation analysis.

Statistical Analyses
All data are expressed as the average value ± the standard deviation (SD) for three replicate treatments. Raw data was summarized and counted in Excel 2010. The significant difference of soil and plant physicochemical indicators among each group was performed by one-way analysis of variance (ANOVA) using SPSS 17.0. Within sample type, α-diversity estimates were calculated by analyzing the Simpson diversity and Shannon diversity in the phyloseq package, as implemented in the tool MicrobiomeAnalyst [42]. The "vegan", "ggplot2" and "ggcor" packages of R 4.2 were used for Pearson and Spearman correlation analysis, and β diversity analysis was based on the Bray-Curtis distance measure [43]. The "vegan" package was used to calculate the number of microorganisms and abundance based on the 16S OTU table. For all samples, we used total sum scaling to calculate the relative abundance and expressed the relative abundance as percentages [44]. The redundancy analysis (RDA) with 999 permutations was carried out with the "vegan" package to analyze the effects of environmental factors on microbial composition [45]. PICRUSt was used to predict the functional diversity of bacterial communities based on the OTU [46]. In addition, structural equation modeling (SEM) was used to evaluate the direct and indirect effects of prominent factors: root microbial diversity and composition, meteorological, plant denote and plant indicators [47]. The SEM was constructed using the "plspm" package in R.

Soil Physicochemical Properties of C. pilosula in Different Seasons
The changes in various indicators of rhizosphere soil and non-rhizosphere soil were compared ( Figure S2 and Table 1 and Table S3). The PRO enzyme activity level in the rhizosphere soil decreased, whereas the enzyme activity level in the non-rhizosphere soil was the highest in autumn (402.51 µg/g/d). The WC and PRO contents in different soils varied significantly, whereas other indices differed, but not significantly. The WC, pH and EC of the soil in different seasons were significantly different (p < 0.05). The contents of AN, TP, AP and OM in the soil showed a significant trend of polarization; that is, the contents in spring and summer were similar, and the contents in autumn and winter were similar (p < 0.001). The content of TN in autumn was the highest, while the contents of TK and AK in spring were significantly higher than those in other seasons. The SUC activity in soil increased significantly with the seasons, with the largest difference in activity occurring between summer and autumn. The activity of URE enzyme decreased significantly with the change in seasons. Based on the above detection results, three environmental indicators reflecting local seasonal changes (monthly average temperature, precipitation and illumination) were used to conduct a correlation analysis of physicochemical factors (Figures in Section 3.2 and Figure S3). The results of the analysis showed that there were significant positive and negative correlations between the physicochemical factors in the soil of different seasons (Figure 1a and (ii) URE, TP and EC were extremely significantly positively correlated (p < 0.001). In addition, there was a significant negative correlation between the two groups (p < 0.01).
Based on the above detection results, three environmental indicators reflecting local seasonal changes (monthly average temperature, precipitation and illumination) were used to conduct a correlation analysis of physicochemical factors (Figures in Section 3.2 and Figure S3). The results of the analysis showed that there were significant positive and negative correlations between the physicochemical factors in the soil of different seasons (Figure 1a,c). The correlation results were mainly divided into two categories: (ⅰ) SUC, CAT, PRO, OM, AP, AN and TN were extremely significantly positively correlated (p < 0.001); and (ⅱ) URE, TP and EC were extremely significantly positively correlated (p < 0.001). In addition, there was a significant negative correlation between the two groups (p < 0.01).  The analysis also revealed that physicochemical factors, such as WC, TK and AK, in the soil were significantly correlated with temperature, precipitation and light (p < 0.05; Figure 1b,d). Combined with the correlation analysis of physicochemical factors and environmental factors, this study hypothesized that soil enzyme activity would also be affected by environmental conditions, but not significantly. As shown in Figure S3, there was no significant correlation between the activity of the PRO enzyme in rhizosphere and non-rhizosphere soils, which indicates that the enzyme's activity is greatly changed during the regulation of C. pilosula. Thus, on the one hand, as a result of the growth regulation of C. pilosula, the physicochemical indicators of the rhizosphere soil changed compared with those of the nonrhizosphere soil, with the changes in PRO and WC being the most obvious. On the other hand, some physicochemical factors, such as WC, TK and AK, are also significantly related to environmental factors. Thus, more physicochemical factors may be significantly correlated with the metabolic regulation between soil and plants.

The Characteristic Changes in Rhizosphere Soil Bacteria and Endophytic Bacteria in C. pilosula during Different Seasons
The α diversity analysis was performed on the sequencing results (Tables S4 and S5), as shown in Figure 2a,b, there were significant differences among the SHCA, SHCS and SHCX in the Shannon diversity index of rhizosphere soil OTUs (p < 0.05), and there were no significant differences among the four groups in the Simpson diversity index of rhizosphere soil OTUs. There were significant differences in the Shannon and Simpson indexes (p < 0.05) of root endogenous OTUs (Figure 2c,d) in different seasons (p < 0.001). The results showed that the abundance and diversity of microorganisms in the rhizosphere soil were the lowest in the summer, and the abundance and diversity of endophytic microorganisms were the lowest in the winter. Simpson index of rhizosphere soil (root endogenous) OTUs. The overall analysis p was asse using an ANOVA, and the analysis between groups was assessed using a t-test. Significant d ence (p < 0.05).
The genus−level Venn diagram ( Figure S4) shows that the rhizosphere soil OTU the spring, summer, autumn and winter were 238, 251, 279 and 276, respectively. Th versity and richness index of the OTUs was the lowest in spring and the highest in tumn. The combined local average temperature and precipitation values were the low in spring (−3.27 °C, 7.03 mm) and the highest in autumn (12.74 °C, 140.62 mm), w indicates that temperature may be an important factor that directly or indirectly aff the significant changes in the OTUs. There was a great correlation between root end The genus−level Venn diagram ( Figure S4) shows that the rhizosphere soil OTUs in the spring, summer, autumn and winter were 238, 251, 279 and 276, respectively. The diversity and richness index of the OTUs was the lowest in spring and the highest in autumn. The combined local average temperature and precipitation values were the lowest in spring (−3.27 • C, 7.03 mm) and the highest in autumn (12.74 • C, 140.62 mm), which indicates that temperature may be an important factor that directly or indirectly affects the significant changes in the OTUs. There was a great correlation between root endogenous OUT's diversity index and seasonal variation. Specifically, it was the smallest in summer and the largest in autumn. It showed that many bacteria forming endogenous OTUs in the roots are more obviously restricted by conditions such as temperature. The endogenous OTUs formed in each of the four seasons (spring, summer, autumn and winter) were 79, 61, 73 and 47, respectively. The diversity and richness of OTUs in the rhizosphere soil were significantly greater than in the roots.
The PCoA results of a β diversity analysis (Figure 3a) revealed that the OTUs of the rhizosphere soil clustered into one category in spring and summer, and clustered into another category in autumn and winter, which indicated that the differences in community structure varied depending on the seasons and were extremely significant (p < 0.01). The PCo1 of OTUs in rhizosphere soil accounted for 53.43%, whereas PCo2 accounted for 14.01%. In addition, as indicated in Figure 3b, there were significant differences in the community structure of the endogenous OTUs in summer and autumn (p < 0.01). Here, PCo1 explained 80.04% of the difference, and PCo2 explained 12.15%. Thus, the above results show that season may be significantly correlated with changes in the species community structure, and random changes in species within a group had little effects on environmental heterogeneity. Additionally, combined with the within group R 2 and p-values, this showed that species in the rhizosphere soil were more correlated with seasons than root endophytes. PCo1 explained 80.04% of the difference, and PCo2 explained 12.15%. Thus, the above results show that season may be significantly correlated with changes in the species community structure, and random changes in species within a group had little effects on environmental heterogeneity. Additionally, combined with the within group R 2 and p-values, this showed that species in the rhizosphere soil were more correlated with seasons than root endophytes. Based on NMDS analysis ( Figure S5), the Shepard fitting results for rhizosphere soil and root endophyte samples showed that the R 2 indices were 0.993 ( Figure S5a) and 0.994 ( Figure S5b), respectively. The two R 2 indices from different samples were relatively close, indicating that the fitting evaluation results were good and could be used for the analysis. The stress analysis results of rhizosphere soil and root endophytic samples in different seasons were 0.084 and 0.078, respectively, indicating that the samples had a high level of reliability. The MDS1 distribution revealed that the samples clustered, indicating that there were obvious differences in the structures of the bacterial communities during different seasons. However, the distribution of MDS2 was scattered and there was no obvious clustering trend, indicating that the time dimension was the main factor leading to the differences in the bacterial community structure of C. pilosula. Based on NMDS analysis ( Figure S5), the Shepard fitting results for rhizosphere soil and root endophyte samples showed that the R 2 indices were 0.993 ( Figure S5a) and 0.994 ( Figure S5b), respectively. The two R 2 indices from different samples were relatively close, indicating that the fitting evaluation results were good and could be used for the analysis. The stress analysis results of rhizosphere soil and root endophytic samples in different seasons were 0.084 and 0.078, respectively, indicating that the samples had a high level of reliability. The MDS1 distribution revealed that the samples clustered, indicating that there were obvious differences in the structures of the bacterial communities during different seasons. However, the distribution of MDS2 was scattered and there was no obvious clustering trend, indicating that the time dimension was the main factor leading to the differences in the bacterial community structure of C. pilosula.

Composition and Correlation of Dominant Bacteria in Rhizosphere Soil and Endophytes of C. pilosula in Different Seasons
The species composition analysis showed that Exiguobacterium was the dominant bacterium in both rhizosphere soil and endophytic samples, and their abundance was the highest in autumn and winter (Figure 4a,b and Figure S6), accounting for about 30% of the total species abundance. The abundance of RB41 [48] in the rhizosphere soil decreased significantly seasonally (p < 0.05), with the largest decrease being 2.99% from spring to summer. The Nitrospira abundance increased significantly seasonally, with the highest increase of 0.54% occurring from summer to autumn. Pseudomonas, Lactococcus, Anoxybacillus and Thermus among the root endophytes significantly decreased by season (p < 0.05). Among them, Pseudomonas had the greatest abundance change, decreasing from 26.3% in summer to 6.3% in autumn. In addition, Acinetobacter, Bacillus and Enterococcus increased significantly with the change in season (p < 0.05). We also analyzed the correlations between dominant bacteria and physicochemical indicators of rhizosphere soil. As shown in Figure 4c,d, the dominant bacteria Exiguobacterium were significantly positively correlated with SUC, CAT, PRO, OM, AP, AN and TN and significantly negatively correlated with URE, TP and EC (p < 0.05). In addition, other dominant bacteria, including Nitrospira, Acinetobacter and Bacillus, showed similar correlations. The correlations between Lysobacter, Pseudomonas, Anoxybacillus and physicochemical indicators were opposite to that of Exiguobacterium. The above results indicated that there was a significant correlation between the composition of the bacterial community in different seasons and the physicochemical indicators of the rhizosphere soil. The changes in bacterial community composition and physical and chemical indexes of the soil changed with the change in season, and this change will affect the growth of C. pilosula in different periods.

Predictions of Bacterial Functions in Different Samples
The community functions of bacteria were predicted using PICRUSt. The prediction results showed that chemoheterotrophy and aerobic chemoheterotrophy were the main functions of bacteria in rhizosphere soil samples (p < 0.01; Figure 5a). The proportions of these functional activities in the summer showed the largest changes, increasing by 3.46% and 3.52%, respectively, compared with in the spring. Additionally, the community functions of soil bacteria were more complex, and the changes during the four seasons were more stable (p > 0.05). The main functions of root endophytic bacteria were invertebrate parasites, chemoheterotrophy and aerobic chemoheterotrophy (Figure 5b). Among them, the abundance of invertebrate parasite function was significantly different in different seasons (p < 0.001). The abundance of Chemoheterotrophy and aerobic chemoheterotrophy We also analyzed the correlations between dominant bacteria and physicochemical indicators of rhizosphere soil. As shown in Figure 4c,d, the dominant bacteria Exiguobacterium were significantly positively correlated with SUC, CAT, PRO, OM, AP, AN and TN and significantly negatively correlated with URE, TP and EC (p < 0.05). In addition, other dominant bacteria, including Nitrospira, Acinetobacter and Bacillus, showed similar correlations. The correlations between Lysobacter, Pseudomonas, Anoxybacillus and physicochemical indicators were opposite to that of Exiguobacterium. The above results indicated that there was a Agronomy 2023, 13, 1545 9 of 17 significant correlation between the composition of the bacterial community in different seasons and the physicochemical indicators of the rhizosphere soil. The changes in bacterial community composition and physical and chemical indexes of the soil changed with the change in season, and this change will affect the growth of C. pilosula in different periods.

Predictions of Bacterial Functions in Different Samples
The community functions of bacteria were predicted using PICRUSt. The prediction results showed that chemoheterotrophy and aerobic chemoheterotrophy were the main functions of bacteria in rhizosphere soil samples (p < 0.01; Figure 5a). The proportions of these functional activities in the summer showed the largest changes, increasing by 3.46% and 3.52%, respectively, compared with in the spring. Additionally, the community functions of soil bacteria were more complex, and the changes during the four seasons were more stable (p > 0.05). The main functions of root endophytic bacteria were invertebrate parasites, chemoheterotrophy and aerobic chemoheterotrophy (Figure 5b). Among them, the abundance of invertebrate parasite function was significantly different in different seasons (p < 0.001). The abundance of Chemoheterotrophy and aerobic chemoheterotrophy were also significantly different among root endophytic microorganisms in different seasons. The metabolic functions of chemoheterotrophy and aerobic chemoheterotrophy functions in the summer were 10.1% and 10.2%, lower than the two functions in the spring, respectively.  Therefore, the community function of bacteria in soil was complex and relatively balanced, and the abundance of each function did not change significantly among the different seasons. The community functions of endophytic bacteria showed opposite characteristics, and the abundance of different functions varied greatly. In addition, the abundance of each function varied significantly in response to seasonal changes, and the main differences were concentrated between summer and autumn.

Correlaton of Bacterial Community Structure and Function with the Growth Indicators of C. pilosula
This paper also studied whether the three growth indicators of C. pilosula, carbohydrate content (sugar content), stress resistance enzymes and root phenotype were driven by seasonal changes. As shown in Tables 2, S6 and S7, through a one−way analysis of variance, the activities of most of the stress-resistance-related enzymes, such as SPS, SOD and RCAT, did not differ significantly (p > 0.05), but the enzyme activities during the different seasons changed significantly. The contents of sugars, such as ReS and Fru, were significantly higher in spring and summer than in autumn and winter. The root lengths and root diameters of C. pilosula were greater in autumn and winter than in spring and summer. The root length increased the fastest in spring and summer, and the root diameter increased the fastest in autumn and winter. We speculate that the roots of C. pilosula mainly grow Therefore, the community function of bacteria in soil was complex and relatively balanced, and the abundance of each function did not change significantly among the different seasons. The community functions of endophytic bacteria showed opposite characteristics, and the abundance of different functions varied greatly. In addition, the abundance of each function varied significantly in response to seasonal changes, and the main differences were concentrated between summer and autumn.

Correlaton of Bacterial Community Structure and Function with the Growth Indicators of C. pilosula
This paper also studied whether the three growth indicators of C. pilosula, carbohydrate content (sugar content), stress resistance enzymes and root phenotype were driven by seasonal changes. As shown in Tables 2 and S6 and S7, through a one−way analysis of variance, the activities of most of the stress-resistance-related enzymes, such as SPS, SOD and RCAT, did not differ significantly (p > 0.05), but the enzyme activities during the different seasons changed significantly. The contents of sugars, such as ReS and Fru, were significantly higher in spring and summer than in autumn and winter. The root lengths and root diameters of C. pilosula were greater in autumn and winter than in spring and summer. The root length increased the fastest in spring and summer, and the root diameter increased the fastest in autumn and winter. We speculate that the roots of C. pilosula mainly grow vertically in spring and summer and grow horizontally in autumn and winter. Consequently, a Mantel correlation test among the bacterial community structure, bacterial community function and C. pilosula growth indicators was performed ( Figure 6 and Figure S6). The results showed that there was a significant positive correlation between RB41, Nitrospira, Exigucbacterium and Resilience index, and between RB41 and plant phenotype in rhizosphere soil. In the root, Exigucbacterium, Acinetobacter, Bacillus, Enierococcus and Annaxybacillus were significantly positively correlated with the resilience index, and Pseudomonas was significantly positively correlated with plant phenotype. A combined functional analysis indicated ( Figure S7) that these bacteria may simultaneously participate in multiple metabolic pathways affecting plant growth regulation. The comparison, shown in Figure 6a,b, indicates that the influence of endophytic bacteria in the root system on each growth index of C. pilosula is more obvious. This is because the endophytic bacteria act more directly on metabolism and other pathways of the C. pilosula root system, and thus, they participate in the growth regulation of C. pilosula more frequently. Consequently, a Mantel correlation test among the bacterial community structure, bacterial community function and C. pilosula growth indicators was performed (Figures 6  and S6). The results showed that there was a significant positive correlation between RB41, Nitrospira, Exigucbacterium and Resilience index, and between RB41 and plant phenotype in rhizosphere soil. In the root, Exigucbacterium, Acinetobacter, Bacillus, Enierococcus and Annaxybacillus were significantly positively correlated with the resilience index, and Pseudomonas was significantly positively correlated with plant phenotype. A combined functional analysis indicated ( Figure S7) that these bacteria may simultaneously participate in multiple metabolic pathways affecting plant growth regulation. The comparison, shown in Figure 6a,b, indicates that the influence of endophytic bacteria in the root system on each growth index of C. pilosula is more obvious. This is because the endophytic bacteria act more directly on metabolism and other pathways of the C. pilosula root system, and thus, they participate in the growth regulation of C. pilosula more frequently.

Correlation of Key Physicochemical Factors Associated with Seasonal Variation with the Growth of C. pilosula
Owing to the variation among components, an RDA analysis between the rhizosphere soil physicochemical indicators and bacteria was conducted, and an SEM model of the key factors affecting the growth of C. pilosula was constructed. The correlation between

Correlation of Key Physicochemical Factors Associated with Seasonal Variation with the Growth of C. pilosula
Owing to the variation among components, an RDA analysis between the rhizosphere soil physicochemical indicators and bacteria was conducted, and an SEM model of the key factors affecting the growth of C. pilosula was constructed. The correlation between the key environmental factors related to seasonal change and the growth regulation of C. pilosula was analyzed. The RDA results showed that 51.49% of the variation in the bacterial community structure in rhizosphere soils could be explained by two canonical axes, with RDA1 and RDA2 accounting for 42.2% and 9.28% of the variation, respectively (Figure 7a). Additionally, 93.34% of the variation in the endophytic bacterial community structure in roots could be explained by RDA, with RDA1 and RDA2 accounting for 92.63% and 0.71% of the variation, respectively (Figure 7b). It also showed that WC, EC, pH, TN, AN, TP and AP were the main factors affecting the bacterial community structure.
3, x FOR PEER REVIEW 12 of 18 7a). Additionally, 93.34% of the variation in the endophytic bacterial community structure in roots could be explained by RDA, with RDA1 and RDA2 accounting for 92.63% and 0.71% of the variation, respectively (Figure 7b). It also showed that WC, EC, pH, TN, AN, TP and AP were the main factors affecting the bacterial community structure.  The SEM results showed that there was a significant negative correlation between meteorological conditions and physicochemical factors of rhizosphere soil, while the latter was positively correlated with physicochemical factors of non-rhizosphere soil, rhizosphere soil bacteria and plant indicators (p < 0.05). Plant growth was obviously promoted by meteorological conditions, rhizosphere soil bacteria and endophytic bacteria (p < 0.05; Figure 7c). As shown in Figure 7d, rhizosphere soil physicochemical factors, rhizosphere soil bacteria and endophytic bacteria were the top three components with the most influence on plant growth.

Discussion
In this study, the physicochemical indexes of rhizosphere soil and non-rhizosphere soil showed the same trend in different seasons, but the overall content was significantly different. ( Figure S2). Similar differences between rhizosphere and non-rhizosphere soils were also shown in previous studies [49,50]. As shown in Figure 1c,d, the WC in the rhizosphere soil was positively correlated with air temperature, precipitation and illumination. Soil moisture may be an important factor affecting the diversity and composition of microorganisms in soil. The higher the soil moisture, the more beneficial it is to the increase in microbial diversity [51]. Temperature, precipitation and light were negatively correlated with TK and AK. The changes in humidity and heat may be the important factors affecting the transformation of potassium. Drought and excessive water will reduce the transformation rate of available potassium [52]. Combined with the changes in soil AK and TK contents, this study shows that the physicochemical properties of soil are restricted by external hot and humid conditions. In addition, due to the lack of tillage or fertilization for a long time, the level of potassium in the soil was insufficient [53,54]. TN is significantly positively correlated with temperature (p < 0.05). The highest average temperature in autumn is 12.7 • C, the highest precipitation is 160.42 mm, and the content of TN is also the highest. Mohamed et al. [55] showed that suitable temperature and higher precipitation are conducive to the increases in TN. The pH value of rhizosphere soil was positively correlated with light, and the soil pH was the highest in summer (p < 0.01). Sheng et al. [56] and Adjout et al. [57] showed that the pH increases along with the illumination. Additionally, it has been shown that EC, pH, AN and OM are main factors that affect soil nitrogen mineralization and nitrification [58][59][60], and this study also showed that they were significantly positively correlated. Based on the consistent changes in AN and TN in the four seasons, this study speculated that the rapid absorption and utilization of TN by plants and the insufficient abundance of nitrifying bacteria in soil led to insufficient transformation of AN [61]. According to Thami et al. [62] and Greenfield et al. [63] the deterioration of the rhizosphere soil environment in autumn may lead to the decrease in PRO activity. Thus, various rhizosphere soil physicochemical indicators are clearly affected by the seasons [64].
The rhizosphere soil, as the main location for material and energy exchange between plants and the outside world, is the area with the most frequent microbial activities [65]. The results of diversity showed that the Shannon index and Simpson index of rhizosphere soil microorganisms increased in different seasons, but the difference was not significant (p > 0.05). In addition, the diversity index of root endophytic bacteria increased significantly (p < 0.01; Figure 2 and Tables S4 and S6). The richness and diversity of endophytic bacterial communities were much lower than those in rhizosphere soil, which was consistent with the results of Hereira et al. [66]. The distribution of the bacterial community structure in rhizosphere soil of C. pilosula was balanced (Figure 4), while the community structure of endophytic bacteria changed greatly (p > 0.05). The abundance of dominant bacteria Pseudomonas decreased gradually with the seasonal variation, while the abundance of Exiguobacterium increased with the seasonal variation (p < 0.001), which indicated that bacteria are driven by seasonal changes [67]. The bacterial community structure is closely related to soil physicochemical indicators [68]. In this study, Nitrospira, Exiguobacterium and Pseudomonas were significantly correlated with TN and AN. Therefore, it is speculated that the changes in the bacterial community structure affect the levels of mineral elements necessary for the growth of C. pilosula, such as nitrogen, phosphorus and potassium [69].
The change in growth index showed that the content of MDA in C. pilosula in summer was three times higher than that in other seasons (Table S6). The content of Pro in C. pilosula was the highest in autumn and the lowest in winter. This shows that the stress resistance of C. pilosula is low in the summer and enhanced in the autumn [70][71][72]. Zhang et al. [73] showed that MDA has antibacterial properties against bacteria such as Staphylococcus and Lactobacillus. Bacteria, such as Exiguobacterium and Nitrospira, play an important role in improving stress resistance and biodegradation in plants [74][75][76]. Therefore, the change in the bacterial community structure driven by seasonal change is closely related to the stress resistance of plants, which is consistent with the conclusion of this study ( Figure 6). Previous studies have shown that changes in the microbial community structure can also affect the function of specific bacteria, such as nitrification, degradation and other metabolic functions of bacteria, which is consistent with the changes in the bacterial community structure [77]. The results of this study showed that the change in microbial community in rhizosphere soil was limited, but the community function of root endophytic bacteria was very different ( Figure 5). The functional changes in these bacteria also affect the growth index of C. pilosula (Figure 7d). Rhizosphere soil bacteria and root endophytic bacteria were greatly affected by seasonal changes, and the latter was more significantly affected. Therefore, the change in bacterial community structure will affect the physicochemical properties of soil, plant resistance and the function of bacteria, and finally, will affect the growth and development of C. pilosula.

Conclusions
The analysis of non-rhizosphere soil, rhizosphere soil and root samples showed that there was a significant correlation between season and soil physical and chemical indexes of C. pilosula, which may further affect the growth, stress resistance and metabolism of C. pilosula. The roots grow mainly vertically in spring and summer and horizontally in autumn and winter. This shows that there is a significant correlation between the growth of spring sowing plants and seasonal changes, and rhizosphere bacteria and root endophytic bacteria play an important role in nutrient supply. This study clarified the correlation between seasonal variation and core bacteria and specific populations of C. pilosula and provided a theoretical basis for further exploration of the quality and yield of medicinal plants improved by microorganisms.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy13061545/s1, Figure S1: The location information of the sampling sites; Figure S2: Seasonal changes in rhizosphere soil and non-rhizosphere soil physicochemical factors; Figure S3: The meteorological conditions of the sample site vary throughout the four seasons; Figure S4: OTUs cluster analysis of C. pilosula at the genus level; Figure S5: The changing of NMDS diversity at the genus level in rhizosphere soil and root endophytic samples; Figure S6: Genus level TOP10 bacteria heat map; Figure S7: Mantel correlation between bacterial function and growth index of C. pilosula; Table S1: PCR amplification system (50 µL system); Table S2: PCR amplification conditions; Table S3: Changes of non-rhizosphere soil physical and chemical factors in different seasons; Table S4: α-diversity index of rhizosphere bacteria; Table S5: α-diversity index of endophytic bacteria in root; Table S6: Resilience indicators of plants in all seasons; Table S7: Phenotypic indicators of plants in all seasons.
Author Contributions: All authors contributed to the study conception and design. Y.W. and J.C. designed and revised the manuscript; T.M. and T.Z. performed the experiments, collected the data and wrote the manuscript; and T.Z., Y.W. and F.L. analyzed the data. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The sequence data information in this study can be found in GenBank of the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/ nuccore, accessed on 13 July 2020), and the access number: PRJNA861447.

Conflicts of Interest:
The authors declare no conflict of interest.