In this study, a mid-channel bar (Taipingkou: 30°30′ N, 112°13′ E) located at the middle reaches of the Yangtze River was selected as to investigate the assembly processes of rhizosphere fungal communities (Figure 1). The mid-channel bar formed in the last 50 years is regarded as a newborn, which has the tendency to move downstream continuously due to constant flow erosion and sediment deposition at the head and tail, respectively (Wang et al., 2015; Wen et al., 2020). As a sandy bar in the middle of the river channel (Wang et al., 2015), Taipingkou is relatively less affected by anthropogenic pollution inputs. However, it is more influenced by the changes in the hydrological regime and sediment deposition caused by the upstream dam operation. The detailed climatic conditions of the study area can be found in the previous article (Ye et al., 2022).

Quadrat delineation and rhizosphere soil collection were performed during a survey on 26 May 2019. Along the mid-channel bar, 5 × 5 m quadrats were setup, which were uniformly distributed in 50 × 50 m sampling grids with a 25-m interval in between (Supplementary Figure 1). A total of 20 quadrats from three elevation-based transects were selected for plant survey and rhizosphere soil collection, including eight quadrats for transect one, seven quadrats for transect two, and five quadrats for transect three (Figure 1). In each sampling quadrat, the above-ground traits, including the number and coverage of each plant were recorded. The species-specific height of each plant was obtained by measuring the five randomly selected individuals. The importance value of each plant was calculated based on its relative height, frequency, and coverage (Ye et al., 2022). The dominant plants in each sampling quadrat were identified according to their importance values.

Five dominant plant species (Phragmites communis, Cynodon dactylon, Hemarthria altissima, Phalaris arundinacea, and Triarrhena sacchariflora) were confirmed across the sampling quadrats (Figure 1 and Supplementary Figure 2). P. communis, C. dactylon, H. altissima, P. arundinacea, and T. sacchariflora was dominant in four, four, seven, twelve and three quadrats, respectively. The values of above-ground traits (i.e., height, richness, coverage, and importance value) of these five plant species are presented in Supplementary Figure 2. Later, the importance value was calculated as the average of the sum of relative height, frequency, and coverage (Ye et al., 2022).

The rhizosphere soil samples were collected according to a previously reported method (Ye et al., 2021). Briefly, five plants of the same species were randomly selected and carefully removed from each quadrat. Soils loosely attached to the roots were gently shaken off. Soils tightly adhered to the roots were regarded as rhizosphere soils, and extracted using saline phosphate buffer. The rhizosphere soils of the same plant species in each quadrat were thoroughly mixed to form a composite sample.

The resulting 30 composite rhizosphere soils were subsequently divided into two parts and either stored at −20°C for DNA extraction or at 4°C for physicochemical analysis. The physicochemical property values of rhizosphere soils are presented in Supplementary Figure 2.

About 0.25 g of the rhizosphere soil from each composition sample was used to extract the DNA following the protocol of the DNeasy PowerSoil Kit (Qiagen, Hilden, Germany). The quality of the DNA was checked by 1% (w/v) agarose gel electrophoresis, and the concentration of the DNA was determined by a NanoDrop Lite Spectrophotometer (Thermo Scientific, Waltham, USA). The extracted DNA was stored in a −80°C refrigerator for subsequent analysis.

The fungal ITS gene was amplified by the polymerase chain reaction (PCR) using primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′)/ITS2R (5′-GCTGCGTTCTTCATCGATGC-3′) (Adams et al., 2013). The PCR products were purified using an AxyPrep DNA Gel Extraction Kit (Axygen, Union City, USA). The purified amplicons were then pooled in equimolar and paired-end sequenced on an Illumina MiSeq PE300 platform (Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China).¹ The raw sequence data were demultiplexed and quality-filtered using Trimmomatic (v. 0.30) (Bolger et al., 2014). The operational taxonomic units (OTUs) were identified at 97% sequence similarity using USEARCH (v. 7.0).² The fungal taxonomy was analyzed via the ribosomal database project (RDP) using the UNITE 8.0 database.³ The sequences in each sample were filtered to an equal sequencing depth based on the fewest sequence reads to compare the samples. Alpha diversity indices, such as Sobs, Shannon diversity, Heip’ evenness, and Phylogenetic diversity, were calculated using Mothur (v. 1.30.1). PCR, sequencing, and data processing were conducted according to a previously reported method (Ye et al., 2021).

A null model approach based on two indices of β-nearest taxon index (βNTI) and raup-crick index (RCI) was introduced to determine the assembly processes of soil fungal communities in different plant species. These two indices were calculated by the “picante” package in R v. 4.0.5 (Wan et al., 2021; Zhao et al., 2022). A βNTI value > 2 or < −2 indicates the dominance of deterministic processes (i.e., variable selection and homogeneous selection) on the assembly of a given microbial community. Similarly, the βNTI value < −2 suggests the effect of homogeneous selection and > 2 variable selection (Stegen et al., 2015). The value of βNTI ranging from −2 to 2 indicates the primary effect of stochastic processes (i.e., dispersal limitation, homogenizing dispersal, and undominated processes). The dispersal limitation is often regarded as the main driver of the community assembly with an RCI value larger than 0.95, while the homogenizing dispersal is identified as the primary contributor with an RCI value smaller than −0.95. Additionally, the RCI value ranging from −0.95 to 0.95 indicates the effect of undominated processes (Stegen et al., 2015).

The potential associations between rhizosphere fungal communities were explored by the co-occurrence network analysis. The network was constructed based on the random matrix theory (RMT) using the molecular ecological network analysis pipeline (MENA)⁴ at an OTU level. Only OTUs present in 10 samples (one–third of the total samples) were retained based on the rarefied fungal OTU abundance table. Subsequently, the relative abundance data of OTUs were uploaded into MENA to conduct network analyses with the default settings and recommended cutoff threshold of 0.81 (Zhou et al., 2011; Deng et al., 2012; Ye et al., 2021). Visualization and layout of the network were performed based on the Fruchterman-Reingold algorithm using Gephi 0.9.2 (Bastian et al., 2009). The topological properties of the constructed network were calculated to describe the associations between the fungal taxa (Ma et al., 2016). Moreover, 1,000 Erdös–Réyni random networks were generated on the pipeline with an identical number of nodes and edges as the constructed network. The potential ecological role of each node in the network was confirmed based on the threshold of within-module connectivity (Zi) and among-module connectivity (Pi) at 2.5 and 0.62, respectively (Deng et al., 2012; Wang et al., 2020; Ye et al., 2021). Furthermore, the nodes highly connected to other nodes in a module were defined as the module hubs with Zi > 0.25 and Pi < 0.62; the nodes connected to different modules were defined as the connectors with Zi < 0.25 and Pi > 0.62; the nodes acting as a module hub and a connector were defined as the network hubs with Zi > 0.25 and Pi > 0.62 (Zhang et al., 2021).

The statistical significances of dominant genera (relative abundance > 1% of the total sequences) and alpha diversity indices (i.e., Sobs, Shannon diversity, Heip’s evenness, and Phylogenetic diversity) of rhizosphere fungal communities among different plant species were determined by one-way analysis of variance (ANOVA) based on the least significant difference (LSD) using SPSS Statistics 20.0 (IBM, USA). The community distribution of rhizosphere fungi was assessed by the principal co-ordinates analysis (PCoA) based on the Bray–Curtis distance using CANOCO 5 software (Ter Braak and Šmilauer, 2012). The significant differences between rhizosphere fungal communities among different plant species were determined by analysis of similarity (ANOSIM) using R v. 3.5.1 with “vegan” package. The relationship between alpha diversity indices and environmental factors and the aboveground plant traits were determined by Spearman’s correlations using SPSS Statistics 20.0 (IBM, USA). The relationships between the fungal community composition and environmental factors, aboveground plant traits were explored by R v. 4.0.5 with “ggcor” package. The Pearson correlation coefficient was conducted in Origin 2022 (OriginLab, USA).

The raw sequencing data have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA)⁵ under the accession number PRJNA903773.

A total of 1,934,139 high-quality sequences were acquired from 30 rhizosphere soils. These sequences were resampled to an even sequence depth of 41,749 across all samples based on the lowest sequence of samples. Subsequently, 4,615 OTUs were identified based on the rarefied sequences with 97% similarity. Among the dominant fungal genera (sequences accounting for more than 1% of the total sequences of all samples), Alternaria and Epicoccum were the two most abundant taxa, accounting for 18.1 and 15.2% of the total sequences, respectively (Figure 2A). The relative abundance of Alternaria in rhizosphere of C. dactylon was significantly higher than that in rhizosphere of H. altissima (one-way ANOVA, P < 0.05), and not significantly different from that of P. communis, P. arundinacea, and T. sacchariflora (one-way ANOVA, P > 0.05; Figure 2A). The relative abundance of Epicoccum in rhizosphere showed no significant difference among different plant species (one-way ANOVA, P > 0.05; Figure 2A).

The alpha diversity indices (i.e., Sobs, Shannon diversity, Heip’s evenness, and Phylogenetic diversity) of rhizosphere fungal communities did not exhibit significant differences among different plant species (one-way ANOVA, P > 0.05; Figure 2B). Exceptionally, significant differences were found between the Sobs and Phylogenetic diversity of H. altissima and P. arundinacea with higher values in the rhizosphere of H. altissima (one-way ANOVA, P < 0.05; Figure 2B). The Sobs and Phylogenetic diversity indices of rhizosphere fungal communities showed significant positive correlations with rhizosphere soil variables of total nitrogen (TN), OM and the aboveground plant trait of richness, but were negatively correlated with soil pH and the aboveground plant trait of importance value (Spearman’s correlations, P < 0.05; Figure 2D).

The PCoA showed an unorganized distribution of rhizosphere fungal communities across different plant species (Figure 2C). The ANOSIM confirmed that there was no significant difference in the rhizosphere fungal communities among different plant species (P > 0.05; Figure 2C). Mantel tests showed that the rhizosphere fungal communities were significantly positively correlated with the rhizosphere soil variables of pH, TN, and total carbon (TC) (P < 0.05; Figure 2D).

In the middle-channel bar, the community assembly of rhizosphere fungi was dominated by the stochastic processes (85.3%), with the contribution of undominated processes and dispersal limitation being 54.5 and 30.1%, respectively (Supplementary Figure 3). Among different plant species, the undominated processes had greater contributions of 52.6, 60.3, 41.9, 62.4, and 47.1% to the fungal community assembly in the rhizosphere of P. communis, C. dactylon, H. altissima, P. arundinacea, and T. sacchariflora (Figure 3). Compared with undominated processes, dispersal limitation contributed most to the fungal community assembly in the rhizosphere of H. altissima (47.8%) and T. sacchariflora (52.9%) (Figure 3). Moreover, the homogeneous selection of deterministic processes exhibited indispensable effects on fungal community assembly in the rhizosphere of P. communis (15.5%), C. dactylon (22.4%), and P. arundinacea (17.8%).

The relationships of the βNTI of rhizosphere fungal communities with soil environmental factors and aboveground plant traits were determined to reveal their impacts on the assembly processes of fungal communities. The changes in rhizosphere soil environmental variables of TC, TN, OM, and pH had significant positive correlations with the pairwise comparisons of βNTI (Figure 4). This result suggested that an increasing discrepancy in these soil environmental variables promoted the effect of stochastic processes on the community assembly of rhizosphere fungi. Moreover, the change in aboveground plant traits of height was significantly positively correlated with βNTI (Figure 4). In contrast, the changes in flooding probability of sampling sites and aboveground plant traits of importance value showed significantly negative correlations with βNTI.

A molecular ecological network composed of 141 nodes with 441 edges (Figure 5A) was constructed based on the RMT to reveal the co-occurrence pattern of the rhizosphere fungal community. The distribution of node degrees for the constructed fungal network showed a scale-free power-law network property (R² = 0.832; Supplementary Figure 4), indicating a non-random structure of the constructed network. Additionally, the values of average path distance (2.667 vs. 2.393) and modularity (0.274 vs. 0.261) were higher than those of the respective randomized network (1,000 Erdös–Réyni), revealing the small-world properties and modular structure of the constructed network. The fungal network was dominated by Ascomycota accounting for 63.8% of the total nodes, followed by unclassified fungi (14.2%) (Figure 5A).

Based on the within-module connectivity (Zi) and among-module connectivity (Pi) of individual taxa within the constructed network, a total of 5, 1, and 51 nodes were identified as the network hubs, module hubs, and connectors, respectively (Figure 5B and Supplementary Table 1). All of these nodes could be regarded as the keystone species in structuring the rhizosphere fungal communities. Nearly half of the keystone species (47.4%) were affiliated to the unclassified fungi at the class level (Supplementary Table 1). Noteworthily, most connectors showed significant correlations with soil TN, OM, and pH, whereas the network hubs and module hubs were not significantly affected by the measured rhizosphere soil environmental variables and aboveground plant traits (Supplementary Table 2).