Chronic fertilization of 37 years alters the phylogenetic structure of soil arbuscular mycorrhizal fungi in Chinese Mollisols

Arbuscular mycorrhizal fungi (AMF) play vital roles in sustaining soil productivity and plant communities. However, adaption and differentiation of AMF in response to commonly used fertilization remain poorly understood. In this study, we showed that the AMF community composition was primarily driven by soil physiochemical changes associated with chronic inorganic and organic fertilization of 37 years in Mollisols. High-throughput sequencing indicated that inorganic fertilizer negatively affected AMF diversity and richness, implying a reduction of mutualism in plant–AMF symbiosis; however, a reverse trend was observed for the application of inorganic fertilizer combined with manure. With regards to AMF community composition, order Glomerales was dominant, but varied significantly among different fertilization treatments. All fertilization treatments decreased family Glomeraceae and genus Funneliformis, while Rhizophagus abundance increased. Plant-growth-promoting-microorganisms of family Claroideoglomeraceae and genus Claroideoglomus were stimulated by manure application, and likely benefited pathogen suppression and phosphorus (P) acquisition. Family Gigasporaceae and genus Gigaspora were negatively correlated with available P in soil. Additionally, redundancy analysis further suggested that soil available P, organic matter and pH were the most important factors in shaping AMF community composition. These results provide strong evidence for niche differentiation of phylogenetically distinct AMF populations under different fertilization regimes. Manure likely contributes to restoration and maintenance of plant–AMF symbiosis, and the balanced fertilization would favor the growth of beneficial AMF communities as one optimized management in support of sustainable agriculture in Mollisols.


Introduction
Arbuscular mycorrhizal fungi (AMF) are ubiquitous in terrestrial ecosystems. They are obligate plant root symbionts found in more than 80% of plant families (Mueller and Bohannan 2015) and they provide multiple benefits for both plant and AMF. For instance, AMF can provide plants with critical nutrients ) and nonnutrient benefits (Walder and van der Heijden 2015), help plants withstand drought , develop soil aggregates (Rillig and Mummey 2006), stimulate photosynthesis (Kaschuk et al. 2009) and increase plant productivity (Maček et al. 2011), thereby influencing ecosystem processes and agricultural sustainability (Rillig 2004). Moreover, AMF confer resistance to plants against many soil-borne pathogens and various nematodes by their colonization of the plant root (Harrier and Watson 2004;Sikes et al. 2009). The plant, in turn, provides photosynthetic carbon (C) as an energy source for AMF, which can alter the fungal community (Kim et al. 2015). In addition to beneficial effects on plants, AMF also enhance ecosystem sustainability by influencing numerous soil properties and structure (Wilson et al. 2009), including soil stability (Mueller and Bohannan 2015), C storage (Treseder and Allen 2000), soil moisture , and nitrogen (N), C and phosphorus (P) cycles (Van Der Heijden et al. 2008). Lower AMF biodiversity can also lead to unsustainable crop production and ecosystem instability (Maček et al. 2011). More importantly, the beneficial effects of AMF decrease under conditions of high fertility (Collins and Foster 2009). Although they are ubiquitous and of ecological importance, AMF have been poorly studied in the important black soils of China.
Black soils, referred to as Mollisols, are one of the most crucial soil resources for crop production (Liu et al. 2012a), due to their key roles in the supply of national staple food (Han et al. 2013). However, excessive cultivation and intensive fertilization have caused serious soil degeneration and substantial losses of soil productivity since the middle of the last century (Singh et al. 2014;Liu et al. 2015a). Intensive inorganic input leads to a shift from fungal-based to more bacterial-based soil food webs (Thiele-Bruhn et al. 2012). Moreover, repeated application of inorganic fertilizer decreases soil pH (Liu et al. 2012b), which in turn can reduce nutrient availability and microbial biomass (Bardgett 2005). Application of organic manure can reduce use of inorganic fertilizers and benefit soil quality by alleviating soil acidification, accumulating soil organic matter (OM), improving soil microbial community structure and reinforcing the selfregulating status of soil systems (Hammesfahr et al. 2011;Insam et al. 2015). The responses of soil bacterial and fungal communities to inorganic fertilizer and manure in black soils were well reported in previous studies (Zhou et al. 2015Ding et al. 2016Ding et al. , 2017, but we know less about the effects on AMF, which work as good indicators of the soil ecosystem and affect the soil nutrition cycle and plant health. With chronic fertilizer input, the changes in soil properties no doubt influence AMF microbial communities and shift the plant-AMF symbiosis towards reduced mutualism or parasitism along a continuum (Verbruggen and Kiers 2010), as they could reduce the benefits provided by these symbionts. Understanding the shift of AMF community composition is critical to assessing the impact of fertilization for determining a more effective strategy in sustainable soil fertility. As previously documented, several mycorrhizal fungi are closely associated with the soil P cycle (Smith et al. 2011), and provide the plant hosts with P in the typically available form of phosphate. Given the relationship between AMF and soil available P (AP), and that soil AP, OM and pH are important contributors to the fungal community (Ding et al. 2017), they may also shape the AMF community. In this study, we investigated the chronic (37-year) effects of fertilization regimes on soil properties and the AMF community, and determined the primary factors driving the shifts in the AMF community. We hypothesized the following: (1) chronic inorganic fertilization would significantly decrease AMF diversity; (2) the shifts in AMF community composition would result from fertilizer-driven changes in soil properties, such as AP, pH and OM; and (3) the application of manure would shift the AMF community towards a better structure and diversity, with an opposite effect for inorganic fertilizer.

Field experiment and soil collection
Long-term experiments can provide realistic and effective conditions for obtaining the exact information required to maintain soil quality by determining changes in soil properties and processes (Li et al. 2008). Our field experiment was conducted at the Scientific Observation Station of Arable Land Conservation and Agriculture Environment, Harbin City, Heilongjiang Province, China. Wheat, soybean and maize had been continuously grown on the field since 1980. The soil at the site was classified a clay loam and the experiment had a completely randomized block design with three replications. Blocks were randomized into plots of 9 × 4 m. More information was provided in our previous studies . We analyzed soil samples from four treatments: CK (no fertilizer), M (manure), NPK (inorganic N, P and potassium (K) fertilizers) and MNPK (inorganic N, P and K fertilizers plus manure). Inorganic N, P and K fertilizers were applied as urea (75 kg hm −2 year −1 ), calcium superphosphate plus ammonium hydrogen phosphate (150 kg hm −2 year −1 ) and potassium sulfate (75 kg hm −2 year −1 ), respectively. The manure was horse manure applied at 18,600 kg hm −2 year −1 .
In June 2016, after 37 years of fertilization application, we collected the rhizosphere soils during florescence of wheat, when the rhizosphere effects were considered to be most pronounced (Cheng et al. 2003). Twenty wheat plants and their roots were collected from the middle of each replication and shaken gently. The remaining adhering soil was carefully collected and mixed thoroughly as a single rhizosphere sample. Thus, a total of 12 soil samples were obtained for further analyses. Each soil sample was passed through 2.0-mm screen and divided into two parts: one was stored at − 80 °C for molecular analysis, and the other used for determining soil chemical properties after air drying.

Soil chemical properties and wheat yields
Soil pH was measured with a pH meter using a 1:1 sample:water extract, after being shaken for 2 h. Soil OM was determined according to the K 2 Cr 2 O 7 -capacitance method (Ciavatta et al. 1991). Total N (TN) was measured using the Kjeldahl method (Strickland and Sollins 1987). Nitrate N (NO 3 − -N) and ammonium N (NH 4 + -N) were quantified by extraction with 2 M KCl, steam distillation and titration (Swift et al. 1996). A modified resin-extraction method was used for AP analysis (Hedley and Stewart 1982). Available K (AK) was analyzed by flame photometry (Shen et al. 2013). Wheat yields were recorded by the Heilongjiang Academy of Agricultural Sciences, after harvest in September 2016.

DNA extraction, PCR and Illumina MiSeq sequencing
A total of 12 soil samples were analyzed. A PowerSoil DNA extraction kit (MOBIO Laboratories Inc., Carlsbad, CA, USA) was used for DNA extraction from 0.25 g of soil of each sample according to the manufacturer's method. Five successive extractions from each soil sample were combined together to obtain enough DNA and to minimize the DNA extraction bias. Purification was performed according to Zhou et al. (2016).
The AMF DNA was selectively amplified from each soil extractions by nested PCR ). The first-round PCR with primers LR1 (5′-GCA TAT CAA TAA GCG GAG GA-3′) and FLR2 (5′-GTC GTT TAA AGC CAT TAC GTC-3′) was carried out in final 50-μL reaction mixtures containing 5 μL of 10 × PCR 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 (Van Tuinen et al. 1998). The PCR protocol used was an initial denaturation at 95 °C for 5 min, 30 cycles of denaturation at 95 °C for 45 s, annealing at 58 °C for 50 s, extension at 72 °C for 45 s and a final extension at 72 °C for 10 min. The primers FLR3 (5′-TTG AAA GGG AAA CGA TTG AAG T-3′) and FLR4 (5′-TAC GTC AAC ATC CTT AAC GAA-3′) were used for the second-round PCR (Gollotte et al. 2004). The FLR3/FLR4 set contained a specific sequencing adapter followed by a 8 mer multiplexing identifier. The reaction mixtures and cycling conditions were the same as for the first round, with modifications of the template using the product of the first amplification. After purification, nested PCR products were quantified using QuantiFluor ™ -ST (Promega, USA) , and then normalized in equimolar ratios and sequenced using the Illumina MiSeq platform. Raw sequencing reads were deposited in the NCBI Sequence Read Archive with accession number SRX3209560.

Bioinformatics and statistical analysis
We used Mothur 1.32 (http://www.mothur.org/) to process the sequencing data as described by Schloss et al. (2009). Barcode sequences were used to assign raw data by samples and low quality sequences were trimmed and excluded to obtain valid tags. Then, operational taxonomic units (OTUs) were defined by clustering at the 97% similarity level using UPARSE (Edgar 2013). Taxonomy was assigned to fungal OTUs according to the NCBI nucleotide sequences database. According to Camenzind et al. (2014) and Williams et al. (2017), the sequences out of Glomeromycota, as well as OTU singletons, were detected by manual BLASTing against the GenBank nonredundant nucleotide database and then removed. As a result, the sequences belonged to Glomeromycota phylum were obtained. The alpha diversity indices of AMF were calculated using Mothur (Schloss et al. 2009). A principal coordinate analysis (PCoA) on the basis of the weighted Fast UniFrac metric was carried out to examine the effect of inorganic and organic fertilizers on AMF community composition and to compare between-sample variations (Marsh et al. 2013). Redundancy analysis (RDA) was performed to visualize the effect of edaphic factors on AMF community structure using CANOCO 5.0 software. SPSS (V19, Chicago, IL, USA) was used for statistical analysis. The differences in soil chemical properties, alpha diversity indices (Chao1, PD whole tree, Shannon) and AMF community abundances among samples were analyzed using one-way analysis of variance. Tukey's procedure was followed for paired comparisons of different treatments. In all tests, P < 0.05 was considered significant.

Soil chemical properties in wheat rhizosphere and wheat yields
Soil chemical properties in the wheat rhizosphere respond differently to fertilization regimes (Table 1). The AP concentrations were higher under the three fertilization treatments than that of CK, with the highest for NPK and MNPK treatments. The NPK treatment significantly increased AK and TN concentrations compared with CK. The NO 3 − -N concentration was significantly higher under the MNPK compared with the NPK treatment. With regards to soil pH, no significant differences were detected between CK and M treatment. Although the application of inorganic N, P and K fertilizer in both NPK and MNPK treatments significantly decreased soil pH, manure significantly alleviated soil acidification. The application of manure also significantly increased soil OM compared with CK. In addition, wheat yields were significantly higher under the three fertilization treatments compared with CK, with NPK and MNPK the most effective.

Alpha diversity of AMF
A total of 381,228 raw tags were obtained from the Illumina MiSeq platform sequence of 12 soil samples. After removal of reads with a low quality score (< 20), ambiguous nucleotides, incomplete barcodes and chimeras, there were 355,011 high quality sequences with a range of read length of 320-360 bp produced, accounting for 93.1% of total sequences. All the high quality sequences were classified into 99 OTUs (97% similarity cutoff ). Numbers of OTUs, Good's coverage, richness and diversity under the different fertilization regimes are shown in Table 2. The Good's coverage values exceeded 99.9%, and the numbers of OTUs under different fertilization regimes were in the range of 62.7-89.0. Moreover, the NPK treatment significantly decreased the richness index (Chao1). The application of manure led to the greatest values of Shannon diversity index, which was lowest for NPK. Compared with CK, the NPK treatment significantly decreased PD whole tree index.

Community composition at different taxonomic levels of AMF
In addition to altering alpha diversity of AMF, chronic fertilization affected their community composition as well. At the order or higher taxonomic levels, AMF community compositions were very similar under the different fertilization regimes. All members of AMF were affiliated with the phylum Glomeromycota, and more than 82% belonged to class Glomeromycetes, accounting for 82.4-91.8% of total sequences. Four different orders were detected: Glomerales, Diversisporales, Archaeosporales and Paraglomerales.
At the family level, the Glomeraceae were dominant, followed by Claroideoglomeraceae, Gigasporaceae, Paraglomeraceae and Archaeosporaceae. The families with significant differences among fertilization regimes are presented in Fig. 1. Abundances of Glomeraceae and Gigasporaceae were lower under all fertilization regimes compared with CK, but Claroideoglomeraceae exhibited the opposite pattern. Furthermore, both the M and MNPK treatments significantly increased the Claroideoglomeraceae abundance, with the highest for MNPK.
The AMF community compositions at genus level in response to chronic fertilization were quite different (Table 3, at least one group with a relative abundance > 1%). Compared with CK, the NPK treatment significantly increased the abundance of Glomus, but this decreased for MNPK. The abundances of Funneliformis were lower under all fertilization regimes compared with CK, with the lowest for NPK. Both inorganic fertilizer and manure significantly increased the abundance of Rhizophagus. The NPK treatment significantly decreased the abundance of Septoglomus, but manure increased it. It should be noted that the application of manure in both M and MNPK treatments had a significant positive effect on Claroideoglomus compared with CK and NPK treatments. In addition, there was a significant decline in Gigaspora abundance for all fertilization regimes compared with CK.

Table 1 Soil properties and wheat yield in different fertilization regimes
Values are means ± standard deviations (n = 3). Values within the same column followed by different letters indicate significant differences (P < 0.05) according to Tukey's multiple comparison

Beta diversity of AMF
The PCoA showed a clear separation among different fertilization regimes (Fig. 2). The axes of PC1, PC2 and PC3 explained 72.78, 12.25, and 6.47% of the variation, respectively. The plots of all fertilization regimes were well distributed along the PC1 axis with an order of CK, M, NPK and MNPK, indicating obvious shifts in AMF community composition. Moreover, AMF community compositions between treatments with or without inorganic fertilizer were quite different, suggesting a strong effect of inorganic inputs.

Relationships between AMF community and edaphic factors
The results of the RDA visualized the effects of edaphic factors on AMF community composition (Fig. 3). All the soil chemical properties selected in the current study accounted for 90.2% of the variation between samples. The first constrained axis explained 65.85% and the second explained 13.56% of the total variance. According to the forward-selection option in CANOCO with the Monte Carlo test, the most important contributors affecting the AMF community were soil AP (F = 11.1, P = 0.002), OM (F = 5.1, P = 0.012) and pH (F = 3.6, P = 0.026), accounting for 58.2, 19.0 and 10.5% of the variation, respectively. The other soil properties affected

Chronic fertilization changed rhizosphere soil pH and OM
Soil properties in the wheat rhizosphere differed significantly among fertilization regimes. Soil pH played a key role in soil quality and crop production. Although soil pH was significantly higher for the MNPK than the NPK treatment, both treatments significantly decreased soil pH compared with CK. This confirmed the positive impact of inorganic fertilizer on soil acidification (Guo et al. 2010), as well as the positive effect of manure on alleviation of soil acidification, likely due to the buffering functions of organic acids, carbonates and bicarbonates (Whalen et al. 2000;Garcia-Gil et al. 2004). Furthermore, manure containing macronutrients obviously led to accumulation of soil OM (Xie et al. 2014).

Chronic fertilization affected AMF richness and diversity
In the current study, chronic fertilization shifted alpha diversity of AMF, probable due to the changes of soil chemical properties. Inorganic P and N fertilizers supplied these soil resources for plants and so, in turn, fewer nutrients would be allocated to underground mycorrhizae (Johnson et al. 2013;Williams et al. 2017), leading to low AMF richness and diversity under NPK regium. Under enriched soil N and P conditions, AMF were found to reduce the numbers of arbuscules and extraradical hyphae (Johnson et al. 2008), and plants invested less in symbiotic development (Smith and Read 2010). This was probably due to less carbohydrates in root exudates when plants experienced high P levels (Bhadalung et al. 2005). Moreover, high N and K concentrations also promote negative effects of P on AMF (Guttay and Dandurand 1989). However, contrasting results of positive effects of N fertilization were reported by Zheng et al. (2014). This might be caused by the P level, because N fertilization was found to increase AMF species richness and diversity in P-limited soils, but caused reductions in P-rich soils (Egerton-Warburton et al. 2007;Cheng et al. 2013). In contrast, organic fertilization was beneficial for AMF richness . As previously documented, OM additions tend to increase AMF colonization of plant roots and hyphal abundance in soils (Gosling et al. 2010), in turn, the hyphae preferentially proliferate in organic substrates in experimental microcosms  and can capture N from organic substrates (Leigh et al. 2009). In addition, multiple lines of evidence suggest a key role for AMF in cycling nutrients via organic sources (Sheldrake et al. 2017).
The AMF functional traits are phylogenetically conserved (Powell et al. 2009), thus, the decreased alpha diversity of AMF probably result in loss or reduction in AMF biological functions . Additionally, the decline in AMF species diversity could increase ecosystem instability and P losses, and result in accumulation of fungal pathogens (Whipps 2004;Verbruggen et al. 2012), while decreasing the exchange of nutrients between plants and AMF (Maček et al. 2011) and so making soil more easily affected by disturbances related to climate change . Plants may allocate biomass to the structures and stimulate the AMF symbiosis to acquire more nutrients given the limited resources or vice versa (Johnson et al. 2013). Consequently, the similar levels of available nutrients with NPK and MNPK treatments probably resulted in similar selection effects on AMF taxa, as well as no significant differences in diversity.

Chronic fertilization shaped AMF community composition
Soil properties play a crucial role in shaping the rhizosphere AMF community . Selection by certain soil nutrients or competition between species can lead to changes in overall AMF abundance (Engelmoer et al. 2014) or even complete exclusion of particular species (Roger et al. 2013), resulting in shifts in AMF community composition. In our study, AMF community composition varied significantly among different (square) NPK, inorganic phosphorus, potassium and nitrogen fertilizer; (cross) MNPK, inorganic P, K and N fertilizer plus manure fertilization regimes, in good agreement with previous findings . The most common AMF family was Glomeraceae, which agreed well with other findings in arable fields and in several other temperate ecosystems Lin et al. 2012). Claroideoglomeraceae was positively related to manure application. The high abundance may have suppressed microbes with known pathogenic traits, as members within Claroideoglomeraceae (i.e., Glomus etunicatum) efficiently inhibited root disease caused by fungal pathogen (Sikes et al. 2009). Furthermore, Claroideoglomeraceae can also promote P uptake . In our previous study (Ding et al. 2017), various microbes besides Glomeraceae with known pathogenic traits also had lower levels following manure application. In comparison with CK, all fertilizer application significantly decreased the abundance of Gigasporaceae, which was negatively correlated with AP in soil ). This could be caused by plants allocating less biomass and energy to roots, thereby negatively affecting the AMF symbiosis under conditions of high AP in soil (Johnson et al. 2013).
At genus level, most taxa belonged to genus Glomus in our study, in accordance with previous findings , as Glomus can survive and propagate more easily due to the ability to colonize via pieces of mycelium or mycorrhizal root fragments (Daniell et al. 2001). Compared with other treatments, the abundance of Glomus was significantly higher in the NPK, which was similar to previous results (Liu et al. 2015b). Additionally, all fertilization treatments decreased Funneliformis abundance, while manure application showed significantly positive effect, probably due to the positive correlation with soil pH (Jansa et al. 2014). Manure application also benefited the abundance of Claroideoglomus, indicating positive effect on facilitating plants to assimilate more P in soil (Wang et al. 2017). The results confirmed other findings that Claroideoglomus was more abundant with organic than conventional farming as well ). In addition, all fertilization treatments significantly increased Rhizophagus abundance, probably due to high levels of nutrients (Thonar et al. 2011). As Jansa et al. (2008) stated, species of R. intraradices within Rhizophagus showed the highest benefit to hosts, indicating a beneficial effect of fertilization, particularly using organic manure. In comparison with CK, all fertilizer treatments with a high level of AP showed lower abundance of Gigaspora. This is likely because sufficient soil nutrients, especially AP, could shift species composition in favor of less efficient mutualists by plants allocating less biomass and energy to roots (Johnson 1993).

Relationships between soil properties and AMF community composition
AMF community compositions were strongly affected by fertilizer-driven changes in soil chemical properties. In turn, shifts in AMF composition could influence plant growth through nutrient turnover, disease incidence and disease suppression, as well as soil quality. For instance, AMF produce large amounts of insoluble glycoprotein, glomalin and polysaccharides, which contribute to aggregate stability (Xie et al. 2015). In the current study, the PCoA revealed a significant correlation between soil fertilization regime and the AMF community. The application of inorganic fertilizer, especially P and N inputs, strongly affected the AMF community (Wang et al. 2011;Chen et al. 2014). In comparison, less disturbance was observed under M treatment, due to the biological processes in the organically managed agriculture system (Harrier and Watson 2004). The results agreed well with previous findings that the AMF community was more similar to those occurring under organically fertilized fields than those under conventional management .
The soil pH and AP were considered important contributors to shifts in the AMF community (Dumbrell et al. 2010;Xu et al. 2017). AP was crucial for AMF, as the P level in soil was directly associated with the allocation of nutrients and energy to underground mycorrhizae (Johnson et al. 2013;Williams et al. 2017). Meanwhile, the AMF played important roles in the P cycle. Most phosphate in soils is in the form of ortho-phosphate, which cannot be directly utilized by plants. This phosphate was absorbed by plant roots or by roots colonized by mycorrhizal fungi (Behie and Bidochka 2014). AMF can secrete enzymes that hydrolyze organic P compounds in soil. In addition, OM was also crucial in shaping the AMF community. As heterotrophs, AMF utilize exogenous C for growth. The soil OM and plant root exudates would have profound influences on AMF community structure. Soil pH was also an important factor that shaped the AMF community, which was in agreement with Xu et al. (2017).