Site description and experimental design
The chosen experimental site was located in the Yining region of Xinjiang Province, China (82°7′E, 43°44′N) [1]. We selected three 10m×5m rectangular plots at the top (altitude:1121m), middle (altitude:1087m), and bottom (altitude:1053m) of a slope in a region containing Ferula sinkiangensis plants in May 2019. We then randomly selected 9 Ferula plants at each of these slope positions and collected rhizosphere soil samples at depths of 0-20 cm, 20-40 cm, and 40-60 cm [55] as per the methods previously described by Riley and Barber [56, 57]. We also collected samples of the surrounding soil within 5 cm of the root system. Next, we pooled together three of the nine plants, and we similarly pooled the rhizosphere soil samples from these three plants at each depth level, as well as the surrounding soil samples. Rhizosphere soil was stored in liquid nitrogen and was sent to Shanghai Maso Bio-Pharmaceutical Techno-logy Co., Ltd. for high-throughput sequencing. Samples of surrounding soil were air-dried and were then used to assess soil physicochemical properties [58, 59]. Soil samples from the top, middle, and bottom slope positions were respectively denoted with the SW, ZW, and XW designations, while samples collected at depths of 0-20cm, 20-40cm, and 40-60cm were respectively designated with the numbers 1, 2, and 3.
Library preparation
Microbial DNA [60] was extracted from 0.5g aliquots of each sample using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, GA, USA) according to the manufacturer’s protocols. Final DNA concentrations and purity were determined by NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, DE, USA), and DNA quality was assessed via 1% agarose gel electrophoresis. The genomic DNA pellet was stored at −30 °C prior to use.
The partial small subunit (SSU) region of the 18S rRNA gene was amplified via nested PCR [61, 62]. AML1F (forward primer) (5’-ATCAACTTTCGATGGTAGGATAGA-3’) and AML2R (reverse primer) (5’-GAACCCAAACACTTTGGTTTCCTTGGTTTCC-3 ') [38, 39] were used as the primers in the first round of amplification using a thermocycler PCR system (GeneAmp 9700, ABI, USA), whereas AMV4.5NF (forward primer) (5’-AAGCTCGTAG-TTGAATTTCG-3') and AMDGR (reverse primer) (5’-CCCAACTATCCCTATTAATCAT-3') [61, 63, 64] were used as the primers in the second amplification. The first PCR reactions were conducted using the following thermocycler settings: 3 min at 95 °C, followed by 32 cycles of 95 °C for 30s, 55 °C for 30 s, and 72 °C for 45 s, followed by a final extension at 72 °C for 10 min. PCR reactions were performed in triplicate, with each reaction being conducted in a 20 μL volume containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase, and 10 ng of template DNA. After amplification, replicates from each sample were pooled and separated via 2% agarose gel electrophoresis. These first-round PCR products were diluted 10-fold and used as templates for the second-round PCR amplification step using the following thermocycler settings: 3 min at 95 °C, followed by 30 cycles of 95 °C for 30s, 55 °C for 30 s, and 72 °C for 45 s, followed by a final extension at 72 °C for 10 min. Reaction mixtures were prepared identically to those for the first round of PCR. The resultant PCR products were purified via 2% agarose gel electrophoresis with the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, CA, USA), and DNA levels were quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer’s instructions.
Illumina MiSeq sequencing
Purified amplicons were pooled in equimolar amounts, and paired-end sequencing (2 × 300) was conducted using an Illumina MiSeq platform (Illumina, CA, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). Sequence read processing was performed using QIIME v.1.9.1.
Processing of sequencing data
Raw fastq files were demultiplexed, quality-filtered using Trimmomatic, and merged via FLASH with the following criteria: (i) The reads were truncated at any site receiving an average quality score <20 over a 50 bp sliding window. (ii) Primers were exactly matched allowing 2 nucleotide mismatching, and reads containing ambiguous bases were removed. (iii) Sequences with >10 bp overlap were merged based upon the overlapping sequence.
Operational taxonomic units (OTUs) were clustered with a 97% similarity cutoff using UPARSE (v.7.1 http://drive5.com/uparse/), and chimeric sequences were identified and removed using UCHIME. The taxonomy of each 18S rRNA gene sequence was analyzed using RDP Classifier (v.2.2 http://sourceforge.net/projects/rdp-classifier/) and the MaarjAM (https://www.maarjam.botany.ut.ee/) 18S rRNA database at a confidence threshold of 70%.
Soil physicochemical properties
Soil properties were assessed as described in a prior study [59, 65]. Analyzed soil physicochemical properties included gravimetric soil water content (by oven drying at 105 °C) [66,67], organic matter content (using the KCr2O7 method) [67], total nitrogen (using the HClO4-H2SO4 digestion method) [68, 69], total phosphorus (using a Mo-Sb colorimetric method) [69, 70], total potassium (as measured via atomic absorption spectrometry) [70, 71], nitrate nitrogen, ammonium nitrogen (as assessed via a 0.01M calcium chloride extraction method using a BRAN+LUEBBE flow analyzer) [68, 69], available phosphorus (using NaH-CO3 extracts analyzed via the Mo-Sb colorimetric method) [69, 71], available potassium (using NH4OAc extracts analyzed via atomic absorption spectrometry), pH (as measured with a Mettler Tolido FiveEasy Plus PH meter) [71], and total dissolvable salts (atomic absorption spectrometry and titration methods) [59, 71].
Data analysis
In the diversity indices, all the rarefaction and diversity analysis for community analysis was performed at the lowest number of reads (14991) per sample. The Mothur (version v.1.30.1 http://www.mothur.org/wiki/Schloss_SOP#Alpha_diversity) was used to calculate Sobs (the observed richness) as well as the Chao1 (the Chao1 estimator), Shannon (the Shannon diversity index), Simpson (the Simpson diversity index) , coverage (Good’s coverage indices), and Phylogenetic diversity (PD), which were respectively used to assess sample richness, diversity, and coverage. Since the data were not normally distributed, Kruskal-Wallis test was used to detect whether there was a significant difference in index values between the groups. The dilution curve was plotted with R (v3.6.1) to count the Alpha diversity index of the corresponding samples of these sequences by randomly selecting 14,991 sequences from the samples. The extracted data volume was taken as the x-coordinate and the Alpha diversity index as the y-coordinate.
The Venn diagram diagram were drawn using R (v3.6.1) to count the number of common and unique species (such as OTU).
The community column were drawn using R (v3.6.1) to count the species composition of different groups (or samples) at each level of classification.
The Circos sample relationship diagram is a visual circle diagram describing the correspondence between samples (or groups) and species, which was drawn by Circos-0.67-7 (http://circos.ca/). The abundance of samples in the group was calculated using mean values. In all samples, species with less than 0.01 abundance were combined as others.
The Unweighted Pair-group Method with Arithmetic Means (UPGMA) clustering was performed as a type of hierarchical clustering method to interpret the distance matrix using average linkage. The UPGMA clustering based on Unweighted Unifrac Distances at the OTU level, which was conducted by QIIME (Version 1.9.1) and drawn by R (v3.6.1).
Principal co-ordinate analysis (PCoA) based on bray-curtis at OUT level analysis the community of AMF by R ( v3.6.1). The Analysis of similarities (ANOSIM) in " vegan " R package is used to examine differences between groups. The test for significance is 999 permutations.
Differences between groups were assessed based upon linear discriminant analysis (LDA) effect size (LEfSe) (http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=lefse_upload). First, the non-parametric factorial Kruskal-Wallis (KW) sum-rank test is applied to detect the significant difference of abundance and find out the group with significant difference. Finally, linear discriminant analysis (LDA) estimate the impact of abundance of each component (species) on the difference effect, the microbes of LDA (>2) are considered to be significant groups of microbes. The multi-group comparison strategy is that species can only be considered as differential species if there are differences in multiple groups.
Since the axis length > 4 in the DCA results, canonical correlation analysis (CCA) was used to test the relationship among environmental factors, samples and microbes. The CCA was estimated using the “vegan” package in R (v3.6.1).
Correlations between soil properties and dominant OTUs were assessed via Spearman correlation analyses, while relationships between rhizosphere composition, soil properties. The Correlation Heatmap was drawn using the “pheatmap” package in R ( v3.6.1).
Since the data were not normally distributed, Kruskal-Wallis test was used to detect whether there was a significant difference in the soil physicochemical factors between the groups. Spearman correlation analyses were also run among the soil physicochemical factors with Alpha diversity index values. The Statistical analysis was carried out with SPSS 19.0 (IBM Inc., Armonk, USA).