Effects of Integrated Rice-Frog Farming on Paddy Field Greenhouse Gas Emissions

Integrated rice-frog farming (IRFF), as a mode of ecological farming, is fundamental in realizing sustainable development in agriculture. Yet its production of greenhouse gas (GHG) emissions remains unclear. Here, a randomized plot field experiment was performed to study the GHG emissions for various farming systems during the rice growing season. The farming systems included: conventional farming (CF), green integrated rice-frog farming (GIRF), and organic integrated rice-frog farming (OIRF). Results indicate that the cumulative methane (CH4) emissions from the whole growth period were divergent for the three farming systems, with OIRF having the highest value and CF having the lowest. For nitrous oxide (N2O) emissions, the order is reversed. IRFF significantly increased the dissolved oxygen (DO), soil redox potential (Eh), total organic carbon (TOC) content, and soil C:N ratio, which is closely related to GHG emissions in rice fields. Additionally, the average emissions of carbon dioxide (CO2) from soils during rice growing seasons ranged from 2312.27 to 2589.62 kg ha−1 and showed no significant difference in the three treatments. Rice yield in the GIRF and OIRF were lower (2.0% and 16.7%) than the control. The CH4 emissions contributed to 83.0–96.8% of global warming potential (GWP). Compared to CF, the treatment of GIRF and OIRF increased the GWP by 41.3% and 98.2% during the whole growing period of rice, respectively. IRFF significantly increased greenhouse gas intensity (GHGI, 0.79 kg CO2-eq ha−1 grain yield), by 91.1% over the control. Compared to the OIRF, GIRF decreased the GHGI by approximately 39.4% (0.59 kg CO2-eq ha−1 grain yield), which was 44.2% higher than that of the control. The results of structural equation model showed that the contribution of fertilization to CH4 emissions in paddy fields was much greater than that of frog activity. Moreover, frog activity could decrease GWP by reducing CH4 emissions from rice fields. And while GIRF showed a slight increase in GHG emissions, it could still be considered as a good strategy for providing an environmentally-friendly option in maintaining crop yield in paddy fields.


Introduction
Rice is the foremost staple food crop for nearly 50% of the current population in the world [1][2][3]. Nevertheless, a recent estimate of cropland GHG emissions indicates that paddy fields account for 48% of the global budget of GHG emissions primarily through discharges of CH 4 and N 2 O [4]. Globally, rice fields are significant sources of atmospheric CH 4 and N 2 O, and they are major contributors to global warming [5][6][7]. They exhibit relative GWP of 28 and 265 times that of CO 2 over a 100-year timescale [8]. Climate change caused by GHG emissions will have a huge impact on agriculture areas. Extreme weather events may result in lower harvestable yields, higher yield variability, and reduction Note: Different field management measures were used in different growing stages of rice. The total number of rice growth days were 131 days. The management of the field was carried out in accordance with the local high-yield pattern. Table 2 adapted from [25].

Sampling and Determination
For this study, soil parameters included TOC, C:N ratio, and soil temperature (T), while the water variables included DO, Eh, pH level, electrical conductivity (Ec), and the water level height of rice field. The DO, Eh, Ec, pH were determined by a multiparameter water quality analyzer (DZS-718L). Soil temperature was measured by a digital thermometer (TP101). The water level height was measured with a graduated scale. TOC was determined by the method of potassium dichromate oxidation. Total N was determined by a Smartchem 200 Discrete Auto Analyzer (Alliance Company, Paris, France). Rice yields per unit area were determined by theoretical method. Rice yields were the product of effective panicles per unit area, average grain number per panicle, The depth of the soil collected in the experiment was 0-20 cm. The soil texture was clay loam, with soil properties as follows: pH 6.9, soil total carbon 21.92 g kg −1 , total N 1.87 g kg −1 , total P 0.62 g kg −1 , total K 11.94 g kg −1 , available N 161.67 mg kg −1 , available P 20.98 mg kg −1 , available K 156.37 mg kg −1 , and soil moisture content 43.33%.
Rice seedlings were transplanted into the field in June and harvested in November. Rice growth stages can be divided into the following stages: pre-transplantation, regreening, tillering, jointing, booting, heading, filling, and maturing.

Experimental Design
The experiment was implemented in the field where continuous rice cultivation has been practiced for 10 years (from 2009 to 2018). Three rice cultivation patterns were carried out in the experiment: (1) Conventional farming (CF): Rice cultivation pattern with full application of chemical fertilizer; (2) Green integrated rice-frog farming (GIRF): Rice-frog co-cropping and applied with 50% chemical and 50% organic fertilizer; and, (3) Organic integrated rice-frog farming (OIRF): Rice-frog co-cropping and applied with 100% organic fertilizer.
In both GIRF and OIRF experimental fields, Chinese milk vetch (Astragalus sinicus L.) was seeded in the paddy fields after the rice harvest and ploughed into the soil the following May as Green manure (a basal fertilizer). The amount of nitrogen fertilizer applied to each treatment was the same at 300 kg N ha −1 [26][27][28]. The frogs used in the experiment was the tiger frog (Rana rugulosa), which has the main advantage of being large in size, strong adaptability, and large capacity of insect-catching. The frogs were released in the fields during the maturing stage and were left to self-sustain. During the rice growth stages, the water management was as follows: flood water layer of 5-8 cm in depth at the early stage, alternate long-term wetting state at mid-season, and field drainage at ten days before harvest (November 4th-14th). Other paddy management measures followed local high-yield field recommendations. Table 1 lists the schedule of nitrogenous fertilization for each treatment [25]. The major agricultural management practices used for the different growth stages of rice (Oryza sativa L.) are shown in Table 2 [25]. Notes: The unit for the fertilizer application rate is kg N ha −1 . The specific fertilization schemes of Chinese milk vetch, rapeseed cake, bulk blending fertilizer, bio-organic fertilizer, and urea were calculated from their total N contents of 0.5%, 5.3%, 26%, 6.2%, and 46%, respectively. No indication that fertilizer was not applied. Table 1 adapted from [25]. Note: Different field management measures were used in different growing stages of rice. The total number of rice growth days were 131 days. The management of the field was carried out in accordance with the local high-yield pattern.

Sampling and Determination
For this study, soil parameters included TOC, C:N ratio, and soil temperature (T), while the water variables included DO, Eh, pH level, electrical conductivity (Ec), and the water level height of rice field. The DO, Eh, Ec, pH were determined by a multiparameter water quality analyzer (DZS-718L). Soil temperature was measured by a digital thermometer (TP101). The water level height was measured with a graduated scale. TOC was determined by the method of potassium dichromate oxidation. Total N was determined by a Smartchem 200 Discrete Auto Analyzer (Alliance Company, Paris, France). Rice yields per unit area were determined by theoretical method. Rice yields were the product of effective panicles per unit area, average grain number per panicle, and weight of each grain. The panicle number, grain number, and grain weight of rice were counted and weighed manually.
Field sampling was carried out manually during the 2018 rice growing season. The CH 4 and N 2 O fluxes were simultaneously measured using the static chamber-gas chromatography (GC) method [29,30]. The sample collection chamber was made from acrylic plexiglass (7 cm thick) material and consisted of a collar, top box, and elevator box. In the center of each experimental plot, a collar with an area of 0.25 m 2 (50 × 50 cm) was permanently inserted in the soil to a depth of 10 cm and maintained in situ over the entire rice growth period. The rice planting density in the steel chamber-base collar was consistent with the experimental plot in accordance with the common practice of local farmers [2]. The top edge of the collar had a groove (6 cm in depth) for filling with water to seal the rim of the chamber. The dimension of the top box, which was sealed at the head, was 50 × 50 × 50 cm, and each top box was equipped with two circulating fans to ensure complete gas mixing. Two holes were cut at the top and middle parts of the top box to determine the temperature and collect gas samples inside the box. The dimensions of the elevator box were the same as those of the top box. It had two covers at both ends that were allowed to open, and a groove at the top to be filled with water to seal the chamber's rim. The elevator box was only used after rice jointing. The sample collection chamber was wrapped with a layer of aluminum foil to minimize air temperature changes inside the chamber during the period of gas sampling [31]. Gas samples were collected once a week from the 10th day after transplanting. Gas samples were taken from 08:00 to 10:00 am to closely resemble the daily average soil temperature and minimize the influence of diurnal variation in CH 4 and N 2 O emissions during the sampling period [32]. Gas samples were collected from each chamber and were placed in pre-evacuated vacuum aluminum film airbags (produced by Dalian Delin Gas Packing Limited Company, Dalian, China) with a capacity of 100 mL at 10 min intervals (0, 10, 20, 30 and 40 min after chamber closure). The exact time of each sampling and the sample number were also recorded [33]. Gas samples in the vacuum air bags were immediately transported to the laboratory for analysis by GC within five days.
The concentrations of CH 4 and N 2 O in the gas samples were determined using a modified GC (Agilent 7890N, Santa Clara, CA, USA), equipped with a flame ionization detector or CH 4 analysis and an electron capture detector for N 2 O analysis [34]. The carrier gases used to carry CH 4 and N 2 O were pure nitrogen (99.99%) and a gas mixture of argon and CH 4 , respectively. The oven and the flame ionization detector were operated at temperatures of 55 • C and 200 • C, respectively. Further details on the principles, techniques, instrument configurations, and operation procedures are discussed by Zheng et al. [35] and Wang et al. [36]. The standard gas was supplied by the National Standard Material Center. In the study, five standard gas samples were utilized to ensure the stability of the gas chromatography. Fluxes were determined based on the slope of the change in mixing ratio for five sequential samples. Sample sets were rejected unless they yielded a linear regression value of r 2 > 0.90 [32]. The GHG emissions flux was calculated by the differences in gas concentrations in accordance with the following equation [27,28]: where F in Equation (1) represents the CH 4 and N 2 O fluxes (mg m −2 h −1 ); represents the CH 4 and N 2 O density at the standard state (kg m −3 ), CH 4 was 0.714 kg m −3 , N 2 O was 1.964 kg m −3 , H is the height of the chamber above the soil/water surface (m); dC/dt is the rate of change in the CH 4 and N 2 O concentration with respect to time (t) in the chamber (mL m −3 h −1 ); and, T is the average air temperature inside the chamber during sampling ( • C).

GWP and GHGI Evaluation
GWP is an index used to evaluate the combined radiative forcing potential of all the GHG, including CO 2 , CH 4 , and N 2 O. In this study, GWP expressed in CO 2 -equivalents (CO 2 -eq) was estimated, taking into account cumulative soil emissions of CH 4 and N 2 O, assuming a 100-year time horizon [8]. The net change in TOC was not measured for this short-term experiment. The CO 2 emissions collected by the dark chamber were not the net flux of the ecosystem, so we were not able to quantify the net CO 2 emissions from the soils [26]. The net effects of CH 4 and N 2 O emissions were expressed in CO 2 -equivalents according to the GWP factors provided by the Intergovernmental Panel on Climate Change (IPCC). The GWP of CH 4 and N 2 O emissions (i.e., CO 2 equivalents [CO 2 -eq]) were calculated using the latest GWP values: 28 for CH 4, and 265 for N 2 O over a 100-year time horizon [8]. The GWP was calculated through the following equation: where GWP in Equation (2) is the global warming potential for CH 4 and N 2 O (kg CO 2 -eq ha −1 ); and, CH 4 and N 2 O are the total CH 4 and N 2 O emissions (kg ha −1 ). The GHGI is the ratio of total warming potential of CH 4 and N 2 O to crop yield [37][38][39], such that: where GWP in Equation (3) is the total warming potential of CH 4 and N 2 O (CO 2 kg ha −1 ); and, Y is the average yield per unit area of the treatment (kg ha −1 ).

Statistical Analysis
All statistical analyses were mainly performed using the SPSS software (version 24.0, IBM, Armonk, NY, USA) and statistical significance was determined at the 0.05 probability level. Origin (version 9.1, OriginLab Co., Northampton, MA, USA) was employed for figure preparation. Differences in CH 4 and N 2 O emissions between treatments were examined using a one-way analysis of variance (ANOVA). Correlation coefficient visibility graph was generated using the 'ggcorrplot' package in software R (version 3.4.2). Redundancy analysis (RDA) was performed to summarize the GHG emissions, which may be explained by the soil and water variables in paddy fields using CANOCO 5.0 software (Microcomputer Power, Ithaca, NY, USA). Finally, the fitting of the structural equation model was accomplished using the Amos 21.0 software (IBM, Armonk, NY, USA). The structural equation model is a multivariate statistical method which can describe the relationship between variables that cannot be directly measured and belongs to the confirmatory model [40].

Environmental Factors in Rice Fields
Statistical analysis reveals that IRFF significantly increased the DO, soil Eh, TOC content, and soil C:N ratio in rice fields. However, the effect of IRFF on soil pH and soil temperature was minimal. In general, the effect of OIRF on soil environment of paddy fields was greater than that of GIRF. Throughout the rice planting season, IRFF significantly increased the DO, soil Eh, TOC content, and soil C:N ratio by 7.95-12.92%, 8.27-9.57%, 25.18-42.39% and 14.86-30.51%, respectively, when compared to CF ( Figure 2). to quantify the net CO2 emissions from the soils [26]. The net effects of CH4 and N2O emissions were expressed in CO2-equivalents according to the GWP factors provided by the Intergovernmental Panel on Climate Change (IPCC). The GWP of CH4 and N2O emissions (i.e., CO2 equivalents [CO2-eq]) were calculated using the latest GWP values: 28 for CH4, and 265 for N2O over a 100-year time horizon [8]. The GWP was calculated through the following equation: where GWP in equation (2) is the global warming potential for CH4 and N2O (kg CO2-eq ha −1 ); and, CH4 and N2O are the total CH4 and N2O emissions (kg ha −1 ). The GHGI is the ratio of total warming potential of CH4 and N2O to crop yield [37][38][39], such that: where GWP in Equation (3) is the total warming potential of CH4 and N2O (CO2 kg ha −1 ); and, Y is the average yield per unit area of the treatment (kg ha −1 ).

Statistical Analysis
All statistical analyses were mainly performed using the SPSS software (version 24.0, IBM, Armonk, NY, USA) and statistical significance was determined at the 0.05 probability level. Origin (version 9.1, OriginLab Co., Northampton, MA, USA) was employed for figure preparation. Differences in CH4 and N2O emissions between treatments were examined using a one-way analysis of variance (ANOVA). Correlation coefficient visibility graph was generated using the 'ggcorrplot' package in software R (version 3.4.2). Redundancy analysis (RDA) was performed to summarize the GHG emissions, which may be explained by the soil and water variables in paddy fields using CANOCO 5.0 software (Microcomputer Power, Ithaca, NY, USA). Finally, the fitting of the structural equation model was accomplished using the Amos 21.0 software (IBM, Armonk, NY, USA). The structural equation model is a multivariate statistical method which can describe the relationship between variables that cannot be directly measured and belongs to the confirmatory model [40].

Environmental Factors in Rice Fields
Statistical analysis reveals that IRFF significantly increased the DO, soil Eh, TOC content, and soil C:N ratio in rice fields. However, the effect of IRFF on soil pH and soil temperature was minimal. In general, the effect of OIRF on soil environment of paddy fields was greater than that of GIRF. Throughout the rice planting season, IRFF significantly increased the DO, soil Eh, TOC content, and soil C:N ratio by 7.95-12.92%, 8.27-9.57%, 25.18-42.39% and 14.86-30.51%, respectively, when compared to CF ( Figure 2). Using Pearson correlation analysis, CH 4 emissions were shown to be negatively correlated with TOC, DO and soil Eh, and were positively correlated with soil C:N ratio and the water level height of rice fields. The N 2 O emissions were positively correlated with TOC and soil temperature and were negatively correlated with soil Eh and Ec. The net CO 2 emissions were significantly affected by soil Ec, Eh, water height, and pH, while TOC and C:N ratio had negligible effect ( Figure 3).

Figure 2.
Effects of IRFF (Integrated rice-frog farming) on soil parameters during the rice growing seasons. * and ** indicate statistical significance at the 0.05 and 0.01 levels, respectively, while ns means not significant. The vertical bars indicate the standard deviation of the means (n = 3 replicates). GIRF: green integrated rice-frog farming; OIRF: organic integrated rice-frog farming.
Using Pearson correlation analysis, CH4 emissions were shown to be negatively correlated with TOC, DO and soil Eh, and were positively correlated with soil C:N ratio and the water level height of rice fields. The N2O emissions were positively correlated with TOC and soil temperature and were negatively correlated with soil Eh and Ec. The net CO2 emissions were significantly affected by soil Ec, Eh, water height, and pH, while TOC and C:N ratio had negligible effect ( Figure  3).

CH4 Emissions
The relationships between the CH4, N2O, CO2 emissions, and soil parameters were assessed using redundancy analysis (RDA) (Figure 4). The first axis explains 99.97% of the total variation in the GHG emissions. The samples from chemical fertilizer treatment are in the third and fourth quadrants, which show that DO and Eh were significantly and negatively correlated with CF. Sample 3 shows that C:N had significant effect on CF. The GIRF data points appear scattered, the DO and Eh were negatively correlated with GIRF, while C:N, TOC, and T were positively correlated. Sample 17 shows that C:N and TOC were the main factors affecting CH4 emissions in totally organic fertilizer treatment of rice.
In the different stages in rice growth, the CH4 emissions fluxes from the three systems were different ( Figure 5). The cumulative CH4 emissions from the entire growth period were also very different, with OIRF being the highest in value and CF being the lowest. At the rice regreening stage, the emissions of CH4 was lower, with values ranging from 1500 to 2000 mg m −2 . There was no significant difference between the different systems. At the rice tillering stage, the CH4 emissions were more than the previous phase. The emissions of CH4 in OIRF were higher than CF and GIRF, possibly due to having lower DO and higher TOC in paddy waters treated with OIRF ( Figure 4). The soil environment of rice fields in OIRF was more favorable in increasing the activity of methanogenic bacteria and in decreasing the activity of CH4-oxidizing bacteria. Thus, the amount of CH4 produced by OIRF is higher than that by CF and GIRF. At the rice jointing stage, CH4 emissions from OIRF treatment increased rapidly, reaching its peak eight and six times higher than those of CF and GIRF, respectively. The sudden increase in CH4 emissions from OIRF treatment was mainly due to the

CH 4 Emissions
The relationships between the CH 4 , N 2 O, CO 2 emissions, and soil parameters were assessed using redundancy analysis (RDA) (Figure 4). The first axis explains 99.97% of the total variation in the GHG emissions. The samples from chemical fertilizer treatment are in the third and fourth quadrants, which show that DO and Eh were significantly and negatively correlated with CF. Sample 3 shows that C:N had significant effect on CF. The GIRF data points appear scattered, the DO and Eh were negatively correlated with GIRF, while C:N, TOC, and T were positively correlated. Sample 17 shows that C:N and TOC were the main factors affecting CH 4 emissions in totally organic fertilizer treatment of rice.
In the different stages in rice growth, the CH 4 emissions fluxes from the three systems were different ( Figure 5). The cumulative CH 4 emissions from the entire growth period were also very different, with OIRF being the highest in value and CF being the lowest. At the rice regreening stage, the emissions of CH 4 was lower, with values ranging from 1500 to 2000 mg m −2 . There was no significant difference between the different systems. At the rice tillering stage, the CH 4 emissions were more than the previous phase. The emissions of CH 4 in OIRF were higher than CF and GIRF, possibly due to having lower DO and higher TOC in paddy waters treated with OIRF ( Figure 4). The soil environment of rice fields in OIRF was more favorable in increasing the activity of methanogenic bacteria and in decreasing the activity of CH 4 -oxidizing bacteria. Thus, the amount of CH 4 produced by OIRF is higher than that by CF and GIRF. At the rice jointing stage, CH 4 emissions from OIRF treatment increased rapidly, reaching its peak eight and six times higher than those of CF and GIRF, respectively. The sudden increase in CH 4 emissions from OIRF treatment was mainly due to the application of organic fertilizer (rapeseed cake) during rice jointing stage. Rapeseed cake contains a large amount of organic carbon, which provides an abundant precursor for the production of CH 4 . Another important reason is that the low content of DO in rice fields also promotes the growth of CH 4 -producing bacteria.
At the rice booting stage, CH 4 emissions also decreased significantly due to a reduction in TOC content and an increase in DO content. The CH 4 emissions of GIRF reached its peak at the booting stage, which may be due to sizable seepage of stored CH 4 as a result of frog movement. During the rice heading stage, CH 4 emissions from the different systems all showed a decreasing trend. Among them, the declines in CF and GIRF were relatively rapid. However, the CH 4 emissions of CF and GIRF were still higher than OIRF, which may be due to the application of urea. At the rice filling stage, rice growth has almost ceased, and CH 4 emissions were at very low levels. During the rice maturing stage, the CH 4 emissions in all systems were very low, mainly caused by lower DO content and Eh. Although the TOC content had also been very low at this stage, the main factor was DO. Concomitantly, the water drainage in paddy fields resulted in the enhancement of soil aeration, reducing the activity of methanogenic bacteria while enhancing CH 4 -oxidizing bacteria.
The DO and Eh were significantly and negatively correlated with CF while C:N was significantly and positively correlated with GIRF. The TOC was significantly and positively related to OIRF. The C:N and TOC had significant positive correlation with CH 4 emissions, while DO and Eh had negative effect [41,42]. Also, Ec and T did not exhibit significant correlation with CH 4 emissions (Figure 4). Based on the entire rice growth cycle, the order of CH 4 emissions was OIRF > GIRF > CF, and they have significant differences. This indicates that the pattern of mixed fertilization, rather than the full application of organic fertilizer, can alleviate CH 4 emissions from paddy fields.
large amount of organic carbon, which provides an abundant precursor for the production of CH4. Another important reason is that the low content of DO in rice fields also promotes the growth of CH4-producing bacteria.
At the rice booting stage, CH4 emissions also decreased significantly due to a reduction in TOC content and an increase in DO content. The CH4 emissions of GIRF reached its peak at the booting stage, which may be due to sizable seepage of stored CH4 as a result of frog movement. During the rice heading stage, CH4 emissions from the different systems all showed a decreasing trend. Among them, the declines in CF and GIRF were relatively rapid. However, the CH4 emissions of CF and GIRF were still higher than OIRF, which may be due to the application of urea. At the rice filling stage, rice growth has almost ceased, and CH4 emissions were at very low levels. During the rice maturing stage, the CH4 emissions in all systems were very low, mainly caused by lower DO content and Eh. Although the TOC content had also been very low at this stage, the main factor was DO. Concomitantly, the water drainage in paddy fields resulted in the enhancement of soil aeration, reducing the activity of methanogenic bacteria while enhancing CH4-oxidizing bacteria.
The DO and Eh were significantly and negatively correlated with CF while C:N was significantly and positively correlated with GIRF. The TOC was significantly and positively related to OIRF. The C:N and TOC had significant positive correlation with CH4 emissions, while DO and Eh had negative effect [41,42]. Also, Ec and T did not exhibit significant correlation with CH4 emissions (Figure 4). Based on the entire rice growth cycle, the order of CH4 emissions was OIRF > GIRF > CF, and they have significant differences. This indicates that the pattern of mixed fertilization, rather than the full application of organic fertilizer, can alleviate CH4 emissions from paddy fields.

N2O Emissions
The pH and T had significant positive correlation with N2O emissions. TOC was also positively associated with N2O emissions, while the correlation between C:N and N2O emissions was very weak. Moreover, Ec, Eh, and DO were negatively associated with N2O emissions (Figure 4). Calculating for the whole growth period of rice, the order of N2O emissions was CF > GIRF > OIRF,

N 2 O Emissions
The pH and T had significant positive correlation with N 2 O emissions. TOC was also positively associated with N 2 O emissions, while the correlation between C:N and N 2 O emissions was very weak. Moreover, Ec, Eh, and DO were negatively associated with N 2 O emissions ( Figure 4). Calculating for the whole growth period of rice, the order of N 2 O emissions was CF > GIRF > OIRF, and they have significant differences.
In the different rice growth stages, the N 2 O emissions fluxes of the three systems were different ( Figure 6). The cumulative N 2 O emissions from the entire growth period were likewise dissimilar with the order of discharge being CF > GIRF > OIRF, antithetical to CH 4 emissions. At the rice regreening stage, the emissions of N 2 O was lower, with values below 5 mg m −2 . This was mainly due to the small biomass of rice. At the rice tillering stage, N 2 O emissions of CF was more than GIRF and OIRF. There was also no significant difference between the two treatments. At the rice jointing stage, N 2 O emissions from CF were significantly higher than those of the other two treatments. One possible reason is that the lower C:N ratio indicates more N in the soil, providing the substrate for nitrification. In addition, less DO in the soil results in enhanced nitrification and promotes N 2 O emissions. At the rice booting stage, N 2 O emissions from CF decreased slightly, the GIRF reached its peak, and N 2 O emissions from OIRF also increased slightly. This may be due to a decrease in N content of CF, while N in both GIRF and OIRF are increased with GIRF having higher N content. Additionally, the frog activity in GIRF and OIRF have increased the DO in the soil, which further promoted the emissions of N 2 O. At the rice heading stage, the N 2 O emissions of CF was relatively high, the GIRF was in decline, and the change in OIRF was not significant. The reason was that in the CF treatment, a large amount of urea had been applied with high N content, resulting in more substrate substances producing N 2 O. At the rice filling stage, the N 2 O emissions in all three treatments were decreased, probably because the rice did not grow nutritionally and the N content in the soil decreased significantly. During the rice maturing stage, the N 2 O emissions in all treatments was very low.

N2O Emissions
The pH and T had significant positive correlation with N2O emissions. TOC was also positively associated with N2O emissions, while the correlation between C:N and N2O emissions was very weak. Moreover, Ec, Eh, and DO were negatively associated with N2O emissions (Figure 4). Calculating for the whole growth period of rice, the order of N2O emissions was CF > GIRF > OIRF, and they have significant differences.
In the different rice growth stages, the N2O emissions fluxes of the three systems were different ( Figure 6). The cumulative N2O emissions from the entire growth period were likewise dissimilar with the order of discharge being CF > GIRF > OIRF, antithetical to CH4 emissions. At the rice regreening stage, the emissions of N2O was lower, with values below 5 mg m −2 . This was mainly due to the small biomass of rice. At the rice tillering stage, N2O emissions of CF was more than GIRF and OIRF. There was also no significant difference between the two treatments. At the rice jointing stage, N2O emissions from CF were significantly higher than those of the other two treatments. One possible reason is that the lower C:N ratio indicates more N in the soil, providing the substrate for nitrification. In addition, less DO in the soil results in enhanced nitrification and promotes N2O emissions. At the rice booting stage, N2O emissions from CF decreased slightly, the GIRF reached its peak, and N2O emissions from OIRF also increased slightly. This may be due to a decrease in N content of CF, while N in both GIRF and OIRF are increased with GIRF having higher N content. Additionally, the frog activity in GIRF and OIRF have increased the DO in the soil, which further promoted the emissions of N2O. At the rice heading stage, the N2O emissions of CF was relatively high, the GIRF was in decline, and the change in OIRF was not significant. The reason was that in the CF treatment, a large amount of urea had been applied with high N content, resulting in more substrate substances producing N2O. At the rice filling stage, the N2O emissions in all three treatments were decreased, probably because the rice did not grow nutritionally and the N content in the soil decreased significantly. During the rice maturing stage, the N2O emissions in all treatments was very low.

CO 2 Emissions
The biomass of rice was measured in the different growth stages by the method of sample plot harvesting. The CO 2 absorbed by rice ecosystem during biomass production was then calculated using the photosynthesis Equation (4). Plants can absorb 264 g CO 2 if they grow 162 g polysaccharide organic matter. That is, the plant can absorb 1.63 g CO 2 for every 1 g of dry matter accumulated by the plant [ The emissions of CO 2 measured in field experiments were a comprehensive flux including photosynthetic fixation, rice respiration, and soil respiration in rice. The index can fully reflect the regulating function of rice ecosystem in the atmosphere, and it shows some negative values. Therefore, it is referred to as the CO 2 absorption flux in paddy fields [43].
The soil temperature was significantly and positively related to the emissions of CO 2 (Figure 4). At the tillering stage of rice, CO 2 emissions reached its peak and then decreased. There was no significant difference among the systems in terms of net cumulative CO 2 emissions. Generally, the average CO 2 emissions from soils during this rice growing period were high and ranged from 2312.27 to 2589.62 kg ha −1 (Figure 7). using the photosynthesis Equation (4). Plants can absorb 264 g CO2 if they grow 162 g polysaccharide organic matter. That is, the plant can absorb 1.63 g CO2 for every 1 g of dry matter accumulated by the plant [43]: 6n CO2 + 6n H2O → n C6H12O6 + 6n O2 → n C6H10O5 (4) 264 108 180 192 162 (4) The emissions of CO2 measured in field experiments were a comprehensive flux including photosynthetic fixation, rice respiration, and soil respiration in rice. The index can fully reflect the regulating function of rice ecosystem in the atmosphere, and it shows some negative values. Therefore, it is referred to as the CO2 absorption flux in paddy fields [43].
The soil temperature was significantly and positively related to the emissions of CO2 (Figure 4). At the tillering stage of rice, CO2 emissions reached its peak and then decreased. There was no significant difference among the systems in terms of net cumulative CO2 emissions. Generally, the average CO2 emissions from soils during this rice growing period were high and ranged from 2312.27 to 2589.62 kg ha −1 (Figure 7).

Rice Yield and GHGI
No growth impairment of rice plants was observed during the cropping period. Rice growth and yield properties were not significantly improved by the GIRF and OIRF systems at rice harvesting stage (Table 3). Rice yields in the GIRF and OIRF were lower (2.0% and 16.7%) than the control. The ripened grain and rice bulk density of the yield in the two IRFF systems were higher than the control. However, the other yield component was reduced by the IRFF. In comparison, rice growth and yield characteristics were not significantly different between OIRF and GIRF treatments ( Table 3). Rice yield was slightly higher (17.7%) in GIRF than the OIRF treatment, but other yield properties and growth characteristics were not significantly different between the two IRFF treatments.

Rice Yield and GHGI
No growth impairment of rice plants was observed during the cropping period. Rice growth and yield properties were not significantly improved by the GIRF and OIRF systems at rice harvesting stage (Table 3). Rice yields in the GIRF and OIRF were lower (2.0% and 16.7%) than the control. The ripened grain and rice bulk density of the yield in the two IRFF systems were higher than the control. However, the other yield component was reduced by the IRFF. In comparison, rice growth and yield characteristics were not significantly different between OIRF and GIRF treatments (Table 3). Rice yield was slightly higher (17.7%) in GIRF than the OIRF treatment, but other yield properties and growth characteristics were not significantly different between the two IRFF treatments. In this study, CH 4 emissions contributed to 83.0-96.8% of GWP; thus, the effect of IRFF on GWP was similar to CH 4 emissions. During the rice growing seasons, IRFF significantly affected the GWP (Figure 8). Compared to CF, the treatment of GIRF and OIRF in the rice growing cycle increased the GWP by 41.3% and 98.2%, respectively. The GWP value of OIRF treatment reached its apex of 4723.63 kg ha −1 at the jointing stage, and subsequently decreased. For the GIRF system, the GWP value reached its peak at the booting stage (2642.57 kg ha −1 ), and then decreased gradually. The GWP value of the control treatment fluctuated slightly in the different rice growth stages. This result indicates that effective measures adopted during rice jointing and booting stages would be beneficial in mitigating GWP from paddy fields. In this study, CH4 emissions contributed to 83.0-96.8% of GWP; thus, the effect of IRFF on GWP was similar to CH4 emissions. During the rice growing seasons, IRFF significantly affected the GWP (Figure 8). Compared to CF, the treatment of GIRF and OIRF in the rice growing cycle increased the GWP by 41.3% and 98.2%, respectively. The GWP value of OIRF treatment reached its apex of 4723.63 kg ha -1 at the jointing stage, and subsequently decreased. For the GIRF system, the GWP value reached its peak at the booting stage (2642.57 kg ha −1 ), and then decreased gradually. The GWP value of the control treatment fluctuated slightly in the different rice growth stages. This result indicates that effective measures adopted during rice jointing and booting stages would be beneficial in mitigating GWP from paddy fields.
GHGI, which indicates net GWP per yield, was approximately 0.41 kg CO2-eq kg −1 grain in the control treatment ( Figure 8). IRFF significantly increased GHGI (0.79 kg CO2-eq ha −1 grain yield), by 91.1% over the control. Compared to the complete bio-organic fertilizer application, the bulk blending fertilizer treatment decreased the GHGI by approximately 39.4% (0.59 kg CO2-eq ha −1 grain yield), which was 44.2% higher than the control.

Structural Equation Modeling
This study is not based only on a single factor treatment, but on the rice ecological integrated production pattern that is now being widely promoted. It is therefore uncertain how much of the GHGI, which indicates net GWP per yield, was approximately 0.41 kg CO 2 -eq kg −1 grain in the control treatment (Figure 8). IRFF significantly increased GHGI (0.79 kg CO 2 -eq ha −1 grain yield), by 91.1% over the control. Compared to the complete bio-organic fertilizer application, the bulk blending fertilizer treatment decreased the GHGI by approximately 39.4% (0.59 kg CO 2 -eq ha −1 grain yield), which was 44.2% higher than the control.

Structural Equation Modeling
This study is not based only on a single factor treatment, but on the rice ecological integrated production pattern that is now being widely promoted. It is therefore uncertain how much of the results are caused by the introduction of frogs and by fertilization. Only by isolating the distinct roles of frogs and fertilization can we fully explain the impact of integrated rice-frog farming in GHG emissions. Since there are too many CO 2 -influencing factors and with the involvement of photosynthesis in rice growth, explaining the emissions mechanism clearly is more complex. And since the contributions of CH 4 to GWP is much greater than that of N 2 O in paddy fields, this study mainly focuses on CH 4 emissions. Using structural equation modeling, the contribution coefficient of complex environmental factors to CH 4 emissions can be estimated.

The Establishment of Conceptual Modeling
Based on the characteristics of CH 4 emissions from paddy fields, a conceptual model of the primary factors influencing CH 4 in paddy fields was established (Figure 9). The model consists of two latent variables (Frog (ξ1) and Fertilization (ξ2)), and five measurable variables (DO (x1), Eh (x2), C:N (x3), TOC (x4), and Methane (y1). Based on the previous analysis of the primary factors, the following hypotheses are given for the conceptual model:  photosynthesis in rice growth, explaining the emissions mechanism clearly is more complex. And since the contributions of CH4 to GWP is much greater than that of N2O in paddy fields, this study mainly focuses on CH4 emissions. Using structural equation modeling, the contribution coefficient of complex environmental factors to CH4 emissions can be estimated.

The Establishment of Conceptual Modeling
Based on the characteristics of CH4 emissions from paddy fields, a conceptual model of the primary factors influencing CH4 in paddy fields was established (Figure 9). The model consists of two latent variables (Frog (ξ1) and Fertilization (ξ2)), and five measurable variables (DO (x1), Eh (x2), C:N (x3), TOC (x4), and Methane (y1). Based on the previous analysis of the primary factors, the following hypotheses are given for the conceptual model:

Model Fitting Index Analysis
For the established conceptual models and assumptions, the first initial model was fitted by using Amos 21.0. After repeated fitting, evaluation, and correction of the model, the final normalized coefficient correction model was obtained (Figure 10). By analyzing the structural equation model with Amos 21.0, better fitting indices can be obtained (Table 4). Based on previous recommendations on the structural equation model, the absolute fitting index, the relative fitting index, and the reduced index had been used in this study. The index includes CMIN/DF, GFI, RMSEA, NFI, TLI, CFI, CLE, AIC, and ECVI. The fitting standard was to be determined. The value of CMIN/DF greater than ten indicates that the model is not ideal; less than five, the model is acceptable, but a value below three would be recommended. The GFI, NFI, TLI, CFI, and IFI values should all be more than 0.9; and the closer the value is to one, the better the effect. Also, the smaller the values of AIC and ECVI, the better the fitting effect [44].
The results of the model fitting index analysis are shown in Table 4. The fitting index of the model was generally acceptable, meeting all fitting index requirements. The relationship model of CH4 emissions in paddy fields obtained from the statistical method was reasonable. The data in the diagram are impacted path coefficients of the modification model for balanced relations between key impact factors ( Figure 10).

Model Fitting Index Analysis
For the established conceptual models and assumptions, the first initial model was fitted by using Amos 21.0. After repeated fitting, evaluation, and correction of the model, the final normalized coefficient correction model was obtained (Figure 10). By analyzing the structural equation model with Amos 21.0, better fitting indices can be obtained (Table 4). Based on previous recommendations on the structural equation model, the absolute fitting index, the relative fitting index, and the reduced index had been used in this study. The index includes CMIN/DF, GFI, RMSEA, NFI, TLI, CFI, CLE, AIC, and ECVI. The fitting standard was to be determined. The value of CMIN/DF greater than ten indicates that the model is not ideal; less than five, the model is acceptable, but a value below three would be recommended. The GFI, NFI, TLI, CFI, and IFI values should all be more than 0.9; and the closer the value is to one, the better the effect. Also, the smaller the values of AIC and ECVI, the better the fitting effect [44].    The results of the model fitting index analysis are shown in Table 4. The fitting index of the model was generally acceptable, meeting all fitting index requirements. The relationship model of CH 4 emissions in paddy fields obtained from the statistical method was reasonable. The data in the diagram are impacted path coefficients of the modification model for balanced relations between key impact factors ( Figure 10).

Model Result Analysis
Fertilization had a positive effect on CH 4 emissions from paddy fields, and the path coefficient was 0.62, indicating it has significant impact in increasing gas discharge ( Figure 10). Frogs had a negative effect on CH 4 emissions from paddy fields. The path coefficient was only −0.37, which was less than fertilization. CH 4 emissions was mainly related to soil fertilization, such that the use and application of fertilizer significantly promotes CH 4 emissions. Frog activities in paddy fields can inhibit the emissions of CH 4 , although their effects was shown to be small. There was a negative correlation between fertilization and frog behavior; the correlation coefficient was −0.81, which was consistent with the direct negative effect of frog behavior on CH 4 emissions ( Figure 10). Therefore, CH 4 emissions from paddy fields can be reduced by modifying the quantity and mode of fertilizer used in the fields.
From the three measurable variables indicative of frog behavior, DO had the strongest contribution to frog behavior, followed by C:N, and then Eh. The path coefficients were 1.59, 0.83, and 0.78, respectively. It is widely recognized that DO has a great effect on soil CH 4 emissions, and it is a measurable variable that affects the relation and trend. The C:N affects the content of methanogenic substrate and is also a measurable variable which contributes greatly to CH 4 production. The influence of Eh on soil CH 4 emissions is theoretically significant; however, Eh is highly influenced by various external factors such as solution temperature, pH, and chemical reaction reversibility. In paddy waters, the complex redox system is formed. Between the two measurable variables of the fertilization-latent variable, C:N had higher contribution to fertilizer application than the TOC. The path coefficients were 1.55 and 0.98, respectively, which indicate that C, N, and TOC have great influence on CH 4 emissions. CH 4 emissions from paddy field is a complex dynamic process, resulting from the interaction between the latent-to-latent variables, latent-to-measurable variables, and measurable-to-measurable variables.
In this study, both OIRF and GIRF systems were assessed on the impact of two principal factors: fertilizer use and introduction of frogs. In both systems, fertilizer and frogs were shown to play major roles in affecting the rate of CH 4 emissions. The use of the structural equation model provides the means to further understand this complex relationship. Thus, the specific contribution of fertilizer and frog behavior to CH 4 emissions were calculated. The results show that the contribution of fertilization to CH 4 emissions in paddy fields is much greater than that of frog activity ( Figure 10).

Effects of Frogs on CH 4 Emissions
"Frog behavior" in this study refers to the activities of frogs in rice fields, including their movement and excretion behaviors. The movement behavior of frogs can disturb the water layer and affect the DO and Eh values in paddy fields. The excretion behavior of frogs primarily involves frog feces production, which affect the TOC and C:N values in rice fields. Thus, DO, Eh, TOC, and C:N become proxy indicators for frog behavior. By using the structural equation model, we have proven that frogs' behavior can reduce CH 4 emissions.
The deeper question now is how do frogs reduce CH 4 emissions? The behavior of frogs in paddy fields is similar to that of ducks, fish, and shrimp. Studies have shown that ducks and fish can affect GHG emissions, control weeds and pests, and minimize diseases in rice fields. The GWP of GHG from integrated rice-duck farming system had been studied by Yuan et al., [45] Xu et al., [46] and Zhan et al., [47]. Their results indicate that while the introduction of ducks in paddy fields can promote N 2 O emissions generated from duck feces, it also increases the concentration of DO in the water layer and reduces CH 4 emissions. Overall, their studies found that the integrated rice-duck farming system decreases the GWP in rice fields. However, Frei et al., [48] Datta et al., [49] and Bhattacharyya et al. [23] concluded that carp production in paddy fields promotes CH 4 diffusion and discharge through the water layer. Fish consume the DO in the water and reduced the Eh, thus increasing CH 4 emissions. As for the introduction of shrimp, previous research have confirmed that the activities of rice shrimp could greatly increase the oxygen content in the soil and water surface [50].
In this study, our statistical analysis indicate that IRFF had significantly increased the content of DO, Eh, TOC and C:N in rice fields. The increase in DO and Eh that helps reduce CH 4 emissions was much higher than the increases in TOC and C:N which promote greater CH 4 emissions. Taken aggregately, the introduction of frogs provides an inhibitory effect on CH 4 emissions in rice fields. However, this study did not dwell in analyzing frog activities in detail. We recommend that future studies investigate specific frog behavior in paddy fields to further explain the effects of frogs on CH 4 emissions more scientifically.

Uncertainty and Prospect
The main limitation of the study is its use of low-frequency measurement in only one rice growing cycle. The soil environment was relatively stable, while we observed unstable meteorological factors frequently. Although we only monitored one rice growing season, we measured many indicators closely related to CH 4 and N 2 O emissions, such as DO, Eh, TOC, soil temperature, air temperature, and water depth of paddy fields. We made sure that we recorded and observed the uncertainty factors affecting GHG emissions, reduced the random factors to a minimum, and increased the reliability and objectivity of the results. It should be noted that CH 4 and N 2 O from paddy soil have a certain level of variability due to the variations in soil attributes and other environmental parameters. The results of GHG emissions may have been influenced by site-specific conditions, and that other locations may not generate similar results. When referring to the results of the study, the local natural meteorological parameters, soil conditions and DO, Eh, pH, TOC during sampling and monitoring should be taken into consideration, as well as the status of instruments for measuring GHG.

Conclusions
The field experiment has shown that GIRF and OIRF increased the GWP by 41.3% and 98.2% respectively for the entire rice growing period. In IRFF (included GIRF and OIRF), CH 4 emissions from rice fields were mainly related to field fertilization, where fertilizer application can significantly promote CH 4 emissions. Although the activity of frogs in paddy fields can inhibit the emissions of CH 4 , their effect was small. Compared to conventional farming, the introduction of frogs into the rice farming system can reduce CH 4 emissions. The result was primarily attributed to the positive effect of frogs' bioturbation on DO and Eh in the water layer. As a whole, although employing GIRF in rice farming showed a slight increase in GHG emissions, it could be considered as a good strategy in paddy fields for improving the agro-ecological environment and maintaining crop yield. Based on this study, the future plan is to (1) carry out more frog experiments and fertilizer experiments, and select the combination mode with the lowest GWP and better ecological benefits. (2) Then converting ecological benefits into economic benefits by adopting afforestation costs and carbon taxes, and select the experimental combination of best net ecosystem economic benefits.