Taxonomic Structure of Rhizosphere Bacterial Communities and Its Association With the Accumulation of Alkaloidal Metabolites in Sophora flavescens

Plant secondary metabolites (SMs) play a crucial role in plant defense against pathogens and adaptation to environmental stresses, some of which are produced from medicinal plants and are the material basis of clinical efficacy and vital indicators for quality evaluation of corresponding medicinal materials. The influence of plant microbiota on plant nutrient uptake, production, and stress tolerance has been revealed, but the associations between plant microbiota and the accumulation of SMs in medicinal plants remain largely unknown. Plant SMs can vary among individuals, which could be partly ascribed to the shift in microbial community associated with the plant host. In the present study, we sampled fine roots and rhizosphere soils of Sophora flavescens grown in four well-separated cities/counties in China and determined the taxonomic composition of rhizosphere bacterial communities using Illumina 16S amplicon sequencing. In addition, the association of the rhizosphere bacterial microbiota with the accumulation of alkaloids in the roots of S. flavescens was analyzed. The results showed that S. flavescens hosted distinct bacterial communities in the rhizosphere across geographic locations and plant ages, also indicating that geographic location was a larger source of variation than plant age. Moreover, redundancy analysis revealed that spatial, climatic (mean annual temperature and precipitation), and edaphic factors (pH and available N and P) were the key drivers that shape the rhizosphere bacterial communities. Furthermore, the results of the Mantel test demonstrated that the rhizosphere bacterial microbiota was remarkably correlated with the contents of oxymatrine, sophoridine, and matrine + oxymatrine in roots. Specific taxa belonging to Actinobacteria and Chloroflexi were identified as potential beneficial bacteria associated with the total accumulation of matrine and oxymatrine by a random forest machine learning algorithm. Finally, the structural equation modeling indicated that the Actinobacteria phylum had a direct effect on the total accumulation of matrine and oxymatrine. The present study addresses the association between the rhizosphere bacterial communities and the accumulation of alkaloids in the medicinal plant S. flavescens. Our findings may provide a basis for the quality improvement and sustainable utilization of this medicinal plant thorough rhizosphere microbiota manipulation.


INTRODUCTION
Plant secondary metabolites (SMs), produced from secondary metabolic pathways in plants, comprise numerous structurally diverse compounds that can be classified into several large families: phenolics, terpenes, steroids, alkaloids, and flavonoids (Pang et al., 2021). Plant SMs play a crucial role in plant defense against pathogens and pests and adaptation to environmental stresses, also acting as symbiotic signals for plants in association with microbes (Pang et al., 2021). In addition, some plant SMs have been used as ingredients or additives of pharmaceuticals and nutraceuticals, making enormous contributions to public health . However, as a material basis of clinical efficacy and vital indicators for quality evaluation of the corresponding medicinal materials, the SMs of medicinal plants can vary among individuals . Genetic, ontogenetic, and environmental factors are recognized as the sources of variation for the intraspecific polymorphisms of plant SMs (Moore et al., 2014).
Terrestrial plants host taxonomically diverse microorganisms, which exert beneficial, harmful, or neutral effects on plant growth and resistance to biotic and abiotic stresses (Müller et al., 2016). Using next-generation sequencing and bioinformatic analysis, the phylogenetic composition and functional diversities of microbial communities associated with various host plants, such as Arabidopsis, several corps, and tree species, have been extensively investigated, revealing associations between plant microbiota and plant nutrient uptake, production, and stress tolerance (Bulgarelli et al., 2012(Bulgarelli et al., , 2015Lundberg et al., 2012;Kembel et al., 2014;Edwards et al., 2015;Jin et al., 2017;Zhang et al., 2018aZhang et al., , 2019Fan et al., 2020;Lee et al., 2021). The community structure and functional capabilities of microbiota associated with medicinal plants have received some research attention recently (Chen et al., 2018;Kui et al., 2021;Liu et al., 2021;Na et al., 2021). It has been proposed that the intraspecific variation in SMs of medicinal plants could be partly ascribed to the shift in the composition of the microbial community associated with the plant host (Huang et al., 2018). However, the associations between plant microbiota and accumulation of SMs in medicinal plants remain largely unknown.
Sophora flavescens, a slow-growing Fabaceae shrub, has been planted as an important medicinal plant in China due to its widespread use in pharmacological drugs (Song et al., 2019). The wide-reaching biological activities of SMs of S. flavescens mainly include antitumor, antiviral, anti-inflammatory, antianaphylactic, antiarrhythmic, and hepatoprotective effects (He et al., 2015). The root of S. flavescens contains extensive bioactive constituents, of which alkaloids are characteristic compounds of the plant and have provoked great interest for their pharmacological effects on cancers, viral hepatitis, and cardiac diseases (He et al., 2015). With the increasing demand for products of medicinal plants and limited arable land, sustainable efforts to improve the production of medicinal plants are required. However, this cannot be achieved by the input of chemical fertilizers and pesticides, which can degrade the quality of medicinal plants (Song et al., 2019). Under these circumstances, the management of favorable microbes is expected to be a potential strategy for facilitating the biosynthesis and accumulation of bioactive components in medicinal plants. Thus, comprehensive studies regarding the association between plant microbiota and bioactive components in medicinal plants are required for the identification of potentially beneficial microbes (Jin et al., 2017).
The plant rhizosphere is a specific microhabitat around the root, where complex soil-microbe-plant interactions occur de Vries et al., 2020). The rhizosphere microbiota, whose genome is referred to as the second genome of plants, contributes to plant phenotype and is vitally important for plant growth and health (Berendsen et al., 2012;Zhang et al., 2017). In addition, the rhizosphere is a hotspot for investigating the interactions between plantassociated microbes and plant SMs (Pang et al., 2021). Therefore, we proposed to address the association between the rhizosphere microbiota and the accumulation of alkaloids in S. flavescens.
In the present research, we sequenced the rhizosphere bacterial communities of 72 S. flavescens individuals grown in four well-separated areas in China. The key drivers shaping the rhizosphere bacterial communities were evaluated. Moreover, the association of the rhizosphere bacterial microbiota with the accumulation of alkaloids in roots of S. flavescens was explored. Potentially beneficial microbes correlated with the total accumulation of matrine and oxymatrine were identified. Finally, the direct and indirect effects of spatial, climatic, and edaphic factors and the preponderant rhizosphere phyla on the accumulation of alkaloids were estimated. We expect that our work will provide a basis for the quality improvement and sustainable utilization of this medicinal plant by rhizosphere microbiota manipulation.

Sample Collection
Rhizosphere and root samples of S. flavescens were collected from 12 fields distributed in four cities/counties in China [Linyuan (LY), Liaoning Province; Changzhi (CZ), Shanxi Province; Dali (DL) and Luonan (LN), Shaanxi Province], which are far-separated with locations between 34 ~ 41° N and 110 ~ 119° E (Supplementary Table S1). The sampling was performed between September and October 2019, roughly corresponding to the period before harvesting. Six individuals from each field were randomly selected as six replicates (n = 6). Taking the influence of growth period on rhizosphere microbiota and alkaloid accumulation into consideration, we sampled 1-3-year-old plants ( Supplementary Table S1). After shaking off soil particles that were loosely attached, firmly adhering soil layers in combination with roots were collected and transported to our laboratory on ice. In the laboratory, the soil particles firmly attached to roots were gently brushed and collected into sterilized tubes as rhizosphere samples. The rhizosphere soil was stored at −80°C until DNA extraction was performed. The roots were washed with distilled water Frontiers in Microbiology | www.frontiersin.org 3 December 2021 | Volume 12 | Article 781316 and oven-dried at 55°C to a constant weight for determination of alkaloid contents.

DNA Extraction, High-Throughput Sequencing, and Sequence Analysis
Total community DNA was extracted from the rhizosphere samples (0.5 g) using a magnetic soil and stool DNA kit (Tiangen Biotech, Beijing, China) following the manufacturer's instructions. The V3-V4 hypervariable regions of the 16S rRNA gene were amplified using the specific primers 341F (5′-CCTAYGGGRBG CASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′; Caporaso et al., 2012) with barcodes. Sequencing libraries were generated using the Illumina TruSeq DNA PCR-free library preparation kit (Illumina, San Diego, United States). The library was sequenced on an Illumina NovaSeq platform at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). Raw reads were processed using the QIIME software package (Caporaso et al., 2010). In detail, the paired-end reads were assigned to each sample based on the unique barcodes, merged, and filtered to obtain high-quality tags. An average of 59,893 qualified reads per sample (ranging from 47,975 to 69,193) was yielded for subsequent analysis. Operational taxonomic units (OTUs) were picked at 97% sequence similarity along with the removal of chimeras (Edgar, 2010;Edgar et al., 2011). Taxonomic assignment was performed against the Silva database (Quast et al., 2013). The nonbacterial OTUs and OTUs whose number of reads was lower than 10 were removed. OTU matrices were rarified to the same depth (42,842 reads per sample). Sequences used in the present research were deposited in the Sequence Read Archive (SRA) dataset of the National Center for Biotechnology Information (NCBI) with the accession number PRJNA780464.

Acquisition of Soil Physicochemical Parameters and Climatic Data
A subset of the rhizosphere soil was air-dried for determination of its physicochemical properties. Soil pH and the contents of organic matter (OM), available N (AN), available P (AP), and available K (AK) were analyzed following the standard testing procedures described by Bao (2000). Soil pH was measured at 1:1 (soil:deionized water, m/v) ratio using a glass electrode pH meter. The OM content was determined by the potassium dichromate (K 2 Cr 2 O 7 ) volumetric method. The AN content was determined by the alkaline hydrolysis diffusion method. AP was extracted with 0.5 M NaHCO 3 solution and quantified by molybdenum antimony chromatometry. AK was extracted with 1 M NH 4 OAC solution and quantified by flame photometry. The air-dried soil was oven-dried at 105°C to constant weight to determine the moisture content. Data of mean annual temperature (MAT) and mean annual precipitation (MAP) were obtained from the Climatologies at High Resolution for the Earth's Land Surface Areas (CHELSA). 1 The soil physiochemical parameters and the climatic data are shown in Supplementary Table S2. 1 http://chelsa-climate.org/

Determination of Alkaloid Contents in Roots
The oven-dried root samples were ground to a fine powder for the determination of the contents of total alkaloids as well as four alkaloids (oxymatrine, oxysophocarpine, sophoridine, and matrine). Total alkaloids were assayed by the acid dye colorimetry method (Sun et al., 2021) with slight modifications. Briefly, pulverized root samples were ultrasound-treated in 0.2% HCl for 40 min to extract the total alkaloids. The extract was centrifuged, and the supernatant was collected and adjusted to pH 10.0. Subsequently, 5 ml of the supernatant was evaporated to dryness, and the residue was reconstituted in 5 ml of 80% ethanol and thoroughly mixed with 0.2 mM bromothymol blue (pH 7.0) and dichloromethane. Then, the dichloromethane layer was collected for colorimetric assay of the total alkaloids at 417 nm. Oxymatrine, oxysophocarpine, sophoridine, and matrine were analyzed by high-performance liquid chromatography (HPLC; Supplementary Figure S1) according to Pharmacopoeia of the People's Republic of China (Chinese Pharmacopoeia Commission, 2020). Powdered root samples (0.1 g) were ultrasound-treated in 0.133 ml of 25% ammonia + 8.33 ml of dichloromethane for 40 min for extraction of alkaloids. After the extract was filtered, 5 ml of the filtrate was evaporated, followed by dissolving the residue in methanol for HPLC analysis. Chromatographic separation was achieved on a Shimpack GIST C18 column (4.6 mm × 250 mm, 5 μm) at 30°C and was detected at an ultraviolet wavelength of 225 nm. A mobile phase consisting of solution A (0.01 M ammonium acetate, pH 8.1) and solution B [acetonitrile: ammonium acetate (0.01 M, pH 8.1), 3:2, v/v] at a flow rate of 1 ml min −1 was used. A gradient (A:B, v/v) program of 0 min (9:1) → 20 min (7:3) → 40 min (6:4) → 50 min (4:6) → 60 min (9:1) → 70 min (9:1) was employed. The content of matrine + oxymatrine was calculated. The alkaloid content data are available in Supplementary Table S3.

Statistical Analysis
The following statistical analyses were performed on the R software platform (R Development Core Team, 2021). The alpha diversity of the bacterial community was calculated for each rhizosphere sample using the species richness and Chao1 and Shannon indices. Significant differences in alpha diversity among the four geographical locations or the three plant ages were determined by Kruskal-Wallis H-tests with the Dunn-Bonferroni post hoc method (P < 0.05). Spearman's rank correlation was used to determine the relationships between the Shannon index and environmental variables (including climatic and edaphic variables). Nonmetric multidimensional scaling (NMDS) and principal coordinate analysis (PCoA), both of which were based on the Bray-Curtis distance matrix, were conducted to explore dissimilarities in the rhizosphere bacterial communities and root alkaloid accumulation among samples, respectively. To assess the significance of variation in the rhizosphere bacterial community composition or root alkaloid accumulation across geographical locations or plant ages, permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis distances was carried out with 999 permutations. In addition, the significant spatial, climatic, and edaphic variables shaping the rhizosphere bacterial communities were evaluated by redundancy analysis (RDA) with a forward selection procedure identifying the significance (Benjamini-Hochberg adjusted P < 0.05). The spatial variables were generated by principal coordinates of neighbor matrices (PCNM) analysis (Borcard et al., 2011; Supplementary Table S2). Based on RDA, variation partitioning analysis (VPA) was performed to measure the relative influence of spatial, climatic, and edaphic factors on the rhizosphere bacterial communities (Peres-Neto et al., 2006). Random forests and structural equation modeling (SEM) were employed for further analyses. Serving as a robust ensemble machine learning method for classification and regression, a random forest model was applied to identify the specific bacteria associated with the total accumulation of matrine and oxymatrine. The relative abundance of bacterial taxa at the OTU level was regressed against the total content of matrine and oxymatrine by the R package "randomForest" (Zhang et al., 2018c). Ranked lists of bacterial OTUs in order of their feature importance in the random forest model were achieved over 100 iterations of the algorithm. The optimal taxa set correlated with the total content of matrine and oxymatrine was generated based on the minimum errors obtained from 10-fold cross-validation with five repeats. The 30 most important OTUs were chosen to display. The SEM, which allows the assessment of multiple interacting factors and potential causal relationships, was applied to evaluate the direct and indirect effects of spatial (the first axis of principal component analysis for all PCNM vectors), climatic (MAT and MAP), and edaphic (pH, OM, AN, AP, and AK) factors and bacterial attributes (the relative abundance of three dominant rhizosphere phyla) on the accumulation of alkaloids (Fan et al., 2020). The SEM was carried out by using the R package "lavaan. "

Diversity and Taxonomic Structure of the Rhizosphere Bacterial Communities
A total of 7,478 non-singleton OTUs were retrieved from 72 rhizosphere samples after normalization of the read numbers. The LY, CZ, DL, and LN samples shared a large proportion of OTUs (6,489 OTUs, 86.77%; Figure 1A). Moreover, 7,407 OTUs (99.05%) were shared by rhizosphere samples of S. flavescens of three ages ( Figure 1B) (Figures 1C,D). The Chao1 index was highest in the CZ samples ( Figure 1D). Moreover, the Shannon index in the CZ samples (10.46 ± 0.08) was higher than that in the LY (9.97 ± 0.53) and DL (9.88 ± 0.29) samples ( Figure 1E). No significant differences in the rhizosphere bacterial diversity were found among the three plant ages (Supplementary Figure S2).
The Spearman's rank correlation test revealed that the rhizosphere bacterial diversity was positively correlated with the MAP and the contents of soil OM and AN, but negatively correlated with the MAT, soil pH, and AK content (Supplementary Table S4).
At the phylum level, Proteobacteria dominated the rhizosphere bacterial communities, followed by Actinobacteria and Acidobacteria (Figures 2A,B). These three most abundant phyla accounted for a large proportion (more than 70%) of the relative abundance of the rhizosphere bacterial communities. Compared with the LY, DL, and LN samples, the CZ samples hosted a higher relative abundance of Actinobacteria (26.07%) but a lower relative abundance of Proteobacteria (39.39%) and Acidobacteria (11.90%). The LY and DL samples presented the highest relative abundances of Proteobacteria (44.49%) and Acidobacteria (19.57%), respectively. We further characterized the rhizosphere bacterial communities at the order level. The results showed that the top 10 most abundant orders included six orders (Burkholderiales, Myxococcales, Rhizobiales, Rhodospirillales, Sphingomonadales, and Xanthomonadales) of the phylum Proteobacteria, two unidentified orders (iii1.15 and RB41) of the phylum Acidobacteria, one order (Actinomycetales) belonging to Actinobacteria and one order (Bacillales) belonging to Firmicutes ( Figure 2C). We found that Actinomycetales dominated the rhizosphere bacterial communities in the CZ and LY samples, whereas the unidentified order iii1.15 (belonging to the phylum Acidobacteria) dominated the rhizosphere compartment in the DL and LN samples. Furthermore, Actinomycetales was more abundant in the CZ (11.03%) and LY (8.41%) samples than in the DL (4.85%) and LN (4.82%) samples. Conversely, the total relative abundance of the two orders of the phylum Actinobacteria in the CZ samples (7.76%) was lower than that in the LY (10.72%), DL (12.95%), and LN (13.15%) samples.

Ecological Factors Shaping the Rhizosphere Bacterial Communities
NMDS ordination plots based on the Bray-Curtis dissimilarity matrix displayed separated clustering of rhizosphere bacterial communities across the geographical locations and plant ages ( Figure 3A). Similarly, PERMANOVA confirmed that both geographical location (P = 0.001) and plant age (P = 0.001) were significant driving factors for the variation in the rhizosphere bacterial community composition. However, geographical location (R 2 = 0.499) had a more pivotal effect than plant age (R 2 = 0.039). The results of RDA showed that two spatial variables (PCNM 1 and 3), two climatic variables (MAT and MAP), and three edaphic variables (pH, AN and AP) were significant drivers (P < 0.05) shaping the rhizosphere bacterial communities (Figure 3B). Among these ecological factors, the MAT (R 2 = 0.185) was the strongest. The results of VPA revealed that the climatic factors (35.0%) explained a larger proportion of variance in the rhizosphere bacterial community composition than the spatial (32.3%) and edaphic (20.1%) factors ( Figure 3C).

Rhizosphere Bacteria Associated With Alkaloid Accumulation in Roots
PCoA grouped with PERMANOVA revealed significant differences in alkaloid accumulation across the geographical locations (P = 0.001; R 2 = 0.298) and plant ages (P = 0.001; R 2 = 0.208; Figure 4). By using the Mantel test, we found that the rhizosphere bacterial microbiota of S. flavescens was remarkably correlated with the contents of oxymatrine (P = 0.001; R = 0.230), sophoridine (P = 0.001; R = 0.351), and matrine + oxymatrine (P = 0.001; R = 0.226) in roots (Supplementary Table S5). Considering that the total content of matrine and oxymatrine in the roots of S. flavescens is used as a quality control marker for this medicinal plant in the Pharmacopoeia of the People's Republic of China (Chinese Pharmacopoeia Commission, 2020), we focused on the association between the rhizosphere bacteria and the total content of these two molecules. A random forest machine learning algorithm was used to establish a model for the association. The 10-fold cross-validation with five repeats generated 935 important bacterial OTUs correlated with the total accumulation of matrine and oxymatrine in roots (R 2 = 0.592; Figure 5A). In order of importance, the top 30 bacterial OTUs are exhibited in Figure 5B according to their relative abundance distribution in roots with different accumulations of matrine + oxymatrine. Within these 30 bacterial markers, OTU 4,135, belonging to Chloroflexi, contributed most to the total accumulation of matrine and oxymatrine. Taxa of the phylum Actinobacteria (10 OTUs, in expectation of OTUs 25,035 and 17,632) dominated the bacteria whose increase in relative abundance favored the total accumulation of matrine and oxymatrine. In contrast, the increase in the relative abundance of certain bacteria presented a negative influence on the total accumulation of matrine and oxymatrine, which included three taxa (OTUs 25,490,13,741 and 14,939) of the phylum Gemmatimonadetes and two taxa (OTUs 1,671 and 13,635) of the phylum Acidobacteria. Additionally, for some bacteria, e.g., OTUs 2,709 and 12,697 (both of which belong to the Proteobacteria phylum), the association of their relative abundances with the total accumulation of matrine and oxymatrine could not be simply defined (Figure 5).
Interactions Between Spatial, Climatic, and Edaphic Factors, Bacterial Attributes, and Alkaloid Accumulation SEM was used to estimate the direct and indirect associations regarding causation between the spatial, climatic, and edaphic factors, the relative abundance of three dominant rhizosphere phyla, and the accumulation of alkaloids. The results indicated the following (Figure 6): (1) The spatial factors directly influenced the MAP.
(2) Soil pH and AK were positively affected by the MAT but negatively affected by the MAP. In addition, soil pH and AN were positively associated with the spatial factors. (3) The climatic factors had positive and negative effects on the relative abundance of Acidobacteria and Actinobacteria, respectively. Besides, the soil AN had positive and negative impacts on the relative abundance of Actinobacteria and Proteobacteria, respectively. In addition, the relative abundance of Acidobacteria was positively influenced by the soil pH. (4) Actinobacteria positively affected the total content of matrine and oxymatrine. Furthermore, the contents of matrine + oxymatrine and oxysophocarpine were negatively influenced by the soil pH. Moreover, the oxysophocarpine content was  positively associated with the spatial factors and the MAT. (5) The sophoridine content was also positively impacted by the MAT.

DISCUSSION
Alkaloids, one of the major classes of plant SMs, are critical bioactive constituents in the medicinal plant S. flavescens (He et al., 2015). Although plant SMs can benefit plant fitness and human health, the accumulation of SMs can vary among individuals (Moore et al., 2014;Pang et al., 2021). In the present study, we also discovered variations in the accumulation of alkaloids among root samples of S. flavescens, with geographic location and plant age being important sources of variation (Figure 4). It is recognized that the intraspecific variation in plant SMs is not only dependent on plant genotype and ontogeny but also influenced by various environmental factors involving biotic and abiotic factors (Moore et al., 2014;Yang et al., 2018). The importance of plant-associated microbes, as vital biotic factors, for explaining plant SM polymorphisms in medicinal plants has been highlighted (Huang et al., 2018). In addition, plant rhizosphere is a hotspot for studying the impacts of plant-associated microbes on the accumulation of SMs in plants. Therefore, we determined the taxonomic structure of rhizosphere bacterial communities of S. flavescens and addressed  the association between the rhizosphere bacterial microbiota and the accumulation of alkaloids in S. flavescens. Plants host a taxonomically diverse microbiota but only a few dominant taxa (Müller et al., 2016). The dominance of the bacterial phyla Proteobacteria, Actinobacteria, Bacteroidetes, Acidobacteria, and Firmicutes was discovered in previous studies (Bulgarelli et al., 2012(Bulgarelli et al., , 2015Lundberg et al., 2012;Zarraonaindia et al., 2015;Fan et al., 2017;Zhang et al., 2019). In the present research, we found that the rhizosphere bacterial communities of the legume S. flavescens were dominated by Proteobacteria (mostly from the Rhizobiales order), followed by Actinobacteria, Acidobacteria, Firmicutes, Gemmatimonadetes, and Bacteroidetes (Figure 2), revealing similarly dominant phyla with previous studies.
Then, we explored the dissimilarities in the rhizosphere bacterial communities of S. flavescens. The results revealed distinct community compositions across the geographic locations and plant ages and indicated that the geographic location accounted for most of the variation (Figure 3A). Previous reports on maize and soybean plants also indicated that the largest proportion of variance in rhizosphere microbiota could be attributed to the geographic location (Peiffer et al., 2013;Zhang et al., 2021). The large spatial distance of sampling employed in the present study may result in the dominant contribution of geographic location in structuring the rhizosphere bacterial communities. The results of VPA showed that the spatial, climatic, and edaphic factors contributed to the variation in the rhizosphere bacterial community composition (Figure 3C), suggesting that dispersal limitation and environmental filtering could account for the turnover of the rhizosphere microbial communities of S. flavescens across the geographic locations (Zhang et al., 2018a(Zhang et al., , 2018b. We also conducted RDA ordination to identify the key spatial, climatic, and edaphic drivers responsible for the variation in the rhizosphere bacterial communities of S. flavescens. The results showed that two spatial variables (PCNM 1 and 3), two climatic variables (MAT and MAP), and three edaphic variables (pH, AN, and AP) were significant in driving bacterial community turnover ( Figure 3B). Besides the geographic location, plant age could influence the assembly of rhizosphere microbiota (Micallef et al., 2009;Marques et al., 2014;Huang et al., 2020;Wang et al., 2020). Plants can select different microbes at the different stages of development, presumably for specific benefits (Chaparro et al., 2014). Overall, in the present study, the geographic location was a larger source of variation for shaping the rhizosphere bacterial communities than the plant age ( Figure 3A). Similar results were obtained with grapevine plants (Manici et al., 2017).
We evaluated the correlation between the rhizosphere bacterial microbiota of S. flavescens and the accumulation of alkaloids in roots of the plant by using the Mantel test. The results elucidated significant correlations between the rhizosphere bacterial communities and the contents of oxymatrine, sophoridine, and matrine + oxymatrine (Supplementary Table S5). Accordingly, we used the random forest model to regress the relative abundance of the rhizosphere bacterial OTUs against the total content of matrine and oxymatrine to identify the potential beneficial bacteria for enhancing the total accumulation of these two molecules in roots of S. flavescens, which is used as a standard for the quality evaluation of this medicinal plant. We found that 59.2% of the variation in the matrine + oxymatrine content could be explained by the random forest model (Figure 5), which reveals the important influence of rhizosphere bacterial microbiota on the accumulation of these two molecules in the plant. Within the top 30 important bacteria, positive markers were mostly from Actinobacteria and Chloroflexi, while negative markers predominantly belonged to Gemmatimonadetes and Acidobacteria.
Actinobacteria members, ubiquitous in terrestrial habitats, are known to produce extensive SMs, many of which are of great importance for medicine and agriculture (Barka et al., 2016;Choudoir et al., 2019). Host plants establishing symbiotic associations with some members of Actinobacteria were reported to have improved resistance to biotic stresses through activation of defense pathways or upregulation of certain SMs (Conn et al., 2008;Kurth et al., 2014). In addition, Yan et al. (2014) found that Streptomyces pactum Act12, a member of Actinobacteria, could improve the growth of hairy roots of Salvia miltiorrhiza and promote the synthesis of tanshinone in hairy roots by upregulating the expression of corresponding biosynthetic genes. As shown in Figure 6, our SEM also indicated a positive and direct effect of Actinobacteria on the total accumulation of matrine and oxymatrine. Acidobacteria is also widespread and abundant in nearly all ecosystems (Kalam et al., 2020). However, owing to difficulties in cultivation, the phylum is rudimentarily explored (Kielak et al., 2016;Kalam et al., 2020). Most of the capabilities and ecological functions of Acidobacteria members, such as decomposition of various biopolymers and plant growth promotion, are predicted from their genomic and metagenomic information (Kalam et al., 2020). Although our knowledge in terms of their influence on plant SM accumulation was limited, the covariance analysis of our SEM revealed a negative correlation between the relative abundance of this phylum and the phylum Actinobacteria (P = 0.002; R = −0.525). This negative correlation could also be indicated by the Spearman's rank correlation (P = 0.000; R = −0.778) and may allude to a competitive relationship between these two phyla, which may account for the negative association between the phylum Acidobacteria and the total accumulation of matrine and oxymatrine in the outcome of the random forest model (Figure 5). The phylum Chloroflexi has been found in anaerobic habitats, such as sediments and hot springs, exerting a role in biogeochemical carbon and nitrogen cycling (Hug et al., 2013;Speirs et al., 2019;Spieck et al., 2020). A previous study by Hu et al. (2018) demonstrated that Chloroflexi bacteria were attracted by root exudates of maize and were involved in host-medicated rhizosphere microbiota assembly to promote the growth and defense capability of the host. The phylum Gemmatimonadetes, containing phototrophic species, is also a nearly unexplored bacterial group (Zeng et al., 2020). The impacts of these microbes on the accumulation of bioactive components in S. flavescens should be illustrated in further studies based on isolation of the identified bacteria and cocultivation of them with the host plant.
Variations in abiotic factors, such as temperature, light, nutrients, and water, could also cause changes in the type and amount of plant SMs, which may be due to the regulation of the biosynthesis of SMs or a consequence of passive dilution or concentration due to the altered plant growth rate (Moore et al., 2014;Isah, 2019). Therefore, we applied SEM to evaluate the direct and Frontiers in Microbiology | www.frontiersin.org 9 December 2021 | Volume 12 | Article 781316 indirect effects of the spatial, climatic, and edaphic factors and the relative abundance of the three dominant rhizosphere phyla on the accumulation of alkaloids in the roots of S. flavescens. The results showed that the spatial, climatic, and edaphic factors had both direct and indirect effects on the accumulation of alkaloids (Figure 6). The direct and indirect (through alteration of the relative abundance of Actinobacteria) effects of edaphic factors on the accumulation of alkaloids implies the impact of comprehensive soil-microbe-root interactions in the rhizosphere on bioactive component accumulation in medicinal plants.
In conclusion, the present study indicated that the accumulation of alkaloids in S. flavescens could be influenced by the rhizosphere bacterial communities of the plant. We identified the potentially beneficial bacteria associated with the total accumulation of matrine and oxymatrine with a random forest machine learning algorithm. In addition, the causal relationship among the spatial, climatic, and edaphic factors, the three most abundant rhizosphere phyla (Proteobacteria, Actinobacteria, and Acidobacteria) and the accumulation of alkaloids was indicated by SEM. We proposed rhizosphere engineering of specific microbial members for the control and improvement of the growth and quality of medicinal plants in the future.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: PRJNA780464 accession is not yet publicly available.