Comprehensive assessment of paddy soil quality under land consolidation: a novel perspective of microbiology

Soil physical and chemical indicators

The soil pH was measured at 1:2.5 (soil:water) by pH meter and the soil water content (SWC) was determined by oven drying at 105 °C for 6 h (Baldoncini et al., 2019; De Paul Obade, 2019). The concentrations of Cu, Zn, Cd, Pb, Ni, Cr, Hg and As in soils were tested by an inductively coupled plasma source mass spectrometer (Agilent 7800; Agilent Technologies, Palo Alto, CA, USA) after being extracted (Wu et al., 2019b). The soil OM, total nitrogen (TN), available phosphorus (AP) and available kalium (AK) were measured by total organic carbon analyzer (BOCS301; Shimazu Enterprise Management (Kyoto) co., LTD, Shimadzu, Japan), automatic kieldahl apparatus (K9860; Shandong Haineng Scientific Instrument co. LTD, Jinan, China), spectrophotometer and flame photometer (Song et al., 2018), respectively.

Microbial diversity

The soil genomic DNA for the PCR amplification was extracted from 0.5 g of triplicate soil samples using the FastDNA SPIN kit (MP Biomedicals, Irvine, CA, USA) according to the manufacturer’s instructions. The quality and concentration of extracted DNA were determined by using the Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA), and the DNA was stored at −20 °C for further use (Ma et al., 2018). The V3–V4 hypervariable regions of bacterial 16SrRNA gene were amplified with the primer pairs 338F and 806R, and the Internal Transcribed Spacer (ITS) regions of the fungal DNA gene were amplified using the primer pairs ITS1F and IT2R (Caban et al., 2018; Nottingham et al., 2018). PCR of bacteria and fungi were all performed to amplify 10 ng of template DNA in a 20 μl reaction system containing 0.8 μl (five μM) of each primer. The PCR reactions were performed by ABI GeneAmp® 9700 with the following program: 95 °C for 3 min; 10 cycles at 95 °C for 30 s, 60 °C for 30 s and 72 °C for 45 s, then a final extension at 72 °C for 10 min and 10 °C for 30 min. Sequencing data were performed on Illumina HiSeq4000 and the raw DNA sequencing data were processed with Quantitative Insights Into Microbial Ecology 2 (QIIME 2) (Jiang et al., 2019; Nottingham et al., 2018). Sequences were clustered into operational taxonomic units (OTUs) at 97% similarity, and the most abundant sequence was selected from each OTU to represent the respective OTU. Finally, the representative sequences were deposited into the NCBI SRA database (accession number: PRJNA532482, PRJNA532470).

The indexes of microbial alpha-diversity were estimated by mothur package (version v.1.30.1), including Shannon and Chao 1 (Barnes et al., 2017; Shao et al., 2019). The student’s test was used to estimate significant differences in diversity indexes between the different groups. And the beta-diversity analyses, including hierarchical clustering, principal component analysis, principal co-ordinates analysis and redundancy analysis, were calculated by QIIME 2 and software R (version 2.1.3) (Jiang et al., 2019).

Furthermore, some bacteria, among the beneficial microbes in the soil, can produce indole acetic acid to increase chlorophyll content, and as a result, increase the photosynthesis rate in the plant. Some other microbes can generate biologically active substances, such as hormones and enzymes, to control soil diseases and accelerate the circulation of nutrients in the soil, thereby improving the crop growth (Antoun, 2012). Therefore, the functions of the beneficial microorganisms play a crucial role in assessing soil quality in the paddy fields. The bacteria such as Azospirillum, Bacillus, Bradyrhizobium, Escherichia-Shigella, Mesorhizobium, Pseudomonas, Rhizobium, Rhodopseudomonas and Streptomyces, as well as fungi such as Glomerales, Mucorales, Orbiliales, Tremellales and Hypocreales are the beneficial microbes in the soil according to the previous researches (Antoun, 2012; Chen et al., 2012; Hayat et al., 2010).

Functional diversity of the bacterial community

16S rRNA gene is a key tool for studying the functional capabilities of microbial communities (Langille et al., 2013). It is necessary to figure out how the information about the soil bacteria helps us predict the extent of effects that soil bacteria may have on farmland ecosystems (Fierer, 2017). The 16S function prediction is a method to standardize the OTU abundance table by the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) software which stores the COG information and KEGG Ortholog (KO) information related to the greengene id. Then, the COG family information and KO information relevant to the OUT were obtained by the greengene id related to each OUT, and the abundance and KO abundance of each COG was calculated. According to the information of the COG database, the description and function of each COG can be analyzed from the eggNOG database, thereby obtaining a functional abundance spectrum. The information of KO, Pathway and EC can be obtained from the KEGG database, then the abundance of each functional category can be calculated according to the OTU abundance. In addition, for Pathway, PICRUSt can be used to obtain three levels of information on metabolic pathways and obtain abundance tables for each level (Malik et al., 2018).

Microbial biomass

The soil microbial biomass carbon (MBC), microbial biomass nitrogen (MBN) and microbial biomass phosphorus (MBP) were measured by chloroform fumigation extraction based on moist paddy soil put forward by Vance (Li et al., 2018a). Briefly, two moist soil samples (five g dry weight equivalent) were inoculated in 150 ml centrifuge tubes. Then, one sample was fumigated with ethanol-free CHCl₃ for 24 h. Both fumigated and non-fumigated soils were extracted with a 20 ml 0.5M K₂SO₄ solution by 30 min horizontal shaking and filtering. The organic carbon and TN in the extracts were determined by a C/N automatic analyzer for catalytic high temperature combustion. Then MBC, MBN and MBP were determined by using conversion factors of 0.45, 0.54 and 0.4, respectively (Murugan et al., 2019).

Minimum data set

The MDS was considered to be a reliable and measurable set of soil assessment indicators that eliminate redundancy and reduce information loss (Wang et al., 2018a).

First, statistical characteristics, variabilities of each indicator and correlations between different indicators are evaluated, then the results are used as a reference to select indicators. Second, the indicators are grouped by PCA in case of data redundancy, and only components with eigenvalues ≥1 are selected. Third, the vector norm value for each indicator is calculated to filter the indicators. Because the vector norm value represents the combined loading of one indicator in all components. The formula is as follows: (1) $$N_{\text{αβ}}=\sqrt{{\displaystyle \sum }_{\text{α}=1}^{\text{β}}\left({\text{λ}}_{\text{β}}{\text{μ}}_{\text{αβ}}^2\right)}$$μ_αβ, λ_β are the loading value and eigenvalue in component β, respectively.

Finally, normal linear transformation is conducted on the norm value and correlation of each indicator to calculate the sum. Indicators whose values are greater than 90% of the maximum are selected for MDS.

Calculation method

Comprehensive assessment of paddy soil quality

Multi-standard quantitation procedures were applied to the comprehensive assessment of PSQ, including soil physicochemical properties, biological indicators, and heavy metal contamination. The indicators of soil physicochemical and biological properties can be integrated into the Eq. (6) based on the results of f(x) in Eqs. (2) and (3). The formula is as follows: (6) $$\text{PB}={\displaystyle \sum }_{i=1}^{n}f{\left(x\right)}_i\times W_i$$

PB is the integrated results of soil physicochemical and biological properties; and f(x)_i and W_i are the score and weight of the indicator i in MDS, respectively. In this study, the weights of the indicators were determined by the common factors from PCA.

Meanwhile, the levels of heavy metal contamination should be changed from Eq. (4) as follows: (7) $$\text{HM}=\{\begin{array}{ll}1,\hfill & \left({\text{GI}}_i<0\right)\hfill \\ 1/\left({\text{GI}}_i+1\right),\hfill & \left({\text{GI}}_i\ge 0\right)\hfill \end{array}$$

HM represents the levels of the standardized environment index, while GI_i is the GI from Eq. (4).

Therefore, the PSQ is illustrated by the results mentioned above, and the formula is as follows: (8) $$\text{PSQ}=\text{MIN}\left\{\sqrt{\frac{\left({\text{PB}}_{\text{ave}}^2+{\text{PB}}_{\text{min}}^2\right)}{2}},\sqrt{\frac{\left({\text{HM}}_{\text{ave}}^2+{\text{HM}}_{\text{min}}^2\right)}{2}}\right\}$$

PB_ave represents the average PB value and PB_min represents the minimum BP values; while HM_ave is the average HM value and HM_min is the minimum HM values.

Soil physiochemical properties

There were 35 fields under land consolidation (LM, LL, DC and AO), and 10 other fields under CT from non-land consolidation area. All the 45 fields had similar paddy management practices and the soils had analogous textures too. The soil physiochemical properties of the samples from the study area are listed in Table 1.

Soil heavy metals

The heavy metal concentrations of the paddy soils and the reference background values of this area are listed in Table 2. On average, the Zn was the main metal, accounting for about 42%, and was ranked by concentrations as follows: Zn > Cr > As > Ni > Cu > Pb > Cd. The concentrations of Cu, Zn, Cd, Ni and As were higher than the background values at all sites, while the concentrations of Pb were lower than the background values at some sites under land consolidation practices of LM and LL. The concentrations of Cr in soils from fields under land consolidation practices of DC and AO were much lower than other fields. The results also suggested the paddy soils in the fields under land consolidation were less polluted by the heavy metals than fields under CT.

Microbial biomass

The concentrations of MBP, MBC and MBN showed an increase after land consolidation, with different land consolidation practices contributing more than 18% of value growth. The concentrations of MBP, MBC and MBN were the highest under the practice of AO, and the concentrations of MBP and MBC were significantly higher than other practices (Table 3).

Microbial diversity

The student’s test revealed that the main factor which significantly affected the diversity of bacterial communities was land consolidation, including LM, LL, DC and AO (p < 0.01). The microbial community richness indices such as Chao 1 and the diversity indices like Shannon were significantly higher in the soils collected from land consolidation area compared with non-land consolidation areas (p < 0.01) (Fig. 1). Meanwhile, most of the alpha-diversity indices were also significantly higher under the AO (p < 0.01).

Similar to the bacterial communities, most of the alpha-diversity indices of the fungal communities were significantly higher in the soil samples from land consolidation area compared with non-land consolidation area (Student’s t-test, p < 0.01). The community richness index—Chao 1, and the community diversity index—Shannon, were significantly higher in soils from land consolidation area (p < 0.01). Meanwhile, the AO also significantly increased the alpha-diversity indices of the fungal communities, except for Shannon. Most of the alpha-diversity indices were also significantly higher in the large-scale field than in small-scale field.

A total of 2,048,688 valid sequences of 16S rDNA were clustered into 11,244 OTUs, which were assigned to 56 phyla and 1,100 genera for bacteria. The OTUs allocated above the phylum level and genus level were approximately 99% and 35%, respectively. Proteobacteria (32.11% of OTUs), Chloroflexi (19.04% of OTUs), Acidobacteria (14.50% of OTUs), Actinobacteria (8.14% of OTUs), Bacteroidetes (5.18% of OTUs) and Nitrospirae (4.31% of OTUs) represented more than 80% of the described bacterial phyla. At genus level, Nitrospira, Sphingomonas, H16, Roseiflexus, Geobacter, Thiobacillus, Haliangium, Anaeromyxobacter, Gaiella and Candidatus_Solibacter were the most abundant genera (Fig. 2).

A total of 2,573,968 valid sequences of ITS were clustered into 5,377 OTUs and assigned to seven phyla and 457 genera for fungi. In addition to unclassified fungi, approximately 95% of the OTUs were assigned to Ascomycota (76.75%), Basidiomycota (13.17%) and Rozellomycota (5.73%). The predominant fungi genera were Gibberella (3.44% of OTUs), Schizothecium (2.92% of OTUs), Pseudeurotium (2.72% of OTUs), Zopfiella (2.65% of OTUs), Fusarium (2.05% of OTUs), Bolbitius (1.92% of OTUs), Guehomyces (1.86% of OTUs), Acremonium (1.41% of OTUs), Mortierella (1.26% of OTUs) and Curvularia (1.07% of OTUs) (Fig. 2).

The relative abundance of beneficial microbes in the soils from land consolidation area was much higher than those from non-land consolidation area (Fig. 3).

Bacterial functional diversity

The 16S rRNA gene is a key tool for studying the functional capabilities of microbial communities (Langille et al., 2013). We obtained all COG and KO annotations by PICRUSt to generate a table of abundances about KO and COG for the soil bacterial genomes that have identifiers in the Greengenes reference tree (Wong et al., 2016). According to the functional profiles of soil samples, amino acid transport and metabolism, signal transduction mechanisms, energy production and conversion and cell wall/membrane/envelope biogenesis were the most abundant functions, especially in the soils under land consolidation (Fig. 4). This suggested that the functions of nitrogen cycling, denitrification and respiration in soil bacteria became more effective after land consolidation (Lopez et al., 2018). Nonetheless, the differences of the functions in the soil bacteria between the soils from land consolidation area were not significant (p < 0.05), which also led to the significant differences of the corresponding bacterial community compositional structures.

Minimum data set

There were five principal components with eigenvalues >1.0 and a cumulative variance of 73.17% (Table 4). The large eigenvector values strongly correlated with each of the five primary PCs were weighted according to the percentage variance demonstrated by the particular PC. The variables with the largest absolute eigenvector values for PC1 were OM and TN. For PC2, the largest absolute eigenvector values were for AK, SB and CB. Then, MBC, BFD, SF and pH were the indicators with the largest absolute eigenvector values in the rest of principal components. The weighting values for these indicators were assigned according to the percentage variance demonstrated by the particular PC. The 25.25% from PC1 was distributed equally between OM and TN, and 18.13% from PC2 was also applied to AK, SB and CB, then the remaining percentage variances were also distributed in the same way. The weighting values were indicated by common factors through dividing each weighting value by the total weighting value of 0.7317. The final calculation of PB formula was showed in Eq. (9). The value of PB and the individual contribution of MDS indicators under different land consolidation practices were exhibited in Fig. 5.

Comprehensive soil quality

The GI was used to assess the degree of heavy metal contamination. Among the eight heavy metals, Cd showed the highest GI values, indicating high level of pollution of Cd. In contrast, the GI values of Pb and As were both under zero, meaning that there was no contamination of Pb and As in the study area. In general, the GI values followed the order of Cd > Cr >Hg >Ni > Cu >Zn > Pb >As. At the same time, there were significant differences in the GI values of heavy metals in different land consolidation practices, except for As and Hg (Fig. 6A). All soils from non-land consolidation area (CT) have the highest GI values (Fig. 6A) and PA values (Fig. 6B), indicating that land consolidation can effectively reduce the pollution level of heavy metals.

Meanwhile, all soils samples from the study area belonged to low soil comprehensive quality level according to the results of comprehensive assessment framework (Fig. 7A). However, the PSQ values of soils under land consolidation were significantly higher than those from non-land consolidation area (Fig. 7B). The soils from CT showed the highest values of PA and the relatively low values of PB, resulting in low level of soil quality. The linear fitting in Fig. 7A was performed to test the relationship among the values of PB, PA and PSQ, and the results showed that the increased values of PB caused soil quality degradation.