High-throughput sequencing reveals rhizosphere fungal community composition and diversity at different growth stages of Populus euphratica in the lower reaches of the Tarim River

Study sites and sampling

The study area was located in the natural P. euphratica forest in the lower reaches of Tarim River basin, Xinjiang province (with the geographical coordinates of 40°28′∼40°55′N, 87°51′∼87°75′E), China. This area belonged to a typical continental extreme arid climate, with the annual average precipitation was less than 50 mm, the evaporation was about 2,960 mm, the annual total solar radiation was 5,692∼6,360 kJ m⁻², and the annual average temperature was 10.5 ∼11.4 °C (Yang & He, 2000), the ecological environment was extremely bad.

In mid-September 2020, soil samples were collected from selected natural P. euphratica forests. According to the classification standard in P. euphratica forest written by Wang, Chen & Li (1995). Four growth stages are selected, which including sapling (A), mature (B), overripe (C), and deadwood (D). The diameter at breast height (DBH) of sapling was about four cm, mature wood was 4∼10 cm, and overripe wood was 30∼70 cm. Three trees with similar growth and no diseases and insect pests were selected for each stage to measure morphological characteristics (Table 1) and collect rhizosphere soil samples (60 cm). At the same time, the bare land without vegetation cover was selected as blank control (CK), in the area and soil samples were collected (same depth), set up three sampling points to take mixed soil samples, and obtain three groups of parallel samples at each place.

The sampling method of rhizosphere soil microorganisms was to dig the soil profile 0.5 m away from the primary root, and start from the fine root (sample’s depth was 60 cm), the soil adhered to the root segment after shaking the fine root was the rhizosphere soil. The collected soil was divided into three parts, one part was placed in a five mL sterile centrifuge tube and stored in a liquid nitrogen tank for the determination of rhizosphere fungal community; other part soil samples were put into sealed bags and brought back to the laboratory, after natural air-drying, they were screened (two mm) for the determination of soil physical and chemical properties; the last part of the soil samples were weighed in aluminum boxes, and drying in an oven (72°/48 h) at the laboratory used for soil moisture measurement.

Soil physicochemical properties

Soil properties were assessed as described in prior studies (Bao, 2008), organic matter content (OM) was assessed using the KCr₂O₇ method, total nitrogen(TN) was assessed using the HClO₄-H₂SO₄ digestion method, total phosphorus (TP) was assessed using a Mo-Sb colorimetric method, total potassium (TK) was measured via atomic absorption spec-trometry, nitrate nitrogen(SNN), ammonium nitrogen(SAN) was assessed via a 0.01 M calcium chloride extraction method using a BRAN+LUEBBE flow analyzer, available phosphorus(AP) wasdetermined by molybdenum antimony anti Colorimetry (sodium bicarbonate extraction), available potassium (AK) was determined by atomic absorption spectrometry (ammonium acetate extraction), pH (as measured with a Mettler Tolido FiveEasy Plus pH meter), total dissolvable salts (TDS) (as assessed via atomic absorption spec-trometry and titration), and conductivity (EC) (measured by Hanna H1 2315 conductivity meter). In addition, the above determination of soil physical and chemical properties was repeated three times.

Soil fungi DNA extraction, PCR amplification and sequencing

Most of the methods adopted here are previously described in Hu, Yesilonis & Szlavecz (2021). Briefly, soil genomic DNA was extracted from rhizosphere soil samples of P. euphratica by cetyltrimethylammonium bromide (CTAB) (Hu, Yang & You, 2010), after that, the purity and concentration of DNA were detected by 2% agarose gel electrophoresis. A proper amount of DNA sample was taken into a centrifuge tube and diluted with sterile water to 1 ng µL⁻¹. Using diluted genomic DNA as a template, ITS1 primers ITS5-1737F (5′-GGAAGTAAAAGTCGTAACAAGG-3′) (Bellemain et al., 2010) and ITS2-2043R (5′-GCTGCGTTCTTCATCGATGC-3′), Phusion^® high-fidelity PCR Master Mix with GC Buffer and efficient high-fidelity enzyme from Biolabs, New England, were selected for PCR (Walters et al., 2016). The PCR product was detected by electrophoresis using 2% agarose gel, and recovered using the gel recovery kit provided by Qiagen company, TruSeq^® DNA PCR-free Sample Preparation Kit (Illumina, USA) was used to construct the library, which was quantitated by Qubit and Q-PCR. Lastly, the constructed library was sequenced and computerized on Illumina HiSeq2500 platform of Beijing Compson Biotechnology Co., Ltd.

Sequence processing and analysis

FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/) (Magoč & Salzberg, 2011) was used to splice the offline data obtained by sequencing to get Raw tags data, through the Qiime (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html) (Gregory et al., 2010) for data quality control, after strict filtering (Bull & Asenjo, 2013), clean tags are obtained, then we do chimera filtering (https://github.com/torognes/vsearch/) (Haas et al., 2011), and ultimately get can be used for further analysis of Effective tags.

Uparse (v7.0.1001, http://www.drive5.com/uparse/) (Edgar, 2013) was used to cluster sequences into Operational taxa (OTUs) with 97% sequence similarity, annotate the OTUs sequences, analyze the species annotation with the BLAST method (http://qiime.org/scripts/assign_taxonomy.html) (Altschul, 1990) in Qiime software and the UNITE (v8.2) database (https://unite.ut.ee/) (Kõljalg et al., 2013), and count the community composition of each sample at each classification level. The MUSCLE (v3.8.31) (http://www.drive5.com/muscle/) (Edgar, 2004) was used for fast multi-sequence alignments to get the phylogenetic relationships of all OTUs sequences, the data of each sample was homogenized, and the data with the least amount of data in the sample was used as the standard for homogenization. The subsequent alpha diversity and beta diversity were analyzed based on the homogenized data.

Qiime (Version 1.9.1) was used to calculate fungal diversity index of soil samples, including microbial richness (Chao1 index and ACE index) and microbial diversity (Shannon index and Simpson index) (Wang et al., 2019).

R (Version 2.15.3) was used to draw the dilution curve, Venn diagram and Principal coordinates analysis (PCoA) diagram. The Analysis of similarities (Adonis)in “vegan” R package was used to examine differences between groups. LEfSe software was used for linear discriminant analysis and effect size analysis with the default filtering value of LDA score set at 4. 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 physicochemical properties and fungal community composition were assessed via Spearman’s correlation analyses.

Statistical analysis (including one-way analysis of variance (ANOVA) and Spearman’s correlation analysis) were carried out with SPSS 22.0 (IBM Inc., Armonk, USA).

Differences in physical and chemical properties of rhizosphere soil

The physical and chemical properties of soil samples were measured (Table 2). The content of SWC in A sample was significantly higher than B and C samples (P < 0.05), however, there was no significant difference in the contents of TP, SNN and SAN in P. euphratica rhizosphere soil at different growth stages (P > 0.05). The content of OM, TN, TK, AK, AP, pH, EC and TDS in D sample were higher than A, B and C samples, among which the contents of TK, AK, pH, EC and TDS were significantly different (P < 0.05). In addition, there was no significant difference between B and C samples except the content of AK (P > 0.05).

Sequencing data and OTU clustering

By sequencing the ITS fragments of soil fungi, 947,049 effective sequences were obtained from soil samples, and 960 OTUs were obtained by clustering with 97% sequence similarity. The dilution curves of all soil samples tended to be flat, indicating that the sequencing depth had basically covered all fungal groups in the samples, which could reflect the real situation of soil fungal community in the rhizosphere of P. euphratica (Fig. 1).

Venn diagram revealed that there were 391 OTUs for A sample, 342 OTUs for B sample, 463 OTUs for C sample, 518 OTUs for D sample, 128 OTUs were shared by A, B, C and D samples (Fig. 2). Besides, the numbers of unique OTUs to each sample were as follows: 75 for A sample, 68 for B sample, 77 for C sample and 191 for D sample, respectively, accounting for 19.18%, 19.88%, 16.63% and 36.87% of the all OTUs.

Differences in fungal diversity

The alpha diversity index of soil fungal community in rhizosphere of P. euphratica at different growth stages was different (Table 3). As shown in Table 3, the Coverage index of each sample was close to 100%, which proved the integrity of the detected samples sequence, indicating that the sequencing results at this level could reflect the true situation of fungal community composition in the measured samples. Shannon index was consistent with Simpson index, the value was the highest in D sample, Chao1 index and ACE index were the highest in B sample. However, there was no significant difference between different growth stages.

Both PCoA (Fig. 3) and Adonis (Table 4) all revealed that there were significant differences in soil fungal community composition between B and C samples (R² = 0.31 P = 0.001), B and D samples (R² = 0.28 P = 0.001). There was no significant difference between A and the other samples.

Differences in fungal community composition at different levels of rhizosphere soil

A total of 10 phyla, 36 classes, 77 orders, 165 families, 275 genera and 353 species were identified through comparative identification of OTUs representative sequences of soil samples (Fig. 4). As shown in Fig. 4, Ascomycota was the dominant phylum in the rhizosphere soil of P. euphratica (average relative abundance of 58.54%), followed by Basidiomycota (15.96%). Compared with phylum classification level, the composition of rhizosphere soil fungal community in different stages differed greatly from class classification level. At the class classification level, Sordariomycetes, Pezizomycetes and Agaricomycetes were the dominant fungi, and the relative abundance of Sordariomycetes was 8.05%∼34.85%, Pezizomycetes (0.21%∼36.09%) and Agaricomycetes (0.35%∼56.94%). At the order classification level, the dominant species were Pezizales, Agaricales, and Hypocreales, accounting for 0.20%∼36.08%, 0.23%∼56.84% and 1.04%∼29.02% of the total sequences in all groups, respectively. The dominant species at the family classification level are Pyronemataceae, Agaricales_fam_Incertae_sedis and unidentified_Sordariales_sp, Pyronemataceae (0.19%∼35.32%), Agaricales_fam_Incertae_sedis (0.001%∼56.43%), unidentified_Sordariales_sp (0.02%∼27.84%). The dominant species at the genus classification level were Geopora, Chondrostereum and unidentified_Sordariales_sp, Geopora (0.08%∼35.26%), Chondrostereum (0.001%∼56.43%), unidentified_Sordariales_sp (0.02%∼27.84%). The dominant species at the species classification level were Chondrostereum_purpureum, Geopora_sepulta and Geopora_sp, Chondrostereum_purpureum (0.001%∼56.43%), Geopora_sepulta (0.01%∼27.38%), Geopora_sp (0.07%∼26.52%). In conclusion, the relative expression abundance of the top ten fungi at each taxonomic level was significantly different at each stage.

Species differences of soil fungal community

LEfSe was used to search for biomarkers, so as to find species with significant differences in abundance between groups. In this study, LEfSe analysis was used to analyze the species abundance data of fungi in rhizospheres soil samples, the rank sum test was used to detect the different species in different groups and LDA score (LDA score = 4) was obtained through LDA. Finally, the evolutionary clade of different species (Fig. 5A) and the histogram of LDA value distribution (Fig. 5B) were drawn, both of them reflected the distribution characteristics of species with different rhizosphere fungal community structure of P. euphratica. A total of 23 biomarkers were obtained, with relatively more in B and C samples (three taxa for A sample, eight taxa for B sample, eight taxa for C sample, and four taxa for D sample). Specifically, Sporobolomyces sp, Sporobolomyces, Fusarium proliferatum were significant in A sample, Thelephoraceae sp, Thelephorales, unidentified, Thelephoraceae, Sordariales, Sordariales sp, unidentified Sordariales sp, unidentified Sordariales sp were significant in B sample, Synchytriales, Synchytriomycetes, Synchytrium endobioticum, Synchytrium, synchytriaceae, Chondrostereum, Chondrostereum purpureum, Agaricales fam Incertae sedis were significant in C sample, Tremellomycetes, Acremonium rutilum, Gibberella intricans, Gibberella were significant in D sample.

Correlation of soil physical and chemical factors and fungal community structure

Canonical correlation analysis (CCA) can reflect the relationship between microflora and environmental factors, and can obtain the important environmental driving factors that affect the distribution of samples (Fig. 6). EC, TDS and AK were the main environmental factors that significantly affected the rhizosphere fungal community of P. euphratica (P < 0.05), and EC was the main driving factor (R ² = 0.704, P <0.01) (Table 5). As shown in Fig. 7, the interpretation amount of the first sorting axis was 12.15%, and that of the second sorting axis was 11.47%. Spearman’ s correlation analysis was used to analyze the correlation between soil factors and the relative abundance of the top 35 species at genus level (Fig. 7). The results showed that EC was significantly positively correlated with Thielavia, Xerombrophila (P <0.05), and negatively correlated with Lecanicillium, unidentified, unidentified Onygenales sp, Fusarium, Geopora (P <0.05). TDS content was positively correlated with Xerombrophila (P < 0.05), and negatively correlated with Lecanicillium, unidentified, unidentified Onygenales sp, Fusarium, Geopora (P < 0.05). AK was positively correlated with Thielavia, Xerombrophila, Didymella, Neomyrmecridium, Aspergillus (P < 0.05), and negatively correlated with unidentified (P < 0.05).