The experimental site location is State-owned Forest seed Farm in Yu’an District, Lu′an City, Anhui Province, China (116°12′E, 31°40’N). With an annual rainfall of 1239.8 mm, the climate of this area is warm and humid subtropical monsoon. The frost free period is nearly 225 days, with an annual average temperature of 15°C, the altitude is 80 ~ 110 m, and the slope is 5° ~ 15°. The soil is mostly clayey yellow soil, with a pH of 5.5 ~ 6.5. The P. massoniana provenance trial forest was established in the spring of 1981. There are 64 provenances in total with a completely randomized block design. The site is divided into six evenly divided districts (Supplementary Figure S1). Each district group contains 64 plots, i.e., 64 provenances, and each provenance has 9 plants as a plot (3 rows, 3 plants per row). The experimental forest is 41 years old and has not been thinned during the period, with a preservation rate of 36.52%. The number of surviving trees in some provenances cannot meet the requirements of statistical analysis due to the growth competition among individual trees, diseases and insect pests in the late growth stage. Eventually, there are 55 provenances that met the requirement of statistical analysis. We use the binary model in the Chinese forestry industry standard (LY/T2263-2014) to estimate the biomass of various sources (Zeng et al., 2010).

M_A: total biomass of aboveground organs for P. massoniana.

M_B: total biomass of the underground part (i.e., root) of P. massoniana.

D: diameter at breast height (DBH), H: height.

We used the proportional function weighting method to control and calculate the biomass of each organ and then summed it. Since the trunk accounts for the largest portion of the whole standing tree, the relative proportion of the trunk biomass to the total biomass was set as 1, and “g” was set as the proportion function for different organs to the total biomass. “M” was set as the biomass of different organs. We obtain the following regression equation:

We multiplied the biomass of different organs estimated by the model by the corresponding carbon content coefficient (Table 1) (Zeng and Xiao, 2011). The carbon storage of each organ was obtained. The carbon storage of each plant from various sources was obtained by summing. Based on the calculation results, we selected the first, twenty-eighth, and fifty-fifth ranked carbon storage among the 55 provenances representing high, medium and low carbon storage levels, respectively. They are Jiangxi Anyuan (AY), Fujian Shaowu (SW), and Zhejiang Shengxian (SX) namely (Table 2).

Four P. massoniana with the same DBH and close to the average value were selected from each provenance. Each tree was set with a 0 ~ 20 cm profile in three directions (120°) to collect fine roots with root diameters less than 2 mm. For the rhizosphere soil, a sterile brush was used to collect the soil within 5 mm of the root surface. Soil mixing of three direction profiles was taken as a repeat, with a total of 4 repeats. Soils were divided into three parts for the determination of different indicators. One part was used directly to determine soil enzyme activity, another part was air-dried for the determination of soil physical and chemical properties, and the last part was frozen at −80°C for microbial group analysis.

Using a glass composite electrode to determine pH value of soil with a water and soil ratio of 1:2.5. The moisture content of the soil was measured by drying it followed by cooling it. Determinations of soil organic carbon (SOC) and available phosphorus (AP) were traditional methods, namely, the potassium dichromate external heating method and molybdenum-antimony-scandium colorimetry (Yeomans and Bremner, 2008). A flow analyzer (Bran+Luebbe, GmbH, Germany) (Bunch and Bernot, 2012) was used to determine the total nitrogen (TN), nitrate nitrogen (NO₃⁻-N) and ammonium nitrogen (NH₄⁺-N) contents of the soil. Nitrite reductase (NIR) is a key enzyme in the nitrogen cycle that degrades nitrite to NO or NH₃. Peroxidase (POD) and polyphenol oxidase (PPO) are key enzymes in the carbon cycle and are mainly associated with the degradation of difficult-to-degrade lignin (Yin et al., 2016). The activities of NIR, PPO, and POD were determined by the 96-well plate microplate method (Saiya-Cork et al., 2002).

DNA was extracted from soil samples by using a DNA Isolation Kit (Bio Teke, Beijing, China). The primers 343F-798R (5′- TACGGRAGGCAGCAG -3′/5′- AGGGTATCTAATCCT- 3′) and ITS1F-ITS2 (5′- CTTGGTCATTTAGAGGAAGTAA-3′/5’-GCTGCGTTCTTCATCGATGC- 3′) were used to analyze the bacterial 16S rRNA V3-V4 region and the ITS genes of the fungi (Chen et al., 2020). PCR amplification and Illumina NovaSeq sequencing were conducted at Shanghai OE Biotech. Co., Ltd. (Shanghai, China).

The Illumina MiSeq platform was used to generate the original paired-end sequences. The project can be accessed in NCBI under BioProject ID PRJNA1018330. Primers for raw data sequences were cut out by cutadapt. Representative sequences and ASV abundance tables were generated by DADA2 using the Qiime2 (Bolyen et al., 2019) default parameters. Representative sequences and ASV abundance tables were generated. Silva database (Version 138) was used to annotate the species database. We normalized each sample to a uniform sequence number to analyze the alpha diversity and principal coordinate analysis (PCoA). The tags used for 16S analysis are 3,232 and the tags used for ITS analysis are 19,217. We used the ggplot2 package in R (Ginestet, 2011) to perform alpha indices, principal coordinate analysis (PCoA) of Bray-Curtis distances and community structure barplot. Redundancy analysis (RDA) and the Mantel test were conducted by the vegan package in R (Dixon, 2003). The taxa exhibiting the most significant intergroup variation were identified by linear discriminant analysis (LDA) effect size (LEfSe). Bacterial functions are predicted using the FAPROTAX databases.¹ FUNGuild databases² were used to annotate fungal functions.

Table 3 shows the soil property variations of the three provenances. The pH of the rhizosphere soil varied from 4.85 to 5.05. The pH of SX provenance soil was significantly higher than that of AY and SW provenances soils (p < 0.05). SOC, NO₃⁻-N, and NH₄⁺-N contents of the three groups were significantly different (p < 0.05), increasing with increasing carbon storage level. We also observed that the TN content of AY was significantly higher than SW. There was no significant change in soil moisture and AP. SW had significantly greater NIR than AY and SX (p < 0.05). Conversely, the POD of AY and SX was significantly higher than SW (p < 0.05). The PPO activity did not change significantly with increasing carbon storage level.

The alpha diversity of bacteria and fungi was estimated by Shannon, Simpson, Chao1, and Ace indices (Figure 1). For bacteria, alpha diversity did not show significant differences among the three provenances. However, the Chao1 and Ace indices generally showed a downward trend with increasing carbon storage level. For fungi, the order for both Chao1 and Ace estimators was AY > SX > SW. AY had significantly greater Chao1 and Ace estimators than SW (p < 0.05). The Simpson and Shannon indices of fungi showed an upward trend with increasing carbon storage level.

For bacteria, the dominant phyla that we identified were Acidobacteriota (44.27%), Proteobacteria (32.57%), Actinobacteriota (13.43%), Desulfobacterota (1.79%), Gemmatimonadota (1.65%), Myxococcota (1.55%), and Bacteroidota (1.43%). Other species, such as Verrucomicrobiota, Firmicutes and Nitrospirota, accounted for only a small part of identified bacteria (Figure 2A) (Supplementary Table S1). For fungal phyla, the most abundant species were Basidiomycota (73.36%) and Ascomycota (24.64%). Zygomycota, Rozellomycota, Chytridiomycota, and Glomeromycota accounted for only a small part of identified fungi (Figure 2B) (Supplementary Table S1).

At the genus level, the dominant bacteria genus that we identified were Subgroup_2 (16.01%), Bryobacter (4.48%), Acidothermus (4.44%), Candidatus (4.10%), Burkholderia (2.96%), Acidibacter (2.19%), Bradyrhizobium (1.95%), Occallatibacter (1.26%), Ellin6067 (1.24%), A21b (1.15%) (Figure 2C). Fungi ASVs mainly comprised Lactarius (16.71%), Rhizopogon (16.53%), Pulveroboletus (3.93%), Tomentella (3.59), Delicatula (3.26%), Sebacina (2.66%), Hydnum (2.55%), Russula (2.02%), Archaeorhizomyces (1.49%) (Figure 2D).

We performed LEfSe analysis to identify the major biomarkers in P. massoniana rhizosphere soil under different carbon storage levels (LDA score > 4) (Figures 3A,B) (Supplementary Table S2). Bacterial taxa mostly included Acidobacteriota, Proteobacteria and Actinobacteriota. Fungi groups were mostly from Basidiomycota and Ascomycota. SX was enriched in five bacterial taxonomic groups (Acidobacteriales, Acidobacteriota, Acidobacteriae, Alphaproteobacteria and Acidobacteriaceae_Subgroup_1) and two fungal taxonomic groups (Boletaceae and Pulveroboletus). The main biomarkers of SW were Actinobacteriota, Thermoleophilia, Gaiellales, Gemmatimonadota, Gemmatimonadales, Boletales, Rhizopogonaceae, Rhizopogon and some uncultured families and groups from Gaiellales, Elsterales, and Gemmatimonadaceae. Taxa with higher abundance in AY were Archaeorhizomycetales, Archaeorhizomycetes, Sebacinaceae, Sebacinales, Athelicaceae, Atheliales, Tremellodendr, Thelephor, and some unidentified families from Archaeorhizomycetales and Atheliaceae.

PCoA based on the Bray-Curtis distance indicated that the community structures of the bacteria and fungi changed at different carbon storage levels (Figures 3C,D). The ANOSIM also supported the result that there were significant differences among the three groups (p = 0.001) (Supplementary Figure S2). For bacteria, there were significant differences between AY vs. SW, SW vs. SX, and AY vs. SX (p = 0.034, 0.02, 0.04, respectively). For fungi, there were significant differences between AY vs. SW (p = 0.026) (Supplementary Figure S2).

The correlation heatmap (Figure 4) revealed that the phyla abundance of bacteria and fungi had no significant correlation with soil AP contents (p > 0.05). For fungi, Basidomycota was significantly negatively correlated with soil SOC, TN, and NO₃⁻-N contents (p < 0.05). While, Ascomycota was on the opposite side which had significant negative correlation with soil SOC, TN, and NO₃⁻-N contents (p < 0.05). The abundance of Rozellomycota had significant negative correlation with soil pH and the activity of POD (p < 0.01). For bacteria, Acidobacteria had significant negative correlation with C/N and NIR, positive correlation with POD and PPO (p < 0.05). Conversely, Gemmatimonadota had significant positive correlation with C/N, NH₄⁺-N and NIR, negative correlation with POD and PPO (p < 0.05). Myxococcota was significantly positively correlated with moisture (p < 0.01). Actinobacteriota was significantly negatively correlated with pH (p < 0.05) and POD (p < 0.01).

The redundancy (RDA) analysis revealed the effects of soil factors on microbial community structures (Figures 5A,B). The first and second RDA axes accounted for 34.99% of the variance in the bacterial communities (Figure 5A). Soil C/N, POD, PPO and NiR activities were the main factors influencing the bacterial communities (p = 0.011, 0.005, 0.02, and 0.001, respectively). The Mantel test showed that SOC (r = 0.263, p = 0.025), C/N (r = 0.379, p = 0.008), NH₄⁺-N (r = 0.432, p = 0.005), and POD (r = 0.568, p = 0.001) had significant correlation with effects on bacterial community structures (Figure 5C).

The first two RDA axes accounted for 42.93% of the variance in fungal communities (Figure 5B). Soil SOC content were the main factor influencing the fungal communities (p = 0.001). Furthermore, the Mantel test suggested that fungal community structures were significantly correlated with SOC (r = 0.425, p = 0.003), TN (r = 0.374, p = 0.005), NH₄⁺-N (r = 0.318, p = 0.025), NO₃⁻-N (r = 0.383, p = 0.007), and POD (r = 0.332, p = 0.015) (Figure 5C). In addition, the mantel test revealed that the carbon storage level of P. massoniana provenance significantly affected the community structure of fungi and bacteria (p < 0.01) (Figure 5C).

The bacterial function was predicted based on FAPROTAX (Figure 6A). Core functions in the bacterial community were chemoheterotrophy (29.61%), aerobic_chemoheterotrophy (28.24), cellulolysis (12.48%), animal_parasites_or_symbionts (7.47%), human_pathogens_all (7.09%), human_pathogens_pneumonia (6.43%), fermentation (2.04%), nitrate_reduction (1.29%), and predatory_or_exoparasitic (1.00%). Compared with AY and SX, SW had significantly higher aromatic_compound_degradation, aerobic_nitrite_oxidation and nitrification (p < 0.05). AY had significantly higher nitrogen_storage than SW and SX (p < 0.05), along with significantly lower nitrogen_respiration, chemoheterotrophy and nitrate_respiration. Fungal communities were divided into different trophic modes, including pathotrophs, saprotrophs, and symbiotrophs, based on the FUNGuild database (Figure 6B). Compared with SW and SX, AY provenance had significantly more lichenization, lichen parasites and plant pathogens (p < 0.05).

Our research suggested that community diversity and structures of soil microorganisms in the rhizosphere soil of P. massoniana provenances changed at different levels of carbon storage. Studies have shown that because the carbon metabolism cycle of fungal mycelia is shorter than that of bacteria (Staddon et al., 2003), the fungi response to the carbon storage effect of plant photosynthesis is more obvious than that of bacteria (Jones et al., 2009), which our research confirms. We did not find significant difference in bacterial alpha diversity of P. massoniana rhizosphere soil at different carbon storage levels. The Chao1 and Ace diversity of fungi (Figure 1) in the P. massoniana rhizosphere soil with a high carbon storage level was significantly higher than that with a low carbon storage level. Researchers (Huang F. et al., 2019) have reported that higher fungal diversity can ensure agroecosystem stability and sustainable crop production, which indicated that the rhizosphere ecosystem of AY P. massoniana provenance representing high carbon storage levels has greater stability than lower carbon storage levels. Our results indicated that under a high carbon storage level, the Chao1 and ACE values of fungi in the P. massoniana rhizosphere soil were significantly increased, but the Simpson index and Shannon index did not change significantly. Thus, the carbon storage level of P. massoniana affected the number of fungal species but not diversity.

The PCoA analysis indicated that the microbial community compositions and structures significantly changed under different levels of carbon storage of P. massoniana provenances. The bacteria that shifts across carbon storage levels were mainly Acidobacteriota, Proteobacteria and Actinobacteriota. Generally, there are extensive studies on Proteobacteria in soil science. Proteobacteria is a diverse bacterium preferring to grow in soil with high effective carbon levels (Fierer et al., 2007). Some Proteobacteria strains are plant growth-promoting bacteria and participate in the nitrogen storage process (Hassen et al., 2011). Actinobacteriota are low-nutrient species that prefer nutrient-poor environments (Liu et al., 2019). These bacteria usually predominate in dry, low nutrition and poor conditions (Maestre et al., 2015). Studies (Eichorst et al., 2018) have shown that Acidobacteriota uses organic nitrogen, inorganic nitrogen, as well as various carbohydrates, as the nitrogen source. It plays a vital role in organic matter decomposition as well as nutrient cycling. The results of FAPROTAX-based online functional prediction (Figure 6C) shows that the carbon storage level of P. massoniana provenance changed the carbon and nitrogen metabolism in the rhizosphere, increased the nitrogen storage effect, and decreased the nitrite respiration, chemical energy heterotrophic and nitrate respiration.

The fungal taxa with high relative abundance were Ascomycota and Basidiomycota in the studied soils (Zhang X. et al., 2020). The PCoA and LEfSe analysis results revealed that fungal community structures significantly shifted across different carbon storage levels, especially for Ascomycota and Basidiomycota. Basidiomycota can effectively degrade lignocellulose organics, participating in litter decomposition (Yelle et al., 2008). Therefore, we inferred that variation in Basidiomycota was directly related to the change in soil SOC content. In addition, we found that Ascomycota in the P. massoniana rhizosphere had more significant competitive advantages at high carbon storage provenances. Ascomycota is the key driving factor of the carbon and nitrogen cycles (Cook, 2015). Ascomycota plays an important role in P. massoniana rhizosphere microorganism response to carbon storage levels through its role in these cycles. Additionally, FUNGuild analysis (Figure 6C) showed that at high carbon storage levels, the lichenized fungi in the P. massoniana rhizosphere increased, which indicated that the nutrient pattern of fungi had changed.

The LDA analysis of bacteria showed that Acidobacteriales, Acidobacteriae, Acidobacteriae_Subgroup_1, Alphaproteobacteria were mainly enriched in SX. Most of these taxa belong to the Acidobacteria phylum. It has been shown that Acidobacteria usually grow in soils with low nutrient levels (Pankratov and Dedysh, 2010). SX is low carbon storage provenance with low content of organic carbon, ammonium nitrogen, nitrate nitrogen in the rhizosphere. It can enrich by acidophilus phylum to decompose cellulose and other carbon compounds to obtain organic carbon and other nutrients required for growth. The biomarkers of SW are Thermoleophilia, Gaiellales, Elsterales, and Gemmatimonadaceae. Compared to the SX group, SW representing medium carbon storage had higher effective nitrogen content, which may have led to the enrichment of Thermoleophilia and Gaiellales (Ramirez, 2011). Besides, Gemmatimonadaceae is plant-promoting bacteria. It protect plant health and promote plant growth by secreting antibiotics that inhibit or kill microorganisms that are harmful to plants, thus indirectly leading to an increase in plant carbon storage (Hui et al., 2020).

The LDA analysis of fungi showed that the biomarkers of SX were Boletaceae and Pulveroboletus, which were suitable for growth between pH of 5–6 (Wu et al., 2014). The pH of SX was 5.05, slightly higher than that of SW and AY. Our hypothesis is that the increase in pH caused by SX may have led to an enrichment of Boletaceae and Pulveroboletus. SW is enriched with more Rhizopogonaceae. This results suggested that medium carbon storage provenance obtained more nitrogen nutrients by altering the biological strategy of the fungi compared to low carbon storage provenance (Zeng et al., 2021). In addition, more diverse fungal taxa were enriched in AY provenance, including Archaeorhizomycetales, Sebacinaceae, Archaeorhizomycetales, Atheliaceae, Tremellodendron, and Thelephora. Archaeorhizomyces has been demonstrated to have the ability to utilize either glucose or cellulose as a carbon source for growth (Menkis et al., 2014). The rhizosphere soil of high carbon storage provenance exhibits elevated levels of organic carbon content, which is influenced by apoplastic material and root secretion. This may result in an increased proportion of decomposers. The presence of Atheliaceae, Sebacinaceae, Tremellodendron, and Thelephora ectomycorrhizal fungi suggested that high carbon storage provenance exhibit greater diversity in nutrient acquisition strategies (Weiss et al., 2004; Smith et al., 2009).

In our research, RDA and Mantel tests revealed that the NH₄⁺-N, PER contents and C/N had the most significant effect (p < 0.01) on bacterial community structures. A potential reason why C:N ratio can significantly affect bacterial communities may be that microorganisms disrupt elemental stoichiometric balance and homeostasis (Bi et al., 2021). Therefore, the C:N ratio is an important marker of changes in microbial communities. We found that NH₄⁺-N had a significant correlation with the bacterial communities in the rhizosphere soil, which indicated that variables associated with nitrogen transformations may play an important role in the changes of bacterial community structure (Zhang et al., 2021a). Some studies believe that for subtropical forest ecosystems with generally acidic soil, compared to NO₃⁻-N, NH₄⁺-N has a more important effect on the bulk soil and soil aggregates, which further affects the microbial community structure (Yan et al., 2017). Studies has shown that increasing the content of NH₄⁺-N will enable more bacterial groups to participate in the recycling of carbon and nitrogen. In addition to the correlation with humus formation (Yan et al., 2017), the POD enzyme also had a significant positive correlation with Actinobaciota, Gemmatimonadota and Rozellomycota. Besieds, RDA and Mantel tests showed that the SOC, TN, and NO₃⁻-N contents had the most significant effect (p < 0.01) on the fungal community structure. As the energy source of microorganisms, the increase in SOC content will significantly affect the microbial metabolic cycle (Sokol and Bradford, 2019). Research has proven that SOC is the direct influencing factor regulating microbial community structure (He et al., 2014). The significant effect of TN on the fungal community structure indicated that TN may limit the carbon utilization by microorganisms (Allen and Schlesinger, 2004). In this study, the nitrate decomposition rate decreased with the increase in P. massoniana carbon storage levels, which indicates that plant carbon storage levels may change the nitrification rate of rhizosphere soil. The significant impact of NO₃⁻-N on the fungal community structure indicated that nitrogen conversion rates of rhizosphere soil played an important role in shaping microbial communities (Wang et al., 2022). Generally, variations in carbon storage levels among P. massoniana provenances lead to differences in the carbon cycle of rhizosphere microorganisms in the P. massoniana forest ecosystem, which also leads to the nutrient preference of soil microorganisms.

Most studies have shown that the AP content in soil is an important driving factor affecting soil microbial diversity, which can improve soil microbial diversity (Tan et al., 2013; Siciliano et al., 2014). However, in our study, soil AP content had no effect on microbial diversity and community structures. Other researchers (Liu and Jiang, 2018) also obtained similar experimental conclusions. This may be related to the small range of AP changes in this study. It has been reported (Jian et al., 2021) that the overall soil comprehensive fertility of P. massoniana forests is poor, and the AP content is low. Long-term low phosphorus stress causes specific microbial communities to be screened out and adapt to this growth condition for a long time (Ky et al., 2020). Microorganisms in low phosphorus-stressed environments were less affected by this than in soil environments with higher AP levels. Therefore, AP was not a driver of microbial community change in the P. massoniana forest ecosystem.

This study confirmed that microbial diversity and community structure were significantly different in the rhizosphere of P. massoniana provenances under different carbon storage levels. For P. massoniana provenances, in the case of small differences in environmental factors, a very important reason for the variation of carbon storage among provenances is due to their own factors, which are mainly reflected in the wood properties of the standing trees and the strength of carbon storage related metabolic pathways (Hou et al., 2016). It has been shown (Yin et al., 2010) that newly fixed plant carbon is allocated to the lower part of the ground for root growth and is continuously lost in the form of rhizosphere sediment. Differences in the amount of carbon storage affect to some extent the type and metabolic rate of plant root metabolites which contains high molecular aromatic compounds and low molecular organic acids, resulting in altered soil sediments (Bi et al., 2021). This will further affect the soil physicochemical properties and microorganisms. Microorganisms are also influenced by the organic matter decomposition, apoplast, nutrient requirements and root turnover (Bi et al., 2022). Variations in organic matter and nutrient content can stimulate microorganisms to secrete more extracellular enzymes that accelerate soil nutrient activation (Luan et al., 2022). In turn, these enzymes can promote nutrient uptake by plant roots. This process can illustrate the relationship between plant–microbe interactions more graphically. A number of researchers have been working on carbon storage in forests from different fields and perspectives (Geddes et al., 2014). From the perspective of plant–microbial interaction, we selected P. massoniana provenances that establish better interaction with environmental microorganisms. While providing excellent varieties for the construction of high carbon storage artificial forests of P. massoniana, we reasonably analyzed the relationship between the carbon storage level of P. massoniana provenances and the environment, which has certain guiding significance for the scientific management of carbon storage forestry in the future. In this study, we believe that AY provenance representing high carbon storage levels have more abundant rhizosphere microbial diversity, which represents a more stable ecological structure of the rhizosphere microenvironment. Therefore, AY provenance is more suitable for building high carbon storage plantations for the purpose of sustainability. In contrast to a series of previous studies (Chen et al., 2022; Zhou et al., 2022) on P. massoniana seedlings, the P. massoniana plantation studied in this paper is an overmature forest. The response of rhizosphere microorganisms to carbon storage capacity and environmental factors at different growth stages is provided with more diversified reference significance and value. Under the new dynamics of carbon neutrality targets, future studies on the interactions between carbon storage and the effects of P. massoniana rhizosphere microenvironments under the coupling of multiple environmental factors (light, water, CO₂ concentration) should also be reinforced.