INTRODUCTION
Cold waterlogged paddy fields (CWPF) are a main type of low-yield paddy soils which are mainly distributed in mountain valleys, low-lying lands of hills, plains and lakes, and downstream of dams in southern China. CWPF have 3 460 000 ha which occupy 15.07% of the total paddy area and 44.2% of the low-yield paddy area in China. Cold waterlogged paddy soils are formed when the groundwater level rises due to the impact of local microclimate and topographical features resulting in water accumulation in the soil surface (Xie et al., 2015). The low productivity of CWPF is attributed to high groundwater table, poor drainage conditions, low soil temperature, poor soil physical characteristics, excessive reducing substances and low contents of available nutrients. With the increasing food demand in China, improvement of productivity of low-yield paddy field is becoming increasingly important. Organic matters in CWPF are easily subjected to accumulate rather than decompose because of reducing conditions. Consequently, CWPF have higher inherent fertility potential. Excess water in field could be rapidly removed after drainage, which could improve the soil aeration condition and accelerate the oxidation of organic matter. Therefore, CWPF can increase their yield through proper drainage. Chinese farmers have employed some techniques for land drainage, such as drainage ditch, mole drainage and tile drainage. Agricultural drainage ditches are essential for removal of surface water and groundwater to allow for crop production in poorly drained agricultural landscapes (Needelman et al., 2007).
Soil microorganisms are important components of agricultural ecosystems as they mediate many processes including cycling of soil nutrients, decomposition of organic matter and formation of soil aggregate (Kamaa et al., 2011). Soil microbial communities are dynamic assemblages of populations with different strategies for responding to changing environmental conditions (DeAngelis et al., 2010). It is well known that agricultural management practices (e.g., fertilization regimes, crop rotations, irrigation pattern) largely determine the structures of microbial communities (Li et al., 2012; Yuan et al., 2013; Zhao et al., 2014). Drainage is an important water management strategy in agriculture for improving the soil physical properties, which could result in shifting the microbial communities. However, relatively few researchers have studied the effect of drainage ditches on the paddy soil microorganisms. Understanding the soil microbial community and its response to the drainage management is important for selecting suitable management practice to improve soil function and increase rice yield. Therefore, experiments examining the effect of the drainage ditch on the microbial communities should help us to determine the suitable drainage practice. In the present study, the composition and structure of soil bacterial and fungal communities were examined in a cold waterlogged paddy soils at the different distances from the drainage ditch on the Illumina Miseq platform (Illumina, San Diego, California, USA).
The objectives of this study were to determine the abundance, taxonomic diversity and compositions of bacterial and fungal community at the different distances; and identify the environmental variables that correlated to changes in the structures of microbial communities.
MATERIALS AND METHODS
Experimental design and soil sampling
The experiment was established in 1979 at Shunchang county (26°42' N, 117°42' E), Fujian Province, China. This region has a subtropical monsoon humid climate with an annual average temperature of 18.5 °C and annual precipitation of 1691 mm. The soil is classified as Gleyi-Stagnic Anthrosols (USDA). To remove excess water in field, a drainage ditch (200 m length, 30 cm width and 1 m depth) was established in the middle of a paddy field. Three distances were selected for soil sampling, which were located at 5, 15 and 25 m from the drainage ditch, parallel to the drainage ditch and marked as "A", "B" and "C", respectively (Figure 1, minor revision from Qiu et al., 2013). Surface soil samples (0-20 cm) were collected in August 2014 (tillering stage during rice cultivation). The amounts of chemical fertilizers were 120 kg N ha-1, 48 kg P ha-1, 84 kg K ha-1 per year in the form of urea, superphosphate and potassium chloride, respectively. Three sampling sites were randomly selected at each distance. At each site, five soil cores (5 cm diameter) were randomly collected and composited. Similarly, two other composite samples at different sites from the same distance were obtained. Therefore, nine composite samples were obtained. The fresh samples were sealed in sterilized plastic bags and transported on ice to the laboratory. Each soil sample was partitioned into three subsamples: one stored at -20 °C for DNA extraction, the second was air-dried and passed through a 2-mm sieve for chemical analysis and the remaining subsample was stored at 4 °C for determining reducing substances.
Soil chemical analysis and soil DNA extraction
Soil organic matter (OM) was determined by the K2Cr2O7 oxidation method. Available N (AN) and available P (AP) were determined by NaOH hydrolyzable and molybdenum methods, respectively. Available K (AK) was determined using atomic absorption spectrometer (AAS, SP-3801, Shanghai-Spectrum, Shanghai, China). Soil redox potential (Eh) was measured at a depth of 5 cm in the field using Eh indicator Spectrum IQ 150 (Spectrum Technologies Inc., Aurora, Illinois, USA). Total reducing substances (TRS), Fe2+ and Mn2+ were measured using the methods described by (Lu (2000).
Soil total DNA was extracted by the SDS-based method (Zhou et al., 1996) and then purified by DNA Fragment Purification Kit (TaKaRa, Dalian, China). DNA concentration and quality were determined using a NanoDrop 2000 (NanoDrop Technologies Inc., Wilmington, Delaware, USA).
Real-time PCR and MiSeq sequencing
Abundances of bacteria and fungi were performed using primers 338F/518R and FF390/FR1, respectively (Prevost-Boure et al., 2011). The PCR mixture was prepared in a total 20 (L using 10 (L SYBR Green I PCR mix (TakaRa, Dalian, China), 0.2 (M of each primer and 2 (L 10-fold diluted extracted DNA. Real-time PCR was carried out in triplicate using an ABI 7500 Real-time PCR Systems (Thermo Fisher Scientific, Waltham, Massachusetts, USA) under the following thermal cycle conditions: 30 s at 95 °C, followed by 40 cycles of 5 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C. The amplification specificity was confirmed by melting curve. Standard curves ranging from 1 × 104 to 1 × 109 copies were prepared by 10-fold serial dilution of known copy numbers of plasmid DNA possessing the genes of interest.
DNA extracted from each soil sample served as a template in 16S rRNA gene and internal transcribed spacer (ITS) region amplification. Bacterial primers 515F (5'-GTGCCAGCMGCCGCGGTAA-3') and 806R (5'-GGACTACHVGGGTWTCTAAT-3') were employed to amplify the V4 hypervariable region of bacterial 16S rRNA gene (Bergmann et al., 2011). The reverse primer was designed to have a 6-bp error-correcting barcodes unique to each sample for the subsequent identification. DNA was amplified following the protocol described by Caporaso et al. (2011). Fungi-specific primers ITS5-1737F (5'-GGAAGTAAAAGTCGTAACAAGG-3') and ITS2-2043R (5'-GCTGCGTTCTTCATCGATGC-3') were used to amplify fungal ITS region (White et al., 1990). Amplicon sequencing was carried out on the Illumina Miseq platform at the Novogene Technology Co. (Beijing, China). Pairs of reads from the original DNA fragments were merged based on the method by Magoc and Salzberg (2011). Sequencing reads was assigned to each sample according to the unique barcode of each sample. Sequences were analyzed with QIIME software package (Quantitative Insights Into Microbial Ecology) and UPARSE pipeline (Caporaso et al., 2010). The reads were first filtered by QIIME quality filters. Then UPARSE pipeline was used to pick operational taxonomic units (OTUs) at 97% similarity. For each OTU, a representative sequence was selected and used to assign taxonomic data using the RDP classifier (Wang et al., 2007). Bacterial and fungal sequences were submitted to Short Read Archive (SRA) under accession numbers SRA307922 and SRA310798, respectively.
Statistical analysis
The community richness index, diversity index, data processing, OTU-based analysis were performed using mothur (Schloss et al., 2009). To compare microbial community structures across all samples based on the OTU and identify environmental factors influencing microbial diversity, principal component analysis (PCA) and redundancy analysis (RDA) were conducted using the Canoco 4.5 (Wageningen University and Research, Wageningen, The Netherlands).
RESULTS
Variations in rice yield and soil characteristics across different distances
The effect of drainage ditch on the paddy soil characteristics and rice yield are listed in Table 1. Rice yield was highly affected by the distance from the drainage ditch. The 5 m-site had the highest rice yield, about 10 160 kg ha-1. The lowest rice yield was observed at the 25 m-site with 8 460 kg ha-1.
With increasing distance from the drainage ditch, AN, AP, AK decreased (P < 0.05). Similar decreasing pattern was also observed in the soil Eh. The minimum TRS, Fe2+ and Mn2+ were observed at the 5 m-site. However, no significant difference in soil OM was observed across all samples. Soil chemical measurements showed that opening a drainage ditch in the cold waterlogged paddy field could improve soil nutrient contents, but the effectiveness of opening a drainage ditch gradually weakened with increasing distance from the drainage ditch. Rice yield was positively correlated with the soil Eh (r = 0.899**), AP (r = 0.723*) and AK (r = 0.854**) and negatively correlated with the total reducing substances (r = -0.891**) and Fe2+ content (r = -0.695*).
Differentiations in the bacterial and fungal abundance
The bacterial 16S rRNA gene abundance ranged from 0.46 × 1010 to 3.34 × 1010 copies g-1 soil (Figure 2a). The highest and lowest bacterial 16S rRNA gene abundances were observed at the 15 m- and 5 m-site, respectively. Fungal 18S rRNA gene abundance rang ed from 1.75 × 109 to 6.12 × 109 copies g-1 soil across different distances (Figure 2b) and 5 m-site had the highest value. The copy number of bacterial 16S rRNA gene was higher than that of fungal 18S rRNA gene at the same distance. Spearman's analysis revealed positive correlations between the bacterial abundance with concentrations of Fe2+ (r = 0.833**), Mn2+ (r = 0.667*). Moreover, there was a positive correlation between fungal abundance and AK (r = 0.683*), and a negative correlation between fungal abundance and TRS was also observed (r = -0.695*).

Vertical bars correspond to standard deviations (n = 3). Different lower-case and uppercase letters indicate significant differences at P < 0.05 and P < 0.01, respectively.
Figure 2 Abundance of soil bacterial 16S rRNA gene and fungal 18S rRNA gene across different distances from the drainage ditch.
Analysis of bacterial communities across different distances
Bacterial diversity and richness are shown in Table 2. The total number of sequences per sample ranged from 40 471 to 63 721 among the treatments. There was a significant difference in the number of OTUs among three distances (P < 0.01). The number of OTUs strikingly decreased with increasing distance from the drainage ditch. Chaol index at the 5 m- and 15 m- sites were higher than the 25 m-site (P < 0.05). There was no significant difference in the Shannon index among all samples.
All sequences were classified into 57 phyla by the UPARSE program. The overall bacterial compositions among nine samples were similar, while the distribution of each phylum varied (Figure 3a). In all samples, Proteobacteria was the most abundant phylum. Acidobacteria, Planctomycetes, Chloroflexi, Crenarchaeota, Verrucomicrobia, Actinobacteria, Nitrospirae, Gemmatimonadetes, Chlorobi, Bacteroidetes, WS3, NC10 were the dominant phyla with the relative abundances above 1% in at least one treatment. The relative abundance of phylum Acidobacteria were higher at the 5 m- and 15 m-sites (20.57% and 19.56%) compared to 25 m-site (15.22%), while phylum Bacteroidetes exhibited the opposite pattern, being higher at the 25 m-site (1.69%) than 5 m- and 15 m- sites (0.80% and 0.96%) (P < 0.05). Phylum Planctomycetes was also higher at the 5 m- and 15 m- sites (10.34% and 8.20%) than 25 m-site (3.99%) (P < 0.01). A higher abundance of phylum WS3 was observed at the 5 m-site (1.98%) than 15 m- and 25 m- sites (1.40% and 1.11%) (P < 0.05). There were no significant differences in the relative abundances of other phyla among all samples. In particular, a higher relative abundance of anaerobic bacteria were observed at the 25 m-site than 5 m- and 15 m- sites (Table 3). For example, one anaerobic phylum Spirochaetes was strikingly higher at the 25 m-site than 5 m- and 15 m- sites, although their relative abundances were low. Geobacter is a genus of Proteobacteria, an anaerobic respiration bacterial species. It was observed that the relative abundance of genus Geobacter was richest at the 25 m-site and lowest at the 5 m-site. Likewise, other anaerobic genera such as Syntrophobacter and Desulfobacterium showed a similar pattern. In contrast, relative abundance of aerobic bacteria such as genus Planctomyces was strikingly richer at the 5 m- and 15 m- sites than the 25 m-site.
To compare the microbiota across three different distances, PCA was conducted based on the bacterial OTUs (Figure 4a). The first two principal components could explain 69.6% of the total variation among all samples. Bacterial community compositions shift greatly between 5 m-, 15 m- and 25 m- sites along the first principal axis demonstrating the community compositions of 5 m-and 15 m- sites were relatively similar. The Bray-Curtis distance was used to compare the relatedness of bacterial communities. Except for B3 (the third replicate of 15 m-site) and C3 (the third replicate of 25 m-site), bacterial communities were partitioned into two clusters: one from treatments of 5 m- and 15 m- sites, and another from treatment of 25 m-site (Figure 5a). More importantly, the replicates of 5 m- and 15 m- sites overlapped each other implicating little shift of microbial communities between 5 m- and 15 m- sites.

a. Relative abundance of the dominant bacterial phyla; b. Relative abundance of the dominant fungal phyla. Vertical bars correspond to standard deviations (n = 3).
Figure 3 Bacterial and fungal community distribution in the surface soil across different distances.
Table 2 Operational taxonomic units (OTU) numbers and diversity indexes of bacteria and fungi across different distances.
Microbial community | Sampling site | Coverage (%) | OTUs | Chao1 | Shannon |
Bacteria | A (5 m) | 95.63 ± 0.26a | 2961 ± 18.71A | 3831.41 ± 169.90a | 9.94 ± 0.06a |
B (15 m) | 95.89 ± 0.16a | 2899 ± 31.75B | 3694.13 ± 160.27a | 9.89 ± 0.10a | |
C (25 m) | 96.33 ± 0.52a | 2807 ± 30.95C | 3354.69 ± 169.58b | 9.99 ± 0.09a | |
Fungi | A (5 m) | 98.77 ± 0.72a | 834 ± 331a | 859.85 ± 8.70ab | 5.16 ± 0.30a |
B (15 m) | 98.26 ± 0.00a | 1021 ± 76a | 928.65 ± 35.09a | 4.60 ± 0.36a | |
C (25 m) | 98.44 ± 0.14a | 955 ± 135a | 801.65 ± 54.51b | 4.58 ± 0.34a |
Data are the means ± SD, n = 3. Values with different lower-case and uppercase letters indicate significant differences at P < 0.05 and P < 0.01 according to LSD test, respectively.
Table 3 Relative abundances of selected bacterial taxa showing significant differences across different distances.
Taxa | A (5 m) | B (15 m) | C (25 m) | |
Spirochaetes (phylum) | 0.0021 ± 0.0001C | 0.0026 ± 0.0001B | 0.0032 ± 0.0000A | Anaerobe |
Planctomyces (genus) | 0.0270 ± 0.0062a | 0.0184 ± 0.0070a | 0.0061 ± 0.0021b | Obligate aerobe |
Geobacter (genus) | 0.0044 ± 0.0006b | 0.0059 ± 0.0004ab | 0.0070 ± 0.0014a | Anaerobe |
Syntrophobacter (genus) | 0.0027 ± 0.0003B | 0.0029 ± 0.0002B | 0.0041 ± 0.0002A | Anaerobe |
Desulfobacterium (genus) | 0.0015 ± 0.0003b | 0.0023 ± 0.0004b | 0.0033 ± 0.0006a | Anaerobe |
Data are the means ± SD, n = 3. Values in the same row with different lower-case and uppercase letters indicate significant differences at P < 0.05 and P < 0.01 according to LSD test, respectively.
Analysis of fungal communities across different distances
Fungal diversity indexes are shown in Table 2. The number of sequences per sample varied from 15 160 to 64 253. Chao1 at the 15 m-site had the highest value and the lowest value was observed at the 25 m-site. There were no significant differences in Shannon index across different distances.
The obtained fungal reads were affiliated to phyla Ascomycota, Zygomycota, Basidiomycota, Un-s-fungal sp WEF9, Chytridiomycota with an average proportion of 37.98%, 29.65%, 9.41%, 1.21%, 0.62%, respectively (Figure 3b). Phylum Zygomycota was dominant at the 5 m- and 15 m- sites while it was strikingly low at the 25 m-site. There were significant differences in the relative abundances of phyla Chytridiomycota and Un-s-fungal sp WEF9 across three different distances (P < 0.01), 25 m-site had the highest value and 5 m-site had the lowest value.
The first two principal components in the PCA targeted fungal OTUs could explain 68.6% of the total variation (Figure 4b). Similar to that of bacteria, a clear separation between 25 m- and 5 m-, 15 m- sites along the first axis could also be observed. Except for A3 (the third replicate of 5 m-site) and C3 (the third replicate of 25 m-site), fungal communities were partitioned into two clusters: one from treatments of 5 m- and 15 m- sites, and another from treatment of 25 m-site. This result was similar to that of bacterial community structure (Figure 5b).

A, B, C represent 5 m-, 15 m- and 25 m- sites, respectively.
Figure 5 Bacterial (a) and fungal (b) hierarchical cluster trees based on the distance matrix that was calculated using Bray-Curtis distance.
Relationships between the microbial communities and environmental factors
RDA was performed to study the effect of environmental variables on microbial abundant phyla. For bacteria, Eh (F = 4.780, P = 0.032) and TRS (F = 4.093, P = 0.048) were the most significant factors in explaining the total variation (Figure 6a). For fungi, AN (F = 4.078, P = 0.020), TRS (F = 3.571, P = 0.036) and Eh (F = 3.255, P = 0.044) were significantly influential (Figure 6b). From above results, it was concluded that Eh and TRS were the most important environmental variables in explaining the variance with respect to bacterial and fungal communities. Specifically, bacterial and fungal communities of the 25 m-site formed a separate group attributed to lower Eh and higher TRS.
DISCUSSION
Previous researches focused on the effect of agricultural management practices such as fertilization regimes, land use, irrigation patterns and water management on microbial community in agricultural land (Li et al., 2012; Chaudhry et al., 2012; Itoh et al., 2013; Dong et al., 2014), however, relatively little research on drainage ditch have been done. Drainage management is an important practice in cold waterlogged paddy fields. Previous research by our team showed that the average underground water level from March 2011 to April 2014 were -16.7, -5.4 and 8.3 cm at 5 m-, 15 m- and 25 m- sites, respectively (Wang et al., 2015). Water is a major factor in rice paddy ecosystems. Soil water content, which has direct and indirect effects on soil oxygen concentration and nutrient availability, is considered an important factor influencing soil microbial composition (Li et al., 2012). In present study, we found that opening a drainage ditch in waterlogged paddy soil has significant impacts on soil characteristics. With increasing distance from the drainage ditch, soil available nutrients and Eh decreased gradually and anaerobic processes such as metal reduction progressed sequentially, as also observed in previous study (Lueders and Friedrich, 2000). Previous research by our team (Wang et al., 2015) showed that 5 m-site is obviously under fluctuating dry-rewetting condition but the water level of 15 m- and 25 m- sites are mostly below the ground level. Soil drying reduces solute, enzyme mobility and the substrate supply to the decomposers (Or et al., 2007). Then the subsequent rewetting condition led to a cascade of responses including mobilizing C physically protected in aggregates, releasing intracellular osmolytes and enhancing metabolic and enzymatic activities, which in turn increase nutrient mineralization (Xiang et al., 2008; Borken and Matzner, 2009). Therefore, in our study, 5 m-site had the higher nutrients than 15 m- and 25 m-sites. The oxidation-reduction potential is associated with the most important soil properties and processes (Tokarz and Urban, 2015). Soil redox potential is of fundamental importance for agriculture as it plays an impact on crop yields. In present study, rice yield decreased sharply when Eh decreased which was also demonstrated by Husson et al. (2001).
Previously, various studies have examined how soil microbes respond to soil water content (Gorden et al., 2008; Xiang et al., 2008). In the present study, opening a drainage ditch in a cold waterlogged paddy soil was selected to study its effect on soil microorganisms. There was no correlation between bacterial abundance and the distance. The copy number of bacterial 16S rRNA gene in 15 m- and 25 m- sites were higher than in 5 m-site, but bacterial abundances were positively correlated with the contents of reducing substrates like Fe2+ and Mn2+. In contrast to bacteria, fungal abundances decreased sharply with increasing distance from the drainage ditch and were negatively correlated with the TRS. Above results demonstrated that fungal growth were easily inhibited under strongly reducing condition. Seo and DeLaune (2010) reported that fungi prefer weakly reducing conditions and exhibit the best development only at Eh > 250 mV.
The sequence analysis revealed that the majority of taxonomic groups detected were presented in all samples but with abundance-dependent shifts across different distances. We have found that among the dominant phyla, Acidobacteria, Planctomycetes, Bacteroidetes and WS3 were rather sensitive to opening a drainage ditch in a cold waterlogged paddy field. Moreover, relative abundances of aerobic bacteria were strikingly higher at 5 m- and 15 m-sites than 25 m-site, in contrast, 25 m-site had higher relative abundances of anaerobic bacteria. The relatively higher abundance of the known Fe reducers Geobacter at 25 m-site may result in much higher Fe2+ concentration. PCA based on soil microbial OTUs at 5 m-, 15 m- and 25 m-sites clearly reflect that microbial community structures at 5 m-and 15 m-sites were more similar as compared to 25 m-site. Excess water at 25 m-site replaces the gaseous phase in the soil pores, and the oxygen supplied from the atmosphere decreases; accordingly, microbial populations switch from aerobic to anaerobic bacteria.
Many studies have shown that environmental factors could shape microbial community structures (Brockett et al., 2012; Liu et al., 2014; 2015). The RDA analysis suggested that Eh and TRS were the common predominant factors affecting microbial community structures. These results were consistent with previous studies, which demonstrated that soil Eh was an important factor in shaping bacterial communities (Husson, 2013).

OM: Organic matter; AN: available N; AP: available P; AK: available K; TRS: total reducing substances; Eh: redox potential.
Figure 6 Redundancy analysis (RDA) illustrating the relationship between microbial abundant phyla and environmental variables: Bacterial community (a) and fungal community (b)
CONCLUSIONS
Our study has documented that opening a drainage ditch in a cold waterlogged paddy field is an efficient agricultural practice in improving soil nutrients and rice yield. The effects of drainage practice on the soil chemical properties, microbial abundance and community structures depend on the distance from the drainage ditch and the effectiveness in our study is about 15 m. Severe anaerobic conditions in soil could result in lower soil nutrients and rice yield.