Influence of 37 Years of Nitrogen and Phosphorus Fertilization on Composition of Rhizosphere Arbuscular Mycorrhizal Fungi Communities in Black Soil of Northeast China

Increased inorganic nitrogen (N) and phosphorus (P) additions expected in the future will endanger the biodiversity and stability of agricultural ecosystems. In this context, a long-term fertilizer experiment (37 years) was set up in the black soil of northeast China. We examined interaction impacts of elevated fertilizer and host selection processes on arbuscular mycorrhizal fungi (AMF) communities in wheat rhizosphere soil using the Illumina MiSeq platform. The soil samples were subjected to five fertilization regimes: no fertilizer (CK) and low N (N1), low N plus low P (N1P1), high N (N2), and high N plus high P (N2P2) fertilizer. Long-term fertilization resulted in a significant shift in rhizosphere soil nutrient concentrations. The N fertilization (N1 and N2) did not significantly change rhizosphere AMF species diversity, but N plus P fertilization (N1P1 and N2P2) decreased it compared with CK. Non-metric multidimensional scaling showed that the rhizosphere AMF communities in CK, N1, N2, N1P1 and N2P2 treatments were distinct from each other. The AMF communities were predominantly composed of Glomeraceae, accounting for 30.0–39.1% of the sequences, and the relative abundance of family Glomeraceae was more abundance in fertilized soils, while family Paraglomeraceae were increased in N1 and N2 compared with CK. Analysis shown that AMF diversity was directly affected by soil C:P ratio but indirectly affected by plant under long-term fertilization. Overall, the results indicated that long-term N and P fertilization regimes changed rhizosphere AMF diversity and community composition, and rhizosphere AMF diversity was both affected by soil C:P ratio and plant.


INTRODUCTION
The interactions between microorganisms belowground and plants aboveground significantly influence ecosystem properties and processes. Arbuscular mycorrhizal fungi (AMF) are mutualistic fungi that form symbiotic relationships with the majority of land plants, including many crops. Within this symbiosis, AMF can provide plants with critical nutrients and water, as the large extra-radical mycelium exploits the soil and transports these nutrients to plants (Hodge and Fitter, 2010). In other cases, AMF can act as protectants against phytopathogens and enhance the sustainability of ecosystems by improving soil structure (Camenzind et al., 2014). The plant, in turn, provides AMF with an energy source, especially carbon (C).
Nutrient depositions by anthropogenic activities have increased in numerous agricultural lands. This elevated nutrient availability would no doubt influence AMF communities, as it could reduce the benefit provided by these symbionts. According to the functional equilibrium model, when soil nutrients become less limited, plants then allocate resources toward compositions of what could acquire other limited resources (Ericsson, 1995) resulting in a decrease in mycorrhizal structures and fine roots (Johnson, 2010). In addition, fungus itself can be nutrient-limited and sensitive to changes in soil properties (Zhou et al., 2016). Studies have showed that elevated concentrations of soil nitrogen (N) as a result of N fertilization decreased AMF biomass, species richness and diversity in some ecosystems (Egerton-Warburton et al., 2007;Verbruggen et al., 2012;Williams et al., 2017). Furthermore, changes to AMF communities caused by N application may be linked to soil P availability (Egerton-Warburton et al., 2007;Cheng et al., 2013). The rhizosphere is considered to be one of the most dynamic interfaces on Earth (Wang et al., 2018). Rhizosphere microbial communities are shaped by interactions between agricultural management and host selection processed (Schmidt et al., 2019), but those studies only focus on the soil or root AMF communities, not on rhizosphere soil. Wang et al. (2019) shown that the root-feeding AMF abundance was increased under elevated N-addition, although other studies found that N decreased its abundance . Those results indicated that the rhizosphere AMF may responsed differently to nutrients addition. However, our understanding of how fertilizer and host selection processes influence on this important mutualistic fungi is still limited, and understanding the shifts in rhizosphere AMF structure and composition following long-term fertilization may have significant implications for using the AMF to increase nutrient availability in soils.
Chinese black soil is distributed in a narrow and long area of ∼6 million ha, and 70% of this is in Heilongjiang Province of Northeast China. The relatively high organic matter and cation exchange capacity and its favorable macronutrient status, as well as a thick (60 cm) mollic epipedon, provides favorable soil structure and conditions for plant growth (Xing et al., 2005). The black soil region contributes to up to 30% of the national staple food using only 20% of the national arable land. Thus, it is important for China's food security. However, extensive agricultural intensification has input large amounts of inorganic fertilizers into this region, resulting in serious degradation of soil physicochemical properties and environmental health since the 1950s (Yin et al., 2015). Zhou et al. (2016) showed that longterm inorganic fertilizer use had significantly decreased soil pH and increased nutrient availability. This nutrient-based alteration significantly changed soil microorganism communities (Yin et al., 2015;Zhou et al., 2016). Our previous study showed that longterm N and P fertilization in black soil reduced fungal diversity and increased ITS gene copy numbers, thus altering the fungal community composition -such shifts were correlated with soil pH (Zhou et al., 2016). However, knowledge of the influences of long-term inorganic fertilizer application on AMF communities and the primary factors driving AMF species diversity in black soils of northeast China is still limited.
To test potential impacts of future nutrient depositions on AMF communities and diversity in black soil, a field experiment using different N and P fertilizer rates was begun in 1980 in Harbin city, northeast China (Wei et al., 2008). In the present study, we investigated the response of rhizosphere AMF species diversity and community composition to different fertilization strategies, and determined the relationships between AMF species diversity, dominant groups and soil parameters associated with these changes. It was hypothesized that longterm N and P fertilizer application would significantly change soil AMF community composition, and AMF species diversity would decrease with the incorporation of N and P fertilization, but increase with N fertilization alone.

Experimental Description
The long-term field experiment was established in 1980 at the Scientific Observation Station of Arable Land Conservation and Agriculture Environment of Heilongjiang Academy of Sciences (45 • 40 N, 126 • 35 E and altitude 151 m), where the mean annual temperature and average annual precipitation of 3.5 • C and 533 mm, respectively. The experimental field contains black soil, which is widespread in the Northeast China. The longterm fertilization experiment has a completely randomized block design with three replicated plots (9 m × 4 m). As previous studies shown that AMF composition and diversity was most affected by N and P fertilizers, but not K fertilizer (Lin et al., 2012;Dueñas et al., 2020), we selected the following treatments: no added fertilizer (CK) and low N (N 1 , 150 kg urea ha −1 y −1 ), low N plus low P (N 1 P 1 , 150 kg urea plus 75 kg P 2 O 5 ha −1 y −1 ), high N (N 2 , 300 kg urea ha −1 y −1 ) and high N plus high P (N 2 P 2 , 300 kg urea plus 150 kg P 2 O 5 ha −1 y −1 ) fertilizer. The N fertilizer was applied as in the form of urea, while the P fertilizer was ammonium hydrogen phosphate and calcium super phosphate. The rates of low N and low P fertilizers in the experimental field follow local customs. Fertilizer treatments were maintained in the same plot location each year. The cropping system was a wheat-maize-soybean rotation, which is one of the main cropping patterns in this area. More detail on this long-term experimental field is given in Wei et al. (2008). The sampling time was during wheat florescence on 29 June 2016.
For each replicated plot, roots were taken from 20 randomly selected plants, loosely adhering soil was shaken off and the tightly adhering soil carefully collected. These rhizosphere soils were pooled to form one composite sample. A total of 15 fresh samples (three replicates × five treatments) were transported to the laboratory on ice, and sieved through a 2-mm mesh to remove plant roots. Each sample was divided into two parts: one stored at room temperature for chemical analysis and the other stored at −80 • C for molecular analysis. The wheat grain yields were measured after harvest.

Chemical Analysis
Soil pH was measured with a glass combination electrode using a soil to water ratio of 1:1. Soil available P (AP) and available potassium (AK) were determined according to Olsen (1954) and Helmke and Sparks (1996), respectively. Nitrate (NO 3 − )-N and ammonium (NH 4 + )-N were extracted with 2M KCl, and determined with a flow injection autoanalyzer (FLA star 5000 Analyzer, Foss, Denmark). Total organic C (TOC) content was determined by wet digestion using a mixture of potassium dichromate and sulfuric acid under heating. Total N (TN) was measured according to Strickland and Sollins (1987).

Total DNA Extraction and Sequencing
We extracted total DNA from 0.25 g of soil using MoBio Power Soil DNA isolation kits according to the manufacturer's protocol (MOBIO Laboratories Inc., Carlsbad, CA, United States). To minimize the DNA extraction bias, we combined six successive DNA lots extracted from the same soil sample and purified using a DNeasy Tissue kit (Qiagen, Valencia, CA, United States).
Amplicons of AMF libraries were produced from each of the 15 soil extracts by nested PCR with first-round PCR primers of LR1 (5 -GCATATCAATAAGCGGAGGA-3 ) and FLR2 (5 -GTCGTTTAAAGCCATTACGTC-3 ) (Van Tuinen et al., 1998), and second-round PCR primers FLR3 (5 -TTGAAAGGGAAACGATTGAAGT-3 ) and (5 -TACGTCAACATCCTTAACGAA-3 ) (Gollotte et al., 2004). The second-round primers were tagged with sequencing adapters followed by an 8-mer multiplexing identifier. The two rounds of PCR (50 µl) both contained 5 µl of 10× Pyrobest Buffer, 4 µl of dNTPs (2.5 mM), 2 µl of each primer (10 µM), 0.75 U of Pyrobest DNA Polymerase and 30 ng of template DNA. The PCR amplification procedure for both rounds was 5 min at 95 • C, followed by 30 cycles of 45 s at 95 • C (denaturation), 50 s at 58 • C (annealing) and 45 s at 72 • C (extension), with a final extension step of 10 min at 72 • C. PCR products from each sample were purified, pooled together in equimolar ratios and sequenced using the Illumina MiSeq platform. Raw sequence data for all the samples were uploaded to the NCBI Sequence Read Archive under accession number SRX3008208.

Bioinformatic Analysis of Sequence Data
Sequence read analysis was carried out using Mothur 1.33 (Schloss et al., 2009). The primers and multiplexing identifier were trimmed, and assignment of samples was based on unique barcodes. Reads with a quality score < 20, with ambiguous nucleotides, lacking a complete barcode or <200 bp were removed and excluded from further analysis. This was followed by checking for chimeras and removal of predicted chimeras. The remaining high-quality sequences were then clustered into operational taxonomic units (OTUs) at 97% identity threshold using UPARSE (Edgar, 2013). The longest sequence from each OTU was selected as the representative sequence, and a manual BLASTing against the GenBank non-redundant nucleotide database was used to detect non-Glomeromycota sequences. The non-Glomeromycota sequences and OTU singletons were removed (Camenzind et al., 2014;Williams et al., 2017). The Chao1 index was used to determine the AMF richness, and the phylogenetic diversity index was used to determine the phylogenetic diversity. The AMF richness and phylogenetic diversity were calculated after subsampling according to the sample with the least sequences in Mothur (Schloss et al., 2009).

Statistical Analysis
Significant differences in soil properties, alpha-diversity and AMF community abundance among samples were determined with a one-way ANOVA and least significance difference (LSD) using SPSS (version 19.1) statistical software (SPSS, Chicago, IL, United States). Pearson's correlation coefficients were used to assess relationships among soil properties, alphadiversity and abundant AMF genera. In all tests, P < 0.05 was considered to be statistically significant. Non-metric multidimensional scaling (NMDS) was performed to examine the effects of fertilization on the AMF community using CANOCO 5.0. Redundancy analysis (RDA) was carried out using CANOCO 5.0 to determine correlations between environmental variables (pH, AP, AK, TOC, and soil C:P ratio) and AMF community composition with the Monte Carlo permutation test (999 permutations).
Structural equation modeling (SEM) was applied to gain a mechanistic understanding of how soil properties and plant mediate alterations in rhizosphere AMF diversity and composition under different fertilizer regimes. The community composition of AMF was obtained by NMDS, and the NMDS1 were used in the subsequent SEM analysis. SEM analysis was performed with the specification of conceptual model of hypothetical relationships (Supplementary Figure S1), assuming that long-term fertilization alters soil properties and plant biomass, which in turn affects AMF diversity, composition and plant. The maximum likelihood estimation method was used to test the data were fitted to the models, and adequate model fits were indicated by the low χ 2 /df (<2), non-significant chisquare test (P > 0.05), a low RMSEA (RESEA < 0.05), and high goodness-of-fit indes (GFI > 0.9) (Zeng et al., 2016;Ning et al., 2020).

Long-Term Fertilization Changed Soil Parameters in Wheat Rhizosphere
After a 37-year application of N and P fertilizer, soil properties in the rhizosphere were significantly altered (P < 0.05, Table 1). Soil contents of TN and NO 3 − -N was not significantly Values are mean ± standard deviation (N = 3). Values within the same column followed by the different letters indicate significant difference (P < 0.05). Fertilizer regimes: CK (without fertilizer), N 1 (150 kg N ha −1 ), N 1 P 1 (150 kg N ha −1 plus 75 kg P ha −1 ), N 2 (300 kg N ha −1 ), and N 2 P 2 (300 kg N ha −1 plus 150 kg P ha −1 ). Soil properties indicated include AP (available P), AK (available K), TOC (total organic carbon), TN (total N), C:N (TOC:TN), C:P (TOC:AP), and WY (wheat yield).
changed by N 1 and N 1 P 1 , but significantly increased by N 2 and N 2 P 2 . However, 37 years of fertilization significantly decreased AK, and soil pH steadily decreased from 6.81 to 5.45 with increasing fertilizer inputs. Soil AP was not significantly changed by N fertilization (N 1 and N 2 ), but significantly increased by N plus P addition. The soil C:P ratio in the N 1 treatment was 43.2% higher than CK, and correspondingly for the N 1 P 1 , N 2 , and N 2 P 2 treatments was 32.8, 79.1, and 92.1% lower. Furthermore, wheat yield was significantly increased by long-term fertilization, but there were no significant differences between N 1 and N 2 or between N 1 P 1 and N 2 P 2 treatments.

Long-Term Fertilization Changed Rhizosphere AMF Community Composition and Diversity
The nested PCR successfully amplified DNA from all 15 samples. After quality control, a total of 409,049 high-quality sequences were obtained. They could be classified with a mean of 27,270 classifiable sequences per sample for use in the subsequent analysis (range 17,625-30,537). All the high-quality sequences were clustered into 220 OTUs; and 6551 sequences (clustered into 60 OTUs) represented non-AMF but this was only 1.6% of all sequences. The Good's coverage values were all above 99.8% with a 97% similarity cutoff, indicating that the current depth of sequencing was sufficient to capture the AMF diversity. The numbers of OTUs in the five fertilizer treatments were in the range of 66-92 (Figure 1). In our study, 99.7% of all sequences matched Glomeromycota, and a total of 7 families were found in this study site. The four families (average abundance > 0.1%) Glomeraceae, Claroideoglomeraceae, Gigasporaceae, and Paraglomeraceae were represented in the sequencing dataset (Figure 2), the remaining three families Diversisporaceae, Archaeosporaceae, and Acaulosporaceae were low abundance in all rhizosphere soils. All soil samples were dominated by Glomeraceae, accounting for 30.0-39.1% of the sequences, followed by Claroideoglomeraceae (2.4-21.9%) and Gigasporaceae (0.03-2.8%). Family Paraglomeraceae was only occasionally detected at a low level. Moreover, 35-54% of sequences did not match any detailed taxonomic at family level. Long-term fertilization significantly (P < 0.05) changed the abundance of families Glomeraceae, Claroideoglomeraceae, Gigasporaceae and Paraglomeraceae. The relative abundance of Glomeraceae in the N 1 , N 1 P 1 , N 2 and N 2 P 2 treatments was 9.2, 7.6, 7.2, and 6.9% higher, respectively, than that of CK. The relative abundance of Paraglomeraceae, which was only detected in N 1 and N 2 treatments, showed a positive relationship with N dosage (Figure 2). At genus level, the abundance of AMF distinctly responded to long-term N and P fertilization. The nine most abundant genera with significant differences are presented in Figures 3, 4. At genus level, Glomus, Claroideoglomus, and Rhizophagus were the dominant genera in wheat rhizosphere soil, occupying 18.26, 15.15, and 14.12% of the data set, respectively. The relative abundance of Glomus was increased by N 1 and N 1 P 1 , but not significantly changed by N 2 and N 2 P 2 , compared with CK. However, the relative abundance of Rhizophagus and Septoglomus was significantly increased and decreased by higher N addition (N 2 and N 2 P 2 ), respectively, but was not significantly changed by lower N addition (N 1 and N 1 P 1 ). Long-term fertilization decreased the relative abundance of Funneliformis and Archaeospora. For N 1 P 1 and N 2 P 2 treatments, the relative abundance of Paraglomus and Diversispora decreased to a greater extent than in the N 1 and N 2 treatments. AMF alpha diversity, including OTUs, AMF richness and phylogenetic diversity, was significant altered in the rhizosphere (Figure 1). The OTUs, phylogenetic diversity and AMF richness were lower in N 1 P 1 and N 2 P 2 treatments compared to CK, but there were no significant differences among CK, N 1 and N 2 treatments. AMF richness in N 1 P 1 and N 2 P 2 treatments was decreased by 17.3 and 20.2%, respectively, compared with CK. Beta-diversity analysis (NMDS) using Bray-Curtis distance suggested that the samples were well separated from each other. NMDS1 generally distributed the AMF communities along with soil pH: the AMF communities with acidic pH were generally in the left portion and those with higher pH to the right (Table 1 and Figure 5A). NMDS2 generally distributed the AMF communities along with soil AP: the AMF communities with high AP concentration were generally in FIGURE 1 | Alpha diversity of detected OTUs, AMF richness and phylogenetic diversity across soils from different fertilizer regimes. Vertical bars represent the standard deviations (N = 3) and the same letters above columns denote no significant difference (P < 0.05, Tukey's test). Fertilizer regimes indicated as described in Table 1. the upper portion and those with low AP concentration in the lower portion.

Relationships Between AMF Community and Environmental Variables
Nine variables were selected by RDA as predictors of species persistence ( Figure 5B). All variables together explained 77.38% of the variation in AMF communities between samples. The first axis explained 37.17% of the data variation, and the second axis explained 24.94%. In Figure 5B, the lengths of the red arrows indicated the relative importance of each variable in explaining the rhizosphere AMF community composition, and soil C:P ratio, pH and AP were the three longest arrows. Indeed, based on this model, soil C:P ratio (F = 6.8, P = 0.002), pH (F = 5.5, P = 0.002), and AP (F = 3.6, P = 0.002) were the three most important contributors to the variation in AMF communities, and individually accounted for 34.2, 20.7, and 11.2% of variation, according to the Monte Carlo test, respectively. The first axis was associated with variables related to soil AP concentration and soil C:P ratio as reflected by the complete set of data. Soil C:P ratio was positive correlated with the rhizosphere AMF composition in CK and N 1 treatments, whereas negative correlated with communities in N 1 P 1 and N 2 P 2 . The oppose trend was found with soil AP. The soil pH changed by long-term fertilization was positive correlated with AMF composition in CK, whereas negative correlated with AMF composition in N 2 and N 2 P 2 .
Pearson's correlation analysis with AMF community composition in genera level and AMF diversity level confirms this relationship with soil C:P ratio, pH and AP. In diversity level, the OTUs, AMF richness and phylogenetic diversity indices were all significantly positively correlated with soil C:P ratio (P < 0.01), but negatively correlated with AP. However, soil pH was only significantly negatively correlated with OTUs and AMF richness (Supplementary Table S1).
In genera level, the soil C:P ratio was significantly positively correlated with relative abundance of Septoglomus and Archaeospora; soil pH was significantly positively correlated with the relative abundance of Funneliformis, Septoglomus, and Archaeospora, but negatively correlated with Rhizophagus; and AP concentration was positively correlated with Rhizophagus, but negatively correlated with Septoglomus and Archaeospora. In addition, soil C:N was also significantly correlated with Glomus, Claroideoglomus, Rhizophagus, and Septoglomus.

Integrated Responses of AMF Structure on Soil Properties and Plant
We also assessed the responses of AMF structure on soil properties and plant by employing a SEM model. The model proved a good fit to the data (χ 2 /df = 0.888; P = 0.447; RESEA = 0.000; GFI = 0.953), and accounted for 96% of the variation in pH, 68% and 83% in NO 3 − -N and C: P ratio, 80% in plant biomass, 87 and 90% in AMF composition and diversity, respectively (Figure 6). Fertilization significantly affected AMF diversity and composition due to decreased soil pH and C: P ratio and changed plant biomass (Figure 6). Decreased in soil C:P ratio and pH directly altered AMF diversity and AMF composition, respectively. In addition, alterations in plant biomass resulted in concomitant shifts in AMF diversity. Those results indicated that rhizosphere AMF structure was directly altered by soil property and indirectly shaped by changing plant biomass. The relationships between those variables were present in Supplementary Table S2, and some variables were not of great significance on their own, but they clearly improved the model fit when incorporated together. FIGURE 2 | Relative abundance of arbuscular mycorrhizal fungal families (average abundance > 0.1%) for different long-term fertilizer treatments. Vertical bars represent the standard deviations (N = 3) and the same letters above columns denote no significant difference (P < 0.05, Tukey's test). Fertilizer regimes indicated as described in Table 1. (I) Diversisppora. Vertical bars represent the standard deviations (N = 3) and the same letters above columns denote no significant difference (P < 0.05, Tukey's test). **The difference is significant at the 0.01 level; * the difference is significant at the 0.05 level according to Tukey's test. Fertilizer regimes indicated as described in Table 1.

DISCUSSION
In this study, we explored the AMF fungal communities in wheat rhizosphere soils. The analysis showed that longterm fertilization had great impacts on rhizosphere AMF species diversity and composition, and the rhizosphere AMF community structure was shaped both by fertilization and plant. The field resource in which the same fertilization regimes were conducted for 37 years enabled us to show that N and P addition produced a strong selective FIGURE 4 | Heat map analysis the relative abundance of rhizosphere AMF genera after long term fertilization. Fertilizer regimes indicated as described in Table 1. pressure against rhizosphere AMF over the long term in this study site.

Long-Term Nutrient Addition Effects on Rhizosphere AMF Species Alpha-Diversity
In our study, N fertilizer did not significantly changed AMF diversity in the wheat rhizosphere, consistent with previous reports (Mueller and Bohannan, 2015;Dueñas et al., 2020). This may reflect the influence of interaction between fertilization and plant on rhizosphere AMF, as the relocation of C to AMF increased with N fertilization. For example, after a large N addition, relatively more plant derived 13 C was allocation to AMF (Williams et al., 2017). In addition, Wang et al. (2019) found that long-term N fertilization increase the quantity of plant root exudates, which also modified the soil properties changes caused by N addition (Wang et al., 2018). Those results indicated that the stability of the AMF diversity after long-term N fertilization was due to both elevated fertilizer and host selection processes.
However, the N 1 P 1 and N 2 P 2 treatments had low AMF diversity in wheat rhizosphere soils than CK, consistent with a study of two grass species (Egerton-Warburton et al., 2007). In terms of OTUs detected and community composition,  changes appeared due to the distinct species induced by those treatments, especially the taxa with lowest abundance (Figures 1-3). The loss of AMF species diversity has been reported to decrease both plant biodiversity and ecosystem productivity, and to increase ecosystem instability (Maček et al., 2011). The present study suggested AMF diversity decrease may result in a less stable agroecosystem and may lead to unsustainable crop production. Egerton-Warburton et al. (2007) and Cheng et al. (2013) already found that changes to AMF diversity caused by N application may be linked to soil P availability; N fertilizer application cause reductions in P-rich soils. This indicated that N plus P addition may have a stronger influence on the rhizosphere AMF species diversity compared with only N addition (Figure 1). Although N fertilization increased plant allocate more C to AMF, N plus P fertilization decreased plant C allocation to AMF (Williams et al., 2017). Considering the most widely cited benefit of AMF is to enhance plant P acquisition (Hodge and Fitter, 2010), and more soil plant-available P reduce plant acquired P from AMF (Williams et al., 2017) and shift AMF functioning from mutualism toward parasitism , some members of the rhizosphere AMF have greater carbon demands or have a less efficient use of carbon may be extinguished. This may explain why N plus P addition decreased rhizosphere AMF diversity compared with CK. In addition, we also found that there was no significance between N 2 and N 2 P 2 treatment on OTUs detected and Phylogenetic diversity (Figure 1). This may be attributing to high N application significantly decreased soil pH and not significantly increased plant biomass compared with low N application (Table 1), as rhizosphere AMF diversity was both shaped by soil pH and plant biomass (Figure 6). After a long time (such as 37 years) of sufficient P available and dramatic changes soil properties like soil pH (caused by long-term N fertilization), the AMF species that cannot acquire enough energy may be extinguished.
Overall, our results suggested that the responses of rhizosphere AMF to elevate fertilizer and plant changes with soil nutrients. When N fertilization increases N supply (P supply is the most limiting resource), the benefit of the symbiosis is enhanced. In this scenario, rhizosphere AMF diversity levels might be maintained. When N and P fertilization increases N and P supply, the benefit provided by this symbiosis is reduced (Dueñas et al., 2020). This may intensify competition among AMF taxa for plant derived C. In this case, a reduction in rhizosphere AMF diversity can be expected.

Long-Term Nutrient Addition Changed Rhizosphere AMF Community Composition
The NMDS analysis demonstrated that AMF communities differed significantly among different fertilization regimes, suggesting that different N and P additions induced a shift in AMF community composition. These results agree with former findings that the abundance of soil nutrients plays a crucial role in shaping the rhizosphere AMF community (Qin et al., 2015;Dueñas et al., 2020). It was reported that nutrients accumulation caused markedly changes in soil chemical properties, species that are insensitive to the disturbances or efficient in acquiring plant C will increase in relative abundance (Koch et al., 2017;Wang et al., 2018). In addition, competition can also lead to the complete exclusion of particular species (Roger et al., 2013) or a decline in overall AMF abundance (Engelmoer et al., 2014). The addition of N plus P fertilization may select for taxa with better ability to hoard P in order to maximize C gains from the host (Whiteside et al., 2019).
In our study, we only found 7 AMF families, similar with the results of Zhu et al. (2020), which found 7 ∼ 8 families in different site of Chinese black soil region. However, Redecker et al. (2013) systematic reviewed of AMF taxa, and summarized 11 AMF families from literatures. The other four families were not detected including Sacculosporaceae, Pacisporaceae, Ambisporaceae, and Geosiphonaceae in our study site. This may be due to those four families were in low abundance, and cannot be detected by Illumina MiSeq platform. For example, Gu et al. (2020) found only less than 0.4% of all sequences belonged to Sacculosporaceae and Pacisporaceae in Chinese black soil. In addition, family Pacisporaceae showed a preferential distribution toward poor nutrient soils, rather than nutrient-rich black soil (Wang et al., 2018;Dudinszky et al., 2019). Cao et al. (2020) and Panneerselvam et al. (2020) found that Ambisporaceae and Geosiphonaceae were more abundant in higher temperature (mean annual temperature 19.1 • C) and elevated CO 2 concentration (550 µmol mol −1 ), respectively. The low temperature, atmospheric CO 2 concentration and high nutrient level may explain why we did not detect those four families in our study site.
In the present study, we observed that long-term N and P fertilization increased the relative abundance of Glomeraceae compared with CK (Figure 2), with the relative Glomeraceae abundance was 39.1%, among them the relative genus Glomus abundance increased by 89.3% compared with CK ( Figure 3A). Previous studies found that the family Glomeraceae could grant better protection from pathogens (Powell et al., 2009), as the majority biomass in the Glomeraceas is found in hyphae growing inside the root and reduce root infection of plant by two soil pathogens (Maherali and Klironomos, 2007). Those results indicated that N and P fertilization shifted the AMF communities toward disease-suppression against root pathogens. Our results also demonstrated that the relative abundance of Paraglomeraceae was higher in N 1 and N 2 treatments than other treatments (Figure 2), similar to observations in some previous studies (Hijri et al., 2006;Chagnon et al., 2013). The Paraglomeraceae with the majority of fungal biomass are located outside the plant were demonstrated to be more efficient in P assimilation (Powell et al., 2009). Those results indicated that long-term N (only) addition changed AMF communities by increasing nutrient assimilation taxa. The results also indicated that long-term N fertilization only may enhance plant P acquisition via AMF communities. In addition, Jiang et al. (2018) found that increased nutrients available would shift mycorrhizal functioning toward parasitism, and that the inefficient in acquiring soil derived C under elevated plant-available N and P levels in soil may be extinguished. However, the links between AMF traits and nutrient requirements or function inferred from a fraction of AMF isolates, and this may prevent us to unequivocally establish whether differential adaptations to nutrient supply are the basis for the patterns reported here.

Factors Affecting Rhizosphere AMF Community Composition
The interactions among microbes, plants and soils play a critical role in ecosystem functioning (Zeng et al., 2016). However, earlier studies about the effects of fertilization on soil AMF communities did not include plant production data, often neglecting the interaction between elevated fertilizer and host selection processes on AMF communities. In our study, we found that N and P fertilization resulted in significant changes in rhizosphere AMF community composition. We further found that rhizosphere AMF composition was directly affected by decreasing soil pH caused by N and P fertilization, while AMF diversity was directly affected by soil C: P ratio and indirectly affected by plant biomass (Figure 6).
In our study, the composition of rhizosphere AMF communities was closely correlated with soil C:P (P < 0.05), as reported by Qin et al. (2015). This could be highly associated with AMF receiving their C supply from their host plants and compensating the plant through enhanced nutrient acquisition, particularly through supply of poorly mobile phosphate ions (Karasawa et al., 2012). In the case of greater nutrient sufficiency in plants in fertilizer treatments, symbiosis with AMF could be less important, reducing C allocation to AMF in the rhizosphere. This could explain why soil C:P was an important factor in shaping AMF communities, and indicated that rhizosphere AMF community was shaped by both soil nutrient and plant selection processes.
The shifts in the relative abundance of species taxonomic groups across different soil C:P values are similar to soil C:P responses observed in other studies. For example, the relative abundance of Gigaspora and Archaeospora had strong positive correlations with high soil C:P in this study ( Table 2), consistent with results of Pagano and Scotti (2010). The relative abundances of Gigaspora and Archaeospora were higher in N 1 and N 2 than N 1 P 1 and N 2 P 2 treatments and this was significantly negatively correlated with AP concentration (P < 0.05). Previous studies demonstrated that P addition could shift species composition in favor of less efficient mutualists, which are competitively superior colonizers (Johnson, 1993). This indicated that the genera Gigaspora and Archaeospora were efficient mutualists with wheat in P-deficient soil in this site.
The soil pH was also an important factor affecting the relative abundance of species taxonomic groups, which was consistent with other reports (Jansa et al., 2014;Qin et al., 2015). In the present study, RDA also indicated that the AMF composition was affected by soil pH (Figure 6). The relative abundance of Funneliformis, Septoglomus, and Archaeospora were positively correlated with soil pH, but Rhizophagus was negatively correlated with it, as similarly reported by Jansa et al. (2014). With the increase in the amount of fertilizer, soil pH decreased from 6.81 to 5.45 (Table 1), and abundance of these AMF taxa changed gradually. The results indicated that these AMF groups were sensitive to soil pH. Rhizophagus was relatively more abundant in N 2 and N 2 P 2 than other treatments (Figures 3, 4), and positively correlated with most available nutrients ( Table 2). This is may be because Rhizophagus can flourish in soils with large amounts of nutrients, and exert a P-uptake function with a large requirement of energy source (Thonar et al., 2011). The available nutrients were higher in N 2 and N 2 P 2 than other fertilizer treatments and CK ( Table 1). The relative abundances of Glomus and Claroideoglomus were significantly positively correlated with soil C:N values, consistent with other reports that soil C:N was also an important factor affecting the soil AMF community (Shi et al., 2014;Treseder et al., 2018;Dueñas et al., 2020).

CONCLUSION
We demonstrated that rhizosphere AMF species diversity was decreased by N plus P fertilization, but not significantly affected by fertilization with N only. Composition of the rhizosphere AMF community was significantly influenced by N and P fertilizers, which appeared to be mediated greatly by soil-plant-microbe interactions. Specifically, the rhizosphere AMF diversity was directly affected by soil C: P ratio and indirectly affected by plant biomass, while rhizosphere AMF community was directly shaped by soil acidification. Our results suggest that rhizosphere AMF community structure was shaped by elevated fertilizer and host selection processes. The 37year inorganic fertilization regimes changed rhizosphere AMF communities with a potential negative impact on AMF transport of nutrients to plants and on the beneficial effect of AMF genera to plant resistance against pathogens, because longterm inorganic fertilization promotes microbes with known pathogenic traits (Zhou et al., 2016). More studies need to be conducted to elucidate the mechanisms that AMF taxa and communities use to cope with the pathogens and nutrients transportation caused by global environmental changes, such as N-or P-deposition. Most importantly, studies about how soilplant-microbe interactions influence on ecosystem functions and progresses remains to be elucidated.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI Sequence Read Archive, SRX3008208.

AUTHOR CONTRIBUTIONS
QW, SW, and JL conceptualized and designed the experiments. MM, XJ, DG, FC, YK, and CC analyzed the data. DW provided the material for this experiment. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank the University of Liège-Gembloux Agro-Bio Tech and more specifically the research platform AgricultureIsLife for the funding of the scientific stay in Belgium that made this manuscript possible.