The field experiment simulating climate change was established in 2010 in Kangbo village (31°30′N, 120°33′E), Changshu Municipality, Jiangsu, China. The area experiences a subtropical monsoon climate with an annual mean temperature (2004–2013) of 16°C and mean precipitation of 1100–1200 mm. The soil is a Gleyic Stagnic Anthrosol formed on a clayey lacustrine deposit which has been cultivated with a rice-wheat rotation for hundreds of years. The basic properties of the topsoil prior to the experiment in 2010 were: soil pH of 7.0, organic carbon of 1.6%, total nitrogen of 1.9 g kg^-1 and bulk density 12 g cm^-3.

The experimental details with elevated CO₂ and warming are presented in the study of Liu et al. (2014). In brief, with a block split-plot design, one field plot was artificially treated with a continuous atmospheric CO₂ concentration enrichment up to 500 ppm (CE) using a liquid CO₂ supply. The crop canopy air of another field plot was warmed by 2°C (WA) above the ambient temperature with infrared heaters. A further field plot was subjected to both CO₂ enrichment and warming (CW), and a final field plot was maintained at an ambient condition as a control (CK). Treatment levels for our investigation were defined according to IPCC (2007): the A2 emission scenario predicts atmospheric CO₂ concentrations to increase by 500 ppm and global mean air temperatures to increase by 2°C. For the elevated CO₂ treatment, pure CO₂ gas from a liquid tank was injected into the plots via perforated pipes surrounding the ring. Sixteen Li-820 CO₂ sensors (Li-COR Inc., Lincoln, NE, USA) were installed and evenly distributed above the canopy to automatically control CO₂ concentrations. CO₂ concentration consistency over the ring was controlled by automatic adjustment to wind direction and velocity; when weed speed was more than 5 m/s, or if it was raining, CO₂ spraying ceased. For the warming treatment, 12 infrared heaters (2000 W, 240 V, 1.65 m long × 0.14 m wide; HS-2420, Kalglo Electronics Co., Inc., Bethlehem, PA, USA) were situated on each ring. The heaters were adjusted every week to maintain a clearance height of 1.2 m above the top of the canopy during the growth period. The air temperature in the experimental plot was monitored by 6 infrared thermometers (Model SI- 121, Apogee instruments Inc., Logan, UT, USA) which were arranged in a hexagonal array (Liu et al., 2015; Cai et al., 2016; Wang et al., 2016). During the rice season, the average CO₂ concentration under elevated CO₂ plots was 515 ± 40 ppm and the increase of canopy air temperature under the warming plots was 1.98 ± 0.2°C. The control plots, surrounded by the same infrastructure, did not receive CO₂ enrichment or any warming treatments. Each treatment was replicated in three rings with the same infrastructure, having 8-m-diameter and covering 50 m² per ring. All the rings were buffered by an adjacent field to avoid treatment cross-over, and the distance between treatment plots was around 28 m.

In the experimental season, rice (Oryza sativa L. cv., Changyou No.5) was transplanted at a density of three seedlings per hill on 20th June, 2013. Plots were treated with local conventional practices, including a soil water regime of flooding during seedling to tillering stages, intermittent irrigation during heading, and drainage for ripening. Urea and ammonium bicarbonate were applied as basal fertilizers at a rate of 150 kg N ha^-2 (120 kg-N ha^-2 as urea and 30 kg-N ha^-2 as ammonium bicarbonate) on 21th June, 2013. Chlorpyrifos was applied as a pesticide at a rate of 800–1000 g ha^-2 at the heading stage. The management practices were consistent across all the treatments.

Rice rhizosphere soils were sampled at the tillering (19th July), heading (4th September) and ripening (24th October) stages in 2013. The rhizosphere of five individual rice plants were randomly collected at a depth of 0–15 cm from each plot, following the procedure described by Butler et al. (2003). The rhizosphere soil (being tightly adhered to the plant roots with about 1 cm thickness) was carefully removed and evenly mixed to form a composite sample. These soil samples were passed through a 2-mm sieve and immediately sealed in a plastic bag before being transferred to the laboratory (within 1 day after sampling). Fresh samples were stored at 4°C and analyzed for soil physico-chemical analyses within 1 week of sampling. A sub-sample of the soil was stored at –20°C prior to DNA extraction; this was undertaken within 1 week of sampling.

Soil microbial biomass carbon (SMBC) was determined by a fumigation-extraction method following Wu et al. (1990). The samples were fumigated with ethanol free chloroform for 24 h at 25°C before being extracted with 0.5 mM K₂SO₄ for 30 min on a shaker; unfumigated samples were also processed using the same method. The extracts were analyzed for extractable C using an automated TOC Analyzer (TOC-500, Japan). A K_EC of 0.45 was used to convert the measured C to SMBC values. Inorganic N (NH₄⁺ -N and NO₃^- -N) was extracted by shaking with 0.5 mol L^-1 K₂SO₄ (1: 5 (w/w) soil: K₂SO₄ solution) for 1 h and then filtering through a 0.45-um-pore-size polysulfone membrane, before colorimetric determination using an automated flow injection analyzer (Skalar Analytical B.V., The Netherlands).

The CH₄ production and oxidation potentials of soil were analyzed with a laboratory incubation method. CH₄ production potentials in the soil samples were determined following the methods of Singh et al. (2012). In summary, 15 g of sample was transferred into a 120 ml glass jar, amended with 25 ml of sterile distilled deionized water and sealed with a butyl rubber stopper. The headspace in the jar was flushed with N₂ for 10 min. Each soil sample was repeated in triplicate and incubated at 28°C in the dark. CH₄ production was analyzed periodically by gas chromatography (Agilent 4890D, USA) equipped with a flame-ionization detector (FID). The CH₄ concentration in the headspace was measured every 24 h for 1 week, and CH₄ production potential was calculated using a linear regression of increased CH₄ concentration with time.

CH₄ oxidation potential was measured following the protocol of Vishwakarma et al. (2010). In summary, 15 g of sample was transferred into gas-tight 120 ml glass jars and incubated at 28°C for 7 days in the dark. All samples were analyzed in triplicate. The headspace in the jars contained 5% v/v methane in air. The CH₄ concentration in the headspace was measured every 24 h for 1 week. CH₄ oxidation potentials were calculated from the initial linear reduction of CH₄ concentration with time and expressed as mg CH₄-C per hour per gram dry weight.

A static closed chamber-GS method was used to monitor CH₄ flux (Zou et al., 2009). Gas samples were collected once a week during the rice growing season. Samples were collected between 08:00 and 10:00 on the collection day and 4 individual gas samples were collected with a syringe at 0, 10, 20, and 30 min after chamber closure. The concentration of CH₄ in a sample was analyzed using a gas chromatograph (Agilent 7890A) equipped with a flame ionization detector (FID). The carrier gas was nitrogen and a flow rate of 40 ml/min was maintained. The oven and FID were operated at 50 and 300°C, respectively.

Total DNA was extracted from 0.35 g of fresh soil with a PowerSoil^TM DNA isolation kit (MoBio, Carlsbad, CA, USA) following the manufacturer’s instructions. DNA quality was assessed on an agarose gel while DNA quantity was determined using a Nanodrop spectrophotometer (Thermo Scientific, DA, USA).

Real-time PCR was performed in a 7500 real-time PCR system (Applied Biosystems, Germany) via fluorometric monitoring with SYBR Green 1 dye. The primer pair mcrAF/mcrAR (Luton et al., 2002) and A189f/mb661r (Horz et al., 2001) were used to quantify methanogenic archaeal mcrA genes and methanotrophic bacterial pmoA genes of all samples, respectively. Each reaction was performed in a 25 μl volume containing 15 ng of DNA, 1 μl of 10 μM of each primer and 12.5 μl of SYBR premix EX Taq^TM (Takara Shuzo, Shinga, Japan). A melting curve analysis was conducted following each assay to confirm specific amplification was not from primer-dimers or other artifacts. A single clone containing the target region was grown in Luria-Bertani media and plasmid DNA was extracted using a plasmid-extraction kit (Takara, Japan). Standard curves were generated using a 10-fold dilution of plasmid DNA, from 10³ to 10⁹ copies of the template. PCR efficiencies were obtained between 98% and 106%, with R² values >0.99. The final methanogenic mcrA gene and methanotrophic pmoA gene copy numbers were calibrated against total DNA amounts and soil water content.

Terminal-restriction fragment length polymorphism was used for analyzing the community structure of methanogens and methanotrophs. Briefly, the functional genes mcrA and pmoA were amplified by PCR using the primer pairs mcrAF/mcrAR and A189f/mb661r, as previously mentioned with the 5′ end of the mcrAF and A189f primers labeled with 6-carboxyfluorescein (6-FAM). All PCRs were performed in duplicate and pooled for subsequent restriction and T-RFLP analysis. PCR products were separated by 1.5% agarose gel, and purified using the PCR solution purification kit (Takara, Dalian, China). Purified PCR products were used in a restriction digest TaqI and MspI (Takara, Dalian, China) for mcrA and pmoA genes as per the manufacturer’s instructions, respectively. Fragment analysis was achieved by capillary electrophoresis (ABI 3100 Genetic Analyzer; Applied Biosystems, Carlsbad, CA, USA) using a GeneScan ROX-labeled GS500 internal size standard. T-RFLP patterns were analyzed using GeneMapper software (Applied Biosystems) by peak height integration of different terminal restriction fragments (T-RFs). The fluorescence intensity (%) represented by a single T-RF was calculated relative to the total fluorescence intensity of all T-RFs. Peaks with heights that were less than 2% of the total peak height were excluded from further analysis to avoid potential noise before calculating relative T-RF abundance.

Based on the obtained T-RFLP results, all soil samples at the ripening stage were chosen to establish clone libraries. Libraries for the functional genes mcrA and pmoA were created by ligating PCR products into pEASY-T3 vectors and being transformed into competent cells Escherichia coli JM109 (Takara, Japan) in accordance with the manufacturer’s instructions. Ninety four methanogenic clones and 102 methanotrophic clones were sequenced. All sequences were checked for chimera by using Bellerophon (Huber et al., 2004) before being grouped into operational taxonomic units (OTUs) using the furthest-neighbor clustering algorithm of the DOTUR software with a 96% threshold. In silico digests with TaqI and MspI were undertaken on the sequences to allow the assignment of phylogenetic identity to individual peaks. The closest relatives of each sequence were checked using a BLAST search within GenBank. The representative sequences recovered in this study have been deposited in the GenBank database under accession numbers KU133526-KU133543 (methanogenic mcrA genes) and KU133544-KU133564 (methanotrophic pmoA genes).

Statistical analysis was performed using SPSS 20.0. One-way ANOVA with Tukey’s HSD test was used to test the difference among the treatments at each growth stage. Repeated measures ANOVA were used to determine the effect of climate change factors and plant growth stage on soil properties, CH₄ production and oxidization rates, and microbial abundance (log₁₀-transformed mcrA and pmoA gene abundances). Principal component analysis (PCA) of the T-RFLP profiles was performed using Minitab v.15 software based on relative fluorescence intensity of T-RFs. The probability level p < 0.05 was considered to be statistically significant.

Soil physico-chemical property data of the rice field soil under the simulated climate change conditions are shown in Table 1. Results showed that soil inorganic N (NH₄⁺ and NO₃^-) did not significantly change under CE, CW, and WA treatments when compared with the control treatment. Soil NH₄⁺ content generally declined with rice growth development across the treatments, ranging from 30.23 mg kg^-1 (CW, tillering) to 7.88 mg kg^-1 (CW, ripening). However, the content of NO₃^- was stable without significant changes with the growth stages. Compared to CK, SMBC significantly increased under elevated CO₂ levels (CE and CW) at all three growth stages, but the WA treatment significantly increased SMBC only at the heading stage. Repeated measures ANOVA showed that elevated CO₂, warming and their combination significantly affected SMBC (p < 0.05), however, the interaction with the growth stage was not significant (Table 2).

The CH₄ production potentials ranged from 141.14 (CK, ripening) to 288.48 (CW, heading) mg CH₄ kg^-1 dw h^-1 across the treatments and growth stages (Supplementary Figure S2). The highest values were recorded at the heading stage before a sharp decline at the ripening stage (Table 1). In contrary, CH₄ oxidation potential declined at the heading stage before increasing to the highest values at the ripening stage (Supplementary Figure S3). Repeated measures ANOVA showed CO₂ enrichment, warming and their interaction resulted in significant effects on CH₄ oxidation potential, while elevated CO₂ only affected CH₄ production potential when crossed with warming (Table 2). Although potential CH₄ production and oxidation potentials significantly differed among the growth stages, the effects of CO₂ enrichment and warming did not depend on the growth stage. CH₄ production potentials significantly increased under CO₂ enrichment (CE and CW) treatments at all three growth stages, but no significant changes were observed under WA treatment or for the control. However, CH₄ oxidation only increased at the ripening stage under warming treatments (33% for CW and 44% for WA).

The seasonal patterns of CH₄ emission profiles were similar between the climate change treatments and the control (Supplementary Figure S1). As the paddy field was waterlogged, CH₄ concentration peaks occurred 20–35 days after transplanting. After 35 days, CH₄ emissions sharply declined and remained at a low rate until harvest. Elevated CO₂ significantly increased CH₄ emissions in this paddy field. Compared to the CK treatment, mean CH₄ emissions increased by 37% and 25% under CE and CW treatments, respectively.

Gene abundance data showed that the abundances of mcrA and pmoA genes generally increased with rice growth development, the highest values being attained at the ripening stage (Figure 1). The abundance of mcrA genes under all treatments ranged from 9.13 × 10⁸ (WA, tillering) to 7.73 × 10⁹ (CE, ripening) copies g^-1 dw. These results were higher than the abundance of pmoA genes, ranging from 1.31 × 10⁸ (CK, tillering) to 4.55 × 10⁸ (CW, ripening) copies g^-1 dw. Repeated measures ANOVA showed that the effects of elevated CO₂ and elevated CO₂ combined with warming were significant (p < 0.05) on the abundance of mcrA genes, but not on the abundance of pmoA genes (Table 2). In this study, no significant changes in the abundance of mcrA and pmoA genes associated with climate change treatments were observed at the tillering stage. However, the abundance of mcrA genes significantly increased under CE and CW treatments at the heading and ripening stages. Compared to CK, the mean abundance of mcrA genes increased by 63% (CE) and 82% (CW) at the heading stage and 98% (CE) and 78% (CW) at the ripening stage, respectively. In contrast, the pmoA gene copy numbers were more stable without significant changes among the climate change treatments at all three growth stages.

The methanogenic and methanotrophic community structures were analyzed using T-RFLP fingerprints. PCA of the T-RFLP profiles at the three growth stages yielded summaries of data, as 51.4% for mcrA genes and 61.8% for pmoA genes of the total variability was explained by PC1 and PC2 (Figure 2). No clear differences in the methanogenic community structure between the treatments and across the growth stages were highlighted by PCA analysis. Figure 2B shows that the methanotrophic community structure under CE, CW, and WA was distinctively separated from the control at each stage.

A total of 12 and 8 T-RFs were obtained from the overall samples for mcrA and pmoA genes, respectively. T-RFLP fingerprinting of mcrA genes revealed that the 268bp, 289bp, 306bp, and 460bp TRFs were dominant in all treatments, while T-RFs of 76bp, 245bp, 437bp, and 510bp were dominant in the methanotrophic pmoA T-RFLP profiles (Figure 3). The main T-RFs in the T-RFLP profiles were identified with in silico restriction of clone sequences. T-RFs related to Methanosarcina were the most predominant (35–45%) across all treatments while T-RFs related to Methanocellales and Methanobacteriales ranged from 19–29% and 13–28%, respectively. The effect of simulated climate change scenarios on the relative abundance of mcrA TRFs was insignificant. However, the relative abundance of the 76bp TRF, related to Methylococcus, decreased under CE and CW treatments at the later growth stages, while that of the 437bp TRF related to Methylocystis increased.

The clone library analysis of mcrA and pmoA yielded a total of 94 and 102 cloned mcrA and pmoA sequences, respectively. Phylogenetic analyses of the mcrA sequences revealed that the methanogenic community in this paddy soil was dominated by members of Methanosarcina, Methanocellales, Methanobacteriales, and Methanomicrobiales (Figure 4). Twenty-one different OTUs of pmoA sequences were confirmed (Figure 5). The phylogenetic pattern of methanotrophic clones indicated that Methylococcus, Methylocaldum, Methylomonas, Methylosarcina, Methylogaea (Type I), and Methylocystis (Type II) were generally dominant during the rice growth period. As shown in Figure 5, Methylococcus (7 OTUs) and Methylocaldum (4 OTUs) sequences were the most dominant in Type I, representing about 31 and 25% of total clones in the pmoA clone library, respectively. To a smaller extent, Type I sequences were affiliated with Methylomonas (9 sequences), Methylosarcina (6 sequences) and Methylogaea (5 sequences). Methylocystis (3 OTUs, 25 sequences), the only groups of Type II methanotrophs, were identified in the rhizosphere.

Methanogenesis favors anoxic conditions with a low redox potential and a neutral pH in submerged soil (Masscheleyn et al., 1993; Orphan et al., 2001). In a laboratory incubation study, Das and Adhya (2012) found that elevated CO₂ significantly increased CH₄ production which was highly associated with decreased soil redox potential and pH in paddy soils. In this study, it was observed that CO₂ enrichment significantly increased CH₄ production potential during the rice growth period, whereas much higher SMBC was found (Table 1). This finding was consistent with results from Das et al. (2011) who observed an increase in SMBC with elevated CO₂ concentrations and temperatures in a rice soil incubation experiment. Das and Adhya (2012) noted that an increase of methanogenic population with elevated CO₂ levels was the most important reason for enhanced CH₄ production under CO₂ enrichment. In a free air CO₂ enrichment (FACE) experiment, Okubo et al. (2015) also found that the number of copies of the mcrA gene significantly increased when CO₂ concentrations were elevated in rice paddy fields. Similarly, our results showed that elevated CO₂ concentrations significantly increased the abundance of methanogens in the rice rhizosphere soil (Figure 1). In this experiment, photosynthesis and total biomass of rice increased with elevated CO₂ (Cai et al., 2016; Wang et al., 2016), resulting in an increase in root exudates and rhizodeposition in the rhizosphere (Bhattacharyya et al., 2013; Okubo et al., 2014). Furthermore, rhizodeposition is regarded as the primary source of CH₄ produced in rice fields (Conrad, 2007). Increased carbon inputs into the rhizosphere soil could also increase microbial biomass and activity (Drissner et al., 2007; French et al., 2009).

Based on T-RFLP and clone sequence analyses, we found that methanogens belonging to Methanosarcina, Methanobacteriales, and Methanocellales (Rice Cluster I) were dominant in this paddy soil. A similar result was found in other studies focusing on methanogenic diversity of paddy soils (Liu et al., 2012; Breidenbach et al., 2015). Methanosarcina, the only species of acetoclastic groups, were present during the rice growth period, suggesting that acetate may be a major substrate for CH₄ production in paddy fields. This observation is consistent with previous investigations showing Methanosarcina as a major species in rice fields under both flooded and drained conditions (Watanabe et al., 2009; Itoh et al., 2013). Based on RNA stable isotope probing analysis, Lu and Conrad (2005) indicated that individuals of the species Methanocellales were the most active for metabolizing rice root exudates, and that they play an important role in CH₄ production in paddy fields. In this study, elevated CO₂ and warming did not significantly alter the relative abundance of methanogenic T-RFs, resulting in minor changes in the composition of the methanogenic community. Thus, it is possible that the increase in CH₄ production with elevated atmospheric CO₂ is driven by the accumulation of substrate input and an increase in the methanogenic population or cell-specific activity without shifting methanogenic community composition.

In our study, a 2°C increase in air temperature did not change the activity, abundance or community composition of methanogens compared with elevated CO₂ treatments. This finding differs from the findings of Das and Adhya (2012); they identified an increase in CH₄ production and methanogenic population in trophic paddy soils associated with a temperature increase from 25–45°C, at intervals of ten degrees. Generally, an increase in temperature stimulates the decomposition of organic matter in submerged soils (von Lützow and Kögel-Knabner, 2009; Karhu et al., 2010) which may lead to higher rates of CH₄ production under anaerobic conditions (Fey and Conrad, 2000). In this study, the 2°C increase in canopy air temperature resulted in a small increase in soil temperature (<1°C, unpublished data), an increase which is in the range of daily/seasonal fluctuations and heterotrophic microbial communities are insensitive to such temperature increases (Malchair et al., 2010). Moreover, the direct effects of warming on soil microbial communities could be confounded by soil moisture content, substrate availability and plant conditions (Allison and Treseder, 2008; Das and Adhya, 2012).

Previous studies indicated that CH₄ oxidation potential increased as soil temperatures increased (Van den Pol-van Dasselaar et al., 1998; Dijkstra et al., 2010). Here, we showed that warming treatments (CW and WA) resulted in significant increases in CH₄ oxidation potential at the ripening stage when the paddy field was drained, whereas no change was observed under CE during the rice growth period. In a ¹³C-labeling study, Kalyuzhanaya et al. (2013) observed that a significantly larger portion of the assimilated carbon (over 62%) was derived from CO₂ in the alphaproteobacterial methanotrophs. Autotrophic methanotrophs, potentially utilizing C-1 compounds for their metabolic activity, could be greatly affected by changing environmental conditions. For example, if an increase of temperature results in the soil becoming drier, CH₄ oxidation may be enhanced (Dijkstra et al., 2010). Furthermore, as well as temperature, the water-logging regime is the other most important factor for CH₄ oxidation in rice fields (Das and Adhya, 2012). In the present study, such an increase of CH₄ oxidation potential under warming treatments was greater at the ripening stage (by 33–45%) when the rice fields were drained than at the tillering and heading stages (by 7–28%) when they were flooded. The growth stage of rice could therefore play an important role in the abundance and activity of methanotrophs during the rice growth period (Eller and Frenzel, 2001; Lüke et al., 2010; Lee et al., 2014). In a study of rice paddy soil, Ho et al. (2011) found that the higher methane oxidation potential related well to the cell-specific activity and population of methanotrophs. Though a nitrite-driven anaerobic CH₄ oxidation was discovered in wetlands (Hu et al., 2014), aerobic methanotrophs in rice field were predominant (Hu and Lu, 2015).

The composition and distribution of methanotrophs were relatively stable across the simulated climate change treatments. In this study, Methylococcus and Methylocaldum (both Type I), and Methylocystis (Type II), were dominant in the rice rhizosphere (Figures 4 and 5). Similar data have been found in rice field soils from Aichi-ken Anjo Research and Extension Center, central Japan (Jia et al., 2007), Gangetic plain of India (Vishwakarma et al., 2010) and National Rice Research Institute of China (Ma et al., 2010). It is generally assumed that Type I methanotrophs to be highly responsive to high substrate resources (Ho et al., 2013; Chen et al., 2014) while Type II methanotrophs are relatively stable (Eller et al., 2005; Krause et al., 2012). As an indication with a fast growth rate, Type I methanotrophs were predominant in occupying niches with abundant resources (Ho et al., 2013), such as in the rhizosphere and on the roots of rice plants (Wu et al., 2009; Chen et al., 2014). In this study, it is possible that increased labile organic compounds under warming are not large enough to affect the population or composition of the Type I methanotrophs.