The experiment was conducted at the Modern Agricultural Plantation of the Ningxia Academy of Agriculture and Forestry, Yinchuan, Ningxia, China (38°28’48” N, 106°13’31” E; altitude 1100 m). The field soil was sand-loam, and a 3-year continuous daylily cropping was implemented for the plantation. The climate of the temperate continental was semi-arid kind with the average yearly temperature (8.5°C), precipitation (200 mm), potential precipitation (2250 mm), and sunshine (2850 h), respectively. The frost-free season lasted for 185 days.

Watermelon (Citrullus lanatus), cabbage (Brassica pekinensis), and kale (Brassica oleracea) were intercropped with daylily. The experimental design included the following 5 treatments: (1) watermelon–daylily intercropping (WD); (2) cabbage–daylily intercropping (CD); (3) kale–daylily intercropping (KD); (4) row mixed intercropping in which watermelon, cabbage, and kale were alternately intercropped with daylily (MI); (5) monoculture of daylily (CK). A completely randomized design was adopted for field experimentation, with 3 plots (each with a surface area of 120 m²) per treatment. Daylily was planted on April 16, 2020, and later, on May 20, 2020, watermelon, cabbage, and kale were intercropped with daylily. The planting densities of different treatments are presented in Figure S1 . Drip irrigation was used for all crops, and the entire cultivation was managed conventionally.

The soil near the roots of daylily were collected from 15 test plots using an auger on August 30, 2020. After the random collection of 5 soil sub-samples (10–20 cm depth) from each experimental plot, all samples were integrated into a whole soil sample. Prior to laboratory transfer, this whole soil sample was sieved through a 2-mm mesh for complete homogenization and elimination of stones, roots, and plant residues. Each sample was divided into 2 aliquots, one of which was stored at an ambient temperature for use in the chemical assay while the other was cryopreserved at –80°C for later use in the biological analysis. The soil samples for chemical analysis were cryopreserved at –20°C after a week of air-drying treatment.

Soil pH was determined using a pH meter (PHS-25, Shanghai, China) at a constant soil/water ratio of 1:2.5. The available nitrogen (AN) concentration in the soil samples was determined using the DigiPREP TKN System (KJELTEC 8400, Foss, Denmark). A UV–Vis spectrophotometer (UH5300, North Points Ruili) was utilized for assessing the available phosphorus (AP) concentration in the soil. The available potassium (AK) in the soil was quantified using an inductively coupled plasma (ICP) spectrometer (Spectro Analytical Instruments, Kleve, Germany). The K₂Cr₂O₇-H₂SO₄ oxidation approach was adopted to assess the organic matter (OM) content in the soil.

Soil enzymes associated with nitrogen, phosphorus, and carbon degradation, including peroxidase (POD), sucrase (SC), and urease (UE), were evaluated for their activity. The soil UE activity was determined as described by Bao (2000) (Bao, 2000) using urea as the substrate. The spectrophotometric approach was adopted for the soil POD quantification in a 96-well microplate, using L-3,4-dihydroxyphenylalanine (L-DOPA) as the substrate (Bach et al., 2013). The soil SC activity was evaluated by determining the glucose discharge from an 8% sucrose solution following 24 h of incubation at 37°C (Chen et al., 2010).

The total genomic DNA was extracted from the soil samples using sterile cotton swabs and the FastDNA spin kit for soil (MP, USA) according to the provided protocol. The extracted genomic DNA was amplified and sequenced by the Gene Denovo company. Li and Wu’s report (2018) was referred to for designing the specific primers ITS4–2409R and ITS3_KYO2F for the ITS2 zone amplification of the ITS rRNA gene, and the specific primers 341F and 806R for the V3-V4 zone amplification of the 16S rRNA gene (Li and Wu, 2018). In order to sequence the purified amplicons, paired-end sequencing (HiSeq 2500, PE250) was performed on the Illumina platform following standard operations.

Quality filtration and sequence fusion were conducted on the raw fastqfiles of 50 microbial samples using FASTP version 0.18.0 and FLSAH version 1.2.11, thereby deriving Tags. The next step was to filter out the low-quality tags following the Tags Quality Control process of QIIME Ver. 1.9.1, which generated high-quality Clean Tags (Edgar, 2013). The chimeras of the tags were recognized and deleted via the UCHIME algorithm by exploiting the reference database, thereby generating valid Clean Tags (Edgar et al., 2011). Finally, using UPARSE version 9.2.64, the Clean Tags were assigned to the identical Operational Taxonomic Units (OTUs) in the similarity setting of tingil (Edgar, 2013). In this process, for every OUT, the tag sequence with the maximum richness was selected as the typical sequence, and for every typical sequence, the classification information was annotated via the Mothur algorithm in the Silva database (Quast et al., 2013).

The “Vegan” R software version 2.5–6 was employed for estimating the alpha- and beta-diversities (Bray–Curtis dissimilarity). SPSS 25.0 (IBM, USA) was used for conducting the one-way analysis of variance (ANOVA) and the least significant difference (LSD) approach (p < 0.05) to determine whether the differences in the inter-sample diversities based on the alpha index were statistically significant. The distribution trends of the soil microorganism communities in different treatments were assessed using the ‘labdsv’ package for the principal coordinates analysis (PCoA). Analysis of similarities (ANOSIM) and permutational multivariate analysis of variance (ADONIS) with 999 permutations were used for determining the significant inter-sample group differences in the beta diversity using Bray–Curtis distance matrices. The taxonomic traits of the microbial species characterizing the inter-treatment differences (https://www.cloudtutu.com) were elucidated through the linear discriminant analysis effect size (LEfSe) assay. In order to simplify the LEfSe process, only the OTUs with a relative richness of >0.01% were selected through OTU table filtration screened. The taxa with significant differential inter-group richness levels were recognized using the factorial Kruskal–Wallis sum-rank test (α = 0.05). Thereafter, the effect size was computed for every discriminative trait based on the logarithmic LDA score (threshold = 2.0). The associations of the environmental parameters with the microbial communities were explored using the “vegan” R software for redundancy analysis (RDA) and the Mantel test. Using the same”vegan” package, the associations of the OTUs with the physicochemical properties and enzymatic activities were visualized through RDA.

Among all treatments, prominently higher levels of AN, AK, AP, and OM (LSD, p < 0.05, Table 1 ) were noted in the intercropping soil systems (WD, CD, KD, and MI) compared to that in the monocropping system (CK). The nutrient levels were higher in WD compared to CK by 12.90% for AN, 15.8% for AP, 35.71% for AK, and 30.95% for OM. The nutrient levels were higher in CD compared to CK by 39.52% for AN, 62.56% for AP, 26.47% for AK, and 34.53% for OM. The nutrient levels were higher in KD compared to CK by 17.21% for AN, 35.90% for AP, 14.90% for AK, and 20.46% for OM. The nutrient levels were higher in MI compared to CK by 31.73% for AN, 3.85% for AP, 2.03% for AK, and 19.08% for OM ( Table 1 ). Among the different intercropping soil systems, CD presented the highest concentrations of soil AN, AP, and OM, while WD presented the highest soil AK concentration ( Table 1 ). The WD, CD, and KD systems were less different from CK in terms of soil pH while presenting significantly increased MI ( Table 1 ). Similarly, soil enzyme activities differed significantly. The UE and SC activities were significantly higher in the intercropping soil systems (WD, CD, KD, and MI) compared to the monocropping soil system (CK) ( Table 1 ). The POD activity was significantly higher in WD and MI compared to CK, while it was significantly lower in CD and KD ( Table 1 ). For yield, the yield of daylily in all intercropping soil systems was higher than that in monocropping systems, with WD having the highest yield (30.46%), followed by MI (18.35%), CD (12.81%) and KD (7.43%) ( Table 1 ).

Approximately 1,665,552 superior-quality V3-V4 sequences and 1,501,148 superior-quality ITS2 sequences were yielded in total in the raw-read sequencing process. The mean read lengths for bacteria and fungi were 456 bp and 357 bp, respectively. The sequences of all samples were clustered into 93,898 bacterial OTUs and 8,372 fungal OTUs, using an identity threshold of 97%. The flat trend of the rarefaction curve ( Figure S2 ) suggested that the desired overall coverage of OTUs was provided by deep sequencing. Figure 1 depicts the calculated alpha diversity (α-diversity) indices of the bacterial and fungal communities for the different intercropping modes based on the Shannon index metric, which incorporates both the abundance and evenness of the microbiota. In the case of the bacterial microflora, the Shannon diversity values were in the following order: KD > CD > MI >WD > CK; in comparison to CK, the Shannon diversity values of CD and KD were pronouncedly higher (LSD, p < 0.05, Figure 1A ), while no significant differences (p > 0.05, Figure 1A ) were observed among WD, MI, and CK. In the case of the fungal microflora, the Shannon diversity values were in the following order: MI > CD > KD > WD > CK; in comparison to CK, the Shannon index of MI was pronouncedly higher (p < 0.05, Figure 1B ), while no significant differences were observed among CD, KD, WD, and CK (p > 0.05, Figure 1B ). The microbiota differences among the different treatments were assessed using the Bray-Curtis distance-based principal coordinate analysis (PCoA). In the case of bacterial microbiota, a sharp separation was observable among CD, WD, KD, MI and CK, while KD and MI were partially overlapped ( Figure 1C ); as demonstrated by the PCoA outcomes, the first 2 axes explained 36.3% and 26.05%, respectively, of the overall variation in the bacterial microflora. In the case of fungal microbiota, a distinct separation was observed among KD, WD, MI, and CK, while CD was overlapped with both KD and WD ( Figure 1D ); the PCoA outcomes revealed that the first 2 axes explained 25.31% and 16.91%, respectively, of the overall variation in the fungal microflora. The results of ANOSIM and ADONIS also indicated significant differences in the bacterial and fungal microbiota in the soil between a minimum of two treatments ( Table S1 ). The Venn diagrams verified that the variation in the microflora among different soil intercropping systems was due to the alterations in the composition of several unique and shared OTUs. Among all the OTUs identified in this study, the bacteria shared across all intercropping soil systems had 2679 OTUs, the WD had 463 unique OTUs, the CD had 481 unique OTUs, the KD had 470 unique OTUs, the MI had 579 unique OTUs, the CK had 527 unique OTUs ( Figure 1E ). Among all the OTUs identified in this study, the fungi shared across all intercropping soil systems had 160 OTUs, the WD had 69 unique OTUs, the CD had 54 unique OTUs, the KD had 63 unique OTUs, the MI had 89 unique OTUs, the CK had 80 unique OTUs ( Figure 1F ).

Figure 2 depicts the outcomes of the sequence analyses conducted at the genus and phylum levels. In the case of bacteria, the richest phylum was Proteobacteria (35.93%), followed by Gemmatimonadetes (16.51%) and Acidobacteria (12.99%), respectively. Among the dominant taxa, the abundances of Proteobacteria, Gemmatimonadetes, Planctomycetes, Actinobacteria, Patescibacteria, and Firmicutes did not change significantly in the intercropping soil systems (CD, WD, KD, and MI) compared to their abundances in the monocropping soil system (CK) ( Table S2 , Figure 2A ). Rokubacteria in KD and Acidobacteria in WD and CD were significantly less abundant compared to CK. A significantly greater abundance of Chloroflexi in WD and Bacteroidetes in MI was observed compared to CK ( Table S2 , Figure 2A ). The richest genera in the different intercropped soils were Sphingomonas (3.94%), Lysobacter (3.33%), RB41 (2.51%), Subgroup_10 (1.49%), Lactobacillus (1.34%), MND1 (1.33%), SWB02 (1.29%), Dongia (0.90%), SM1A02 (0.80%), and Arenimonas (0.67%) ( Figure 2C , Table S2 ). Sphingomonas, MND1, SWB02, Dongia, and Arenimonas were not significantly altered in the intercropping systems (WD, CD, KD, and MI) compared to the monocropping system (CK) ( Table S2 ). A significantly greater abundance of Lysobacter in CD, KD, and MI and that of Lactobacillus in WD was observed compared to CK. RB41 in WD and CD, Subgroup_10 in WD, and SM1A02 in MI were significantly less abundant compared to CK ( Table S2 ).

In the case of fungi, the richest phylum was Ascomycota (89.05%), followed by Mortierellomycota (6.11%) and Basidiomycota (3.88%). Among the dominant taxa, the abundances of Ascomycota, Chytridiomycota, Glomeromycota, Rozellomycota, and Aphelidiomycota did not change significantly in the 4 investigated intercropping systems compared to those in the monocropping soil system (CK) ( Table S3 , Figure 2B ). A significantly greater abundance of Mortierellomycota in WD and CD, Basidiomycota in KD and MI, and Mucoromycota in CD was observed compared to CK ( Table S3 , Figure 2B ). The richest genera in the different intercropping soil systems were Fusarium (10.65%), Pseudaleuria (6.92%), Mortierella (6.11%), Kotlabaea (5.42%), Schizothecium (5.06%), Coprinellus (2.11%), Peziza (1.88%), Madurella (1.46%), Acaulium (1.13%), and Microascus (0.92%) ( Figure 2D , Table S3 ). Kotlabaea and Schizothecium were prominently richer, while Fusarium and Peziza were pronouncedly less abundant in the 4 investigated intercropping systems compared to their abundances in the monocropping soil system (CK); the differences in the abundance of Microascus were not significant between the two kinds of systems ( Table S3 ). Significantly greater abundances of Mortierella in WD, CD, and KD, Pseudaleuria in KD, Coprinellus in KD and MI, and Madurella and Acaulium in WD were observed compared to CK. Significantly lower abundances of Pseudaleuria in WD and CD and Madurella in CD, KD, and MI were observed compared to CK ( Table S3 ).

Further, the effects of the different intercropping patterns on soil microbial community composition were determined using the linear discriminant analysis effect size analysis (LEfSe). LEfSe results revealed seventy-seven bacterial taxa that differed significantly in terms of relative richness among the different intercropped soil systems. As depicted in Figure 3A , 12 bacterial taxa in WD, 18 bacterial taxa in CD, 13 bacterial taxa in KD, 11 bacterial taxa in MI, and 23 bacterial taxa in CK were identified as biomarkers. As revealed by LEfSe, 60 fungal taxa differed significantly in terms of relative richness among the different intercropped soil systems. Among the identified biomarkers were one fungal taxon in WD, 4 fungal taxa in CD, 15 fungal taxa in KD, 26 fungal taxa in MI, and 14 fungal taxa in CK ( Figure 3B ).

Redundancy analysis (RDA) revealed a relationship between the microbial community structure, soil properties (soil chemistry and enzymatic activity) and daylily yield. In the case of bacteria, AK, AP, AN, and OM were closely associated with WD and CD; UE, pH, and SC were closely related to KD and MI; POD was closely related to CK; Yield was closely related to WD and CD, and was also affected by AK, AP, AN and OM ( Figure 4A ). The ANOVA-like permutation test revealed that OM, UE, SC, AN, AK, Yield, and AP were significantly associated with the bacterial communities ( Table S4 , p-value < 0.01), while pH and POD did not exhibit a significant association ( Table S4 , p-value > 0.05). In the case of fungi, AP, AK, and OM were closely related to WD, CD, and KD; AN, UE, and pH were closely related to MI; POD and SC were closely related to CK; Yield was closely related to WD, CD and KD, and was also affected by AK, AP, AN and OM ( Figure 4B ). The ANOVA-like permutation test revealed that AP, AK, and SC were significantly associated with the fungal communities ( Table S4 , p-value < 0.05), while pH, OM, AN, UE, POD, and Yield did not exhibit a significant association ( Table S4 , p-value > 0.05).

The co-occurrence patterns between the microbiome and soil properties were investigated using network analysis. As depicted in Figures 4C, D , overall, the module was more densely connected in the case of bacteria (average degree = 1.72) compared to fungi (average degree = 1.28). Several microbial taxa were associated significantly with various soil parameters, particularly AN, which exhibited high connectivity. For instance, Blastocatella was linked positively to SC and negatively to AN, AP, and OM, Gemmatimonas was associated positively with SC and POD and negatively with OM and AK ( Figure 4C ), Fusicolla was associated positively with AN, AP, and UE, and Stachybotrys was associated negatively with OM, AN, AP, and UE ( Figure 4D ).