This study was conducted in the Yellow River Delta, Shandong Province, China (118°33′E–119°20′E, 37°35′N–38°12′N) in 2019. This region has a temperate monsoon climate with an annual average temperature of 12.1°C and 196 frost-free days. The average annual precipitation is 530–630 mm, of which 70% is concentrated from June to August. The soil types are latent regional fluvial and saline. Five typical P. australis habitats were selected for soil sampling. Detailed information about habitat location is shown in Figure 1 and Supplementary Table S1. Surface soil sampling and plant community surveys were conducted in August, when P. australis was the most vigorous. At each of the three repeat sampling locations, five subsamples of topsoil (0–15 cm) were collected after removing the topsoil plant litter and thoroughly mixed to obtain a representative topsoil sample. We collected 15 samples, which were immediately stored in an ice box and transported to the laboratory. Each homogenised sample was divided into two parts: one was stored at −20°C for further microbial DNA extraction, and the other was air-dried to analyse its physicochemical properties.

The air-dried soil samples were filtered through a 2 mm mesh screen to remove visible plant roots and pebbles. The pH value of the water extract of each sample was determined using a distilled water ratio of 1:5 (w/v) and a pH metre (FiveEasy Plus^™, Mettler Toledo, Switzerland). The same water extract was used to determine the electrical conductivity (EC; FiveEasy Plus^™, Mettler Toledo, Switzerland). Soil organic matter (OM) content was determined using the potassium dichromate volumetric method. Based on the method provided by Zhu et al. (2021a), the total nitrogen (TN) and total phosphorus (TP) contents of the soil samples were determined using the Kjeldahl method and molybdenum antimony colorimetry, respectively. The soil extract extracted with a 2 M KCl solution was used to determine nitrate nitrogen and ammonium nitrogen using a continuous flow analyser (San⁺⁺ SYSTEM, Skalar, Netherlands). Available potassium (AK) was extracted by 1 M ammonium acetate solution and subsequently measured based on the atomic absorption spectrum.

The total DNA of the soil sample was extracted from a 0.500 g soil sample using the DNA Soil Kit (TIANGEN BIOTECH, Beijing, China) according to the manufacturer’s instructions. The quality and concentration of DNA were determined using a NanoDrop^® ND-1000 Spectrophotometre (Thermo Fisher Scientific Inc., Waltham, MA, United States) and 1.0% agarose gel electrophoresis. The DNA was stored at −20°C until further processing. Using the Illumina MiSeq platform by Majorbio Biotech. Co., Ltd., bacterial and fungal communities were assessed using 16S rRNA gene and internal transcribed spacer (ITS) high-throughput sequencing, respectively. Non-conserved 16S V3V4 regions were amplified using primers 341F (5-ACTCCTACGGGAGGCAGCAG-3) and 806R (5-GGACTACHVGGGTWTCTAAT-3), whereas ITS1F (5-CTTGGTCATTTAGAGGAAGTAA-3) and (5-GCTGCGTTCTTCATCGATGC-3) were used to amplify the ITS regions. To improve the quality of downstream analysis, the QIIME pipeline (version 1.7.0) was used to perform sequence stitching and quality trimming of raw data from Illumina MiSeq sequencing (Caporaso et al., 2010). After quality filtration, clean reads were clustered to generate operational taxonomic units (OTUs) with a 97% similarity identity cut-off using Uparse (version 7.0.1090). The most abundant sequence in each OTU, taken as a representative sequence, was selected to annotate taxonomic level using RDP classifier (version 2.11) based on the Silva database (Version 138). Unidentified and non-microbial sequences were defined as ‘unclassified’ for subsequent analysis.

According to Xu et al. (2021), based on the relative abundance and occurrence frequency of OTUs in samples, they were divided into six classifications: always abundant taxa (AAT) with relative abundance >1% in all samples, conditionally abundant taxa (CAT) with relative abundance >1% in some samples and never <0.01%, always rare taxa (ART) with relative abundance <0.01% in all samples, conditionally rare taxa (CRT) with relatice abundance <1% in all samples and <0.01% in some samples, moderate taxa (MT) with relative abundance between 0.01 and 1% in all samples, and conditionally rare and abundant taxa (CRAT) with relative abundance ranging from rare (< 0.01%) to abundant (>1%). Among these classifications, AAT, CAT, and CRAT were considered abundant taxa, whereas ART and CAT were regarded as rare abundant taxa.

Based on the R (Version 4.0.2) software package ‘Hmis’, the correlation matrix of microbial community was obtained by calculating the Spearman correlation coefficient between bacterial and fungal OTUs detected in more than eight samples. The false discovery rate (FDR) method was used for multiple testing corrections of the generated p-values (Yoav and Yosef, 1995). Significantly related OTUs (p < 0.05, r > 0.8) were used to construct a co-occurrence network graph. Co-occurrence of the microbial community was visualised using Gephi (version 0.9.2). Furthermore, the network characteristics were calculated using Gephi. The topological role of each node was assessed using the threshold values of Within-module connectivity (Zi) and Among-module connectivity (Pi), as proposed by Zhu et al. (2021a). Moreover, based on topological properties of sub-networks calculated by the ‘subgraph’ function, we evaluated the correlations between soil physicochemical properties and topological properties of the network.

Statistical analysis and figures were carried out in R (Version 4.0.2). One-way analysis of variance (ANOVA) was conducted using the ‘agricolae’ package to assess the discrepancy of microbial alpha diversity and soil physicochemical properties. Based on the ‘vegan’ package, principal coordinate analysis (PCoA) and permutational multivariate analysis of variance (PERMANOVA) were conducted to visualise and digitise the dissimilarity in microbial community composition and function, respectively. Subsequently, the histograms and heatmaps were performed using the ‘ggplot2’ and ‘pheatmap’ packages, respectively. Procrustes analysis was applied to explore the correlations between the microbial community and soil physicochemical properties (TN, TP, OM, NH 4 + − N , NO 3 − − N , AK, pH, and EC) in the ‘vegan’ package. A chord diagram showing the distribution of links among the bacterial and fungal phyla in co-occurrence network was constructed using the ‘circlize’ package. Canonical correspondence analysis (CCA) was performed using the ‘vegan’ package. In addition, based on the Majorbio I-Sanger Cloud Platform¹ designed for bioinformation analysis, we predicted bacterial and fungal community functions using PICRUST2 and FUNGuild tools. PICRUST2 classifies bacterial community functions according to the metabolic pathways of bacterial genes in KEGG (Douglas et al., 2020). FUNGuild is an open annotation tool for parsing fungal community datasets by ecological guild (Nguyen et al., 2015). Similarly, linear discriminant analysis (LDA) combined with linear discriminant analysis effect size (LEfSe) was conducted to identify the biomarker microorganisms based on the cloud platform.

Phragmites australis is distributed in different habitats with different physicochemical properties in the Yellow River Delta. As shown in Table 1, most soil physicochemical indices exhibited significant differences among the five soil samples, except NH 4 + − N (p > 0.05). As the dominant stress factors affecting plant community in saline-alkali land, the values of pH and EC of soil in different P. australis ranged from 7.65 to 8.07 and 1.27 to 5.36 mS/cm, respectively. Among the soil nutrient indices, OM, TN, NH 4 + − N, NO 3 − − N, and TP ranged from 17.24 to 58.99 mg/g, 70.00 to 169.79 mg/kg, 3.08 to 3.81 mg/kg, 1.14 to 5.41 mg/kg and 1.01 to 1.36 mg/g, respectively.

In total, 512,759 high-quality bacterial reads and 615,503 high-quality fungal reads were obtained, after removing potential chimeras. Based on 97% sequence similarity, the bacterial and fungal sequences were clustered into 11,299 and 1,635 OTUs, respectively. The rarefaction curve of bacterial and fungal OTUs reached the saturation platform, which suggested that most of the bacterial and fungal taxa in the soil samples were detected (Supplementary Figure S1). The number of bacterial and fungal OTUs detected in each P. australis habitat ranged from 4,335 to 5,506 and 348 to 633, respectively, and the number of shared bacterial and fungal OTUs among P. australis habitats was 965 and 48, respectively (Figures 2A,B).

Alpha diversity indices (Shannon index and Chao1) are important indicators reflecting the distribution patterns of species. The observed species, Shannon and Chao1 indices of bacterial community among five habitats ranged from 845.33 to 1.302, 4.72 to 5.68, and 988.70 to 1532.15, respectively. Regarding the fungal community observed species, Shannon and Chao1 indices ranged from 150.33 to 280.67, 2.26 to 3.57 and 165.72 to 301.54, respectively (Table 2). There were no significant differences in the observed species numbers, Shannon index, or Chao1 index of both the bacterial and fungal communities among the five P. australis habitats (ANOVA, p > 0.05).

A PCoA based on Bray–Curtis distance metrics was constructed to depict the differences in microbial communities among the P. australis habitats (Figures 2C,D). Coordinate axes 1 and 2 explained 16.9 and 14.5% of the variation in the bacterial community and 20.9 and 11.1% of the variation in the fungal community, respectively. This showed that there were clear bacterial and fungal distribution patterns in different P. australis habitats. PERMANOVA indicated that the bacterial communities of the five P. australis habitats were significantly different from each other (R² = 0.429, p < 0.05), while the whole fungal community showed no significant difference (R² = 0.314, p > 0.05), indicating that different habitats of P. australis have a significant influence on the bacterial community but not on the fungal community.

After classification analysis of 16S rRNA and ITS genes sequences, a total of 57 bacterial phyla and 9 fungal phyla were identified. The five dominant bacterial phyla (Proteobacteria, Chloroflexi, Bacteroidetes, Acidobacteria, and Gemmatimonadetes) accounted for 76.8–90.3% of the identified bacterial sequences (Figure 3A). Among them, Proteobacteria had the highest relative abundance in the P. australis soil, followed by Chloroflexi or Bacteroidetes. Bacteroidetes with the second highest relative abundance in S1 and S4, while Chloroflexi had the second highest relative abundance in S2, S3, and S5. For fungal community, dominant Ascomycota and Basidiomycota, representing 16.5–98.7% of identified fungal sequences within 15 soil samples (Figure 3A). Among them, Ascomycota was found to have the highest relative abundance in all the soil samples, with relative abundances ranging from 11.8 to 96.3% (Figure 3B).

At the family level, we focused on the responses of bacterial and fungal families with a top 30 relative abundance in P. australis habitats (Figures 3C,D). The top 30 bacterial families accounted for 30.0–65.4% of the total sequences, in which 15 families belong to Proteobacteria, 5 families belong to Bacteroidetes, 4 families belong to Acidobacteria, and 4 families belong to Chloroflexi. In terms of average relative abundance, the most predominant bacterial family was Methylophagaceae, with an average value of 4.7%, followed by norank_class__Gemmatimonadetes (4.0%), A4b (3.8%), norank_order__SBR1031 (3.2%), Flavobacteriaceae (3.0%), Nitrincolaceae (2.5%), Desulfobulbaceae (2.4%), and Kiloniellaceae (2.0%). Similarly, the top 30 fungal families accounted for 78.1–99.7% of the total sequences, in which 18 families belong to Ascomycota and 8 families belong to Basidiomycota. The most predominant family was unclassified_kingdom_fungi (23.3%), followed by unclassified_class__Sordariomycetes (17.9%), Pleosporaceae (5.9%), Aspergillus (5.6%), unclassified_order_Agaricales (5.3%), Glomeraceae (4.4%), and Cladosporiaceae (4.0%).

LEfSe tools were used to investigate biomarkers that caused differences among different P. australis habitats (Supplementary Figure S2). The results showed that 21 bacterial clades exhibited significant differences in different P. australis habitats with an LDA threshold of 3.5, whereas the number of fungal clades was 20 (Supplementary Figures S2C,D). In each P. australis habitat, the significantly enriched bacterial and fungal taxa varied. Among them, the number of bacteria and fungi enriched in the S5 plot was the largest, at 9 and 12, respectively. Moreover, no enriched fungal taxa were detected in S4.

The functional characteristics of the soil bacterial and fungal communities in different P. australis habitats were predicted using PICRUSt2 and FUNGuild, respectively (Figure 4). In bacterial community function, six categories at level 1 containing metabolism (65.7–68.7%), genetic information processing (10.6–11.6%), environmental information process (7.1–8.5%), cellular processes (6.5–8.2%), human diseases (4.0–4.8%), and organismal systems (2.0–2.2%) were identified (Figure 4C). At level 2, 45 pathways were identified, of which most abundant functional categories were global and overview maps (12.2–12.9%), carbohydrate metabolism (11.9–12.9%), amino acid metabolism (10.4–11.7%), energy metabolism (6.8–7.6%), and metabolism of cofactors and vitamins (6.3–6.6%). PCoA based on the Bray–Curtis distance revealed that different P. australis habitats significantly affect soil bacterial community function (PERMANOVA, p < 0.05; Figure 4A).

The FUNGuild results showed a total of 68 guilds belonging to 8 trophic modes: saprotrophs, saprotroph-symbiotrophs, pathotroph-saprotrophs, pathotrophs, symbiotrophs, pathotroph-saprotroph-symbiotroph, pathotroph-symbiotroph, pathogen-saprotroph-symbiotroph, and pathogen-saprotroph-symbiotroph (Figure 4D). Generally, the information on fungal species data is far less than that of bacterial data, resulting in a large amount of fungal OTU classification information that has not been identified. Among the identified fungal guilds, undefined saprotrophs were the most predominant guild (mean relative abundance accounting for 13.5% of the total relative abundance), followed by animal pathogen-endophyte-plant pathogen-wood saprotrophs (5.6%), arbuscular mycorrhizal (4.8%), animal pathogen-endophyte-lichen parasite-plant pathogen-wood saprotrophs (4.3%), and animal-pathogen-plant pathogen-undefined saprotrophs (2.2%), among others. Similar to the discrepancy in microbial community composition, the function of the soil fungal community was not sensitive to P. australis habitats (PERMANOVA, p > 0.05; Figure 4B).

To identify potential interactions between bacteria and fungi in the soil of P. australis, a co-occurrence network of microbial communities was constructed (Figure 5). Several important topological properties of the microbial community co-occurrence network were listed in Table 3. The network exhibited non-random co-occurrence patterns, as indicated by the R² of the power-law higher than 0.7. Moreover, the values of the average clustering coefficient (ACC) and average path length (APL) were much higher than those of random network, indicating that the network has ‘small world’ properties. In the soil microbial co-occurrence network of P. australis, bacteria and fungi mainly showed positive synergistic interactions, as indicated by the 95.17% positive edges (Supplementary Table S2). In terms of the co-occurrence network composition, 86.6% of the nodes were affiliated with Proteobacteria (46.4%), followed by Chloroflexi (14.4%), Bacteroidetes (9.5%), Gemmatimonadtes (8.2%), and Acidobacteria (8.1%) (Figure 5B). Correspondingly, the distribution of links in the co-occurrence network conformed to the distribution of nodes. Bacteria-involved links accounted for 97.3%, whereas fungi–fungi links only accounted for 0.5%. Among these links, those derived from Proteobacteria and Chloroflexi with other phyla were predominant, accounting for 69.3 and 27.1%, respectively (Figure 5B).

The topological properties of individual nodes were determined based on Zi-Pi analysis (Figure 5C). According to the within-module connectivity (Zi) and among-module connectivity (Pi) of individual nodes, nodes were divided into four classifications: peripherals (Zi < 2.5, Pi < 0.62), connectors (Zi < 2.5, Pi > 0.62), module hubs (Zi > 2.5, Pi < 0.62), and network hubs (Zi > 2.5, Pi > 0.62) [Olesen et al., 2007]. Among these four classifications, connectors, module hubs, and network hubs were identified as generalists, and were also considered to be keystone taxa that maintain the balance of co-occurrence networks (Montoya et al., 2006). In the present study, only 2.3% of the nodes were affiliated with generalists in the co-occurrence network. Among these 24 generalists, 23 nodes were affiliated with bacteria, and 1 node was affiliated with fungi. It is worth noting that most bacterial keystones were rare taxa, whereas the only fungal keystone was abundant rare taxa (Supplementary Table S3).

To better understand the causes of changes in microbial communities, we conducted a Procrustes analysis of the relationships between the microbial community composition and soil physicochemical properties (Figure 6). The results revealed that bacterial community composition was significantly correlated with soil characteristics (M² = 0.731, p < 0.05), whereas the fungal community composition was not significantly affected by soil characteristics (M² = 0.893, p > 0.05). The measured soil physicochemical properties explained 23.14% of the variation in the bacterial community and 19.70% of the variation in the fungal community, as indicated by CCA (Figure 6). Among these physicochemical properties, AK, OM, and NH 4 + − N played important roles in affecting the bacterial community composition, while NO 3 − − N and OM played vital roles in shaping the fungal community. Pearson correlation analysis showed that the Shannon indices of bacterial and fungal communities were significantly negatively correlated with NH 4 + − N and NO 3 − − N , respectively (p < 0.05; Supplementary Figure S3). However, analysis of microbial community function and soil physicochemical properties revealed that bacterial (M² = 0.623, p < 0.05) and fungal (M² = 0.665, p < 0.05) community functions were significantly associated with soil physicochemical properties. Moreover, among all soil properties detected, NO 3 − − N was the only soil physicochemical property that significantly affected the topological properties of the co-occurrence network of the microbial community, but the network modularity was not affected (Figure 7).

Generally, microbial community composition in the Yellow River Delta is affected by changes in soil properties and vegetation characteristics caused by human disturbance and habitat heterogeneity (Bickel and Or, 2020; Chi et al., 2021; Wu et al., 2021). Surprisingly, neither bacterial nor fungal ɑ-diversity, such as the Shannon and Chao1 indices, were affected by the habitats with different soil properties (Tables 1 and 2). Similarly, other research conducted in the Yellow River Delta showed the bacterial chao1 index had no significant difference between the supratidal and intertidal habitats of the Suaeda salsa community (Liu et al., 2020). In the present study, although soil TN and different forms of N ( NO 3 − − N and NH 4 + − N ) significantly related to the diversity of bacteria and fungi, differences in these nutrients between habitats did not cause significant differences in microbial diversity (Supplementary Figure S3). In terms of the composition of microbial communities, our results showed that bacterial communities were significantly affected by different P. australis habitats, whereas fungal communities were not (Figures 2C,D). This result can be explained by the relationships between soil physicochemical properties and microbial communities, as indicated by Procrustes analysis (Figures 5A,B), that is, soil physicochemical properties drove variation in bacterial communities. Liu et al. (2020) highlighted that compared with bacterial diversity, bacterial community composition was more sensitive and responds more quickly to a variety of soil physicochemical properties. Another study evaluating global surface soil microbial communities also supported our conclusion that soil physicochemical properties had a greater impact on bacteria than fungi, as geographic distance played an important role in fungal community construction (Bahram et al., 2018). Furthermore, Sweeney et al. (2021) revealed that plant phylogenetic relatedness and plant traits were important factors that determine rhizosphere fungal community composition. The number of samples in this study may also be a factor contributing to the fungal community composition being unaffected by P. australis habitat. Compared to the low intra-group variation in bacterial communities, the greater intra-group variation in fungal communities may result in fungal communities that were unaffected by P. australis habitat. Therefore, consideration should be given in subsequent studies to increasing the number of samples in order to reduce the potential limitations posed by the small number of samples. Combined with the above results, we suggested that plant types in a community were more likely to cause differences in the soil fungal community than soil physicochemical properties.

In terms of microbial community composition, the most dominant bacterial phylum was Proteobacteria, which has also been reported in previous studies conducted in the Yellow River Delta (Gao et al., 2021; Wu et al., 2021). Various substrate utilisation capabilities, energy acquisition methods, strong environmental adaptability, and rapid growth have been considered to help Proteobacteria become the predominant bacterial phylum in the Yellow River Delta (Bernhard et al., 2005; Mukhopadhya et al., 2012). Moreover, other dominant bacterial phyla have also been reported to play key roles in elemental cycling and energy flow, such as Acidobacteria and Bacteroidetes, which participate in the nitrogen cycle; Acidobacteria and Chloroflexi are the dominant decomposers of OM (Eichorst et al., 2018; Fang et al., 2018; Zhai et al., 2018). The composition of the fungal community was simpler than that of the complex bacterial community. Owing to the strong environmental adaptability and reproductive ability of Ascomycota, its relative abundance was an absolute advantage in each sample, which has also been reported in a variety of habitats, such as pioneer forests and tropical oligotrophic peatlands (Zhong et al., 2018; Morrison et al., 2021).

At the family level, the bacterial community was mainly composed of Methlophagaceae, Desulfobulbaceae, Flavobacteriaceae, Nitrincolaceae, and Marinobacteraceae, which were widely distributed in many coastal wetland soils (Zhang et al., 2012; Risgaard-Petersen et al., 2015; Chen et al., 2020; Yang et al., 2021). Desufobulbaceae, regarded as ‘cable bacteria’, was involved in the formation of a centimetre-deep suboxic zone in marine sediment through their role in the electrical couple of sulfide oxidation and oxygen reduction (Risgaard-Petersen et al., 2015). Methylotroph Methlophagaceae, an obligate that uses C1-compounds contaning methanol and methylamine produced by methane monooxygenase bacteria, had a relatively high relative abundance in all samples, indicating the adaptation of the bacteria to the unique environment of the coastal wetland (Neufeld et al., 2008; Moussard et al., 2009). Moreover, the high content of NO 3 − − N in S1 enabled Nitrincolaceae with strong heterotrophic nitrate reduction to ammonium ability to be recruited in large amounts, which was beneficial to the growth of plants (Yang et al., 2021). Thus, the significant differences among functional bacteria caused by different habitats significantly affected bacterial community function (Supplementary Figure S2; Figure 4). As the amount of data in the fungal classification database is much less than that of bacteria, multiple fungi with high relative abundances have not been effectively classified at the family level, such as unclassified_kingdom_Fungi, unclassified _class_Sordariomycetes, and unclassified_phylum_Ascomycota (Figure 3D). Among fungal families with good classification information, salt tolerance, and adaptation to anaerobic environments, fungal phyla Pleosporaceae, Aspergillaceae, and Cladosporiaceae were found to have high relative abundances in all habitats, being widely distributed despite the different degrees of saline-alkali soil in the Yellow River Delta (Yang and Sun, 2020). Moreover, Glomeraceae, which is involved in mycorrhizal formation, was largely enriched in S1, S2, and S4, but was rarely detected in other habitats, indicating the existence of mycorrhiza in the root system of P. australis; however, it was affected by different habitats (Ban et al., 2017).

As ecosystem decomposers, bacteria and fungi directly participate in biogeochemical cycles and energy fluxes, and changes in their compositions and functions directly and indirectly affect ecosystem function (Treseder and Lennon, 2015; Laforest-Lapointe et al., 2017; Wagg et al., 2019). The results of PERMANOVA and Procrustes analysis revealed that bacterial community function was affected by habitat (p < 0.05) and soil physicochemical factors (M² = 0.623, p < 0.05) (Figure 4A). Soil physicochemical properties can significantly affect the function of the bacterial community by altering the composition of the bacterial community or influencing the activity of bacteria, such as the nitrogen fixation capacity of Trichodesmium spp., which was affected by the content of nitrogen (Mulholland and Capone, 2000). Among the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway groups at level 1, metabolism was the most predominant process, especially carbohydrate metabolism, amino acid metabolism, cofactor, and vitamin metabolism, which were the top three in total relative abundance, indicating high C and N turnover rates in coastal wetlands, which was also reported by Zhang et al. (2016) and Liu et al. (2020). Carbohydrate metabolism, involving the degradation of hemicellulose and cellulose, can release various inorganic substrates and readily available substances to promote microbial and plant growth (López-González et al., 2015). It is not surprising that the prevalence of amino acid metabolism was due to the developed root system containing rhizomes of P. australis, which can release root exudates containing amino acids into the soil (Coops and Velde, 1995; Weidenhamer et al., 2013). Microbial metabolism of amino acids not only enhances their growth and development but also promotes plant growth by releasing inorganic nitrogen. Moreover, vitamins can be involved in carboxylation reactions, fatty acid metabolism, and amino acid metabolism as enzyme cofactors, and their metabolism can maintain metabolic homeostasis (Petri et al., 2019). Overall, the high relative abundance of metabolic pathways suggested that bacteria play a key role in the elemental cycle and biotransformation of OM.

The FUNGuild prediction results showed that the soil fungal community function of P. australis was not significantly different among the five habitats but was significantly affected by soil physicochemical properties. Considering that the effect of habitat on microbial community function is mainly through biotic (e.g., plant community) and abiotic factors (e.g., environmental factors), we speculate that the effect of the P. australis community on fungal community function was stronger than the effect of the environment on fungal community function (Waters et al., 2017). Saprophytic fungi are the main degradators of plant and animal residues, and their high relative abundance in the soil of P. australis habitats indocated that fungi play key roles in the degradation of OM and elemental cycling. As another relatively abundant fungal trophic mode, arbuscular mycorrhizal fungi (AMF) can form symbionts with roots to promote nutrient absorption and stress tolerance of plants (Berruti et al., 2016; Mathur et al., 2019). However, AMF were enriched only in S1, S2, and S4, indicating that, compared with 80% common infection, the proportion of P. australis roots and AMF forming symbiosis was still relatively low (Figure 4D). Moreover, a variety of plant and animal pathogens were classified in all habitats with high relative abundance, indicating that wetland habitats of P. australis were repositories of fungal pathogens (Chen et al., 2020). Similarly, the large intra-group variation caused by the small number of samples could be a potential reason for the lack of significant effects of P. australis habitat on fungal community function. In general, FUNGuild is a functional classification system based on existing studies and is limited by the amount of information in the fungal comparison database; as such, a more realistic understanding of fungal community function requires further exploration.

Soil microorganisms do not exist as individual species but form complex ecological networks through positive, negative, and neutral interactions, including but not limited to competition, commensalism, and mutualism (Faust et al., 2015; Graham et al., 2016). Menezes et al. (2017) highlighted the importance of a comprehensive understanding of the interactions between fungi and bacteria to improve soil ecosystem services. In the present study, we constructed intra- and inter-kingdom co-occurrence networks based on bacterial and fungal OTUs. The results showed that intra-kingdom interactions (bacteria, bacteria, fungi, and fungi) form more easily than inter-kingdom (bacteria–fungi) interactions, indicating that microorganisms were more susceptible to the influence of intra-kingdom microorganisms (Duan et al., 2021). Moreover, considering that the positive correlations between microbes were determined by their positive ecological interactions, such as commensalism and mutualism, it is not surprising that most interactions were positive because the degradation of OM requires the cooperation of various microorganisms, especially recalcitrant OM containing cellulose, hemicellulose, and lignin (Wang et al., 2016). Recently, microbial co-occurrence networks that were more inclined toward intra-kingdom symbiosis have been reported in a variety of habitats, such as the rhizosphere and mountainous areas (Duan et al., 2021; Qian et al., 2021).

A keystone microbial taxa is a highly connected taxonomic group in the co-occurrence network that significantly influence microbial community structure and function irrespective of their abundance across time and space (Banerjee et al., 2018). According to the Zi-Pi classification principle of keystone taxa and the identification principle of abundant/rare OTUs, we found far more bacteria than fungi in the identified keystone taxa, indicating that bacteria play a more critical role in maintaining the structural stability of the co-occurrence network (Zhu et al., 2021a). Furthermore, most of these bacterial keystone taxa were rare species, whereas the only fungal keystone identified was an abundant species. Indeed, a growing number of studies have highlighted the key role that rare species play in maintaining microbial community structure (Power et al., 1996; Xu et al., 2021; Zhu et al., 2021b). Generally, rare taxa of bacteria are highly active, as indicated by higher 16S rRNA/rRNA genes, and can increase rapidly under changing environments (Campbell et al., 2011). In addition, the lack of sensitivity of rare species to abiotic changes caused by the fact that these rare species are neither opportunistic nor take advantage of the reduction of dominant species may be the reason why rare species become the majority of keystones (Kurm et al., 2019). Notably, most rare bacterial species and abundant fungal species as keystone taxa that played important roles in maintaining the stability of intra-kingdom microbial co-occurrence networks have also been reported in cadmium-contaminated soil, but whether or not this phenomenon was random requires further study (Xu et al., 2021).

An increasing number of studies have focused on the relationships among soil physicochemical properties, ecological functions, and topological properties of microbial co-occurrence networks. Ma et al. (2016) revealed the importance of geographic distance, climatic factors, and soil physicochemical properties in network topological features. Zhang et al. (2018) emphasised that soil physicochemical properties, especially pH, available nitrogen, and available Mg, were significantly related to the topological properties of the microbial network. Moreover, the multifunctionality of the soil ecosystem was positively correlated with the network complexity index of the number of edges and nodes (Qiu et al., 2021). In this study, Spearman’s correlation analysis was conducted to explore the impact of soil physicochemical properties on the key topological properties of the network. Among the detected soil physicochemical properties, NO 3 − − N was the only one that could affect the topological properties of the microbial network (Figure 7). This was consistent with the high amino acid metabolism capacity of soil bacterial communities that dominate the co-occurrence network in different P. australis habitats (Figure 4C). Dissimilatory nitrate reduction to ammonium bacteria, such as Nitrincolaceae, was widely distributed in wetlands and can effectively transform the existing form of nitrogen, and thus affects the interactions between microorganisms (Zhao et al., 2020). Furthermore, as the only network topological feature unaffected by soil physicochemical properties, modularity may be more susceptible to the influence of phylogenetic relatedness, as indicated by the classification information of microorganisms with the same modularity being more similar (Hu et al., 2017).