Long-term combined application of manure and chemical fertilizer sustained higher nutrient status and rhizospheric bacterial diversity in reddish paddy soil of Central South China

Bacteria, as the key component of soil ecosystems, participate in nutrient cycling and organic matter decomposition. However, how fertilization regime affects the rhizospheric bacterial community of reddish paddy soil remains unclear. Here, a long-term fertilization experiment initiated in 1982 was employed to explore the impacts of different fertilization regimes on physicochemical properties and bacterial communities of reddish paddy rhizospheric soil in Central South China by sequencing the 16S rRNA gene. The results showed that long-term fertilization improved the soil nutrient status and shaped the distinct rhizospheric bacterial communities. Particularly, chemical NPK fertilizers application significantly declined the richness of the bacterial community by 7.32%, whereas the application of manure alone or combined with chemical NPK fertilizers significantly increased the biodiversity of the bacterial community by 1.45%, 1.87% compared with no fertilization, respectively. Moreover, LEfSe indicated that application of chemical NPK fertilizers significantly enhanced the abundances of Verrucomicrobia and Nitrospiraceae, while manure significantly increased the abundances of Deltaproteobacteria and Myxococcales, but the most abundant Actinobacteria and Planctomycetes were detected in the treatment that combined application of manure and chemical NPK fertilizers. Furthermore, canonical correspondence analysis (CCA) and the Mantel test clarified that exchangeable Mg2+ (E-Mg2+), soil organic carbon (SOC) and alkali-hydrolyzable nitrogen (AN) are the key driving factors for shaping bacterial communities in the rhizosphere. Our results suggested that long-term balanced using of manure and chemical fertilizers not only increased organic material pools and nutrient availability but also enhanced the biodiversity of the rhizospheric bacterial community and the abundance of Actinobacteria, which contribute to the sustainable development of agro-ecosystems.

nutrient content and nutrient release rate of organic fertilizers make them less likely to meet crop requirements, what's more, compared with chemical fertilizers, the higher amount and inconvenient usage of organic fertilizers resulted in a less prevalence. Nonetheless, combined using of organic and chemical fertilizers has been proven as an effective approach to improve soil fertility, crop yields and environmental quality compared with the single use of either of them [10][11][12] .
Microorganisms are the key components of the soil ecosystems and involve in nutrient cycling, energy flow, and organic matter degradation 13 . Bacteria, as the most abundant and diverse group of soil microorganisms 14 , which participate in biogeochemical nutrient cycling and organic matter decomposition processes 2 and maintain ecological balance 15 . Clearly understanding the responsive shifts of bacterial community structure and composition to different fertilization regimes is vital for developing better fertilization practices and for further improving soil fertility and function 2 . High-throughput sequencing technologies had provided important ways to determine microbial information in recent years, and many studies had generated more valuable soil microbial information through experiments involving long-term fertilization practices. However, most had focused on dryland soil agro-ecosystems [16][17][18][19] , just a few studies had reported microbial community characteristics in paddy soils. For example, the study by Daquiado et al. 20 was performed in silty mixed mesic Typic Haplaquepts and specific climates, while the observation of Chen et al. 2 focused on the effects of different ratios of organic and inorganic fertilizers, both of them viewed on bulk soil bacteria. Actually, the microbial community in the rhizosphere had more important impacts on the growth and health of plants than that in bulk soil 21,22 . Therefore, exploring the rhizospheric microbial community is vital for understanding the interactions among soil, microbes and the host plants under long-term fertilization, which can help identify better fertilizer regimes and plant breeding methods for heavily fertilized cropland 23 .
In this study, a long-term field experiment with four different fertilization regimes in a rice-rice cropping rotation system, established in the Red Soil Experimental Station of the Chinese Academy of Agricultural Sciences since 1982, was conducted to investigate the influences of fertilization regime on the rhizospheric bacterial community structure and composition in the Quaternary reddish paddy soil of Central South China. The objectives were (i) to identify the effects of different fertilization regimes on the rhizospheric bacterial diversity and richness of reddish paddy soil; (ii) to explore the relative abundances of dominant bacterial phyla under different fertilization regimes; (iii) to screen specific biomarkers of different fertilization regimes under different bacterial taxonomic levels; and (iv) to clarify the key driving factors for shaping bacterial communities in the reddish paddy soil of Central South China.
General analyses of the sequencing data. Across all soil samples, a total of 985,898 sequence reads were successfully elicited. After removing short and low-quality reads, singletons, replicates and chimeras, 601,336 sequences, ranging from 41,702 to 57,123 per sample, were retained. Based on 97% similarity, a total of 5585 OTUs, ranging from 2279 to 2616 per sample, were obtained across all samples (Table S1). Among the total sequences, ~99.4% were classified as bacteria, with 62 phyla, 146 classes, 185 orders and 303 families.
A rarefaction analysis showed that the number of OTUs observed for 16S did not reach saturation (Fig. S1), which indicated that the sequencing capability was not large enough to capture the complete biodiversity of these communities, as the curves did not reach a plateau by increasing sample size, a result similar to several previous  16,20 . However, the data were sufficient to show differences among the treatments and to reveal the influences of different long-term fertilization regimes on the soil bacterial community.
Richness and biodiversity of the bacterial community in the rhizosphere. The coverage indices of all treatments were more than 0.98, which also indicated that the sequencing capability was large enough to capture most of the bacterial community characteristics of each treatment (

Effects of different fertilization regimes on bacterial community structures in the rhizosphere.
The unweighted pair group method using arithmetic averages (UPGMA) was used to compare the phylogenetic relationships of bacterial communities based on the weighted UniFrac similarity index of OTU abundance, and the result is presented in Fig. S2. In general, the 12 investigated communities (three replicates for each of the 4 treatments) were divided into two major groups: For NF and NPKM treatments, the three replicate samples of each treatment were closest and grouped together, the NF and NPKM were further grouped together; for M and NPK treatments, a total of the 6 samples were roughly grouped together, and the two treatments were not separated clearly. Furthermore, two-dimensional principal coordinate analysis (PCoA) based on the weighted UniFrac distance clearly revealed that the soil bacterial communities varied among the different fertilization treatments (Fig. 1).  The four different fertilization treatments are clearly separated and distributed in four quadrants along the first and second components (PCoA1, PCoA2), and the first two axes can explain ~63.17% of the observed variance in total, which indicates that the soil bacterial community structures of the four treatments may maintain great differences, though more information is needed.

Biomarkers of different treatments subjected to different long-term fertilization regimes.
To identify the specific bacterial taxa associated with different fertilization levels, we compared the bacterial communities in NF, NPK, M and NPKM treatments using linear discriminant analysis (LDA) effect size (LEfSe). The cladogram representative of the structure of the microbial communities and their predominant bacteria is shown in Fig. 4. The greatest differences (LDA ≥ 4) in taxa among the four communities are displayed. A total of 9 biomarkers were screened under different taxonomic levels across the 4 treatments. The relative abundances of the class Deltaproteobacteria and order Myxococcales were dramatically higher in M treatment than in other 3 treatments, and they were considered as the biomarkers of M treatment. Similarly, the phylum Acidobacteria, class Betaproteobacteria and order Sc-I-84 were considered as the biomarkers of NF treatment, the phylum Verrucomicrobia and family Nitrospiraceae were regarded as the biomarkers of NPK treatment, and the phyla Actinobacteria and Planctomycetes were deemed as the biomarkers of NPKM treatment.

Relationships between soil properties and bacterial communities.
At the phylum level, canonical correspondence analysis (CCA) was employed to detect the relationships between soil properties and bacterial communities based on the relative abundances of species, which revealed that the first and second CCA components were able to explain 78.63% of the total bacterial variation (Fig. 5). The first component (CCA1), explaining 50.79% of the total variation of bacterial phyla, separated the M and NPKM treatments from NPK and NF treatments, respectively. The second component (CCA2), explaining 27.84% of the total variation of bacterial phyla, separated the NF treatment from NPK treatment. Moreover, CCA identified that different environmental variables had different impacts on the overall bacterial community, with an order as follows: E-Mg 2+ > SOC, AK, AN, E-Ca 2+ , SOC/TN> pH, AP. Furthermore, the Mantel test was employed to investigate Spearman's correlated relationships between the relative abundances of bacterial phyla and the soil properties, which indicated that a series of properties, including SOC, AN, AP, and E-Mg 2+ , were significantly (p < 0.01) correlated with the bacterial community composition in reddish paddy soil. Moreover, the pH, AK, and E-Ca 2+ were also significantly (p < 0.05) correlated with  the bacterial community composition, but the SOC/TN was observed to slightly affect the bacterial community composition.

Discussion
In this study, long-term application of composted manure either alone or combined with chemical NPK fertilizers significantly enhanced the soil pH value, whereas single application of chemical NPK fertilizers significantly declined the pH value compared with no fertilization treatment (Table 1), which was similar with previous studies 2,24 . Above those indicated that different response mechanisms of pH to fertilization regimes, guessing that manure application increased soil pH due to the liming effect of added organic matter and carbonates 25 , but single application of chemical NPK fertilizers decreased the pH for the nitrification processes that converting mineral fertilizer ammonium to nitrate 26 . Moreover, long-term fertilization (NPK, M, NPKM) increased the soil nutrient status in different degrees, among them, combined application of manure and chemical NPK fertilizers can effectively improve the concentrations of SOC, AN, AP, AK, and E-Ca 2+ (Table 1), which in accordance with other long-term experiments in various ecosystems 18,20,27 . All these indicated that long-term balanced using of manure and chemical fertilizers increased organic material pools and nutrient availability and further improved the physical environment of soil and the yield of rice grain 2 . The biodiversity and richness of the microbial community are considered to be critical to the integrity, function, and long-term sustainability of soil ecosystems, but they are usually diminished by agricultural perturbation 28 . In this study, long-term application of chemical NPK fertilizers resulted in a significantly lower rhizospheric bacterial richness index (chao1) ( Table 2), the similar results were observed in bulk and rhizospheric soil by long-term fertilization experiments based on 454 pyrosequencing 2 or Illumina Hiseq sequencing 17 technologies, which mainly resulted from the lower pH subjected to long-term chemical fertilizers input 17 . Additionally, long-term application of composted manure either alone or together with chemical NPK fertilizers facilitated a higher soil bacterial biodiversity, which was significantly distinct from that shaped by long-term application of chemical NPK fertilizers or no fertilization ( Table 2). These results supported previous reports from other long-term fertilizer trials 26,27 , probably due to the effects of manure application on soil physicochemical and biological properties, especially for soil pH 17 and microbial biomass carbon 26 . All these results indicated that multiple bacterial groups dwell in the soil with long-term manure input, which could help for developing more balanced microbiome that allows rapid recovery from environmental stresses or natural perturbations, and then suppressing pathogens from flourishing 29 .
Any environmental change could alter the soil microbial community structure to some extent 30 . The PCoA based on the weighted UniFrac algorithm clearly demonstrated the variations among different fertilization regimes (Fig. 1). These four treatments were definitely separated either along the first component (PCoA1) or the second component (PCoA2) with respect to the bacterial communities. Above results were consistent with the observations by many long-term fertilization studies based on 454 pyrosequencing or Illumina Miseq sequencing methods 29,31 .
Analysis for phylum abundance indicated that Proteobacteria, Acidobacteria, Chloroflexi and Nitrospirae were the top 4 dominant phyla across all treatments or samples (Figs 2, S3), which was roughly in accordance with previous studies 2, 32 . Analysis of variance and LEfSe jointly identified that the phyla Actinobacteria and Planctomycetes were most abundant in NPKM treatment (Figs 3, 4). Actinobacteria play key roles in organic matter decomposition and the humus formation process 33,34 , what's more, which can produce variety of antibiotics 35 to protect roots from pathogenic microorganisms 17 . Similar abundant tendency of Planctomycetes among treatments was reported by Chen et al. 24 based on long-term fertilization trial, which accounting for seven genomes with distinct functional and transcriptional activities 36 . Meanwhile, the most abundant phylum Verrucomicrobia and family Nitrospiraceae observed in NPK treatment (Figs 3, 4). Similarly, Zhou et al. 37 indicated that Verrucomicrobia was more abundant in long-term higher chemical fertilizer input than that in lower or no chemical fertilization, which involved in organic carbon utilization for its abundant transcripts 36 , while Nitrospiraceae was a potential keystone species for iron uptaking and nitrite oxidation 38 . Moreover, the most abundant Acidobacteria was detected in NF treatment (Figs 3, 4), supported previous studies 9, 19,37 , which belong to the oligotrophic (or K-selected) group 39 , its function keeps obscure for the majority of them are ungrouped 40 . Additionally, the most abundant phylum Proteobacteria, class Deltaproteobacteria and order Myxococcales observed in M treatment (Figs 3, 4). Deltaproteobacteria, as one class of phylum Proteobacteria, the abundance of which tremendously increased probably leading to more abundant Proteobacteria. In fact, Deltaproteobacteria comprises a branch that predominantly composed of the genus Haliangium, including the relatives Haliangium luteum 41 and Haliangium ochraceum 42 , which can produce haliangicin, a novel antifungal metabolite that can inhibit the growth of broad-spectrum fungi 43 . Myxococcales, as another abundant species in M treatment, had been reported to antagonize soil-borne plant pathogens 44 . Previous studies had indicated that soil environmental factors affect the microbial community structure 28,45,46 . For example, Tu et al. 16 confirmed that the available Ca and available Mg play important roles in shaping soil bacterial community composition in Paulownia fortunei plantations based on long-term fertilizer treatments. Wei et al. 18 screened the SOC and TN as the main drivers that driven microbial community structure from six environmental variables (pH, SOC, TN, TP, IN and AK) based on 35 years of manure and chemical fertilizer application experiment. However, Wang et al. 17 indicated that soil pH, organic matter and available P concentrations were the most important factors in shaping bacterial communities in the maize rhizosphere. In this study, CCA revealed that E-Mg 2+ , SOC, AK, AN, E-Ca 2+ and SOC/TN were the key factors that driving the distribution and composition of the soil bacterial community (Fig. 5), while the Mantel test screened the SOC, AN, AP and E-Mg 2+ as the key factors (Table 3). Different from previous studies, we concluded that E-Mg 2+ , SOC and AN were the key driving factors in shaping soil bacterial communities based on joint analysis of CCA and Mantel test. Guessing the above various results mainly caused by different soil parent materials, crop species, specific climates, etc. In spite of this, the carbon and nitrogen, as the two mainly important energy materials of soil microorganisms, is vital for shaping bacterial communities, which still were broadly supported by previous studies 15,18,47 .

Conclusions
Long-term various fertilization regimes shaped different bacterial community structures in the rhizosphere of reddish paddy soil. Interestingly, balanced using of manure and chemical fertilizers not only increased the organic material pools and nutrient availability but also enhanced the biodiversity of the rhizospheric bacterial community. Moreover, balanced fertilization regime significantly increased the abundance of Actinobacteria, which involve in organic matter decomposition and the humus formation process, more importantly, it can produce various antibiotics to antagonize pathogenic microbes. Above advantages contribute to the sustainable development of agro-ecosystems. Additionally, CCA and the Mantel test clarified that E-Mg 2+ , SOC and AN are the key driving factors for shaping bacterial communities in the reddish paddy rhizospheric soil of Central South China.

Materials and Methods
Natural profile of experimental site. The long-term fertilization experiment was set up in 1982 at the Red Soil Experimental Station of the Chinese Academy of Agricultural Sciences, Qiyang County, Hunan Province, China (latitude, 26°45′N; longitude, 111°52′E; elevation, 120 m), which has a typical subtropical monsoon climate, with an average annual temperature of 18 °C, annual precipitation of 1250 mm, annual evaporation of 1470 mm, and 300 frost-free days. The paddy soil is derived from Quaternary red clay, the experiment mimicked the typical early rice-late rice double cropping system which is widely used in the hilly area of central south China, and the equivalent fertilization treatments were carried out in early and late rice, respectively. The initial properties of the plough layer (0∼20 cm) soil are listed as follows: pH, 5.97; soil organic matter, 1.98%; total nitrogen, 1.5 g kg −1 ; alkali-hydrolyzable nitrogen, 158.0 mg kg −1 ; total phosphorus, 0.5 g kg −1 ; available phosphorus, 9.6 mg kg −1 ; total potassium 14.2 g kg −1 ; and available potassium, 65.9 mg kg −1 .
Experimental design. The completely randomized block experiment included four treatments: no fertilizer application (NF); application of chemical NPK fertilizers (NPK); application of composted manure (M); and the combined application of composted manure and chemical NPK fertilizers (NPKM). Each treatment was replicated three times for a total of twelve plots, and each experimental plot was 27 m 2 (1.8 m × 15 m). These plots were managed with separate irrigation inlets and drainage outlets. Chemical NPK fertilizers were applied as urea (46% N), calcium superphosphate (12% P 2 O 5 ) and potassium chloride (60% K 2 O); manure was cow dung from the farm nearby and was composted completely before application, the muti-year average N, P 2 O 5 and K 2 O content of which was 0.32%, 0.25% and 0.15%, respectively. All these fertilizers were used as basal fertilizers, and the application rates of the treatment plots are listed in Table 4.

Soil sampling and physicochemical analysis.
To minimize the interference of recent fertilization, maximize the differences among different fertilization regimes and avoid interference from harvesting, sampling was performed prior to late rice harvesting in October 2016. A total of five plants were excavated randomly in each plot, and the rhizospheric soil was collected from each of the 5 plants according to the description by Han et al. 48 and then pooled to form a composite sample. All samples were sealed in sterile plastic bags, packed on ice and transported to the laboratory immediately, and they were then sieved through 2-mm meshes and thoroughly homogenized after removing plant residues and gravel. Each sample was then divided into two parts: one was air-dried for analysis of soil properties, and the other was stored at −80 °C for DNA extraction. The selected soil properties are listed in Table 2 and were measured following the description by Bao 49 . Briefly, soil pH was determined using a pH meter (FE20-FiveEasy TM pH, Mettler Toledo, Germany) with a 1:5 ratio of soil to distilled water. Soil organic carbon (SOC) was determined using the vitriol acid-potassium dichromate oxidation method. Alkali-hydrolyzable N (AN) was measured by the alkali-hydrolysis and diffusion method. Available P (AP) was extracted with sodium bicarbonate and determined using the molybdenum-blue spectrophotometric method. Available K (AK) was extracted with ammonium acetate, and concentrations were determined using a flame photometer. Exchangeable Ca 2+ (E-Ca 2+ ) and Mg 2+ (E-Mg 2+ ) were extracted with ammonium acetate, and concentrations were determined using an atomic absorption spectrophotometer. DNA extraction, PCR and sequencing. For each sample, total genomic DNA was extracted from 0.5 g of fresh soil using the PowerSoil DNA Isolation Kit (MoBio Laboratories Inc., USA) according to the manufacturer's instructions, and three successive DNA extractions of each sample were pooled before PCR to minimize DNA extraction bias. The DNA concentration and purity were monitored on 1% agarose gels, which were then stored at −40 °C for subsequent analysis.
An aliquot (50 ng) of DNA from each sample was used as a template for amplification 50 , and the V4-V5 hypervariable regions of bacterial 16S rRNAs (Escherichia coli positions 515-907) 51 were amplified using the selected primers 515 F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 907 R (5′-CCGTCAATTCCTTTGAGTTT-3′), which could minimize the overestimation of bacterial diversity due to intragenomic heterogeneity within 16S rRNA genes 52 . The forward and reverse primers were tagged with adapter, pad and linker sequences. Each barcode sequence was added to the reverse primer for the pooling of multiple samples in one run of HiSeq sequencing. All primers were synthesized by Sangon Biotech. co., Ltd. (Shanghai, China).
PCR amplification of each sample was performed in triplicate using a Gene Amp PCR-System 9700 (Applied Biosystems, Foster City, CA, USA) in a total volume of 30 μl, which contained 15 μl of Master Mix (2×), 0.5 μM the final concentrations of the forward and reverse primers, 10 ng of template DNA and double-distilled water up to 30 μl. The PCR was executed under the following conditions: initial denaturation at 94 °C for 5 min, followed by 30 cycles at 94 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 30 s, with a final extension at 72 °C for 10 min. The triplicate PCR products were mixed equally after being purified with a Qiagen Gel Extraction kit (Qiagen, Germany), and then they were sent to Novogene Bioinformatics Technology Co., Ltd. (Beijing, China) for further sequencing based on the Illumina HiSeq2500 platform. The complete sequencing data were submitted to the NCBI Short Read Archive under Accession No. SRP135963.
Processing of sequencing data. The raw sequences, after the barcodes and primers were trimmed, were assembled for each sample according to the unique barcode using QIIME 53 . The split sequences of each sample  were merged using FLASH 54 , and low-quality sequences were discarded using QIIME. The retained sequences were analysed according to the UPARSE pipeline, using USEARCH and Perl scripts to generate an OTU table and to pick representative sequences 55 . Briefly, sequences with a quality score lower than 0.5, a length shorter than 200 nt or singletons were discarded, and the retained sequences were then assigned to operational taxonomic units (OTUs) based on 97% similarity after filtering chimeras. The representative sequences of each OTU were aligned and classified using the SILVA reference database for bacterial sequences 56 .
Bioinformatic and statistical analysis. Rarefaction was performed to compare the relative levels of OTU richness across all soil samples at an OTU cut-off of 0.03. To correct the sampling effects, a randomly selected subset of 15,494 sequences for each sample were chosen for further bacterial community analysis. The OTU table was transformed into a suitable input file for further alpha and beta analyses using MOTHUR 57 . Species richness and biodiversity were estimated by the Chao1 estimator (Chao1), Shannon diversity index (Shannon) and Good's nonparametric coverage (Coverage). Principal coordinate analysis (PCoA) and cluster analysis (unweighted pair group method using arithmetic averages, UPGMA) based on the calculated weighted UniFrac distance were used to explore the differences in the bacterial community structure among soil samples. Based on the Kruskal-Wallis (KW) sum-rank testing, linear discriminant analysis effect size (LEfSe) was performed to identify significantly different species (biomarkers) of bacterial taxa among groups, and then linear discriminant analysis (LDA) was performed to estimate the effect size of each biomarker 58 . Canonical correspondence analysis (CCA) and the Mantel test were performed to analyse the relationships between soil properties and dominant bacterial phyla, followed by 999 and 9,999 permutations to test the significance, respectively. An analysis of variance was performed using IBM SPSS 20.0 (SPSS Inc., USA).