The experimental field site was located in Changzhai Village, Changshun County, Guizhou Province, China (26°01′25″ N, 106°30′55″ E; elevation, 1,004 m). The strata are mainly light-colored limestone of the Permian Qixia Formation and Maokou Formation. The soil type is yellow soil according to the Chinese soil classification system (Gong, 1999) and Orthic Acrisols according to the World Reference Base (WRB) soil classification system (USS Working Group WRB, 2015). The region has a subtropical monsoon humid climate. The annual average temperature is 15.1°C, the annual average precipitation is 1396.7 mm, the annual sunshine duration is 1202.1 h, and the annual frost-free period is 275 days.

Four treatments were set up in the field experiment: blank control group (CK), no fertilization; CF group, compound CF was applied; OF group, ordinary OF was applied; OMF group, potassium-containing OMF was applied. Four parallel plots (5 m × 5 m) were set for each treatment group (total, 16 plots). The randomized complete block design was used to divide the sample area blocks (Supplementary Figure S1). On 7 May 2019, P. frutescens was planted and fertilized. The P. frutescens seeds were the 1^st generation hybrids cultivated by the Oil Materials Research Institute of Guizhou Academy of Agricultural Sciences (China).

The fertilizer specifications used were as follows:

The raw materials for OF fermentation were mushroom residue, distiller’s grains, straw, and chicken manure, which were mixed according to a mass ratio of 1:1:1:2 and EM microbial agent was added for fermentation (1 kg bacterial agent was added for every 10 T substrate). The mixture was then fermented in a fermentation tank for 30 days with periodic stirring. The total nutrient content (N + P₂O₅ + K₂O) was 7.2%, among which the content ratio of N:P₂O₅:K₂O was 0.9%:2.9%:3.4%, organic matter was 48.3%, and pH was 7.7.

Potassium-containing OMF was made by mixing the raw fermentation materials of the OF and potassium-containing rock powder, which contained 76% potassium feldspar, and passing the mixture through a 2-mm sieve. In accordance with a mass ratio of 3:1, the chemical composition was as follows: Al₂O₃, 17.11%; SiO₂, 54.06%; K₂O, 9.09%; CaO, 1.9%; Fe₂O₃, 6.15%; and MgO, 3.41% (Sun et al., 2019). Then, EM microbial agent was added for fermentation (1 kg bacterial agent per 10 T substrate). The raw materials were thoroughly mixed and fermented for 30 days with periodic stirring. The total nutrient content of OMF was 6.4%, with a content ratio of N: P₂O₅: K₂O of 1.6%:1.1%:3.7%, the content of organic matter was 56.4%, and the pH was 7.6. OF and OMF were produced by Guizhou Guifu Ecological Fertilizer Co., LTD. (China).

The compound CF was produced by Guizhou Xiyang Fertilizer Co., LTD. (China) and had a total nutrient content ≥45%, N:P₂O₅:K₂O content ratio was 15%:15%:15%. Urea was produced by Guizhou Chitianhua Tongzi Chemical Co., LTD. (China) with a total N ≥ 46.4%. Phosphate fertilizer was produced by Fuda Phosphorus Chemical Co., LTD. (China) with P₂O₅ ≥ 12%. The proportions of N, P, and K in OF and OMF groups were balanced by urea and phosphate fertilizer. All treatment groups had fertilizer applied with a content ratio of N:P₂O₅:K₂O of 15:15:15, and application amounts are shown in Table 1.

Rhizosphere soil sampling was conducted on 19 August 2019, which was the 105^th day of P. frutescens growth. Rhizosphere soil samples were arbitrarily collected from eight P. frutescens plants in each plot and mixed into one sample. A total of 16 mixed soil samples were collected. When collecting rhizosphere soil, the whole plant was first dug up, and the scattered soil at the root was gently shaken off. The remaining soil attached to the root system was considered the rhizosphere soil (Sun et al., 2021). The soil samples were frozen and transported to the laboratory on dry ice, passed through a 2-mm sterile sieve, and the plant roots were removed. Each sample was divided into two parts: one part was naturally air-dried to determine physicochemical properties, and the other was stored in a − 80°C freezer for DNA extraction.

Soil pH was analyzed by vibrating slurry with a water:soil ratio of 2.5:1 (v/w) and determined using a pH meter (Mettler-Toledo FE28, Switzerland; Marcos et al., 2019). Soil total organic carbon (TOC), total organic nitrogen (TON), total carbon (TC), and total nitrogen (TN) were determined using an elemental analyzer (Vario MACRO Cube, Germany; Liu et al., 2014). Total phosphorus (TP) and total potassium (TK) were determined by sodium hydroxide melting flame spectrophotometry (Ren et al., 2016). Available phosphorus (AP) was determined by the NaHCO₃ method, and available potassium (AK) was determined by ammonium acetate extraction–flame spectrophotometry (Lu, 1999).

On 21 September 2019, P. frutescens were harvested on the 137th day of growth. The P. frutescens yield was measured after harvest. The main indicators were plant height, biomass per plant, number of stem nodes, number of effective branches at one time, length of effective branches at one time, number of branch angles, number of panicles per plant, length of the main panicle, number of fruits in a single row of the main panicle, number of fruits on the main ear, and number of grains in the 10 main ears. After the perilla seeds were harvested, seed quality was inferred by determining fatty acid content (e.g., palmitic acid content, stearic acid content, oleic acid content, linoleic acid content, and linolenic acid content), total lipid content, and crude protein content. The fatty acid content was determined by gas chromatography–mass spectrometry (GCMS-QP2010, Schimadzu, Japan; Javier et al., 2018; Ko et al., 2018). Fat content was determined by sequential Soxhlet extraction according to ISO 659:2009 (Kourimska et al., 2018). The crude protein content was determined by the Kjeldahl method, the nitrogen concentration of the sample was calculated with a conversion factor (6.25), and the total nitrogen and protein mass were determined to obtain the crude protein content (Yaldiz and Camlica, 2020).

Total soil DNA was extracted from 0.5 g soil according to the manufacturer’s instructions for the E.Z.N.A.® Soil DNA Kit (Omega Bio-Tek, USA). The DNA extraction quality was detected using 1% agarose gel electrophoresis, and the DNA concentration and purity were determined using a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific Co., LTD., USA). Using the extracted DNA as a template, the V3–V4 region of the bacterial 16 s rRNA gene was amplified using the primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′; Xu et al., 2016). The fungal ITS region was amplified using the primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2R (5′-GCTGCGTTCTTCATCGATGC-3′; Adams et al., 2013). PCR amplification conditions and high-throughput sequencing were conducted as described in Li et al. (2021). Sequencing was performed on Illumina’s MiSeq PE300 platform (Shanghai Majorbio Bio-pharm Technology Co., LTD., China). The raw sequence data reported in this paper were deposited in the NCBI SRA database (serial numbers: bacteria, PRJNA836163; fungi, PRJNA836186).

Trimmomatic (version 0.33²) was used for quality control of the original sequences (Bolger et al., 2014), and FLASH (version 1.2.11³) was used for splicing (Magoc and Salzberg, 2011). The splicing sequence data were analyzed using UPARSE (version 7.1;⁴ Edgar, 2013), and sequences with a similarity of ≥97% were assigned to the same operational taxonomic unit (OTU). After quality control and concatenation of the original sequences of all samples, 1,088,672 and 976,526 high-quality sequences of bacteria and fungi were obtained, respectively. At the 97% sequence similarity level, the sequences clustered into 5,685 and 2,585 OTUs, respectively. For each representative sequence, the SILVA (bacteria⁵) and UNITE (fungi⁶) databases were used to annotate taxonomic information (Fan et al., 2020; Kang et al., 2022).

Mean value, standard deviation, and variance analysis of soil physicochemical properties and Perilla yield and quality were analyzed using Microsoft Excel 2010 and SPSS Statistics (version 20.0, IBM, USA). Differences between mean values were determined by one-way ANOVA and LSD post-hoc test (p < 0.05). All bioinformatic analyses were performed using R (version 3.6.1; https://cran.r-project.org/bin/windows/base/old/3.6.1/).

The alpha diversity (Sobs, Chao, and Shannon indices) of Perilla rhizosphere soil microbial communities was estimated based on OTUs. All indices were calculated by the “vegan” (Dixon, 2003) and “picante” (Kembel et al., 2010) packages in R (version 3.6.3). The linear discriminant analysis (LDA) effect size (LEfSe) method was used to assess potential bacterial and fungal biomarkers (from phylum to genus) within soil microbiomes that were specifically enriched under different fertilization management types based on p < 0.05 and an LDS score > 4.0 (Segata et al., 2011). Principal coordinates analysis (PCoA) was performed to calculate the gradient of compositional changes for bacterial and fungal microbial communities (based on Weighted-Unifrac distance matrix) using the ggplot2 package (Lozupone and Knight, 2005). Differences in bacterial and fungal communities between different samples were analyzed by Adonis test.

After variance inflation factor (VIF) analysis, pH, TC, TOC, TN, TK, AP, and AK with a VIF threshold less than 10 were selected for redundancy analysis (RDA) between environmental factors and soil microbial communities. Variance partitioning analysis (VPA) was used to quantitatively evaluate the individual and common explainability of environmental factor variables for microbial community differences. In addition, FAPROTAX (Louca et al., 2016) and FUNGuild (Nguyen et al., 2016) were used to analyze the ecological functions of soil bacteria and fungi, respectively. Kruskal–Wallis H test was used to test the significance of differences between groups.

To study the effect of different fertilization management techniques on the relationship of soil microbial communities, soil bacterial and fungal communities were combined based on fertilization management technique, and a soil microbial co-occurrence network based on genus classification was constructed. The co-occurrence network was constructed with genera that had a relative abundance greater than 0.1% based on random matrix theory (Deng et al., 2012). To simplify the networks for better visualization, a Spearman’s correlation between the two genera was considered statistically significant if the Spearman’s correlation coefficient (r) was >0.6 and the p value was <0.05. Moreover, p values were adjusted using the Benjamini–Hochberg FDR method (Benjamini and Hochberg, 1995). Spearman’s correlation and network attributes were calculated using the WGCNA, Psych, Igraph, and fdrci packages in R (version 3.6.3) (Csardi and Nepusz, 2006). The network attributes included the number of edges, average clustering coefficient, average degree, modularity, average path length, graph density, and betweenness centrality. Higher numbers of nodes and edges, graph density, average degree, and lower average path length indicate a more complex and connected network (Ma et al., 2016; Jiao et al., 2021). The higher the betweenness centrality value of microbial species, the greater the critical role of the species in the network. The Fruchterman–Reingold layout algorithm was used in the interactive platform Gephi (version 0.9.2⁷) for network visualization and network topology parameter calculation (Bastian et al., 2009). The network stability was evaluated by removing the nodes in the static network to estimate the speed of robustness decline, and the network robustness was evaluated by the natural connectivity of the nodes (Yuan et al., 2021; Zhu et al., 2022).

A structural equation model (SEM) was used to identify the direct and indirect effects of soil physicochemical properties (such as pH, TC, AK, and AP) on bacterial and fungal diversity (Shannon index), and P. frutescens yield (biomass per plant) and quality (linoleic acid content). To reduce SEM complexity, the representative indices of soil physicochemical properties were calculated by PCoA (Sun et al., 2021). All variables were standardized using Z-transformation (mean = 0, standard deviation = 1; Du et al., 2022).

The theoretical model assumed that: (1) soil pH has a direct impact on soil nutrient content, microbial community diversity, and P. frutescens quality or yield, (2) the soil available phosphorus, potassium, and total carbon had direct or indirect effects on the P. frutescens quality and yield and the soil microbial community, and (3) the P. frutescens yield has a direct effect on the soil microbial community. Model fitting was performed using root mean square error of approximation (RMSEA), probability level p value, Bentler comparative fit index (CFI), maximum likelihood goodness of fit (χ²), and degrees of freedom (df) tests (Du et al., 2022). The SEM was constructed using Amos Graphics (version 24.0, IBM Corp., USA; Liu, L. et al., 2019).

Compared with the CK group, the OMF group significantly increased soil TC and TK contents (p < 0.05; Table 2). Fertilization also significantly increased soil AK content (p < 0.05; Table 2). In addition, there were no significant differences in the physicochemical properties between the different fertilization treatment groups. The measurement results of different fertilization treatment groups indicated that short-term fertilization treatments may not significantly improve soil physicochemical properties. However, numerical analysis demonstrated that OF and OMF treatments tended to produce better results than no fertilization and CF treatment.

The analysis of P. frutescens yield during the harvest period showed that, compared with the CK group, the CF, OF, and OMF groups all significantly improved the yield indicators, including the number of stem nodes, number of effective branches at one time, and number of branch angles. However, there were no significant differences among the three fertilization treatments (Table 3). In addition, the number of panicles per plant and length of the main panicle in the OMF group were significantly greater than those in the other treatment groups.

Analysis of P. frutescens seed quality showed that the contents of oleic acid, linoleic acid, α-linolenic acid, and fat in the OF and OMF groups were significantly higher than those in the CK and CF groups. Additionally, the OMF group had the highest contents of oleic acid, α-linolenic acid, and fat (Figure 1). The crude protein content of the OMF group was significantly higher than that of the other treatment groups (p < 0.05). This indicated that the application of OFs, especially OMF, helped improve P. frutescens quality.

Alpha diversity analysis showed that different fertilization treatments did not significantly affect soil bacterial community richness and diversity, whereas OMF significantly increased soil fungal community richness (Figure 2). PCoA showed that there was no clear distinction between samples from different treatment groups on the PC1 and PC2 axes (Figures 3A,B); this indicated that there was no significant difference in soil bacterial community and fungal community composition among different fertilizer treatment groups in the short term.

Fertilization management did not cause significant changes in phylum-level bacterial composition (p > 0.05, ANOVA analysis; Figure 3C). Proteobacteria, Actinobacteria, Acidobacteria, and Chloroflexi were the main dominant bacterial phyla in each fertilization treatment group. Fertilization management also did not significantly affect the fungal community composition at the phylum level (p > 0.05, ANOVA; Figure 3D). Ascomycota, Mortierellomycota, Basidiomycota, and Glomeromycota were the dominant fungal phyla in each treatment group. The taxonomic analysis of the dominant bacterial genera with relative abundance >1% showed that different short-term fertilization treatments had no obvious effect on most dominant bacterial genera (Supplementary Table S1). However, the relative abundances of Bacillus and Candidatus_Solibacter significantly increased after fertilization, especially in the OF and OMF groups; the relative abundance of Bryobacter significantly decreased in the OF and OMF groups. The relative abundances of Clonostachys and Gonytrichum in the OMF group were significantly lower than those in the other groups, and the relative abundance of Penicillium in the OF group was significantly lower than that in the other treatment groups. In addition, there were no significant differences in other dominant fungal genera among the treatment groups (Supplementary Table S2).

The LEFSe results showed that Bryobacter was significantly enriched in the CK group; the nitrogen-fixing bacteria Saccharimonadia and iron-reducing bacteria Desulfobacca were significantly enriched in the CF group; Actinomadura, Nakamurellaceae and Nakamurella were significantly enriched in the OF group; and Iamia and Iamiaceae was significantly enriched in the OMF group (Figure 4A). The LEFSe analysis of the fungal community showed (Figure 4B) that Cyphellaceae, Piskurozymaceae, Filobasidiales, Solicoccozyma, and Funneliformis were mainly enriched in the CK group, and Mycosphaerellaceae, Didymosphaeriaceae, Paraphaeosphaeria, and Pseudopithomyces were mainly enriched in the CF group. Helotiales, Glomerales, Cladorrhinum, and Pseudallescheria in the OF group, and Phaeosphaeriaceae, Bulleribasidiaceae, and Sodiomyces were mainly enriched in the OMF group.

Analysis of the topological indicators of the bacterial community co-occurrence network revealed that the OMF group had the highest increase in the number of network edges and average degree (Table 4). The OF group had the highest modularity index, followed by the OMF group. This finding indicated that the application of OF and OMF could enhance the association of soil bacterial communities because the degree of modularity was higher. The robustness analysis results showed that the robustness of OF and OMF treatment groups were higher than that of CK and OF groups, and OMF groups was highest (Figure 5C).

Analysis of the network topology indicators of the fungal community showed that, compared with other treatment groups, the OMF group significantly increased the total number of edges in the network, number of positive correlation edges, number of nodes, graph density, and average degree (Table 4). The OMF group had the highest modularity index; this demonstrated that OMF can improve the association of fungal communities, and OMF can make the fungal community structure more modular and more stable. The robustness analysis results showed that the robustness of OMF and CF groups had little difference, but was higher than that of CK and OF groups (Figure 5D).

Analysis of the scale proportions of the top three modules in each group revealed that the proportion of the top three modules of bacterial and fungal communities was greatest in the OF group followed by the OMF group (Supplementary Figure S2). This indicated that the bacterial communities and the fungal communities of the OF and OMF groups were more closely related.

The main nodes in the bacterial community network of each treatment group belonged to Proteobacteria, Actinobacteria, Acidobacteria, Chlorobacteria, Bacteroidetes, and Firmicutes; the main nodes in the fungal community network belonged to Ascomycota, Basidiomycota, and Glomeromycota. This finding indicated that these bacterial and fungal phyla were keystone microbiota in all treatment groups (Figure 5). According to the analysis of the betweenness centrality values (Supplementary Table S3), the top 10 genus-level species were different in each treatment group, which indicated that the genus-level microbiota that played a key role in the co-occurrence network differed among treatment groups.

According to FAPROTAX functional analysis, chemoheterotrophy, aerobic chemoheterotrophy, nitrification, aerobic ammonia oxidation, and nitrogen fixation were the top five bacterial community functions in each treatment group (Figure 6A). Among the top 50 functions, only human pathogens all significantly differed among different treatments, and the OF and OMF treatments reduced the proportion of human pathogenic bacteria in soil compared with CK and CF treatments. There were no significant differences in other bacterial community functions among different treatment groups.

FUNGuild functional analysis indicated that Saprotroph, Saprotroph–Symbiotroph, Pathotroph–Saprotroph–Symbiotroph, Pathotroph, and Symbiotroph were the top five trophic modes of the fungal community (Figure 6B). Analysis of the relative abundance of arbuscular mycorrhizal fungi (AMF) in symbiotic trophic fungi showed that the AMF abundance was significantly reduced in the CF group (Supplementary Figure S3A); this indicated that CFs inhibited AMF growth. F_unclassified_o_Paraglomerales, Glomeraceae, Diversisporales_fam_Incertae_sedis, f_unclassified_o_GS24, and Paraglomeraceae were the dominant families in AMF, but there were significant differences in their relative abundances among different fertilization treatments groups (Supplementary Figure S3B).

RDA indicated that the selected environmental factors explained 38.20% of the total change in bacterial communities (Figure 7A) and 31.32% of total changes in fungal communities (Figure 7B). The results of RDA showed that pH (ANOVA, p = 0.021) and AP (ANOVA, p = 0.024) were the main factors affecting the soil bacterial community composition, and pH (ANOVA, p = 0.034) and AK (ANOVA, p = 0.011) were the main factors affecting the soil fungal community composition (Figure 7).

The SEM fit the measured data well (bacteria, χ²/df = 0.658, p = 0.764, CFI = 1.000, RMSEA = 0.000; fungi, χ²/df = 1.493, p = 0.135, CFI = 0.920, RMSEA = 0.181); this indicated high consistency between the hypothesized model and the observed data (Figure 8). The SEM showed that soil physicochemical properties accounted for 68% of P. frutescens quality, among which AK and TC were significantly positively correlated with P. frutescens quality, and pH and AP were significantly negatively correlated with P. frutescens quality (Figure 8). AK and AP were significantly positively correlated with bacterial diversity, and P. frutescens yield was significantly negatively correlated with bacterial diversity; these factors accounted for 68% of bacterial diversity (Figure 8A). AP was significantly positively correlated with and explained 61% of P. frutescens yield. TC and P. frutescens yield were directly and significantly positively correlated with fungal community diversity, whereas AK was directly and significantly negatively correlated with fungal community diversity (Figure 8B).

Compared with the CK group, OMF application significantly increased the TC and TK contents of the soil (Table 2). The experiment proved that although short-term application of OMF could not significantly improve most of the soil physicochemical properties, it could guarantee the same effect as the same amount of inorganic fertilizer. Numerical analysis showed that OF especially OMF application had a trend of improving soil fertility compared with no fertilization and CF application, which was consistent with previous continuous fertilization results (Zhao et al., 2016; Du et al., 2022).

The yield test results of P. frutescens at harvest showed that the number of panicles per plant and length of the main panicle were significantly greater in the OMF group compared with the other treatment groups (Table 3); this indicated that OMF increased P. frutescens yield to a certain extent. The quality inspection results of P. frutescens showed that OF and OMF treatments significantly increased the contents of unsaturated fatty acids, total fat, and total protein in P. frutescens, and the OMF treatment values were higher (Figure 1). These results were similar to those of a previous study that investigated the effect of OMF on Purslane growth, which also showed that OMF treatment increased unsaturated fatty acid content (Yang et al., 2020). The SEM results also demonstrated that soil AK content was significantly positively correlated with P. frutescens quality (Figure 8). Therefore, OMF application was beneficial for improving P. frutescens quality.

Microorganisms are the driver of soil fertility changes; they can directly indicate soil quality and play an important role in plant growth and crop yield (Fan et al., 2020; Finkel et al., 2020). In this study, short-term fertilization treatments did not significantly affect alpha diversity (Figure 2). Similar studies also showed that soil microbial alpha diversity was stable and not easily affected by agricultural management practices (Coller et al., 2019; Gui et al., 2021; Kang et al., 2022). How fertilization management affects soil microbial diversity depends on soil properties (Mendes et al., 2015), such as soil pH, which is generally considered a decisive factor underlying microbial diversity (Bissett et al., 2011).

Actinobacteria and Firmicutes are generally considered to be beneficial microorganisms for plants (Yang et al., 2017). Actinobacteria can control plant bacterial diseases by producing various antibiotics, secreting cell wall-degrading enzymes, and inducing host resistance (Conn et al., 2008; Chater et al., 2010; Liu et al., 2012). LEFSe analysis showed that the treatment of OF and OMF promoted the enrichment of some bacteria belong to Actinobacteria, including Actinomadura, Nakamurellaceae, Nakamurella, Iamiaceae, and Iamia (Figure 4). Actinomadura can produce several antibiotics that inhibit the growth of soil pathogens and reduce the occurrence of crop diseases and insect pests (Li et al., 2022). Bacillus and Candidatus_Solibacter were also significantly increased in OF and OMF groups (Supplementary Table S1). Bacillus is commonly formulated as biocontrol agents because they secrete antibiotics or antimicrobial proteins (Ahimou et al., 2000; Weller et al., 2002; Moyne et al., 2004), and improve soil fertility by increasing soil mineral nutrient availability (Chen et al., 2016). Additionally, Candidatus_Solibacter is a bacterium that decomposes organic matter (Rime et al., 2015).

LEFSe analysis showed that the main enriched fungal species (e.g., Glomerales, Cladorrhinum, and Pseudallescheria) in the OF groups belonged to the phyla Ascomycota and Glomeromycota (Figure 4B). Among them, Pseudallescheria is a biocontrol fungus; it is an important natural enemy of some plant parasitic nematodes that can parasitize eggs and infect larvae and females, and can significantly reduce the damage of plant nematode diseases such as those caused by root-knot, cyst, and stem nematodes of various crops (Wang et al., 1997; Ko et al., 2010; Zhu et al., 2020). Cladorrhinum is an effective biocontrol fungus for controlling the soil-borne Rhizoctonia solani pathogen (Liu, H. et al., 2019). Ascomycota is a key driver of the degradation of organic residues in soil (Richardson, 2009; Ma et al., 2013); therefore, the Ascomycota abundance may increase with increasing organic matter content (Du et al., 2022). Glomeromycota can undergo symbiosis with terrestrial plants to form arbuscular mycorrhizae, and this symbiosis can help plants absorb inorganic salts in soil, especially phosphorus (Smith and Read, 2008; Calaca and Bustamante, 2022). Our results showed that, OF treatments promoted increase in the number of beneficial fungi in the karst soil.

In conclusion, short-term fertilization treatments affected soil microbial communities, and OF and OMF had advantages over CF. However, compared with previous long-term experimental results, some differences in this study were not significant. Therefore, extending the fertilization period and intensifying fertilization may produce significant fertilization effects (Kox et al., 2020).

There is a complex association network among soil microbial communities, and they do not exist alone. When soil microbial community composition changes because of fertilization management, the microbial co-occurrence network also changes (Kang et al., 2022). The OF and OMF groups had higher modularity indices, which indicated that OF and OMF improved the soil microbial community connectivity and made the community connected more closely (Figure 5; Supplementary Figure S2; Table 4). This result is generally consistent with those of other studies on OF application (Ling et al., 2016; Wang et al., 2017; Liu et al., 2020; Kang et al., 2022).

Complex networks with higher connectivity are more tolerant of environmental disturbances than simple networks with lower connectivity (Santolini and Barabási, 2018). In this study, the network connectivity of both bacterial and fungal communities was highest in the OMF group (Table 4); this indicated that OMF treatment resulted in higher anti-interference ability of soil microorganisms. As the core components of soil organic matter degradation, bacteria and fungi usually form different functional groups and change the interaction between their ecological networks because of the decomposition or utilization of organic and inorganic nutrients (Wang et al., 2017; Dai et al., 2018; Samaddar et al., 2019), and they tend to maintain a complex network structure (Kang et al., 2022). The robustness analysis showed that the application of OF, especially OMF, could improve the stability of bacterial and fungal network structures (Figures 5C,D).

Keystone microbial groups play an important role in maintaining ecosystem homeostasis (Banerjee et al., 2018; Fan et al., 2020). According to the degree values of the network nodes, it was found that the key bacterial and fungal phyla in the nodes were not significantly different among the treatment groups (Figure 5); this was consistent with the composition and distribution of dominant species in the community (Figures 3C,D). At the phylum level, microorganisms had strong stability and were not easily affected by fertilization management. The genera with the highest betweenness centrality scores are generally considered keystone taxa (González et al., 2010; Vick-Majors et al., 2014). In this study, the keystone genera differed among fertilization treatment groups. The keystone genus of the bacterial community in the OMF group was Pseudomonas, and the keystone genus of the fungal community was Metarhizium (Supplementary Table S3), both of which are recognized as biocontrol microorganisms. Pseudomonas can adsorb heavy metals in soil and promote plant growth (Costa-Gutierrez et al., 2020; Ghorbanzadeh et al., 2020; Wu et al., 2022), whereas Metarhizium can kill plant pests and mitigate plant diseases (Riguetti Zanardo Botelho et al., 2019; Gebremariam et al., 2021; González-Pérez et al., 2022). In conclusion, this study demonstrated that the application of organic fertilizer, especially OMF, can not only enhance the connectivity, cohesiveness and stability of microbial community network structure, but also increase the abundance of beneficial microorganisms.

The RDA results showed that soil pH was the most important factor affecting bacterial and fungal communities (Figure 7). This results was consistent with the previous study (Cho et al., 2016; Gu et al., 2019). Changes in soil pH can alter soil structure, fertility, and vegetation communities, thereby directly or indirectly affecting soil microbial community composition (Lauber et al., 2009; Qi et al., 2018). In addition, in this study, AP was also the main factor affecting bacterial communities, and AK was the main factor affecting fungal communities (Figure 7). Fertilizers may be absorbed and used by plants after entering the soil, or they may remain in the soil, leading to changes in the composition of bacterial communities (Sun et al., 2016) and fungal communities (Zhou et al., 2016). Therefore, different fertilization treatments affected soil microbial community composition by mediating the effects of soil physicochemical properties, especially pH, and AP and AK.

SEM results showed that both soil TC and AK content was positively correlated to the quality (linoleic acid content) of P. frutescens (Figure 8), indicated that improve soil carbon and available potassium content could improve the quality of P. frutescens. There was a direct and significant positive correlation between P. frutescens yield and AP (Figure 8); this indicated that P. frutescens yield was mainly affected by the AP content in soil, and increasing the AP content can improve P. frutescens yield. Soil bacterial diversity was significantly positively correlated with soil AK (Figure 8A), but fungal diversity was significantly negatively correlated with soil AK (Figure 8B). Some previous studies showed that soil bacteria and fungi exhibited different patterns in response to fertilization treatments (Álvarez-Martín et al., 2016; Ai et al., 2018), and our study also supported this conclusion. Fungi are generally considered more closely related to plants and they are able to provide nutrients to plants in a symbiotic relationship (Chen et al., 2017), whereas bacteria are more affected by soil properties and environmental factors (Singh et al., 2008; Delgado-Baquerizo et al., 2016, 2018); our study obtained similar results (Figures 7, 8). In this study, we found that soil bacterial diversity was significantly negatively correlated with P. frutescens yield, and fungal diversity was significantly positively correlated with P. frutescens yield. This result may be opposite to many previous long-term fertilization studies. Because, in this study, short-term application of OF and OMF significantly increased P. frutescens yield (biomass per plant), but there was no significant change in bacterial diversity (Shannon index), and the value of OMF group even showed a downward trend (Figure 2). The application of organic fertilizer may stimulate the rapid growth of some dominant bacteria and beneficial bacteria in soil in a short time, while the abundance of some oligotrophic microorganisms that are not adapted to the existence of organic fertilizer will decline, resulting in the decrease of bacterial diversity. Through the correlation analysis among soil physicochemical properties, microbial communities and P. frutescens agronomic efficiency, it can be concluded that the application of OF or OMFs can improve the physicochemical properties of soil, especially the contents of total carbon, available potassium and available phosphorus, which can promote the quality and yield of P. frutescens.