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

Background Populus euphratica is one of the most ancient and primitive tree species of Populus spp and plays an important role in maintaining the ecological balance in desert areas. To decipher the diversity, community structure, and relationship between rhizosphere fungi and environmental factors at different growth stages of P. euphratica demands an in-depth investigation. Methods In this study, P. euphratica at different growth stages (young, medium, overripe, and decline periods) was selected as the research object, based on the determination of the physicochemical properties of its rhizosphere soil, the fungal community structure and diversity of P. euphratica and their correlation with soil physicochemical properties were comprehensively analyzed through high-throughput sequencing technology (internal transcribed spacer (ITS)) and bioinformatics analysis methods. Results According to the analysis of OTU annotation results, the rhizosphere soil fungal communities identified in Populus euphratica were categorized into10 phyla, 36 classes, 77 orders, 165 families, 275 genera and 353 species. The alpha diversity analysis showed that there was no obvious change between the different growth stages, while beta diversity analysis showed that there were significantly differences in the composition of rhizosphere soil fungal communities between mature and overripe trees (R2 = 0.31, P = 0.001), mature and deadwood (R2 = 0.28, P = 0.001). Ascomycota and Basidiomycota were dominant phyla in the rhizosphere fungal community and the dominant genera were Geopora, Chondrostereum and unidentified_Sordariales_sp. The relative abundance of the top ten fungi at each classification level differed greatly in different stages. Canonical correspondence analysis (CCA) and Spearman’s correlation analysis showed that conductivity (EC) was the main soil factor affecting the composition of Populus euphratica rhizosphere soil fungal community (P < 0.01), followed by total dissolvable salts (TDS) and available potassium (AK) (P < 0.05). Conclusions Our data revealed that the rhizosphere fungal communities at the different growth stages of P. euphratica have differences, conductivity (EC) was the key factor driving rhizosphere fungi diversity and community structure, followed by total dissolvable salts (TDS) and available potassium (AK).


INTRODUCTION
Soil microorganism is one of the key components of desert soil microecosystem, which plays an important role in soil nutrient cycle and vegetation nutrient supply and recovery (Neilson et al., 2012;Bull & Asenjo, 2013). Rhizosphere soil fungi are an important part of soil-plant ecosystem, and 'their community structure closely affects the growth and development of plants. The diversity and ecological function of rhizosphere soil fungi have become a hot topic in the field of soil microecology at home and abroad (Wu, Lin & Lin, 2014). Previous studies on soil fungal communities are mainly based on the traditional plate culture method, which has proved the dependence of plants on fungal communities (Costa et al., 2006). However, the number of soil fungi obtained by this method is very small, and many fungi cannot be isolated and cultured directly (Wang et al., 2021a;Wang et al., 2021b), so that it cannot fully reflect the composition of the community structure.
With the development of science and technology, high-throughput sequencing technology has gradually become the main method to study soil microorganisms Luo et al., 2020). In the current research, it was found that there is a close correlation between rhizosphere microorganisms and plants, and their interaction mechanism was complicated (Morgan, Bending & White, 2005). The composition of rhizosphere microbial community will be affected by vegetation type (Sinha et al., 2008), soil type (Lu et al., 2011;Acharya et al., 2021), human factors and other factors (Li et al., 2015;Huang et al., 2021;Yin, Li & Du, 2021). For the same plant, different planting methods (Durrer et al., 2021), different development stages, even genetic background and other factors will lead to the change of the rhizosphere microbial community (Marschner et al., 2001;Dang et al., 2020).
Populus euphratica is one of the most ancient and primitive tree species of Populus spp. It is a unique desert forest tree species, having the characteristics of drought resistance, saline alkali resistance, heat resistance, wind and sand resistance, and plays an important role in maintaining the ecological balance in desert areas (Nekoa et al., 2018). China has the largest distribution range and the largest number of P. euphratica species in the world. More than 90% of P. euphratica forests in China are concentrated in Xinjiang region, and mainly in the lower reaches of Tarim River and many downsteam in the southern edge of Tarim basin (Wang, 1996). However, some researchers indicated that P. euphratica population regeneration in the lower reaches of Tarim River showed a decline type, the proportion of young plants in the population decreased significantly or even lacked, population was mostly over mature forest plants and the overall performance of the decline trend (Zhou et al., 2018). Currently, most studies of P. euphratica were focus on heteromorphic leaves (Li et al., 2020a;Li et al., 2020b), photosynthetic physiological characteristics (Wang et al., 2014), water use efficiency (Zhou et al., 2019), population structure (Miao et al., 2020), etc, while few studies on the relationship between plant and rhizosphere soil fungal community composition and diversity. Therefore, the goal of our work was to (1) analyze the composition and diversity of rhizosphere soil fungi of P. euphratica at different growth stages in the lower reaches of the Tarim River; (2) explore the dominant fungi in the P. euphratica rhizosphere at different stages, and the change of rhizosphere soil physical and chemical properties; (3) elucidate the correlation between fungal community composition and environmental factors. This study will provide scientific basis for the study of rhizosphere microorganisms and population rejuvenation of P. euphratica and the interaction between plants and microorganisms in arid areas.

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 −2 , 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 2 O 7 method, total nitrogen(TN) was assessed using the HClO 4 -H 2 SO 4 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 −1 . 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 R 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 R 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.
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 2 = 0.31 P = 0.001), B and D samples (R 2 = 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%,  The dominant species at the species classification level were Chon-drostereum_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

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 2 = 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

Figure 5 Cladograms (A) and LDA value distribution histogram (B) in samples.
In cladograms (A), the circle radiating from inside to outside represents the taxonomic level from the Phylum to the species. Each small circle at a different taxonomic level represents a taxonomic at that level, and the diameter of the small circle is proportionate to the relative abundance of species. The figure shows the species with LDA Score greater than the set value (default setting is 4) (B), that is, species with significant differences in different groups. The length of the histogram represents the size of the influence of species with significant differences. The English letters in the figure is the group name (A, B, C and D: bare soil, sapling, mature, overripe and deadwood, respectively). The length of the arrow line represents the degree of correlation between a certain environmental factor and community and species distribution, and the longer the arrow, the greater the correlation. When the angle between the environmental factors is acute, it means that there is a positive correlation between the two environmental factors, while when the angle is obtuse, there is a negative correlation. (SWC, soil water content; OM, organic matter; TN, total nitrogen; TP, total phosphorus; TK, total potassium; AK, available potassium; AP, available phosphorus; SNN, nitrate nitrogen; SAN, ammonium nitrogen; pH, hydrogen ion concentration; EC, electrical conductivity; TDS, total salt; respectively; A, B, C, D and CK: sapling, mature, overripe, deadwood and bare soil). Full-size DOI: 10.7717/peerj.13552/ fig-6 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).

DISCUSSION
The lower reaches of the Tarim River is located in an extremely arid climate area, with a harsh ecological environment and a very fragile ecosystem (Zhao et al., 2015). Water and salt content are the key factors limiting plant growth and development in this habitat, our study found that SWC reached the highest in the root of young P. euphratica, followed by the deadwood, and was significantly higher than mature, overripe and bare land (P < 0.05), TDS content reached the highest in the rhizosphere of deadwood, which was three times of that in other growth stages, the high SWC in the rhizosphere soil of saplings helps Figure 7 Heat maps of Spearman's correlation analysis. Ordinate is the information of environmental factors, and abscissa is the information of species at the genera level of taxonomy. The correlation coefficient r of Spearman is between −1 and 1, r < 0 is negative correlation, r > 0 is positive correlation, and the mark * is significance test p < 0.

Notes.
r 2 reflects the relationship between soil physical and chemical properties and fungal community structure, the P values are the correlation coefficients. ** P < 0.01. *P < 0.05. them to translocate, grow and develop smoothly in this environment. With the increase of tree age and the weakening of growth potential, the capacity of plants to absorb soil salt decreased, and the salt accumulation significantly increased the salt content in the rhizosphere soil of the deadwood, it was consistent with the results of Guan's study (Guan et al., 2020). In addition, the accumulation of salt and some nutrients in the rhizosphere soil may be the result of the comprehensive action of roots, water, microorganisms and other factors. P. euphratica is a salt-tolerant plant, during its growth and development, it can selectively absorb some salt ions, especially potassium ions, so as to resist the stress of saline soil environment. Therefore, EC and AK in the rhizosphere soil of sapling, mature and overripe were significantly lower than deadwood.
In this study, Shannon index and Simpson index were used to calculate community diversity, ACE and Chao1 index were used to calculate community richness (Bokulich et al., 2013). Through the calculation of the above indexes, there was no significant difference between different groups. However, in terms of fungal community structure, there were significant differences between mature and overripe, mature and deadwood. According to previous studies, with the growth and development of plants, the metabolic activities of plant roots, the nutrients, water and ventilation in the environment around the roots have different changes, and the diversity of rhizosphere microorganisms and community structure also change (Qiu et al., 2016;Zhao, Zhou & Ren, 2020). Deadwood root's SWC, OM and part of the soil nutrient content values were significantly higher than other stages, the number up to 518 OTUs, alpha diversity index is relatively high, studies have shown that SWC, OM, available nutrient content is higher, can create favorable conditions for the growth of fungi, which can protect soil fungi and enhance their community abundance (Zhang et al., 2021).
Ascomycota and Basidiomycota were the dominant phyla, but the abundance of the two fungi was different at different developmental stages of P. euphratica. Among them, Ascomycota had the highest relative abundance in sapling (81.20%), followed by mature (72.08%), and the lowest in overripe (38.74%), Basidiomycota had the highest relative abundance in overripe (57.25%), and the relative abundance in other periods was low (range 1.41%∼2.69%). Many studies have shown that Ascomycota and Basidiomycota were the dominant fungi in plant rhizosphere (Chen et al., 2021). For example, in the rhizosphere of Picea asperata (Liu et al., 2021), Castanopsis hystrix and Pinus massoniana (Wang et al., 2021a;Wang et al., 2021b), the relative abundance of both fungi were higher than other fungi. Ascomycota was the dominant fungi in soil, most of which were saprophytic fungi (Paungfoo-Lonhienne et al., 2015), they could degrade the organic matter such as lignin and keratin in soil (Beimforde et al., 2014), and have a rapid evolution rate in various soil ecosystems (Wang & Guo, 2016). In addition, Ascomycota might have the ability to adapt to saline alkali or relatively arid soil environment, which made it the main dominant fungal community in P. euphratica rhizosphere soil under the harsh environment in the lower reaches of Tarim River. Basidiomycota was mostly saprophytic or parasitic fungi, as an important decomposer in the soil (Yelle et al., 2008), it played an important role in the nutrient cycle of P. euphratica rhizosphere soil.
Compared with phylum classification level, dominant flora differed greatly among different groups at genus level. Chondrostereum, as a biomarker in C sample, had the highest relative abundance (56.43%) in C sample, but less than 0.01% in B and D samples. Geopora had the highest relative abundance in A sample (35.26%), but less than 0.1% in D sample. Unidentified_Sordariales_sp, as the biomarker in B sample, had the highest relative abundance in B sample (27.84%), but the average abundance in other periods is less than 0.1%. Otherwise, the abundance of the above dominant fungi in CK was less than 1%. In conclusion, dominant fungal communities changed significantly in different stages. with the growth and development of P. euphratica, its roots had different selective enrichment effects on specific fungi in the soil, which further indicated that these fungi may be closely related to the growth and development of P. euphratica. Moreover, previous studies have shown that the dominant fungal community changes dynamically with the growth and development of plants (Li et al., 2020a;Li et al., 2020b).
The microbial community structure in the rhizosphere of plants was affected by various biological and abiotic factors. The species and growth stage of plants determine which microorganisms can enrich in the rhizosphere. On the other hand, the physical and chemical properties of the soil in this area have a more macroscopic influence on the microbial community (Chu et al., 2020). The results showed that EC, TDS and AK had significant effects on soil fungal community, and EC was the most important factor, this is reflected in the Spearman's correlation analysis between dominant fungi and soil physical and chemical factors. EC was significantly correlated with seven dominant fungi species, of which it was significantly negatively correlated with four dominant fungi species. Some studies showed that the main soil factor affecting the metabolic characteristics of P. euphratica rhizosphere fungal community was EC, and EC and AK were negatively correlated with fungal metabolic activity (Wang et al., 2017), which were consistent with our study. Due to the scarcity of precipitation and high daily evaporation, the study area belonged to the harsh arid and salt-alkali environment. TDS not only had a certain influence on plant growth, but also could directly inhibit the activity of microorganisms (Wang et al., 2009), and the influence on soil fungal community should not be ignored.

CONCLUSIONS
Based on high-throughput sequencing technology, we studied the differences in the composition and structure of soil fungal community in the rhizosphere of P. euphratica at four development stages. The results showed that the dominant phyla of rhizosphere fungi were Ascomycota and Basidiomycota. The dominant fungal genera were Chondrostereum, unidentified_Sordariales_sp, and Geopora. The relative abundance of the top 10 fungi at each classification level varied greatly in different growth stages. The study on the relationship between environmental factors and fungal community showed that EC was the main soil factor affecting the composition of rhizosphere fungal community, followed by TDS and AK.