The experimental site was established in 2004 at the National Field Scientific Observation and Research Station for Dry Crop Farming Systems in the Inner Mongolia Autonomous Region, China (41° 08′ 22.8″ N and 111° 17′ 43.6″ E). The area features a mid-temperate continental monsoon climate at an altitude of 1,570 m, annual rainfall of 250–400 mm, potential evaporation of 1,848.3 mm, and more than 80% of the rain falls between June and September. The area exhibits an average annual temperature of 1.5–3.7°C and a frost-free period of 90–120 d. The ecosystem represents a typical semi-arid agro-pastoral ecotone. The cropping system is a one-season crop and the soil type is chestnut soil. The initial soil physical-chemical properties in 2004 (in the 0–20 cm layer) were shown in Table 1. The precipitation (195.5–414.9 mm, average 300.08 mm) and the average temperature change (14.6–17.9°C, average 16.2°C) during the crop growth period (from May to September) during the 2004–2019 experiment.

Different fertilization treatments were established within the long-term fertilization experimental site including chemical fertilization alone (CF), organic manure alone (M), organic manure with chemical fertilization (CFM), and a non-fertilization control (CK). All of the treatments were established in a randomized complete block design with three replicates and plot sizes of 50 m² (5 m × 10 m). Fertilizer dosages and ratios were designed and recommended by the International Plant Nutrition Institute (IPNI) based on soil testing results. Urea was used for nitrogen nutrients, calcium superphosphate for phosphate, potassium chloride for potash, and sheep dung for organic manure. Its water content is 60–65%, organic matter content is 36–45%, and nutrient content (%) is N-P₂O₅-K₂O = 0.57-0.23-1.04. Sheep dung was added at 7,500 kg hm⁻² in 2004–2005 and 2007-2012, and 15,000 kg hm⁻² in 2006 and 2013-2019, the specific fertilization amounts are shown in Table 2, including the N, P₂O₅, and K₂O components within sheep dung. Potato, rape, and oat were used as rotation crops, with one crop planted every year. In 2019, the potato was planted. Crops were generally planted in early May and harvested in late September. During the 16 years (2004–2019), organic manure was evenly spread on the soil surface each year before sowing and artificially poured into the soil. Chemical fertilization was applied when sowing and all fertilizers were only applied once. The experiment was conducted under rainfed conditions, while other field management measures including pest control followed standard tillage measures.

A total of 60 soil samples (four treatments × five random points × three replicates) (2 cm diameter and 0–20 cm depth) were collected from each plot after potato harvest in 2019 in this experiment, and then totally mixed the soils from five random points in each plot. Finally, and the 12 soil samples were obtained, the 12 soil samples were divided into two portions after removing impurities, with one taken to the laboratory in a self-sealing bag, followed by natural airy-drying and filtering with a 2 mm sieve to determine soil nutrient contents. The other portion was sealed in sterile bags and immediately taken back to the laboratory for storage at −80°C for microbial community analyses.

A total of 36 soil samples (four treatments × three random points × three replicates) for the surface soil (0–20 cm) bulk density calculation were taken using cutting rings (volume 100 cm³, inner diameter 5.05 cm, 200–300 g) on three randomly selected points in each plot.

Soil pH was measured with a compound electrode (PE-10, Sartorius, Germany) (water: soil = 2.5:1). Other parameters were measured as previously described including soil bulk density (BD) (O’connell, 1975), organic matter (OM) (Nelson and Sommers, 1982), total nitrogen (TN) (Sparks et al., 1996), available nitrogen (AN), available phosphorus (AP), and available potassium (AK) (Lu 2000).

Total genomic DNA was extracted from 0.5 g of composite soil samples using the CTAB extraction method following the manufacturer’s protocol (MP Biomedicals, Illkirch, France). DNA concentration and purity were assessed with 1% agarose gel electrophoresis. Fungal ITS sequences were amplified with the primer ITS5F (5′GGA​AGT​AAA​AGT​CGT​AAC​AAG​G3′) and ITS2R (5′TCC​TCC​GCT​TAT​TGA​TAT​GC-3′) (White et al., 1990) while the V4 hypervariable region of bacterial 16S rRNA genes were amplified using the primer pair 515F (5′-GTGCCAGCMGCCGCCGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Berthrong et al., 2013). Illumina MiSeq sequencing of the amplicons was then conducted at Beijing Novogene Technology Co., Ltd. (Novogene, Beijing, China).

PCR reactions were conducted in 30 μL reactions including 15 μl of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM of forward and reverse primers, and about 10 ng of template DNA. The PCR conditions included an initial denaturation at 98°C for 1 min, followed by 30 cycles with denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 30 s, all followed by a final extension at 72°C for 5 min. PCR amplification success was evaluated with gel electrophoresis using 1x loading buffer (containing SYBR green) and PCR products, with a 2% agarose gel. PCR products were mixed in equimolar ratios and the product pools were purified with the GeneJET^TM Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using the Ion Plus Fragment Library Kit (Thermo Scientific) following manufacturer’s recommendations. Sequencing library quality was assessed with the Qubit@ 2.0 Fluorometer (Thermo Scientific). The library was then sequenced on an Ion S5TM XL platform by generating 400 bp/600 bp single-end reads (Sheng et al., 2019; Wang et al., 2021).

Single-end reads were assigned to samples based on their unique barcodes, followed by truncation via removing barcode and primer sequences. Quality filtering of raw reads was performed to obtain high-quality clean reads using the Cutadapt (Martin, 2011) (V1.9.1, http://cutadapt.readthedocs.io/en/stable/) quality control pipeline. The reads were compared against the SILVA reference database (https://www.arb-silva.de/) (Quast et al., 2013) using the UCHIME algorithm (http://www.drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011) to detect chimeric sequences, followed by chimera removal (Haas et al., 2011). Sequence clustering was conducted with the UPARSE software package (v.7.0.1001, http://drive5.com/uparse/) (Edgar, 2013), wherein sequences with ≥97% nucleotide similarity was assigned to the same operational taxonomic units (OTUs). Representative sequences for each OTU were compared against the SILVA Database (https://www.arb-silva.de/) (Quast et al., 2013) using the Mothur classification algorithm to taxonomically annotate the OTUs. To investigate the phylogenetic relationships among OTUs and differences in dominant species among different sample groups, multiple sequence alignments were conducted using the MUSCLE software program (v.3.8.31, http://www.drive5.com/muscle/) (Edgar, 2004).

The raw sequences were submitted to the Sequence Read Archive (SRA) database under the BioProject database of the National Center for Biotechnology Information (NCBI) platform (project identification numbers PRJNA785382 (https://submit.ncbi.nlm.nih.gov/subs/bioproject/SUB10748118/overview) and PRJNA787050 (https://submit.ncbi.nlm.nih.gov/subs/bioproject/SUB10779674/overview) for bacterial and fungal sequences, respectively).

Soil physical and chemical properties, relative abundances of dominant bacterial and fungal taxa (average relative abundances ≥1%), and microbial diversity were analyzed using the least significant difference (LSD) at 5% level by SPSS program (v.24.0, IBM, United States). Permutational multivariate analysis of variance based on Bray–Curtis dissimilarity matrices using the “adonis” function in R and Principal coordinates analysis (PCoA) were used to test and visualize cultivar on the microbial communities using the R packages “vegan”, “pairwise”, “ade4” and “ggplot2”. Taxa with differential abundances among treatments were identified using linear discriminant analysis (LDA) effect size (LEFSe) analyses. Differentially enriched taxa (e.g., biomarkers) were defined based on discriminative characteristics of LDA score thresholds (log10 value) > 3.0 and p values <0.05 for between-group factorial Kruskal-Wallis tests. The relationships among soil properties and dominant soil microbe were examined using Pearson correlation analysis by SPSS program. Redundancy analysis (RDA) was used to evaluate changes in soil microbial community composition in association with physico-chemical properties and to determine the primary factors that drove these changes. RDA was performed using the CANOCO software program (v. 5.0, Ithaca, NY). Partial least squares path models (PLS-PMs) were used to model the relationships among soil physico-chemical properties in association with bacterial and fungal diversity and community composition. The path coefficients and the coefficients of determination (R²) in the path models were estimated in R (4.0.3) using the plspm package (1,000 bootstraps) (Sheik et al., 2012; McMurdie et al., 2013).

Sixteen years of fertilization led to significant changes in soil physico-chemical properties (Table 3). Soil BD decreased in M and CFM compared with that in CF and CK. Soil TN, AN, AP, and AK contents had the greatest increases in CFM. Soil OM content increased significantly in M and CFM, but decreased in CF compared with CK. The CF treatment significantly reduced soil pH compared with that in CK.

The number of bacterial operational taxonomic units (OTUs) in M was significantly higher than that in CK, whereas the Chao1 index in CF was lower than that in the other treatments. Fungal OTU numbers were higher in M and CFM than in CK, the Chao1 index was significantly higher in M, whereas it was not significantly different in CF and CFM. The Shannon index of bacterial and fungal communities were not significantly different among treatments (Table 4).

Variation in structure of soil bacterial and fungal communities in different treatments was investigated using principal coordinates analysis (PCoA) (Figure 1). Treatments M and CFM were clustered together, indicating that structure of those microbial communities was similar, and they were also significantly separated from CF and CK along the PCo1 axis.

In bacterial community, the dominant phyla across all samples included Proteobacteria (28.6–35.4%), Actinobacteria (25.4–33.5%), Acidobacteria (11.9–14.8%), Bacteroidetes (7.1–12.1%), with Gemmatimonadetes, Chloroflexi, Verrucomicrobia, Thaumarchaeota, Firmicutes, and Nitrospirae accounting for 13.5–16.8%. At the phylum level, 2.6%–3.3% of the bacterial reads were unclassified (Figure 2A). Abundances of Proteobacteria were significantly higher in M and CFM than in CK and CF, whereas abundances of Actinobacteria and Chloroflexi were significantly lower than those in CK and CF. Abundance of Acidobacterial was highest in CK, whereas abundances of Bacteroidetes and Gemmatimonadetes were significantly higher in CFM than in CK (Supplementary Table A1). Responses of the 20 most abundant genera of bacteria indicated that abundances of seven dominant genera were significantly different among fertilization treatments. Abundances of Sphingomonas, unidentified_Acidobacteria, and Bryobacter were higher in M than in other treatments. In addition, abundances of Lysobacter and Iamia were higher in CFM, whereas abundances of Nocardioides and Pseudarthrobacter were higher in CF communities than in other treatments. Blastococcus and Gaiella abundances were lower in CF, M, and CFM than in CK (Supplementary Table A2). Cluster analysis indicated that dominant genera in CFM were similar to those in M and CF but were significantly different compared with those in CK (Figure 2C).

In fungal communities, Ascomycota was the dominant phylum (69.1–72.6%), followed by Mortierellomycota (3.7–14.4%), Basidiomycota (2.3–7.5%), and Mucoromycota (0.5–3.4%). In addition, Chytridiomycota, Glomeromycota, Olpidiomycota, Aphelidiomycota, Rozellomycota and Monoblepharomycota accounted for 0.4–0.9% of all fungal communities. At the phylum level, 12.9%–18.0% of all fungal reads could not be classified (Figure 2B). Abundances of Mortierellomycota and Basidiomycota were significantly higher in CF than in the other treatments, whereas abundances of Ascomycota were significantly lower in CFM than in CK, CF, and M (Supplementary Table A3). Only six of the dominant genera had abundances that were significantly different among fertilization treatments. Compared with CK, Mortierella was significantly more abundant in CF, M, and CFM, whereas Acremonium, Alternaria, and Dactylonectria were significantly less abundant in CF. In addition, Fusarium was significantly more abundant in CF, whereas Rhizopus was less abundant in CF and CFM compared with CK (Supplementary Table A4). Cluster analysis indicated that dominant genera in CFM and CF were similar, whereas genera in CK were significantly different from those in the other treatments (Figure 2D).

In a LEfSe analysis, LDA values >3.5 indicated significant enrichment of microbial taxa, and those responsive taxa were considered “biomarkers” for different treatments. Forty-five taxa, inclusive of all taxonomic levels, were identified that best separated soil bacterial communities among the four fertilization treatments. Twenty-two enriched taxa were identified in CFM, fourteen in CK, seven in M, and only two in CF (Supplementary Figure A1). The class Gammaproteobacteria was one of the most predominant biomarkers in CFM, whereas relatively high abundances identified the families Sphingobacteriaceae and Pyrinomonadaceae as biomarkers in M and CK. The CF treatment included only a single biomarker, the class Chloroflexia. Fifty significantly enriched fungal taxa (LDA score >3.5) distinguished soil fungal communities among different treatments, including 12 taxa in M, 18 in CK, 10 in CFM, and 10 in CF (Supplementary Figure A2). The genus Mortierella, order Hypocreales, and the class Sordariomycetes were identified as biomarkers in CFM, M, and CK, respectively. In addition, abundances of the family Cladosporlaceae were significantly higher in CF than in other treatments.

Correlation analysis revealed that soil chemical properties affected microbial community composition (Table 5). In bacterial communities, Proteobacteria was significantly positively correlated with soil OM and TN, whereas Actinobacteria was significantly negatively correlated with soil TN, AN, and pH. In fungal communities, Ascomycota was significantly negatively correlated with soil AN, AP, AK, and pH. Redundancy analysis (RDA) revealed that soil chemical properties were associated with phylum-level composition of bacterial and fungal communities in CFM. Composition of the bacterial communities in CF and composition of fungal communities at the phylum-level in M and at the genus-level in CK were greatly affected by soil pH and BD. Consequently, changes in soil pH induced by fertilization were reflected in responses in microbial abundances (Figure 3).

PLS-PM analysis was used to evaluate relations among fertilization, soil physico-chemical properties, and microbial community composition and diversity (Figure 4). The model constructed indicated that differences in fertilization types (CK→CF→M→CFM) had significant and direct positive effects on OM (1.052), nutrient contents (a combination of N, P, and K; 0.895), and especially fungal communities (1.527), in addition to significant and direct negative effects on soil BD (−0.610) (Supplementary Table A5). However, differences in fertilization did not directly affect on bacterial communities. In addition, soil pH directly and strongly affected composition of bacterial and fungal communities, whereas soil OM and nutrient contents did not significantly affect composition of microbial communities.

Proper fertilization is a critical practice to stabilize soil productivity, improve soil quality, and promote the sustainable development of farmland ecosystems. Soil BD is a direct physicochemical indicator of soil quality conditions and is one of the most sensitive soil properties to long-term fertilization (Lin et al., 2018). In CK in this study, soil BD increased by 5.56% compared with that in 2004. By contrast, soil BD significantly decreased in CFM and M, which could be because organic manure supplements provided a cementing agent for soils, thereby resulting in changes to soil structure and promoting crop root growth. The combination of root secretions and mineral colloids can help improve soil BD (Qiao et al., 2020). Soil OM and nutrient (TN, AN, AP, and AK) contents were higher in M and CFM than CK. Humic acids in soil organic C pools can increase plant root vigor and increase root respiration along with nutrient uptake, thereby facilitating accumulation of soil nutrients (Furukawa et al., 2014). The CF treatment significantly decreased soil pH, with a clear trend toward soil acidification. Low pH due to proton deficiency caused by chemical fertilizer application also reduced soil nutrients (Tian and Niu, 2015). In M and CFM, the application of organic manure increased levels of alkaline cations in soils and helped to neutralize soil acidity. Thus, M and CFM treatments maintained a high soil pH, and protonation of organic anions to form central molecules was the main mechanism to resist soil acidification (Shi et al., 2019). Alternatively, the experimental area had arid chestnut soils with high CaCO₂ contents, which could have reduced the risk of soil acidification from chemical fertilizer application (Sarfraz et al., 2017). In addition, the pH in both M and CFM soils decreased slightly compared with the initial soil pH (Table 1), and it was inferred that the mechanism of the slight decrease in soil pH might be related to the high initial pH of the alkaline soils. Most studies conclude that chemical fertilizer application does not directly increase soil nutrient contents, but primarily increases accumulation of soil nutrients by increasing plant root stubble, root systems, and root secretions while stimulating the decomposition of pre-existing soil nutrients. Thus, although the CF treatment also increased soil nutrients, overall increases were limited (Knoblauch et al., 2017). The CF treatment reduced soil OM content because of long-term deficiencies in soil organic matter inputs. In this study, CFM was the treatment that most effectively improved soil properties, while also providing a supplementary input of in-season nutrient sources to soils. However, the risk of further soil acidification associated with the long-term application of chemical fertilizers is also a concern and reinforces the importance of using organic manure–chemical fertilizer to improve soil chemical properties.

Microorganisms have significant roles in forming soil fertility, improving ecological characteristics, and preventing soil borne diseases. As a consequence, abundances and diversity of microorganisms affect the conversion of soil nutrients and the modulation of biological effectiveness, which are critical signs of soil health (Ji et al., 2014). According to the PCoA changes in fertilization had relatively large effects on composition of microbial communities. Proteobacteria, Bacteroidetes, and Gemmatimonadetes are copiotrophic taxa that can effectively deplete labile C pools in soils, and they also have high nutrient requirements (Ghosh et al., 2016). In this study, Proteobacteria, Bacteroidetes, and Gemmatimonadetes were significantly more abundant in CFM and M, which also had higher OM and nutrient contents than those in CK. Organic manure application favors the growth of copiotrophs, whereas high abundances of copiotrophic taxa also indicate high soil fertility. The Proteobacteria comprises many taxa involved in soil N cycling, including Sphingomonas and Lysobacter, which varied in abundance similarly to that of Proteobacteria overall. Thus, organic manure application can improve soil N fixation capacity. Acidobacteria and Chloroflexi are oligotrophic taxa that consume recalcitrant organic C pools and have high nutrient affinities, despite relatively slow growth rates (Delgado-Baquerizo et al., 2017). In this study, Acidobacteria and Chloroflexi were most abundant in CK with low OM contents and soil nutrients and therefore overall low fertility because of prolonged lack of fertilization. Such conditions can drive bacterial communities toward dominance by oligotrophs. Unidentified_Acidobacteria and Bryobacter are members of Acidobacteria, and abundances of those taxa varied similarly to that of Acidobacteria overall. Actinobacteria are mostly pathogenic bacteria, and in this study, abundances were highest in CK and CF. Thus, long-term chemical fertilization and the lack of fertilization can lead to soil nutrient imbalances and increased abundances of pathogenic bacteria, thereby increasing the risk of soil bacterial diseases that endanger healthy soil environments. Blastococcus, Nocardioides, Gaiella, Pseudarthrobacter, and Iamia are members of Actinobacteria and abundances of those taxa varied similarly to that of Actinobacteria overall, except for Iamia. Lysobacter and Sphingomonas are associated with antagonistic activities toward plant pathogens (Khan et al., 2014; Liu et al., 2020), and in this study, those taxa had high abundances in CFM, suggesting that application of organic manure with chemical fertilizer improved soil microenvironments.

Ascomycota was dominant in fungal communities in this study. Soil fungi use OM as energy sources and nutrients, participate in OM mineralization, decompose OM from crop residues or organic manure that enter soils, and supply plants with nutrients (Detheridge et al., 2016). However, most fungi are also plant pathogens. In this study, 16 of the 20 most abundant fungal genera were members of Ascomycota, and abundances of those genera were significantly lower in CFM than in the other treatments. To some extent, those results reflected a reduced incidence of fungal soilborne diseases. Fusarium, Gibberella, and Dactylonectria can cause plant stem rot or spikelet rot (Larkin and Griffin, 2007; Cannon et al., 2012), whereas Paramyrothecium and Plectosphaerella can cause plant root rot (Carlucci et al., 2012). Paramyrothecium abundances were significantly lower in CFM than in other treatments. Furthermore, Alternaria, which can cause potato brown spot and potato early blight, was most abundant in CK and CF (Xu et al., 2022). Cladosporium is a saprophytic taxon but is also an important plant pathogen that mainly infects plant leaves, branches, and fruits, resulting in poor quality and reduced yields (Bensch et al., 2012). Cladosporium abundance was significantly higher in CF than in the other treatments. The genus Mortierella in Mortierellomycota is important in improving nutrient uptake efficiency and protecting crops from adverse conditions by improving the availability of soil P, K, and ferric iron (Ozimek and Hanaka, 2021). Vishniacozyma is a member of Basidiomycota, which had high abundances in CF, although the ecological functions of the genus remain unknown. Long-term fertilization increased relative abundances of fungi in communities, which might lead to soil fungalization. Overall, the results of this study suggest that long-term fertilization affects the structure of microbial communities, with chemical fertilization increasing pathogenic fungal loads and thereby increasing the risk of disease infestation. However, combining organic manure and chemical fertilization can reduce fungal pathogen loads and increase abundances of beneficial fungal genera.

Soil physical and chemical properties greatly affect the composition of soil microbial communities. In this study, dominant microbial taxa were consistently positively correlated with soil OM and nutrient contents in CFM. Those results were likely because the combined manure–chemical fertilization application increased soil OM and nutrient contents, which are necessary nutrients for microbial growth. Thus, high OM content contributes to improving soil quality and maintaining soil health (Mishra et al., 2018). High soil fertility results in increased abundances of copiotrophic populations. In this study, Proteobacteria was significantly positively correlated with soil OM and TN, as were microbial taxa that could inhibit some pathogens (Raoul des Essarts et al., 2016), including the fungus Mortierella, which was highly abundant in CFM. Soil fungal diversity and abundances of dominant fungi were significantly higher in CF than in other treatments, although soil OM content was lowest in that treatment. Fungi primarily use C sources and nutrients for growth and reproduction, which are obtained from the decomposition of organic matter. This is particularly evident for large fungi that can consume large molecular-weight organic matter. Such interactions might explain the decrease in OM resulting from chemical fertilization. The dominant microbial taxa in CK were negatively correlated with soil nutrient contents, suggesting that poor fertility induces explosive growth of fungal pathogens, as indicated by significant negative correlations of Ascomycota with soil AN, AP, and AK. The PLS-PM analysis also indicated that fertilization treatments directly affected fungal communities, in contrast to bacterial communities, indicating that fungi were more sensitive to fertilization treatments than bacteria (Cassman et al., 2016; Ai et al., 2018). In addition, the PLS-PM analysis indicated that the soil pH was primary factor that significantly affected soil microbial community composition. Those results are in contrast to those of many studies suggesting that nutrient changes due to fertilization directly drive shifts in microbial communities. Xun et al. (2015) proposed the“rugby ball model”to explain this phenomenon. In that model, near-neutral soils (characterized by the central part of the rugby ball) are strongly influenced by soil nutrients that also affect microbial community transformations. The radius of the central part of the rugby ball is larger and represents the strength and direction of the driving effects of nutrient indicators. When soil acidification or alkalization increases (characterized by the two ends of the rugby ball), soil pH has an increasing effect on microbial community transformations, whereas the effects of soil nutrients are gradually minimized (as represented by movement towards the ends of the ball, with the cross-sectional radius of the rugby ball). In this study, different fertilization treatments resulted in significant changes in soil pH in dryland soils. Changes in pH were especially evident with long-term chemical fertilization alone, with increases in soil pH in that treatment increasing the risk of soil acidification. In this case, the negative effects of soil acidification outweighed improvement in soil fertility, as indicated by the decrease in soil OM content. By contrast, the combined use of organic manure and chemical fertilization significantly improved soil pH while also supplementing soil nutrients, thereby promoting and realizing the full positive effects on soil nutrients.

In this study, different fertilization types affected soil physico-chemical properties and microbial community characteristics. However, only the overall changes in soil microbial communities were evaluated, and microorganisms in unique functional groups were not investigated. Thus, differences in soil metabolic pathways and functions associated with unique functional groups of microorganisms necessitate further research in order to determine whether differences in biomass, diversity, and community structure under different fertilization practices lead to changes in soil fertility, health status, and ecosystem function.