Drivers and assemblies of soil eukaryotic microbes among different soil habitat types in a semi-arid mountain in China

Study area

The area south of the Taihang Mountains (112°28′–112°30′E, 35°01′–35°03′N) is a semi-arid region in China. The study area has a continental monsoon climate in most parts, with an average annual sunshine rate of 54%, average air temperature of 14.3 °C, and annual precipitation of 440–860 mm. The most common soil is Ustalf (United States Department of Agriculture (USDA) classification). Local vegetation restoration began in 1958. The main land-use methods are grass soil, forest soil, shrub soil, and farmland. Most of the forest soil habitat types are artificial forests, including Platycladus orientalis (L.) Franco, Quercus variabilis Bl., and Robinia pseudoacacia L. Lespedeza bicolor Turcz., Vitex negundo L., and Ziziphus jujuba Mill. var. as the major shrub types. Artemisia princeps H. Lév., Setaria viridis (L.) Beauv., and Arthraxon hispidus (Thunb.) Makino are the dominant herbaceous plants in the natural grass soil. All farmland species were planted by local residents and include Zea mays L., Triticum aestivum L., Ipomoea batatas L., Brassica campestris L., and Lycopersicon esculentum Mill.

Sampling design

Based on the status of the local vegetation resources (Zhao, 2007), we established study plots, which represented the ecological systems classified as the forest, shrub, grass, and farmland soil habitat types in October 2017. Samples were collected from the different soil habitat types (each type had three plots), and 12 independent plots (10 × 10 m) were set up. Soil samples were collected at a soil depth of 5–10 cm and were taken along the “S” model from five sampling points, then pooled together as one large sample for each plot. After cleaning up the impurities, 200 g soil in the O horizon was collected at a depth of zero to five cm. Root samples in the vicinity of the soil samples were also collected with a soil auger (inner diameter = 4 cm) at a depth of 5–10 cm. We carefully collected the soil samples and placed them in freezing boxes for storage at −70 °C until further soil geochemical and molecular biological analyses.

Topography, soil properties, and fine root data

Topographic data were taken from each plot of the different soil habitat types. GPS was used to record the longitude, latitude, and elevation of all plots. A clinometer was used to analyze the slope. The parameter settings followed Harms et al. (2001) and Valencia et al. (2004).

The composition of soil organic matter (SOM) was measured by the potassium dichromate volumetric method (Bao, 2000). Water content was measured according to a soil agricultural chemistry analytical method (Bao, 2000). A soil auger were used to collect the O horizon soil (five samples in a diagonal line for each plot), and the O-elements (O horizon soil elements) were measured by inductively coupled plasma mass spectrometry (Agilent Technologies, Palo Alto, CA, USA). A soil thermograph (TZS-1W; Zhejiang Tuopu Instrument, Zhejiang, China) was used to record the temperature of each soil sample.

Fine root biomass was determined by rinsing the surface of the root samples with clean water. Then, the roots were dried in the shade. Vernier calipers were used to sift out root diameters >2 mm in all root samples. The fine roots of each sample were weighed fresh and biomass was determined by the formula (biomass = fine root fresh weight × 10⁴/π*(inner diameter of earth drill/2)²).

Soil DNA extraction and miseq sequencing step

Soil total DNA was extracted from 50 mg of soil in each plot using the Fast DNA Isolation Kit according to the manufacturer’s instructions (Qbiogene, Heidelberg, Germany). After checking DNA concentration and purity, the extracted soil DNA was stored at −20 °C for further use. The DNA extraction methods for the eukaryotic soil microbes were described by Wu et al. (2011) and Bates et al. (2013).

The V4 region was defined by the primer set forward (5′-CCAGCASCYGCGGTAATTCC-3′) and reverse (5′-ACTTTCGTTCTTGATYRA-3′) to amplify the sequences (18S rRNA gene fragment). This primer set has been successfully used for eukaryotic microbes in several studies (Guillou et al., 2013; Logares et al., 2014). We used a two-step polymerase chain reaction (PCR) to construct the data library. First, specific primers were used to amplify the target fragment, and the target fragment was recycled with a Gel Recovery Kit (ASJ0013 Gel extraction AxyPrep DNA; Axygen Scientific Inc., Union City, CA, USA). Then, the recovery product was used as the template for the second PCR amplification. The purpose of this method was to add the necessary sequencing joints, barcodes, and sequencing primers to the ends of the segment for sequencing on the Illumina platform. More detailed steps of the two-step PCR are shown in Table S1 and S2. The libraries were sequenced using the PE300 and a MiSeq v3 Reagent Kit (Tiny Gene Co., Shanghai, China).

Bioinformatics method

Mothur v. 1.33.3 and UPARSE (version v8.1.1756) software was used to analyze the sequence libraries following the methods of Schloss et al. (2009). The original FASTQ files were denoised and split by the barcode. The portable executable reads for each sample were run in Trimmomatic (version 0.35), using the parameters (maxambig = 0, length = 200–580, which protected the longer correct fragments through a higher threshold, maxhomop = 8) and avoided low quality base pairs. Moreover, the sequence data were aggregated with 97% homology to operational taxonomic units (OTUs) (Edgar, 2013), and the NCBI database was used to BLAST each OTU representative sequence. The OTUs were also determined based on NCBI (the highest score in BLAST was classified as the matched species, and incorrect OTU information was corrected or removed).

Statistical procedures

Statistical analyses of alpha diversity (including OTU richness, Coverage, Chao 1, and Shannon’s indices) were performed in Mothur v. 1.33.3 software following Schloss et al. (2009). The topography, soil properties, and fine root data were analyzed by one-way analysis of variance using SPSS v. 19.0 software (SPSS Inc., Chicago, IL, USA). The eukaryotic community tree with the four soil habitat types was prepared using MEGA software (version 5.0) (Fig. S1). We used statistical software to analyze the dissimilarity of the dominant soil eukaryotic microbial groups among each soil habitat type, and a Venn diagram of shared OTUs was also prepared using the “VennDiagram” library; the indicator species analysis was prepared using the “indicspecies” library (Ri386 3.4.4). The network was generated with Cytoscape v. 2.8 software (Smoot et al., 2011). A principal component analysis was performed using Ri386. 3.4.4 (Table S3). Distance-based redundancy analysis (db-RDA) is an efficient method to detect soil eukaryotic microbial group-environmental correlations regarding the response of the soil eukaryotic community to environmental variables. Thus, the Monte Carlo permutation and db-RDA tests were charted using Ri386. 3.4.4 (“vegan” library), and the graph was optimized using Canoco software (Windows 4.5 package) (Braak & Smilauer, 2002). All raw sequences have been deposited into the NCBI database under accession number SRP148713.

Soil eukaryotic microbial data analyses

We obtained a total of 362,717 sequences and 1,056 OTUs (at 3% evolutionary distance) from the entire dataset, and the number of OTUs ranged from 396 to 603 per soil sample. The average coverage estimator value was 98.3% (range, 97.5–98.7%), indicating that the sequence data sufficiently covered the diversity of the eukaryotic population. Based on the Shannon’s index analysis, the eukaryotic diversity in the study region ranged from 1.17 to 3.5 (Table 1). Tukey’s HSD test was used to evaluate the significance of the results. Eukaryotic richness in the forest soil habitat type was significantly higher than the richness of the farmland, shrub, and grass soil habitat types. The diversity of eukaryotes in the forest and shrub soil habitat types was significantly higher than that in the other habitat types; farmland had the lowest eukaryotic diversity (Table 2). The indicator species analysis revealed soil eukaryotic microbe characteristic (P < 0.05) were different among the soil habitat types; a total of 98 indicators were detected in forest soil, rhizarian taxa (35 indicators) had a dominant position; 49 indicators were detected in farmland soil, metazoa taxa (19 indicators) comprised the main group; 74 indicators were found in shrub soil (25 indicators in metazoa taxa), and 48 indicators were observed in grass soil (35 indicators in rhizarian taxa) (Table S4). More detailed information about the OTUs in all soil samples is shown in Table S5.

Soil eukaryotic microbial composition south of the Taihang Mountains

The species composition of soil eukaryotic microbes remarkably differed among the different soil habitat types. Fungal (mainly Ascomycota, Basidiomycota, Chytridiomycota, Cryptomycota, and Glomeromycota), Cercozoa, Chlorophyta, Amoebozoa, Stramenopiles, Alveolata (mainly Ciliophora, Apicomplexa, and Dinophyceae), and Metazoa (mainly Nematoda) were the dominant soil eukaryotic microbial groups. OTUs that could not be identified were assigned as unclassified (Fig. 1). Within the soil fungi groups, Cryptomycota occupied the largest component, particularly in farmland (11.3–13.1%). Nematoda and Arthropoda dominated the Metazoa in the study area. Nematoda was extremely abundant among the soil metazoan groups (12.7–31.5%). Eucoccidiorida was comprised mainly of order Apicomplexa sequences. Heterocapsaceae was the main group of the Dinophyceae sequences.

Similarities and differences in the soil eukaryotic microbes among the different soil habitat types

We applied Venn and network diagrams to reveal the differences in the OTU distributions among the different soil habitat types. The Venn diagrams showed that the number of soil habitat type OTUs ranged from 552 (farmland soil) to 814 (grass soil), and only 270 of 1,056 OTUs were shared by the four soil habitat types. Moreover, there were 32 shared OTUs in the soil metazoan groups (200 OTUs), 23 shared OTUs in soil fungi (126 OTUs), 30 common OTUs in Alveolata (105 OTUs), and 90 OTUs existed in all of the rhizarian communities (268 OTUs) (Fig. 2). Network analyses revealed that the dominant OTUs were different among the four soil habitat types. Forest soil was dominated by OTU_7 (Metazoa), OTU_9 (unclassified), OTU_10 (Alveolata), and OTU_12 (Metazoa). Shrub soil was dominated by OTU_3 (Alveolata), OTU_17 (Metazoa), and OTU_22 (Fungi). Farmland soil was dominated by OTU_31 (Metazoa), whereas grass soil was dominated by OTU_20 (unclassified) (Fig. 3). Furthermore, network analyses provided more insight into these shared OTUs (OTUs existed in different land types). The results show that soil habitat types with high species diversity (forest and shrub habitat types) shared larger OTUs (a higher number of sequences) than in the low species diversity habitat types (grass and farmland habitat types). The number of sequences of shared OTUs in the grass type was usually smaller than the others. The OTUs with few sequences (<200) was the main portion of shared OTUs among the four soil types, suggesting that the eukaryotic microbes with few sequences can more easily adapt to different soil habitats (Fig. 3).

Effect of environmental factors on the soil eukaryotic community

The db-RDA of all soil habitat types and eukaryotic communities mainly showed the variation in the composition of the soil eukaryotic groups, which was explained by soil properties (water content, organic matter, O-elements, and soil temperature,) and partly by plant factors (fine root biomass) and geographical position (elevation and slope) (Fig. 4). These environmental factors explained 81.06% of the variation in the soil eukaryotic community, while 18.94% of the variation was not explained in the ordination diagram. TOM 1 (r² = 0.7311, P < 0.01), elevation (r² = 0.6032, P < 0.01), water content (r² = 0.4706, P < 0.05), and SOM (r² = 0.5092, P < 0.05) were the most important factors (Table 3), whereas the effects of soil temperature, slope, and fine root biomass were less influential.