Effects of wine-cap Stropharia cultivation on soil nutrients and bacterial communities in forestlands of northern China

Experimental site

The experimental forestland was an area of 20 × 150 m, located in Liying Town, Jining City, Shandong Province (116°37′E, 35°30′N, 43 m above sea level). The nursery stock is made up by 7- to 10-year-old trees of horse chestnut (Aesculus chinensis Bunge), which were planted with 2 m spaces between plants in rows and 5 m between rows to achieve a canopy density of 0.7. This location is considered a warm temperate, semi-humid monsoon climate characterized by hot, rainy summers and cold, dry winters, with an annual average temperature of 13.2–14.1 °C. The highest temperature in July exceeded 27 °C, and the annual average temperatures above 10 °C accumulated to 4060.7 °C (growing degree days). The annual precipitation is 650–700 mm, with rainfall from May to August accounting for more than 65% of the total rainfall for the whole year. The soil type was non-calcareous cinnamon tide with a clay loam texture. All these data were obtained from the Jining Soil and Fertilizer Workstations, China (1990).

Sample plots and S. rugosoannulata cultivation

The experimental forestland was divided into five 20 × 30 m grids, which were marked as Y000, Y010, Y011, Y001, and Y101 respectively. Among them, Y010, Y011, Y001, and Y101 were used for mushroom cultivation with different regimes, and each of them was divided into three, 10 × 20 m plots for independent replicates; Y000 was used as a no-cultivation control. The cultivation year of each grid is shown in Table 1.

The cultivation of S. rugosoannulata began in 2013 and was performed every November as described by Gong et al. (2016). The basic materials included 48.9% rice husk, 30% corncobs crushed into particles with a diameter of 0.5 cm–1 cm, 20% sawdust, which was a mixture that contained a variety of hardwood chips, 1% soil acquired from each plot before cultivation and 0.1% lime. These materials were mixed, stacking fermentation was performed, and then distributed onto the sample plots between the plant rows with a thickness of approximately 25 cm. The S. rugosoannulata spawn was divided into blocks of approximately 3 cm in length and inoculated into the fermented material using superimposed square planting. Then, 3 cm of the forest surface soil was sprinkled onto the surface of the fungal bed. The fungal bed was vented and kept moist by a 2–3 cm cover of straw under black plastic film. A micro-spray system was installed in each plot, and the ditch between the cultivation beds drained into a stagnant water well. By April of the next year, fruiting had begun, and by late June, the harvest was complete. The soil was subjected to rotary tillage in November, i.e., the material rotting stage (MRS).

Sample collection and measurements of soil properties

A five-point sampling method was used to collect soil samples in October 2016. The surface organic materials of Y011, Y001, and Y101, and 1 cm of the surface soil of Y000 and Y010 were removed to distinguish the effect of the raw organic materials added from the mushroom cultivation. Five soil cores (5 cm diameter) were collected from each plot with a depth of 30 cm, fully pooled and then sifted using a 2-mm sieve. Subsequently, each soil sample was divided evenly into two portions: one was air dried and used for soil nutrient measurements, and the other was stored at −20 °C before soil DNA extraction.

The soil properties were measured in the Shandong Provincial Key Laboratory of Soil Erosion and Ecological Restoration (Tai’an, Shandong, China). The soil organic matter (OM) content was determined with the potassium dichromate external heating method (Ciavatta et al., 1991). The total nitrogen (TN) content was determined by the dichromate oxidization method (Bremner, 1965). The total phosphorus (TP) content was determined by molybdenum-blue colorimetry after digestion by HF-HClO4 (Jackson, 1958). The alkaline hydrolysable nitrogen (AN) content was determined using the alkaline-hydrolysable diffusion method (Xiong et al., 2008). The available phosphorus (AP) was extracted with sodium bicarbonate and determined using the molybdenum-blue method (Olsen, 1954). The available potassium (AK) was extracted by ammonium acetate and then determined by flame photometry (Carson, 1980). The soil pH was determined according to the international standard with a soil/water ratio of 1:5 (ISO 10390: 2005). The soil field capacity (FC) was measured using the laboratory Wilcox method (Duan et al., 2010).

Soil DNA extraction and polymerase chain reaction (PCR) amplification

The hexadecyl trimethyl ammonium bromide (CTAB) method was used for the soil DNA extraction (Zhou, Bruns & Tiedje, 1996), and the purity and concentration of genomic DNA was monitored by 1% agarose gel electrophoresis. DNA was diluted to 1 ng/µL using sterile water for the PCR. The forward specific primer 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) (Turner et al., 1999) and reverse specific primer 907R (5′-CCGTCAATTCMTTTRAGTTT-3′) (Lane, 1991) were employed to amplify the V4–V5 region of 16S RNA. PCR-based amplifications were performed using Phusion® High-Fidelity PCR Master Mix with GC Buffer and high-fidelity DNA polymerase (New England Biolabs, Ipswich, MA, USA) following an amplification programme of 1 cycle at 98 °C for 1 min, 30 cycles composed of three steps for each cycle (98 °C for 10 s, 50 °C for 30 s, and 72 °C for 30 s), and a final elongation step of 72 °C for 5 min.

Equal volumes of 1× loading buffer (containing SYBR green) and PCR products were mixed and electrophoresed on a 2% agarose gel. Samples with bright main bands between 400 and 450 bp were selected for further experimentation. The PCR products were mixed in equidensity ratios and then purified with a Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany). The library was constructed using TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, USA), and the library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 systems. The library was sequenced on an Illumina HiSeq 2500 platform at Novogene Bioinformatics Technology Co., Ltd., Beijing, China, and 250 bp paired-end reads were generated. All paired-end reads were deposited in Sequence Read Archive (SRA), BioProject: PRJNA453134.

Bioinformatic analysis

Paired-end reads were assigned to samples based on their unique barcodes and truncated by trimming the barcode and primer sequences. After that, paired-end reads were merged using FLASH (V1.2.7; Magoč & Salzberg, 2011) to obtain raw tags. The raw tags were then subjected to quality filtering using QIIME V1.7.0 (Caporaso et al., 2010) to obtain high-quality clean tags (Bokulich et al., 2013). Default settings (r = 3; p = 0.75 total read length; q = 3; n = 0; Sun et al., 2014) was used for quality filtering. These clean tags were compared with the reference database (Gold database, http://drive5.com/uchime/uchime_download.html) using the UCHIME algorithm (Edgar et al., 2011) to detect and remove chimaera sequences (Haas et al., 2011). Thus, we obtained effective tags. Uparse software (v7.0.1001; Edgar, 2013) was used to assign sequences with more than 97% similarity to an operational taxonomic unit (OTU). Representative sequences that showed the highest frequency for each OTU were screened for further taxonomic assignment. The Mothur method with a threshold of 0.8–1 was selected in QIIME (Version 1.7.0), and the SSU rRNA database (Quast et al., 2013) in SILVA (Wang et al., 2007) was used for taxonomic assignment. To obtain the phylogenetic relationships among different OTUs, multiple sequence alignments were conducted using MUSCLE software (Version 3.8.31; Edgar, 2004). The phylogenetic tree for each sample plot was visualized using GraPhlAn (Asnicar et al., 2015).

The OTU abundance data were rarefied using a standard sequence number corresponding to the sample with the fewest sequences. Subsequent analyses of the alpha diversity and beta diversity were performed based on the rarefied output data. The alpha diversity indices, including Good’s coverage estimator and the Shannon and Simpson diversity indices, were calculated using QIIME (Version 1.7.0). The differences in taxonomic composition were evaluated using a beta diversity analysis. The methods of principal component analysis (PCA), principal co-ordinates analysis (PCoA) and non-metric multi-dimensional scaling (NMDS) were used to illustrate the clustering of different samples. PCA was calculated in the R packages FactoMineR (Lê, Josse & Husson, 2008) and ggplot2 packages (Wickham, 2010), and the Hellinger transformation method (Rao, 1995) was used for PCA. PCoA of the weighted and unweighted UniFrac distances was calculated in the R package “ape” (Lozupone & Knight, 2005). An NMDS of the weighted and unweighted UniFrac distances was calculated according to Peck (2010). A canonical correspondence analysis (CCA) calculated using the R package “vegan” (Oksanen et al., 2007) was used to visualize the relationship between edaphic factors and the bacterial community structure in each sample plot. Prior to performing the CCA, we filtered out the intercorrelated environmental factors that affected sample distribution by using a variance inflation factor (VIF) analysis (Gross, 2003).

Statistical analysis

The soil chemical concentration, dominant taxa and alpha diversity indices were measured, and a one-way analysis of variance (ANOVA) was performed to determine whether differences existed among treatment means at a significance level of α = 0.05. Multiple comparisons were conducted for significant effects using the Tukey’s test at α = 0.05, and FDR of Benjamini & Hochberg (1995) was used for Tukey’s test. These statistical analyses were implemented using the Statistical Program for Social Sciences SPSS (Version 22; IBM, USA).

The linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011) was used to identify significantly different taxa among groups using the LEfSe software with a default LDA score value of 4. An analysis of molecular variance (AMOVA, Excoffier, Smouse & Quattro, 1992), analysis of similarities (ANOSIM, Clarke, 1993) and permutational multivariate analysis of variance (PERMANOVA or ADONIS, Anderson, 2001) was used to determine differences in the microbial community structure between the groups using the amova function in Mothur software (https://www.mothur.org/). Correlations among edaphic factors with estimated diversity levels were tested for significance via Spearman’s correlations (Algina & Keselman, 1999) performed in R ( R Core Team, 2013).

Soil properties

As shown in Table 2, the cultivation of S. rugosoannulata in forestland changed the soil field capacity, pH, organic matter, total nitrogen, total phosphorus, alkaline hydrolysable nitrogen, available phosphorus and available potassium contents. The ANOVA showed that the soil organic matter and available phosphorus increased significantly in all cultivating regimes of S. rugosoannulata compared with the no-cultivation control. The soil nutrient concentrations in the one-year interval cultivation regime were the highest among all grids. Additionally, the soil total phosphorus and alkaline hydrolysable nitrogen in the fallow for one year after prior cultivation regime, the soil total phosphorus, alkaline hydrolysable nitrogen, and available potassium in the two years of continuous cultivation regime and the soil total nitrogen and alkaline hydrolysable nitrogen in the current-year cultivation regime decreased significantly compared with those of the control. In addition, the soil field capacity and pH in all cultivating regimes were not significantly changed (P < 0.05).

Bacterial community composition

After removing potential chimaeras, a total of 1,127,888 high-quality V4–V5 16S rDNA sequences were analysed across the five grids. These sequences were assigned to 8,751 OTUs. The number of OTUs in the grids ranged from 4,756 to 5,011 (Table S1).

The phylogenetic relationship of different OTUs of each sample was illustrated in Figs. S3–S17. The top ten most abundant phyla represented 94% of the sequences (Fig. 1A), of which, Proteobacteria and Acidobacteria were the most dominant phyla in all groups, representing 55–61% of the total sequences. Among them, only Acidobacteria, Planctomycetes and Gemmatimonadetes showed significant changes in relative abundance between forestlands with one of the cultivation regimes and the no-cultivation control (Fig. 1A, Table S2). Acidobacteria in the current-year cultivation regime was the most abundant, and was significantly more abundant than in the control. The abundance of Planctomycetes in forestland with cultivation was greater than that of the no-cultivation control and was greatest in the one-year interval cultivation regime. Significantly fewer Gemmatimonadetes were observed in the one-year interval cultivation regime than in the other sample groups.

The LEfSe analysis was used to identify the specific bacterial groups in the soil from forestland with the cultivation regimes and in the no-cultivation control. Major differences were observed in the bacterial groups among the samples. Notably, Acidobacteria was the most common group in the soil with cultivation (Figs. 1B–1E). The most frequently observed differences were between the one-year interval cultivation regime and the no-cultivation control (Table S3).

Bacterial α-diversity

Before performing the α-diversity analysis, the OTU abundance data were normalized with a cutoff value of 59,458. In all samples, the Good’s coverage values reached 0.98 (Table S1), indicating that the normalised sequencing data was sufficient to capture the bacterial diversity. The Shannon and Simpson indexes were calculated to evaluate the bacterial diversity (Table S1), and no significant difference was observed among the five grids (p = 0.05), even though slightly lower Shannon and Simpson values were observed in the forestlands with cultivation than in the no-cultivation control.

OTU-level bacterial β-diversity analysis

The PCA of the bacterial community construction in different samples is shown in Fig. 2A. The five treatments were clearly distinguished in the PCA. The first two principal components, PC1 and PC2, best reflected the differences between these treatments and represented variations of 12.35% and 10.04% in the bacterial community, respectively. Within the PC1 axis, the one-year interval cultivation regime was distinct from the other regimes. Within the PC2 axis, the no-cultivation control was distinct. Similarly, the weighted Unifrac-based analysis of PCoA and NMDS (Fig. S1) and unweighted Unifrac-based analysis of PCoA and NMDS (Fig. S2) all showed that the one-year interval cultivation regime was distinct from the other regimes. These data indicated that the bacterial community composition in the one-year interval cultivation regime was relatively distinct from other regimes. However, significant differences in the bacterial community composition were not found between the no-cultivation control and the one-year interval cultivation regime via the AMOVA (Fs = 4.05, P = 0.074), the ANOSIM (R = 1, P = 0.1) or the ADONIS (R2 = 0.52, P = 0.1) (Table S4).

The VIF analysis suggested that the soil field capacity, pH, organic matter, total phosphorus, alkaline hydrolysable nitrogen were the uncorrelated edaphic factors in the CCA that could represent the relationship between soil physicochemical properties and bacterial community composition. Based on this model, a total of 60.39% of the variance was explained by CCA1 (24.46%) and CCA2 (19.49%), which were the first two constrained axes of the CCA (Fig. 2B). The CCA suggested that the organic matter, total phosphorus, and alkaline hydrolysable nitrogen were the determinants among the edaphic factors. The correlation analysis showed that only organic matter was significantly associated with the soil bacterial composition (r =  − 0.579, P = 0.024) and the Shannon index (r =  − 0.571, P = 0.026).

Effects of cultivating S. rugosoannulata under nursery stock shade on soil properties

Recently, growing S. rugosoannulata under nursery stock shade has been considered a win-win agricultural practice that can improve the quality of nursery stock soil in China. In this study, higher organic matter and available phosphorus content were observed in forestlands with cultivation of S. rugosoannulata compared with the no-cultivation control. However, other soil nutrients did not increase consistently in the forestlands with cultivation; some significantly decreased (Table 2). These results were unexpected because a positive correlation has been observed between the organic matter content and soil fertility (Sadikhani, Sohrabi & Fard, 2014). However, the acute angles between the arrow line representing organic matter content and the arrow lines representing other nutrients in the CCA (TN, TP, AN, AP and AK; Fig. 2B) indicate that organic matter content was positively correlated with the other nutrients in this study also. These results are consistent with those of other studies (Sihi et al., 2017; Zhou et al., 2017).

Cultivation in temperate climates usually results in a significant loss of mineralized organic N in soil (Tiessen, Cuevas & Chacon, 1994). Although farming S. rugosoannulata in forestland is a type of agricultural practice, it is distinct from traditional crop cultivation. The decrease in alkaline hydrolysable nitrogen content in forestland under cultivation (except for the one-year interval cultivation regime) indicated that the following cultivation regimes resulted in nitrogen loss: fallow for 1 year after prior cultivation regime (Y010), two years of continuous cultivation regime (Y011) and current-year cultivation regime (Y001). The one-year interval cultivation regime (Y101) effectively suppressed the nitrogen loss and significantly increased the alkaline hydrolysable nitrogen content. In addition, the one-year interval cultivation regime performed well in maintaining soil fertility and had the highest soil nutrient content (Table 2). In contrast, the two years of continuous cultivation regime resulted in a loss of soil nutrients with a significant decrease in total phosphorus, alkaline hydrolysable nitrogen, and available potassium content (Table 2). The current-year cultivation regime resulted in a significant decrease in total nitrogen and alkaline hydrolysable nitrogen content.

Soil use and management, such as less intensive management, can cause the loss of phosphorus (Leinweber et al., 1999). Herein, total phosphorus loss was found in the fallow for one year after prior cultivation regime, two years of continuous cultivation regime, and current-year cultivation regime. This loss may have resulted from the use of a large amount of water during the fruiting stage of S. rugosoannulata. In contrast, a significant increase in total phosphorus was observed in the one-year interval cultivation regime, indicating that the amount of phosphorus increase was greater than the amount of phosphorus lost in the one-year interval cultivation regime. We hypothesize that spent mushroom compost left after the harvesting of fruiting bodies would add a certain amount of macronutrients like nitrogen, phosphorous and potassium (NPK) (Kim et al., 2011). However, this hypothesis needs to be further tested in the future.

It has been suggested a higher organic matter content may lower the soil pH (Hodges, 1996) and increase the water content at field capacity (Hudson, 1994; Tale & Ingole, 2015). Inconsistent with this, the two years of continuous cultivation regime and the current-year cultivation regime showed an increase of organic matter content, together with a decrease in the field capacity and an increase in the pH. The decreased field capacity may be related to disturbances from farming practices that disrupt the aggregates in the soil structure (Dong, 2017), whereas the increased pH may have resulted from the application of quicklime (Moir & Moot, 2010) on the soil surface before S. rugosoannulata cultivation.

Effects of cultivating S. rugosoannulata under nursery stock shade on the soil bacterial community composition

Using high-throughput sequencing analyses, we observed a consistently higher abundance of Acidobacteria and a consistently lower abundance of Actinobacteria and Firmicutes in the forestlands with cultivation (Table S1). However, reports (Ramirez, Craine & Noah, 2012; Zhou et al., 2017) showed that the abundance of Acidobacteria was reduced with increased nutrient inputs because of the oligotrophic properties of these organisms, and the abundance of Actinobacteria and Firmicutes was increased with increased nutrient inputs because of the copiotrophic properties of these organisms. These inconsistencies may be ascribed to the increased organic matter in forestlands with cultivation that creates an oligotrophic soil environment due to its ability to slowly release nutrients (Tiessen, Cuevas & Chacon, 1994). In addition, we observed a consistently increased abundance of Planctomycetes and Bacteroidetes and a consistently reduced abundance of Chloroflexi and Nitrospirae in our study (Table S2). The further study is needed to reveal the underlying mechanisms for this phenomenon.

Interactions between soil fungi and bacteria are common in nature. For example, fungus-released compounds may impact bacterial selection (Warmink, Nazir & Elsas, 2009; Nazir et al., 2010). During the cultivation cycle of S. rugosoannulata, high-density hyphae are observed in the culture substrate for long periods and are even found in the spent mushroom compost. Additionally, the soil contains a considerable amount of hyphae. Therefore, the changes in bacterial communities in the soil after the incorporation of spent mushroom compost would be consistent with changes in bacterial communities in soil environments that surround the dense fungal hyphae, such as soil microhabitats, i.e., hyphospheres or mycospheres (Johansson, Paul & Finlay, 2004; Nazir et al., 2010), that more or less are densely permeated by the fungal hyphae. In our study, a decrease in bacterial diversity was found in forestlands with cultivation (Table S1), which was consistent with reports showing that the bacterial community diversity is lower than that of bulk soil (Warmink, Nazir & Elsas, 2009; Halsey, De & Andreote, 2016). Additionally, the selection of bacteria by the hyphae of S. rugosoannulata may represent a factor that contributes to the emergence of specific bacterial groups, such as Acidobacteria and Subgroup_6 which also belongs to Acidobacteria in forestlands (Figs. 1B–1E). However, it is possible that spent mushroom compost could be more influential on the soil nutrients and bacterial communities because the hyphae of the wine-cap Stropharia disappear along with the deposition of spent mushroom compost (data not published). Further study is needed to understand the impacts of spent mushroom compost and fungal hyphae on soil texture and microbial communities.