The field experiment was conducted in Yaozhou District, Tongchuan City, Shanxi Province, China (34°48’–35°19’N, 108°34’–109°06’E) in the spring of 2022. The study area exhibits a characteristic temperate continental climate, characterized by an average annual temperature of 10.3°C, a precipitation level of 554.5 mm and a frost-free period of 217 d. Soil physicochemical characteristics were as follows: organic matter 0.22×10⁻⁴ mg/kg, total nitrogen 1246 mg/kg, available nitrogen 0.85×10⁻⁵ mg/g, available phosphorus 10.7762 mg/kg, and available potassium 0.1831×10⁻³ mg/kg.

A single-factor 4-level randomized block design was used to arrange the field experiment, with a plot area of 4 m × 4 m = 16 m². Each fertilization treatment was set up in three replicates, a total of 4 × 3 = 12 plots, and a protection area of 350 m² was set up around the experimental plot. The plant samples used in this study were one-year-old S. miltiorrhiza seedlings. Four fertilization treatments including base fertilizer (F1), foliar fertilizer (F2), base fertilizer +foliar fertilizer (F3) and blank control (CK) were applied. The base fertilizer was a compound fertilizer of nitrogen, phosphorus, and potassium from STANLEY from the local market, applied at a rate of 37.5 kg per acre. For foliar fertilizer, STANLEY potassium dihydrogen phosphate foliar fertilizer was used, with a spraying rate of 50–60 g per acre. Before cultivating seedlings, base fertilizer was applied to the soil, and then starting from the vigorous growth period of S. miltiorrhiza seedlings, spray foliar fertilizer three times every 15 d. Three ridges of 120 cm in width and 30 cm in height were set in each test plot, and a drainage ditch of 25 cm in width was set around the furrow. Two rows of S. miltiorrhiza were planted in each ridge with a row spacing of 15 cm×30 cm.

In March 2022, the experiment was planted and plant samples were collected at the end of the growing season in November 2022. Three healthy plants were randomly selected from each plot, and approximately 10 g of samples were collected from the leaves, stems, and roots using sterile surgical blades. Root samples for rhizosphere microbial analysis were placed in a 4°C cooler and subsequently brought back to the laboratory. Rhizosphere soil microbial analysis samples were from soil samples 10–15 cm away from the roots, and 200 g of soil sample was collected from each plot and placed in a sealed plastic bag. All soil samples were transported to the laboratory in insulated containers. Before experimenting, all samples were sieved (<2 mm mesh) to remove rocks, coarse roots and other litter. Soil samples used for enzyme analysis were stored in a 4°C refrigerator. Other soil subsamples were air-dried and used for determination of soil physicochemical properties. A total of 12 leaf, stem, root samples and 12 soil samples were collected. The samples for high-throughput sequencing were stored in a −80°C freezer for preservation.

Five grams of plant sample (leaves, stems, or roots) were weighed and put into a 50 mL centrifuge tube, with 50 mL of 0.1 M potassium phosphate buffer (PPB, pH=8.0) added. Plant sample in tubes were washed with 1 min sonication and 10 s vortex, and repeated. Then the samples were transferred to new tubes with 50 mL of 0.1M PPB and washed again. The suspension from two washes was mixed and filtered through a 0.2 µm membrane. The filter membranes with epiphytes were snap frozen in liquid nitrogen and stored at −80°C in the refrigerator for subsequent DNA extraction (Bodenhausen et al., 2013).

The genomic DNA from epiphytic microorganisms was extracted from filter membranes using the FastDNA^® Spin Kit for Soil (MP Biomedicals, USA) according to the user’s manual. The DNA purity and concentration were measured with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA), and DNA integrity was examined using 1% agarose gel electrophoresis and stored at −20°C in a refrigerator for subsequent experiments.

The ABI GeneAmp^® 9700 PCR thermal cycler (ABI, USA) was used to amplify the 16S V3-V4 region (5’-GTGCCAGCMGCCGCGGTAA-3’ and 806R, 5’-GGACTACHVGGGTWTCTAAT-3’) of epiphytic bacteria and the ITS1 region (ITS1F/ITS2, 5’-CTTGGTCATTTAGAGGAAGTAA-3’ and 5’-GCTGCGTTCTTCATCGATGC-3’) of epiphytic fungi (Tamaki et al., 2011; Tepper and Gaynor, 2015). The PCR reactions were performed in a 20 μL system, including 2 μL of 10×buffer, 2 μL of 2.5 mM dNTP, 0.8 μL of each 5 μM primer, 0.2 μL of Taq polymerase, 0.2 μL of BSA, 10 ng of template DNA, and topped up to 20 μL with ddH₂O. The amplification of the bacterial 16S V3-V4 region was carried out under the following conditions: denaturation at 95°C for 3 min; 30 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 45 s; and a final extension at 72°C for 10 min. For the fungal ITS1 region, PCR amplification was performed under the same reaction system and conditions, but with 35 cycles. Each amplification was repeated three times. The PCR products were recovered on a 2% agarose gel, and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). Subsequently, the recovered PCR products were quantified using the Quantus™ Fluorometer (Promega, USA). The purified amplicons were mixed in equal amounts, and a library was constructed using the NEXTFLEX^® Rapid DNA-Seq Kit. The final sequencing was conducted using the Illumina MiSeq PE300 platform at Shanghai Majorbio Bio-pharm Technology Co., Ltd. The raw data have been deposited in the NCBI SRA database (PRJNA983712, PRJNA982701).

Quality control for the raw reads was performed using software tools fastp (version 0.19.6) and FLASH (version 1.2. 11) with following steps: (1) filter bases with a mass value of less than 20 and read content containing N bases, set a window of 50 bp, and truncate bases if the average mass value in the window was less than 20, finally, the reading below 50 bp was filtered after quality control; (2) pairs of reads were spliced according to the overlap between PE reads, and the minimum length of overlap was 10 bp; (3) a maximum mismatch ratio of 0.2 was allowed in the overlap region of the spliced sequences, and inconsistent sequences were removed; (4) samples were demultiplexed based on the barcode. Quality-control splicing sequences were clustered into operational taxonomies (OTUs) based on 97% similarity using the UPARSE software (version 7.1). All sequences with mitochondrial and chloroplast annotations were removed. Taxonomic placement of epiphytic bacteria and fungi was annotated according to Silva 16S rRNA gene database (V 138) and UNITE database (version 8.0), respectively, based on an RDP classifier (version 2.11) with a 70% confidence threshold. The community composition of each sample was analyzed at various species taxonomic levels. Bacterial function prediction was based on the PICRUST2 database (Douglas et al., 2020), and bacterial function prediction was based on the FUNGuild database (Nguyen et al., 2016).

The weight of roots and stems was determined using an analytical balance. The lengths of the roots and shoots were ascertained using a tape measure. Plant water content was determined using the drying method. As described previously, the NBT photochemical reduction method was selected to test superoxide dismutase (SOD) activity (Cheng et al., 2015). Soluble protein content was determined using the Bradford assay (Kielkopf et al., 2020). Abscisic acid (ABA, ng/mL), cytokinin (CTA, ng/mL), gibberellin (GA, pmol/mL), Indole-3-acetic acid (IAA, nmol/L), and nitrate reductase (NR, µg/g.h) were measured using enzyme-linked immunosorbent assay (ELISA) detection kits.

The medicinal ingredients of S. miltiorrhiza referred to in the 2020 edition of Chinese Pharmacopoeia I, and the lipid soluble chromatographic conditions were as follows: Silica gel bonded with octadecyl silane was used as the filler, acetonitrile was employed as mobile phase A, and a 0.05% phosphoric acid solution served as mobile phase B. The gradient elution conditions were indvided four parts: 0–6 min (A was 61% and B was 39%), 6–20 min (A changed from 61% to 90% while B changed from 39% to 10%), 20–20.5 min (A changed from 90% to 61% and B changed from 10% to 39%) and 20.5–25 min (A was 61% and B was 39%). The flow rate was maintained at 1 mL/min. The fat-soluble components determined included cryptotanshione, Tanshinone I, and Tanshinone IIA. Water soluble chromatographic conditions were as follows: using a C₁₈ chromatographic column, acetonitrile was employed as mobile phase A, and a 0.05% phosphoric acid solution served as mobile phase B. The four parts of gradient elution for 0–15 min, 15–30 min, 30–40 min and 20.5–25 min were A of 17–23%, 23–25%, 25–90%, and 90% and B of 83–77%, 77–75%, 75–10%, and 10%, respectively. The column temperature was maintained at 30°C with a flow rate of 1 mL/min at a detection wavelength of 286 nm. The water-soluble components determined included salvianolic acid B and rosmarinic acid.

Soil pH was measured in a 1:2.5 (w/w) soil: water suspension with pH 3000 (STEP Systems Gmbh, Germany). Soil organic carbon (SOC) was determined using the combustion loss method (Heiri et al., 2001). The available phosphorus in soil (AP) was measured using the sodium bicarbonate extraction-molybdenum antimony colorimetric method (Olsen et al., 1954), with 0.5 mol/L sodium bicarbonate solution used for extracting the available phosphorus, which reacted with molybdenum antimony to generate phosphomolybdenum blue. The available potassium (OP) was determined using the tetraphenylboron method (Chen et al., 2019), with 1 mol/L NaNO₃ solution used for extracting soil K⁺, which reacted with tetraphenylboron in a weakly alkaline medium to produce a barely soluble white precipitate. The measurement of ammonium nitrogen (NH₄⁺-N), nitrate nitrogen (NO₃⁻-N), total phosphorus (TP), and total nitrogen (TN) was performed using the Smartchem 200 analyzer (Alliance, France) (Xie et al., 2017). The activity of sucrase (SC) was determined using the method of Guan (1986). The sucrase in soil catalyzes the hydrolysis of sucrose into reducing sugar, which reacts with 3,5-dinitrosalicylic acid to produce orange 3-amino-5-nitrosalicylic acid under the boiling condition, and the depth of color is positively correlated with the content of reducing sugar. The measurement of nitrate reductase (NR) activity was carried out using the sulfanilamide diazotization colorimetric method (Zhao et al., 2009), where the reaction of nitrite with sulfanilamide and α-naphthylamine under acidic conditions produces a red compound. The urease (URE) activity was measured using the improved Hoffmann and Teicher colorimetric method (Kandeler and Gerber, 1988), with urea as the substrate, and indophenol, which is generated by the reaction of enzyme products and phenol-sodium hypochlorite, is used for analyzing the activity of urease. The activity of alkaline phosphatase (ALP) was measured using the method described by Tarafdar and Marschner (1994), with disodium phenyl phosphate (pNPP) as the substrate. The substrate is hydrolyzed by soil acid phosphatase to produce yellow p-nitrophenol (pNP), and the amount of pNP produced is directly proportional to the absorbance of the yellow solution, which can be used for quantitative analysis. The data analysis was performed using SPSS software version 24.0 (IBM SPSS Statistics, USA).

Alpha diversity (Shannon, Chao1, goods_coverage, and pielou_e indices) and beta diversity analysis (NMDS) were computed using the QIIME2 software. Statistical analysis and plotting for Mantel tests were performed using three R packages (dplyr, linkKEt, and ggplot2; http://www.R-project.org). Bivariate correlations between plant growth parameters, soil physicochemical properties, and microbial communities were analyzed using SPSS 24.0 software, with Tukey’s test used to compare means (P<0.05). Variance partition analysis (VPA) was used to investigate the effects of growth parameters, soil properties, and plant hormones on the microbial abundance in different parts of Salvia miltiorrhiza. Pheatmap and Vegan were used for statistical analysis and data visualization. Network parameter analyses were conducted using six R packages (igraph, psych, Hmisc, vegan, dplyr, and reshape2). The βNTI value was calculated to assess community assembly processes, using the formula below (Equation 1):

The regression curves of OTUs observed at the OTU level for S. miltiorrhiza epiphytic bacteria and fungi gradually flattened, indicating that the amount of sequencing data for epibiotic bacteria and fungi was reasonable and the OTU depth met the requirements for diversity analysis (Supplementary Figure S1). A total of 4256 OTUs of epiphytic bacteria were obtained, belonging to 361 species, 544 genera, 146 orders, 251 families, 53 classes and 27 phyla. In addition, 1702 OTUs of epiphytic fungi were identified, belonging to 299 species, 226 genera, 123 families, 56 orders, 23 classes, and 8 phyla.

At the phylum level (Supplementary Figure S2A-a), Actinomycetes and Firmicutes were the dominant groups of epiphytic bacteria in S. miltiorrhiza leaves, with the abundance of Actinomycetes (13.25%, 8.42%) and firmicutes (13.25%, 8.42%), respectively. The relative representation of Firmicutes (10.22%) and Bacteroidetes within the root system exhibited a reduction, whereas the prevalence of Proteobacteria (75.77%) experienced an elevation. Interestingly, Acidobacteriota expression decreased under both fertilization treatments (Supplementary Table S1). At the genus level (Supplementary Figure S2A-b), a decline in the prevalence of Ralstonia and Massilia was observed across both fertilization treatments. The concentration of Hymenobacter (8.21%, 9.26%), Rhizobium (7.53%), and Sphingomonas (8.97%) also increased. The proportion of Pseudomonas was elevated to 19.14% and 13.82% in the group F2 in leaves, heightened in the groups F2 and F3 in stems, rose in the group F1 in roots, and reduced in the other treatment groups. At the phylum and genus level (Supplementary Figures S2A-c, d), fertilization treatments had a lesser effect on altering the abundance of the epibiotic fungal community in S. miltiorrhiza.

After the fertilization treatment, the unique bacterial OTUs of the three parts of S. miltiorrhiza were lower than those of the blank treatment except for the YF1 treatment (Supplementary Figure S2B-a), and the number of bacterial OTUs was curtailed after the fertilization treatment. The most common bacterial OTUs were found in roots. The number of unique fungal OTUs in the three parts of S. miltiorrhiza increased as a whole except for a few (Supplementary Figure S2B-d). The combination of foliar fertilizer plus base fertilizer showed the greatest increase in OTUs, doubling the amount on average in all three sites. The leaves had the most common fungal OTUs, and the number of bacterial and fungal OTUs decreased under the combination of leaf fertilizer and basal fertilizer.

When evaluating the α-diversity of epibiotic bacteria in different parts of S. miltiorrhiza based on the OTU level we found that the Good’s coverage index values were all close to 1, indicating that the sequencing depth was qualified. This further indicated a negligible likelihood of undetected sequences across the analyzed samples (Supplementary Table S2). The Shannon index and Chao1 index of bacteria and fungi were tested by T-test. After fertilization, no significant change was found the diversity of epiphytic bacteria and fungi from in the three parts of S. miltiorrhiza, but there was a substantial variation in their abundance (Figure 1). Different fertilization treatments had varying effects on abundance. For example, in the leaf section, the concentration of bacteria and fungi was significantly higher in groups F2 and F3 than in group CK. This pattern was also observed in the root section, while only bacteria showed a difference in the stem section. Under fertilization treatment, the diversity of epiphytic bacteria and fungi (excluding stem-associated bacteria) had undergone alterations across groups, albeit not apparently, namely there was a reduction in bacterial diversity and an increase in fungal diversity. Excluding stem bacteria, pronounced fluctuations in microbial community proportions were evident, with the F3 treatment cohort exhibiting the most marked alterations, typified by diminished bacterial and heightened fungal levels.

The NMDS and ANOSIM tests indicated that there was no salient dissemblance in the bacterial and fungal community structure on different parts of S. miltiorrhiza, and the same phenomenon was observed among treatment groups. (Supplementary Table S3; Figure 2). Bacterial community similarity analyses revealed divergent structure compositions in leaves between F2 and F3 cohorts versus groups CK and F1. Analogous differentiation patterns were noted across stem-based bacterial populations. Stem-based bacterial community structures exhibited variability between groups F2 and CK. In the roots, all treatments showed divergent community structures (Figures 2A–C). For fungal communities, the community structure of the F3 treatment group was distinct from the CK, F2 and F1 treatments. The fungal community composition in the F2 treatment diverged from the control (CK). In stem, the communities within the F3 treatment were distinct compared to those in the CK, F2, and F1. Community structure in the F1 exhibited variations from the CK. The community structures of epiphytic fungi and bacteria across the three parts of rhizome and leaves were different but these changes were not prominent.

Based on the zero βNTI system diversity model analysis, the results showed that in S. miltiorrhiza (root, stem, and leaf), at different positions of the surface, bacterial and fungal microbial communities were built in a largely deterministic process (defined as β|NTI| > 2) (Dini-Andreote et al., 2015) (Figure 3). This choice was mainly homogeneous. In fungal communities, all three fertilization methods converted the microbial construction process into a deterministic one. Among the bacterial communities, group F1 had the greatest impact on the construction process compared to the other two treatments, and roots were the most sensitive to the construction process.

A co-occurrence network map was constructed with the measured OTUs to illustrate the general symbiosis model of bacteria (Figure 4) and fungi (Figure 5) under different fertilization treatments for the key species group of epiphytic bacteria. In the leaf bacterial network, Bacteroidota and Hymenobacter were significantly increased in the F1 compared with the CK. Firmicutes increased in F2 and F3, and the composition of the major genus-level nodes in this phylum was no longer dominated by Lactobacillus and Faecalibacterium, but switched to other genera (Figure 4A). Moreover, the modularity index and the number of nodes increased after fertilization, and the proportion of positive correlation decreased, which might make the community structure more stable. In the stem bacterial network, compared to the CK treatment, the F1 had more nodes and edges, but the modularity index was reduced. It can also be seen from Figure 4B that most nodes were clustered together, and the distribution was not uniform. For groups F3 and F2, the node distribution was more uniform, meaning that even though the number of nodes and edges decreased, the modularity index was increased, thereby suggesting that the bacterial community in the stem may be more stable after the application of foliate fertilizer. In the root bacterial network, the community structure of roots was more stable than that of leaves and stems, and the modularity index of all treatment groups was greater than 0.4. The positive and negative relationships among species and the node distribution were also more uniform (Figure 4C). There were no significant differences between the treatment groups, but the number of nodes and edges in the network decreased after fertilization.

In the leaf fungal network, compared with CK treatment, Basidiomycota was the dominant phylum in groups F1 and F2, and Cryptococcus and Rhodotorula were the leading genera in groups F1 and F2. In the CK, Dioszegia, Moesziomyces and Tilletiopsis were the main genera. Ascomycota was the principal phylum in CK and F3, while Alternaria and Epicoccum were the major genera in CK and F3, and Mycosphaerella was added in F3 (Figure 5A). Fertilization treatment gradually reduced the modularity index of leaf epiphytic fungi and the stability of the network. In the stem fungal network, compared to the CK treatment, the F2 treatment group was dominated by Basidiomycota, and the species at the genus level had changed greatly. The F2 treatment group was mainly controlled by Vishniacozyma and Cryptococcus. Moreover, Tilletiopsis was the primary species in the CK treatment. Ascomycota was the prevailing phylum in groups CK, F1, and F3. Setophaeosphaeria and Epicoccum were the superior species in CK. The genera Plectosphaerella, Knufia, and Selenophoma were prevalent in the F1, and Alternaria and Paraphoma were predominant genera in the F3 (Figure 5B). Importantly, fertilizer treatment also changed the composition of fungal communities in the stems. In the root fungal network, compared to group CK, the distribution of nodes in F1 was more uniform, the number of nodes and edges in F2 was reduced, and the modularity index of the two treatments was greater than 0.4. Additionally, Mucoromycota and Mortierellomycota were added as new phyla to the F3. But the modularity index was less than 0.4 (Figure 5C).

PICRUST2 database was used to predict the function of the bacterial community under different fertilization treatments (Supplementary Figure S3-a). It is remarking that the protein function was higher in leaves than in other parts. For instance, the protein function of glycosyltransferases involved in cell wall biosynthesis, as well as the acyl coenzyme A dehydrogenase associated with the alkylation reaction protein AidB, were significantly enhanced. These functions play a pivotal role in plant development and stress response.

The function of epibiotic under different fertilization treatments was determined by the The FUNGuild database, and the result was plotted in Supplementary Figure S3-b. The predicted fungal functions mainly include three categories: pathogenic, symbiotic, and saprophytic. It is subdivided into animal pathogens, arbuscular mycorrhizal fungi, ectomycorrhizal fungi, lichenized fungi, mycoparasites, plant pathogens, undefined saprotrophs, and wood saprotrophs. Animal pathogens, plant pathogens, and undefined saprophytic fungi were saliently enhanced.

The results revealed significant disparities in the growth parameters of the leaves, stems, and roots following fertilization. (Figures 6A, B). There was a reduction in the levels of abscisic acid (ABA), superoxide dismutase (SOD) and nitrate reductase (NR), while the levels of Cryptotanshinone (CTS), rosmarinic acid (RosA) and salvianolic acid B (SalB) increased, The gibberellin (GA), soluble protein (SP), and tanshinone IIA (TSN-SS) content decreased in the stems, while the content of tanshinone I (TI) and RosA increased. The concentration of GA, ABA and CTS decreased in the roots, but the content of liposoluble medicinal ingredients increased. The comparison among different treatment groups revealed that the impact of group F3 was the most pronounced.

Analysis of the ANOVA (Figure 6C) showed that leaf growth parameters and hormones together explained 18% of the variation in the epibiotic bacterial community but without the fungal variation. For the stem, growth parameters and hormones accounted for 14% and 45% of the bacterial community variation, respectively, and were the main factors affecting the bacterial community. 11% of the variation in fungal communities can be attributed to growth parameters, but the combined effects of growth parameters and hormones can only cause 9% of changes. In terms of roots, growth parameters and hormones were the main factors affecting the bacterial community, and the variation in the bacterial community can reach 20% and 41%, respectively. Variations in the fungal community are predominantly governed by growth parameters and hormonal activity, contributing to 32% and 38% of the modifications, correspondingly.

Many studies have shown that reasonable combined application of fertilizers can improve crop yield and quality. For example, for tea, fertilizers elevated the content of amino acids and tea polyphenols while reducing the phenol-to-ammonia ratio (Liu et al., 2023a). Similarly, cucumber yield and vitamin C content were improved with fertilizer application (Wang et al., 2023). In apples, fertilizers increased yield, soluble sugar and vitamin C content, resulting in a higher sugar–acid ratio (Wang et al., 2022b). In this study, plant growth parameters, biomass and the levels of active components were increased, indicating that fertilization could improve the quality of S. miltiorrhiza. Importantly, the effect of the F3 treatment was the most significant.

Recent investigations have shown that the structure of microbial communities is sensitive to fertilization (Lozupone et al., 2012). That is, nutrients can evidently affect the community changes of key microbial species and ultimately regulate the construction of microbial communities (Schmidt et al., 2014). These keystone species are closely related to the nutrient cycling process of carbon, nitrogen, and phosphorus in soil and play a considerable role in improving crop productivity (Li et al., 2017b). The analyzed results from the βNTI model showed that the composition of bacterial and fungal communities in different parts of S. miltiorrhiza (roots, stems, leaves) was a deterministic process. As an abiotic factor, fertilization determined the microbial community, which had a significant effect on the composition of bacterial and fungal communities in S. miltiorrhiza. Jiao and Lu (2020) analyzed soil fungal communities in farmland, forest, wetland, grassland, and desert ecosystems in the Hexi Corridor of China and found that the communities of rare taxa were mainly regulated by deterministic processes of homogeneous selection. Similarly, Xiong et al. (2021) studied fungal community structure in soil, root surfaces and leaf surfaces under different fertilization practices in maize/wheat rotation and maize/barley rotation field systems in different places and found that fungal community structure was mainly affected by deterministic processes. Although only basal fertilizer was applied, the microbial community construction in the shoot was still a deterministic process, indicating a strong connection between the belowground section of plants and the shoot.

Recently, the concept of microbial-root-shoot was also proposed (Almario et al., 2017). In the fungal communities, all three fertilization methods made the microbial construction process a deterministic one in our study. Among the bacterial communities, the F1 treatment had the greatest impact on the construction process compared to the other two treatments, and the root was the most sensitive to the construction process. The different responses of these bacteria and fungi may stem from changes in community stability, as well as varying sensitivities to fertilizers among fungi and bacteria. The results showed that the selection process of the construction of the epiphytic communities of S. miltiorrhiza was distinctly affected by fertilization. Fertilization had a filtering effect on the species of epiphytic communities of S. miltiorrhiza, which directly affected the survival of the epiphytic communities. Moreover, the selection was mainly homogeneous, and the similar abiotic environments such as fertilization interact with the surface microorganisms of S. miltiorrhiza to form the selection effect. Combined current studies show that fertilization can evidently affect microbial community construction, markedly impact the composition and abundance of key species (Lin et al., 2019), and drive changes in microbial community structure and ecosystem function (Fan et al., 2019; Han et al., 2022).

It is proved that the fertilization can significantly change plant microbial community composition and diversity (Sabir et al., 2021; Guan et al., 2022; Zhang et al., 2022). An increase in nutrient availability correlates with accelerated microbial growth rates and enhanced substrate metabolism (McCann et al., 1998). Fertilizers are even shown to impact light competition after eutrophication (Hautier et al., 2009), which leads to a reduction in plant microbial diversity. In this study, the abundance of the dominant taxa and genus groups of the epimicroorganisms in different ecological niches of S. miltiorrhiza decreased after fertilization, and the changes of fungi were reduced compared to bacterial communities, which may be due to the different establishment methods of fungal and bacterial community structure. Among them, in comparison to the other treatment groups, the effect of F3 was significantly stronger. These results indicated that fertilization greatly affected the composition and abundance of vital species in the epimicrobial community of S. miltiorrhiza (Lin et al., 2019). Data from diversity index of α, NMDS, and ANOSIM tests showed that the diversity of bacterial and fungal flora in different parts of S. miltiorrhiza was altered but there was no meaningful difference between the groups, while the abundance of bacterial flora was markedly different. The findings indicated that fertilization does not exert a substantial impact on the diversity of S. miltiorrhiza epiflora. The lack of this effect may be attributed to the complexity of the epiphytic environment, where there are multiple factors that affect microbial colonization, and fertilization is not the main determining factor (Bringel and Couée, 2015). Despite the absence of statistical significance, modifications in certain pivotal species were observed, warranting additional investigation into the alterations of these key species.

Microbes live in complex communities (Ley et al., 2006; Fuhrman, 2009; Faust and Raes, 2012). In addition to abiotic factors affecting microbial diversity, interactions among organisms also play an important role in the overall composition, stability, and biodiversity of microbial ecosystems (May, 1988; Wardle, 2006; Mougi and Kondoh, 2012). They can compete for resources, inhibit the growth of other communities by producing antibiotics, or support each other by cross-feeding (Ives and Carpenter, 2007; Ptacnik et al., 2008; Pande et al., 2014). In general, higher biodiversity often (but not always) implies greater ecosystem stability (Hector et al., 1999). The structure and diversity of microbial communities are sensitive to external environmental factors such as fertilization and irrigation (Lozupone et al., 2012; Bhattacharyya et al., 2017; Chen et al., 2020; Dangi et al., 2020; Khan et al., 2021). That is, nutrients affect the community changes of key microbial species and ultimately regulate the construction of microbial communities. These keystone species are closely related to the nutrient cycling process of carbon, nitrogen, and phosphorus in soil, and have a significant effect on improving crop productivity (Li et al., 2017b).

The node, edge, node type, positive correlation, and modularity index of the co-occurrence network were all increased after different fertilization treatments, indicating that community structure became more stable. The root and stem network had an increased abundance of Proteobacteria, a nutrition-sensitive phylum that include a variety of pathogenic bacteria (Dai et al., 2018). Acidobacteria species decreased in leaf sites, and Acidobacteria abundance was negatively correlated with nutrient content (Fierer et al., 2007; Banerjee et al., 2016), which can promote the dissolution of soil inorganic P, and their abundance was significantly correlated with soil available P (Flieder et al., 2021; Wu et al., 2021). Increased root Rhizobia, microorganisms that favor crop–microbe interactions, usually perform beneficial functions for the host by providing a variety of nutrients and metabolites (Erlacher et al., 2015; Gazdag et al., 2018; Chang et al., 2022). Ascomycota was the most abundant phylum in the fungal network, and fertilization promoted the growth of Ascomycota (Guo et al., 2020). After fertilization, Cryptococcus, Rhodotorula, Apicophorus, and Alternaria were widely distributed in the fungal network structure. The structure of group F3 was more stable among different treatments, mainly origin the combination of two kinds of fertilization resulted in more sufficient and comprehensive fertility across the different parts of the plant. The results showed that fertilization treatment can change the composition and structure of key communities in bacterial and fungal community structures of S. miltiorrhiza (Feng et al., 2018). This may be due to the difference in the type of fertilizer applied, which can notably influence microbial abundance and network complexity (Fan et al., 2019), thereby increasing the abundance of keystone species and radically changing the structure of the epiphytic microbial communities of S. miltiorrhiza.

In the Mantel correlation analysis of growth parameters and microbial community, no matter which site, the significance of microbial community and growth parameters in F3 was greater than that in other treatment groups, indicating that F3 treatment group had a more noteworthy effect on the epiphytic communities of S. miltiorrhiza. The medicinal components in the roots were markedly positively correlated with Sphingobium, Mortierella, Fusarium, and Epicoccun. The medicinal components of the stem were evidently positively correlated with Novosphingobium and Dioszegia, while the medicinal components of the leaves were clearly positively correlated with Blautia and Genolevuria. The variance decomposition showed that the hormones ABA, CTK, and GA were the main factors responsible for the changes in the epiphytic communities of S. miltiorrhiza. Studies have shown that Sphingomonas can improve plant resistance to drought, salt and alkali, and heavy metals, and can synthesize some plant hormones to promote plant growth (Haney et al., 2015). The Fusarium can produce plant-stimulating hormones (gibberellin), which increases crop yield. These keystone species possess the capability to promote plant nutrient uptake, and enhance plant tolerance to environmental stressors (Trivedi et al., 2020; Zhang et al., 2020).

In the Mantel correlation analysis showed that the influence of environmental factors on the rhizosphere bacteria of S. miltiorrhiza surpasses that on fungi, and these factors exert a stronger impact on the roots compared to other plant parts. Apparently, within the various treatment groups, the F3 treatment demonstrated a more pronounced influence on environmental factors in comparison to the other treatment groups. The phyla Proteobacteria, Actinobacteria, and Bacteroidetes showed a positive correlation with NO₃⁻-N and a negative correlation with NH₄⁺-N (Liu et al., 2023b). In addition to fertilization, these indicators are influenced by various other factors, such as soil pH (Wang et al., 2017), root exudate metabolites (Wu et al., 2017), and microbial diversity levels (Xun et al., 2019). Spatially different biotic and abiotic environmental factors may lead to different selection pressures between bacteria and fungi (Li et al., 2017b; Saleem et al., 2018). Keystone species significantly affected the process of community construction, likely because they exhibited high connectivity in community microbial networks and were good predictors of community deviation and turnover (Herren and McMahon, 2017, Herren and McMahon, 2018).

Analysis of the epiphytic communities in niches showed that fertilization significantly affected some key species in the epiphytic communities of S. miltiorrhiza. In the bacterial co-occurrence network, Bacteroidota and Firmicutes were the dominant phyla, and other indices in the network increased, indicating that the bacterial community network was more stable after fertilization. In the fungal co-occurrence network, Basidiomycota, Mucoromycota, and Mortierellomycota were the dominant phyla, while other indices in the network decreased, showing that the fungal community network was no longer stable after fertilization. In the functional prediction, the protein functions of epiphytic bacterial communities were markedly enhanced. The function of pathogenic and saprophytic fungi was ​elevated, and that of commensal fungi was reduced. Mantel correlation analysis revealed heightened coefficients across all plant segments under the F3 treatment compared to other treatment cohorts, denoting F3 exerted a more pronounced impact on S. miltiorrhiza epiphytic populations. According to the variance decomposition, ABA, CTK, and GA were the main factors responsible for the changes in the epiphytic communities of S. miltiorrhiza. Studies have shown that Sphingomonas can improve the resistance of plants to drought, Salinization, heavy metals, and can synthesize some plant hormones to promote plant growth. Fusariums can produce plant stimulating hormone (GA), which can improve crop yield. These key species could enhance plant nutrient absorption, boost plant disease resistance, and improve plant stress tolerance (Trivedi et al., 2020; Zhang et al., 2020). Our work found that fertilization can select dominant and stable key species in the community, and then stimulate plants through hormones secreted by these key species to improve the yield and quality of S. miltiorrhiza. The effect of F3 treatment on the epiphytic communities was significant and gave the highest quality of S. miltiorrhiza. Therefore, we determined the best fertilization method is the combination of base fertilizer and foliar fertilizer.