Peanut-based intercropping systems altered soil bacterial communities, potential functions, and crop yield

Experimental site and design

The test site was in 2019 at Anjiazhuang Town, Feicheng City, Shandong Province, China (116°46′E, 35°57′N). The site has a continental monsoon climate with four distinctive seasons, sufficient light, and a warm climate. The mean annual temperature was 12.9 °C and the annual rainfall was 660 mm. The frost-free period was 200 d. And the annual sunshine hours totaled 2607 h. Physicochemical properties of the surface soil layer (0–20 cm) were the following: organic matter, 11.6 g kg⁻¹; total nitrogen, 1.0 g kg⁻¹; alkaline nitrogen, 75.9 mg kg⁻¹; available phosphorus, 28.0 mg kg⁻¹; available potassium, 168.6 mg kg⁻¹; pH, 7.0. The climate data was collected from a small, self-built weather station.

The field experiment was carried out using a randomized complete block design. In this experiment, sole millet (Setaria italica, labeled Millet), sole sorghum (Sorghum bicolor, labeled Sorghum), and sole peanut (Arachis hypogaea, labeled Peanut), and two intercropping patterns of sorghum/peanut intercropping of 3:4 (Fig. S1A, labeled Sorg/Pean) and millet/peanut intercropping of 4:4 (Fig. S1B, labeled Mill/Pean) were set. The field area of each plot was 0.45 ha, with three plot replicates per treatment. The tested cultivars of millet, sorghum, and peanut were “Jigu 20”, “Jiliang 3”, and “Huayu 36”, respectively. The compound fertilizer (N/P₂O₅/K₂O = 15:15:15) of 750 kg ha⁻¹ was broadcast on the field as basal fertilizer before sowing. The millet, sorghum and peanut are all sown simultaneously in mid-May. Millet and sorghum are harvested in early September, and peanuts in mid-September. The field management was implemented according to the practices of local farmers.

Sample collection and indices calculation

The soil samples were obtained during the harvest period of crops on 12 September 2019. A soil sampler was used to collect 0–20 cm of soil from each plot. Soil sampling was carried out at five points following the five-point sampling method. These soil samples were then blended into one sample. Thus, a total of 15 samples (five treatments with three replicates) were obtained.

After removing the impurities such as sand, stone, and root, fresh soil samples were transferred into aseptic cryopreservation tubes and stored in a −80 °C refrigerator to analyze soil bacterial community structure. Moreover, the remaining samples were dried naturally and used to determine physical and chemical properties.

At maturity, two ridges × 4 m (3.2 m²), two rows × 4 m (6.4 m²), and three rows × 4 m (7.2 m²) per plot were randomly harvested to examine the yields of peanut, millet, and sorghum, respectively.

The land equivalent ratio (LER) was calculated as follows (Mead & Willey, 1980): (1.1)${\displaystyle \mathrm{LER}={\mathrm{pLER}}_{\mathrm{A}}+{\mathrm{pLER}}_{\mathrm{B}}=\frac{{\mathrm{Y}}_{\mathrm{iA}}}{{\mathrm{Y}}_{\mathrm{sB}}}+\frac{{\mathrm{Y}}_{\mathrm{iB}}}{{\mathrm{Y}}_{\mathrm{sB}}}}$

where pLER_A and pLER_B indicate partial LER for crop A and crop B, Y_iA and Y_sB are the yield of crop A in intercropping and sole-cropping, and Y_iB and Y_sB are the yield of crop B in intercropping and sole-cropping, respectively. The LER value >1 indicates that intercropping system has yield advantage.

Soil physical and chemical properties

The total nitrogen (TN) content was measured using the micro-Kjeldahl procedure (Cao et al., 2017). The available nitrogen (AN), available phosphorus (AP), and available potassium (AK) content were measured by alkaline diffusion, NaHCO₃ leaching molybdenum-antimony reverse absorption spectrophotometry and ammonium acetate atomic absorption spectrophotometry, respectively (Cao et al., 2017). The soil organic carbon (SOC) was measured by potassium dichromate oxidation (Wang et al., 2015). And the pH was determined in a suspension with a soil: water ratio of 1:2.5 according to the potentiometric method (Li et al., 2016).

Soil DNA extraction, PCR amplification and sequencing

Total DNA was extracted from the tissue using E.Z.N.A.® soil DNA kit (Omega Bio-tek, Norcross, GA, US) according to the manufacturer’s instructions. The DNA quality was examined by 1% agarose gel electrophoresis. The concentration and purity of extracted DNA were determined using the NanoDrop2000.

The primer set of V338f (5′-ACTCCTACGGGAGGCAGCAG- 3′) and V806R (5′-GGACTACHVGGGTWTCTAAT-3′) was used to amplify the V3–V4 region of the bacterial 16S rRNA (Tang et al., 2020). The PCR products from the same sample were mixed and recovered on a 2% agarose gel. The recovered products were purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). The recovered products were detected by 2% agarose gel electrophoresis and quantified by Quantus™ Fluorometer (Promega, Madison, WI, USA). The sequencing library was constructed using the NEXTFLEX Rapid DNA-Seq Kit following manufacturer’s recommendations. Sequencing libraries were mixed with PhiX. Reads were removed during data demultiplexing (Demultiplexing refers to splitting multiplexed reads from different or the same lane based on the index, generates the sample’s corresponding fastq file). The Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) Illumina Miseq PE300 platform was used for high-throughput sequencing. Negative controls were used when performing amplification and were not used during sequencing.

Sequence analysis

The paired-end data from the sequencing were demultiplexed according to the barcode to remove the adapters and barcodes (each barcode was only used for a single sample). The adapter sequences at the 3′ and 5′ ends of the reads were quality clipped using the FASTP version 0.20.0 (https://github.com/OpenGene/fastp). A standard procedure was performed for quality control using DADA2 in Qiime2 version 2019.7.0 (https://qiime2.org/) to obtain the amplicon sequence variants (ASVs) (Bolyen et al., 2019; Wang et al., 2023). We removed any sequences having (i) an average base mass <20, (ii) a length <50, and (iii) N bases. The DADA2 plugin in QIIME2 was used for sequence quality trimming, denoising, merging, and chimera detection (Callahan et al., 2016). The consensus sequences for the ASVs were classified with a Naive Bayes classifier trained against the SILVA 16S rRNA gene reference (release 132) database.

Statistical and bioinformatics analysis

Statistical analysis was carried out using R (version 3.3.1; R Core Team, 2016). Alpha diversity indices (Shannon and Chao1 indices) were calculated by Mothur (version 1.30.). Nonmetric multidimensional scaling (NMDS) and analysis of similarity (ANOSIM) based on the Bray–Curtis similarity matrix were used to test the differences in the compositions of the bacterial communities among different cropping systems. Linear discriminant analysis (LDA) effect size (LEfSe), according to the cutoff value of LDA >3.5, was used to identify certain specific biomarkers between different groups (Segata et al., 2011). Redundancy analysis (RDA) was carried out to investigate the effect of soil physicochemical properties on bacterial abundance at the bacterial phylum and genus levels. The data were analyzed using SPSS 20 (IBM, USA). One-way analysis of variance (ANOVA) and least significance difference (LSD) test were performed to determine the differences of soil factors and microbial diversity in different cropping treatemnts.

Grain yield and LER

The total grain yield in the grain sorghum/peanut intercropping system was 8075 kg ha⁻¹, which increased by 29.2% and 10.3%, compared with the Peanut and Sorghum, respectively (Fig. 1A). The yield of Mill/Pean was slightly higher than Peanut but significantly higher than Millet. Additionally, the LER values of two intercropping systems were all greater than one, and the higher LER value was obtained in Sorg/Pean (Fig. 1B).

Soil physicochemical properties

Soil physicochemical properties of all samples are shown in Table 1. Intercropping practices changed the soil characteristics. Compared with sole millet, intercropping with peanut significantly increased the AN and TN content by 9.9% and 4.9%, respectively (p < 0.05). Compared with sole sorghum, intercropping with peanut significantly increased the AN and TN content by 16.0% and 10.3%, respectively (p < 0.05). The pH showed no significant difference among the five cropping systems. The AP was positively correlation with SOC, while the SOC, AP, and AK were negatively correlated with TN (Fig. S2).

Bacterial community diversity

A total of 4157 ASVs were detected in all samples. The number of sequencing reads after denoising and removing chimeras for each sample was shown in Fig. S3. Venn diagrams of observed ASVs in each treatment are shown in Fig. S4. The five treatments shared 683 ASVs (17%). Each treatment has its own unique ASVs. Sorghum had the most (23%), whereas Sorg/Pean had the least (16%). There were no statistically significant differences in bacterial community richness (Chao1) or alpha diversity indices (Shannon) among Sorghum, Millet, Peanut, Sorg/Pean, and Mill/Pean (Fig. S5).

Nonmetric multidimensional scaling (NMDS) analysis was carried out using the Bray–Curtis dissimilarities to reveal the soil bacteria beta diversity of five cropping systems. The results showed that the taxa of soil bacterial communities were distinct among five cropping systems (Stress = 0.04, Fig. 2A). Furthermore, the analysis of similarities (ANOSIM) also confirmed that the bacterial community differed among five cropping systems ( R = 0.858, p = 0.001, Fig. 2B).

Bacterial community structure

The taxonomic distribution of the phylum and the genus is presented in Fig. 3. At the phylum level, all treatments possessed similar dominant phyla of microbial communities. Actinobacteria (27%), Proteobacteria (23%), Acidobacteria (20%), and Chloroflexi (10%) were the most abundant phyla, followed by Firmicutes (3.8%), Gemmatimonadota (3.1%), Bacteroidota (3.0%) and Myxococcota (2.5%). Methylomirabilota and Nitrospirae were detected at low relative abundances (<2.0%). However, the abundance of dominant phyla varied among treatments. Sorg/Pean enhanced the relative abundance of phyla Actinobacteriota and Chloroflexi while reducing the relative abundance of Acidobacteriota, Bacteroidota, and Myxococcota compared to the corresponding monocultures. Mill/Pean enhanced the relative abundance of Proteobacteria, Acidobacteriota, and Nitrospirota while decreasing the relative abundance of Chloroflexi and Firmicutes. At the genus level, the relative abundances of the top ten dominant genera accounted for 17.4% of total sample sequences. The top ten genera were Sphingomonas (3.0%), RB41 (2.8%), Bacillus (2.5%), Gaiella (2.0%), Nocardioides (1.4%), Microvirga (1.3%), Streptomyces (1.3%), MND1 (1.2%), Nitrospira (1.2%) and Lysobacter (0.8%).

Taxonomic composition

Linear discriminant analysis (LDA) effect size (LEfSe) showed biomarkers of soil bacterial lineages in five cropping systems (Fig. 4). Cropping systems specifically affected the bacterial community from the phylum to the genus level. For example, Sorg/Pean had high abundances of Streptomyces and Microvirga. Mill/Pean soil increased the proportions of Kribbella and Reyranella. Candidatus_Solibacter, Nitrospira, and Sphingomonas were significantly enriched in Sorghum soil. Enrichment of Nocardioides and Gaiella was significant in Millet soil. In addition, Rubrobacter and Pontibacter were enriched in Peanut soil.

Correlation between soil bacterial groups and physicochemical properties

Redundancy analysis (RDA) (Fig. 5) and Pearson correlation analysis (Table 2) were carried out to determine the relation of the soil factors to microbial communities in soils. The pH and TN were positively correlated with each other, and the SOC, AK, and AP were positively correlated with each other. The SOC had the most significant impact on bacterial community composition. Other soil properties such as AN, AP, and TN were also important factors affecting the bacterial community structure.

At the phylum level, Actinobacteriota, Proteobacteria, Chloroflexi, Gemmatimonadota, Bacteroidota, and Nitrospirota among the top ten bacterial phyla were highly correlated with SOC. At the genus level, genera Sphingomonas, Nitrospira, and Lysobacter showed a positive correlation with SOC, while the opposite result was observed with the Bacillus.

Potential metabolic pathways

To understand the potential function of the soil microbial community, metabolic properties and pathways in different cropping systems were predicted by Tax4fun. Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG), 7834 KEGG orthologues (KO) were found across all samples, mainly belonging to metabolism, genetic information processing, environmental information processing, and environmental information processing, cellular processes, and human diseases. In second hierarchy level bacterial functions, the top abundant functional pathways (relative abundance >1%) are shown in Fig. S6. Ten pathways were observed for metabolism, three for genetic information processing, two for environmental information processing, two for cellular processes, and one for human diseases. Compared to other pathways, carbohydrate metabolism, amino acid metabolism, and membrane transport were overrepresented.

PLS-PM analysis

The result of partial least squares path modeling analysis (PLS-PM) showed that beta diversity and soil properties had significant impacts on yield (Fig. 6). Among them, soil bacterial community (beta diversity, phylum abundance, and genus abundance) had a direct impact on yield and indirectly influenced yield via changing soil properties.