Diversity of rhizosphere and endophytic fungi in Atractylodes macrocephala during continuous cropping

Study site and experimental design

Samples were collected from a grassland located in a hilly area in Lin’an County, Zhejiang Province, China, with a typical mid-subtropical monsoon climate (30°14′N, 119°25′E, 155 m elevation). The study area was divided into five zones according to different cropping practices (no A. macrocephala planted, A. macrocephala planted for 1 year, A. macrocephala planted for 2 years, 1-year fallow, and 2-year fallow). Each zone consisted of a 20 × 20 m plot, which was managed as uniformly as possible. Due to limited land access, only one plot was used for each zone. However, we randomly collected 15 soil cores from the plow layer (0–30 cm in depth), thoroughly pooled them as a composite sample and divided the pooled samples into three replicates to make the samples more representative (Hu et al., 2017; Tan et al., 2017a). Fallow soil zones were left unplanted for the specified number of years after an initial planting and harvesting of A. macrocephala. For the fallow zones, three soil samples were randomly selected, with each soil sample composed of 15 soil cores randomly collected from the plow layer (0–30 cm in depth).

Plants were sampled by inspecting them for root rot symptoms, and were then placed either into the group with or without symptoms. Symptoms were recorded regarding the general appearance of the entire plant, and more specifically regarding the health of the roots, the stems, and the foliage. In the absence of overt symptoms in above-ground parts, roots were checked for necrotic flecking, dead, broken roots, and water-soaked lesions (Sims et al., 2019). For the 1-year and 2-year A. macrocephala planted zones, 15 plant samples for each plant type (1-year healthy plant, 2-year healthy plant, and 2-year diseased plant) were collected and put into sterile plastic bags, placed on ice and transported to the laboratory immediately. The soil which was loosely attached to the roots was removed by gentle shaking. Tightly bound soil was considered rhizospheric and was obtained by firmly shaking the root system in a sterile plastic bag (Kobayashi et al., 2015; Tan et al., 2017b). The rhizospheric soil samples were thoroughly mixed, and randomly divided into three subsamples, as described previously (Yu et al., 2019). The plants were washed with tap water and then rinsed three times with distilled water. They were then separated into leaf, stem, tuber, and root tissues. For surface disinfection, samples from each plant tissue were successively immersed in 75% ethanol for 1.5 min, a fresh 2.5% sodium hypochlorite solution for 3 min, and 75% ethanol for 30 s, then washed with distilled water three times (Zhu et al., 2018a). Each plant tissue was thoroughly mixed to make a composite sample, and randomly divided into three subsamples. In total, there were 18 soil samples (3 fallow soils and 3 rhizospheric soils, with 3 replicates each) and 36 plant samples (3 plant types, with 4 tissue types and 3 replicates each). Collected soils were sieved through a <4 mm filter and stored at 4 °C prior to biochemical characterization, or frozen at −80 °C less than two weeks before DNA extraction.

The soil and plant conditions were as follows: 1-year fallow soil (FS_1), 2-year fallow soil (FS_2), 1-year healthy rhizospheric soil (RS_H1), 2-year healthy rhizospheric soil (RS_H2), 2-year root-rot diseased rhizospheric soil (RS_D), blank control soil with no A. macrocephala planted (S_CK), 1-year-old healthy root/stem/leaf/tuber (R/S/L/T_H1), 2-year-old healthy root/stem/leaf/tuber (R/S/L/T_H2), 2-year-old root-rot diseased root/stem/leaf/tuber (R/S/L/T_D) (Fig. S1).

Soil physicochemical properties

The chemical and physical properties of all soil samples were analyzed. Soil pH was measured by a pH meter and soil electrical conductivity by an electric conductivity bridge, according to the manufacturers’ protocols (Satapute et al., 2019). The moisture content of the soil was determined after dry weights were obtained by drying the samples at 105 °C. The soil organic matter content was measured using Yang’s method (Yang et al., 2019). We used a nickel crucible to digest soil total K with sodium hydroxide at 750 °C, and extracted soil available K using 1 mol/L ammonium acetate (Zhao et al., 2014). All K contents were determined by an atomic absorption spectrophotometer (AAnalyst 400, PerkinElmer, US) (Yao et al., 2017). Soil available P and total P and were determined by the ammonium molybdate method (Murphy & Riley, 1962). Soil total N was determined using the semi-micro Kjeldahl method, according to Nelson and Sommers (Nelson & Sommers, 1980). Soil hydrolysis N was determined by the alkali-hydrolyzed distillation method described by Zhang et al. (2018). Heavy metals Pb, Cu, Hg, As, Cr and Cd were detected using flame atomic absorption spectrometry (Analytikjena AAS vario6, Germany) according to the method of Ye et al. (2011). All assays were performed in triplicate.

DNA extraction, PCR amplification, and sequencing

Genomic DNA was isolated from the rhizosphere, fallow soil, blank control soil samples and plant organs using the E.Z.N.A.^® soil DNA Kit (Omega Bio-tek, USA) according to the manufacturer’s instructions. The DNA extractions for each sample were conducted in triplicate. The quantity and purity of the DNA samples were assayed by a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, USA). The DNA samples were also run on a 1% agarose gel by electrophoresis to further assess their quality and integrity. The DNA samples were diluted to 1 ng/µL with sterile water and stored at −80 °C until further processing. The ITS1 region was amplified from fungal genomes with primers ITS1-F (CTTGGTCATTTAGAGGAAGTAA) and ITS2-F (GCTGCGTTCTTCATCGATGC) in a thermocycler PCR system (Eppendorf Mastercycler pro, Germany). The A adaptor-containing end was the sequencing end, and a 10 bp barcode was added between the adaptor and the forward primer sequence to allow for dis-crimination of the samples during sequencing (Bai et al., 2015). The PCR mixtures contain 5 × TransStart FastPfu buffer 4 µL, 2.5 mM dNTPs 2 µL , forward primer (5 µM) 0.8 µL, reverse primer (5 µM) 0.8 µL, TransStart FastPfu DNA Polymerase 0.4 µL, template DNA 10 ng, and finally ddH₂O up to 20 µL. The reaction program was: initial denaturation at 95 °C for 3 min, followed by 27 cycles of denaturing at 95 °C for 30 s, annealing at 55 °C for 30 s and extension at 72 ° C for 45 s, and single extension at 72 °C for 10 min, and end at 4 °C. The PCR reactions for each DNA sample were conducted in triplicate, resulting in nine total PCR products for each plant/soil sample (3 DNA replicates, with 3 PCR replicates each). The resulting amplicons were separated by electrophoresis using a 2% agarose gel and extracted with the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, USA). The nine purified PCR products were then pooled before paired-end sequencing (2 ×300) on an Illumina MiSeq platform (Illumina, USA) (Saenz et al., 2019).

Bioinformaitics analysis

The raw ITS rDNA gene sequencing reads were demultiplexed, quality filtered by Trimmomatic and merged by FLASH with the following criteria: (1) Reads were removed if any site had an average quality score <20 over a 50 bp sliding window. Reads containing Ns or with a length of less than 50 were also removed. (2) The remaining pair-end reads were assembled into larger contigs according to their overlaps, with a minimum overlap length of 10 bp. (3) The maximum mismatch ratio allowed in the overlap area of the merged sequence was 0.2. (4) The directionality of reads was corrected based on their barcodes and primer sequences, with no mismatches allowed in the barcode and 2 mismatches allowed in the primers (FLASH and Trimmomatic).

We assigned sequences with ≥ 97% similarity to the same Operational Taxonomic Units (OTUs) in UPARSE (version 7.1, Table 1) and then selected a representative OTU for each community in order to perform taxonomic annotation using a Ribosomal Database Project classifier. The OTU clustering process and criteria were as follows: (1) In order to reduce the amount of redundant calculation in the process of analysis, non-repetitive sequences were obtained from optimized data (http://drive5.com/usearch/manual/dereplication.html). (2) Singletons were removed (http://drive5.com/usearch/manual/singletons.html). (3) With 97% clustering, an OTU sequence was required to be at least 3% different from all other OTUs, and OTU representative sequence were the most abundant sequences in their neighborhood. During the clustering process, chimeric sequences were discarded. (4) All optimized sequences were mapped to OTU representative sequences, and sequences with a similarity >97% were selected to generate an OTU table (Usearch, vsesion 7.0, http://drive5.com/uparse/). Alpha within samples was analyzed by in-house Perl scripts.

Next, we rarified the OTU table and calculated Observed Species, Shannon’s index, and Chao1to calculate Alpha Diversity (Tan et al., 2017a). QIIME and the compute_core_microbiome.py script were used to identify OTUs that occurred within 95% of the A. macrocephala or soil samples (Mahoney, Yin & Hulbert, 2017). A species was defined as dominant if P_i > 1/S, where S represents species richness, and P_i is the relative abundance of a species i, defined as the number of competing species present in the community (Rivera-Orduna et al., 2011; Wu et al., 2013).

Statistical analysis

One-way ANOVA of soil physical and chemical properties as well as the changes in dominant fungal genera of rhizospheric and fallow soil samples were performed using SPSS (version 19.0; SPSS, Chicago, IL, USA). Significance was calculated by Duncan’s mean test (P < 0.01). Mean values of replicates are expressed as mean ± standard error (SE). False discovery rate (FDR) was used to adjust the significance value in multiple comparisons.

A heatmap was drawn to show the relative abundances of the 25 predominant genera in each sample. Principle component analysis (PCA) was conducted using all OTUs relative abundances in the plant and soil samples.

We performed a canonical correspondence analysis (CCA) using normalized OTU abundance and soil physicochemical data using the vegan package implemented in R (https://www.r-project.org/) (Dixon, 2003). Significant correlations between soil properties and the fungal communities at OTU level were determined via Spearman’s correlations using >0.8 or <−0.8 as the threshold, with a cutoff of P <0.01 (Faust et al., 2012). All the data analysis for correlations was finished in online ‘i-sanger’ (http://www.i-sanger.com/) developed by Majorbio Bio-Pharm Technology Co. Ltd (Lin et al., 2018).

Soil properties

Continuous cropping practice noticeably decreased soil pH compared with the S_CK and fallow soil. Electrical conductivity, total P, total K, hydrolysis N, available P, and available K contents were much higher in the continuous cropping soil relative to the S_CK, while organic matter and total N contents were higher in the S_CK than in continuous cropping soil samples. In addition, the continuous cropping soil exhibited much higher Pb, As, Cd, Cr, and Cu contents, while S_CK had higher Hg levels (Table 2).

Composition and alpha diversity of the fungal communities

A total of 484,844 and 251,911 reads were obtained from A. macrocephala plant and soil samples, respectively. A total of 48,550 reads mapped to unidentified phylum, while 271,914 reads mapped to unidentified genera. As indicated by the rarefaction curves, all samples reached the saturation phase with a satisfactory level of confidence and a Good’s coverage index ≥ 99.59% (Table 1, Fig. 1). Sequences from plant samples clustered into 2,545 OTUs using a 3% dissimilatory threshold. The core fungi of plant samples consisted of 276 OTUs, each of which had a relative abundance of greater than 1% and were found in 95% of plant samples, and included Pleosporales (with 73 OTUs), Hypocreales (45), Capnodiales (81), Helotiales (33), Tremellales (19), Cantharellales (5) and Sporidiobolales (20) (Fig. 2A). Samples T_D and T_H2 contained the highest levels of Hypocreales (93.43%) and Pleosporales (72.57%) sequences, respectively (Fig. 2C).

Sequences from soil samples clustered into 3,596 OTUs, and their core fungal composition (431 OTUs) included Eurotiales (63), Hypocreales (104), Pleosporales (65), Chaetothyriales (51), Mortierellales (27), Sordariales (44), Helotiales (42), Geminibasidiales (3), Chaetosphaeriales (12) and Capnodiales (20) (Fig. 2B). Samples S_CK and FS_1 contained the highest levels of Eurotiales (40.08%) and Hypocreales (37.04%), respectively (Fig. 2C). At the phylum level, Zygomycota, Ascomycota, and Basidiomycota were the dominant fungal phyla in both continuous cropping and fallow samples. During continuous cropping, rhizospheric soil showed a significant increase of Talaromyces, Fusarium, Cotylidia and Acremonium, and significant decrease of Penicillium (P < 0.01) (Fig. 2D). A higher abundance of the genera Cotylidia, Phialophora, and Fusarium in rhizospheric fungi was observed in root-rot diseased plants than the S_CK (Fig. 2D). Additionally, rhizospheric soil fungal alpha-diversity decreased with long-term A. macrocephala cropping. The fungal diversity index was the highest in blank control soil, and the fungal diversity of the fallow soil samples was also higher than that in samples from continuous cropping. Moreover, root-rot disease decreased the endophytic fungal diversity of A. macrocephala (Table 1).

Structure of the fungal communities

We performed a hierarchical clustering analysis on the top 25 most abundant fungal genera-across 18 plant and soil samples (Fig. 3A). The analysis showed that A. macrocephala aerial samples clustered into one category (except for T_D), while all rhizospheric and fallow soil samples clustered into another category. R_D, T_H1, R_H1, and T_H2 endophytic samples clustered together, but were separate from the R_H2 sample. PCA analysis based on the OTU composition revealed obvious variations in fungal communities among the 18 plant and soil samples (Fig. 3B). The first two axes (PC1 and PC2) explained 41.59% and 19.86% of the total variance in the fungal OTUs of the plant and soil samples, respectively, in A. macrocephala continuous cropping fields.

Correlation between soil properties and the fungal communities

We then removed the redundant variables, and the eight remaining environmental characteristics were subjected to detrended correspondence analysis (DCA). The DCA showed that the responses of fungal community composition to soil properties fit a single-peak model (Length = 6.32). Therefore, we further analyzed how fungal communities associated with soil physicochemical factors based on CCA (Fig. 4). Soil physicochemical factors Hg, electrical conductivity, available P, Cr, pH, and hydrolysis N positively correlated with each other, while total N and organic matter positively correlated with each other. Moreover, pH (p = 0.001), hydrolysis N (p = 0.001), electrical conductivity (p = 0.001), and Hg (p = 0.001) had a significant correlation with fungal community structure (Fig. 4). Spearman’s correlation analysis revealed that the relative abundance of the genus Acremonium was significantly negatively corelated with soil electrical conductivity and Hg content, while Penicillifer genus abundance had a negative correlation with soil Cr content (Tables S1, S2).