*Article* **Altered Organic Matter Chemical Functional Groups and Bacterial Community Composition Promote Crop Yield under Integrated Soil–Crop Management System**

**Qi Li 1,2, Amit Kumar 3, Zhenwei Song 4, Qiang Gao 1, Yakov Kuzyakov 5, Jing Tian 2,\* and Fusuo Zhang <sup>2</sup>**


**Abstract:** Sustainable agricultural production is essential to ensure an adequate food supply, and optimal farm management is critical to improve soil quality and the sustainability of agroecosystems. Integrated soil–crop management based on crop models and nutrient management designs has proven useful in increasing yields. However, studies on its effects on the chemical composition of soil organic carbon (SOC) and microbial community composition, as well as their linkage with crop yield, are lacking. Here, we investigated the changes in SOC content, its chemical functional groups, and bacterial communities, as well as their association with crop yield under different farmland management based on four farmland management field trials over 12 years (i.e., FP: farmer practice; IP: improved farmer practice; HY: high-yield system; and ISSM: integrated soil–crop system management). The crop yield increased by 4.1–9.4% and SOC content increased by 15–87% in ISSM compared to other farmland management systems. The increased proportion of Methoxy C and O-alkyl C functional groups with a low ratio of Alkyl C/O-alkyl C, but high Aliphatic C/Aromatic C in ISSM hints toward slow SOC decomposition and high soil C quality. The relative abundances of r-strategists (e.g., Firmicutes, Myxobacteria, and Bacteroidetes) was highest under the ISSM. Cooccurrence network analysis revealed highly complex bacterial communities under ISSM, with greater positive links with labile SOC functional groups. The soil fertility index was the main factor fueling crop yields, as it increased with the relative abundance of r-strategists and SOC content. Our results indicated that crop yield advantages in ISSM were linked to the high C quality and shifts in bacterial composition toward *r*-strategists by mediating nutrient cycling and soil fertility, thereby contributing to sustainability in cropping systems.

**Keywords:** SOC functional group; bacterial community structure; crop yield; soil fertility; integrated soil–crop management system

#### **1. Introduction**

Soil organic matter (SOM) supports multiple ecosystem functions and is closely linked to consistently high crop yields and sustainable agriculture [1,2]. Therefore, maintaining and enhancing SOM content in agricultural fields is important, especially within the context of food security and climate change mitigation [3,4]. Although agricultural land covers 38% of the Earth's land surface, large areas of agricultural land suffer from medium to strong degradation [5]. Therefore, improving soil organic carbon (SOC) is an efficient way

**Citation:** Li, Q.; Kumar, A.; Song, Z.; Gao, Q.; Kuzyakov, Y.; Tian, J.; Zhang, F. Altered Organic Matter Chemical Functional Groups and Bacterial Community Composition Promote Crop Yield under Integrated Soil–Crop Management System. *Agriculture* **2023**, *13*, 134. https:// doi.org/10.3390/agriculture13010134

Academic Editor: Ryusuke Hatano

Received: 18 November 2022 Revised: 30 December 2022 Accepted: 31 December 2022 Published: 4 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to increase C sinks, reduce greenhouse gas (CO2) emissions, and improve crop productivity in carbon-depleted agricultural soils [6,7].

Inappropriate management has depleted 25–75% of SOC in global cropland [8], resulting in a growing interest in optimizing agricultural management. Current optimization techniques include judicious use of fertilizer, adding manure, or returning straw to the field [9–13]. Straw returns maximize the use of natural resources, improve soil structure, and increase the SOM content, providing an ideal environment for crop plantations [14,15]. Organic fertilizers are rich in organic matter and beneficial microorganisms [16]. Many field trials have shown that long-term organic fertilization can increase the conversion of organic and inorganic materials in the soil, which is conducive to improving soil fertility and promoting crop growth and development. [12,17–19].

SOM is chemically diverse and consists of pools with various availability and turnover rates [20,21]. The chemical composition of SOM, which defines its decomposability and stability, is receiving attention in the context of research on soil fertility and quality [22,23]. Intensive management has a direct and indirect impact on the quality and quantity of organic C entering agricultural soils. For instance, organic fertilization generally increases Alkyl C and aromatic C, while decreasing O-alkyl C in soil [23,24]. On the contrary, combined mineral NPK fertilizers (combined with organic manure) increase the O-alkyl group levels but decrease the levels of Alkyl and aromatic groups and the ratio of Alkyl-C/O-alkyl-C (A/O-A) groups [25,26]. The inconsistent compositions of SOM chemical functional groups could be attributable to site-specific soil conditions, input types, climate, and the complexity of the microbial decomposition processes [27–30].

SOC is considered an overarching edaphic factor that shapes bacterial diversity [12,31,32]. This is because heterotrophic soil microorganisms rely on increased SOC to obtain their nutrients [33–35]. After long-term organic matter input, greater SOC is observed, and it is accompanied by greater microbial activity [36,37]. Specifically, when SOC increases, there tends to be a shift in microbial community composition; the community moves toward having more *r*-strategists compared to unfertilized plots [38]. Beyond determining life history strategies, SOC also determines which bacterial phyla are present in the soil [31,32]; this may result from the ability of the microbial groups to utilize SOC pools [39]. Longterm organic amendments (e.g., sewage sludge and chicken manure) increase the dissolved organic carbon content and support an indicative shift of bacterial community to r-strategist taxa such as Proteobacteria and Actinobacteria [40]. This is because the r-strategists have high nutrient requirements and can maximize their reproductive outputs when resources are abundant [31]. The application of organic fertilizers mitigates the negative effects of mineral fertilization on microbial diversity, and the shift of microbial communities toward *r*-strategists is owing to the efficient input of C that characterizes organic fertilizers [36–38]. All these studies suggest that SOC quality is vital for adjusting microbial community composition and life history.

To increase crop production while minimizing impact on the environment, the integrated soil crop system management (ISSM) strategy was recently introduced. The strategy is based on a Hybrid-Maize simulation model and a seasonal root zone nitrogen management strategy (IRNM) to determine the most appropriate combination of planting date, crop density, and fertilizer application rate at the trial site [41]. ISSM has been shown to increase yields (maize, rice, and wheat) by an average of 10.8–11.5% from 2005–2015 [41–43].

Northeast China is considered the "first granary", with grain production accounting for a quarter of the country's total grain production [44]. SOC levels in Northeast China have declined by 22% over the past three decades [2,45], and the decrease is mainly due to intensive cropping and inefficient fertilization [46,47]. Therefore, it is crucial to identify appropriate soil management practices that are conducive to soil C sequestration and sustainable development of agricultural ecosystems in this region. In this study, we aimed to (i) compare the changes in total SOC content, chemical functional groups, and bacterial community composition across various farmland management systems; and (ii) explore changes in SOC quantity and quality in relation to bacterial communities and their impact on crop yield. We hypothesize that (1) compared to farmer practice (FP), ISSM will increase SOC, especially the labile functional pools, because of the direct input of nutrients associated with this strategy [48]; (2) changes in SOC quantity and quality in the ISSM system are important factors shaping microbial community composition, with increasing r-strategist bacteria ascribed to preferentially exploit labile organic compounds [49,50]; (3) changes in SOC quantity and quality impact crop yield via regulating bacterial community composition and enhancing soil fertility [18].

#### **2. Materials and Methods**

#### *2.1. Study Site*

The field experiment was initiated in April 2009 using a maize monoculture, in Zhuling Gong County (43◦12 N, 124◦66 E, 206 m elevation), Ji Lin Province, Northeast China. The site has a temperate continental monsoon climate, an annual precipitation of 595 mm, and an annual mean temperature of 6.3 ◦C. The soil in this region is classified as Phaeozem (equivalent to Hapudoll as per USDA Soil Taxonomy). At the start of the experiment, the topsoil (0–20 cm) had the following basic characteristics: pH, 6.7; SOC, 16 g kg−1; total N (TN), 1.6 g kg−1; available P (AP), 7 mg kg−1; and available K (AK), 151 mg kg−1. Four farmland management systems were established: (1) farmer practice (FP)—set up based on a survey of various production management systems that local farmers are accustomed to applying; (2) improved farmer practice (IP)—a model based on farmer practice by optimizing the combination of existing techniques (e.g., planting density and fertilizer application) to increase crop yields by 15–20%; (3) high-yield system (HY)—a management model intended to increase crop yield by 30%, regardless of the associated environmental costs; and (4) integrated soil–crop system management (ISSM)—a sustainable integrated program to increase crop yield and resource use efficiency by optimizing planting density, fertilizer application, and tillage systems with the aim of achieving an ultrahigh yield while lowering resource and environmental costs.

The field trial design was a randomized complete block design with four farmland management treatments and four replications, for a total of 16 plots. The size of each plot was 6 m × 23 m. For each farmland management system, details of the fertilization rate, planting density, and straw management are shown in Table 1. Maize (ZhongDan99 cultivar) was sown in late April and harvested in early October. The straw was returned to the field based on the amount of maize root and shoot (leaf and stem) residues after above-ground harvesting in each plot, representing 30% of the total dry matter weight.


**Table 1.** The detailed field management information for four farmland management systems.

FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. Straw management: amount of straw returned to the field.

#### *2.2. Soil Sampling and Analysis of Physicochemical Properties*

Five topsoil cores (0–20 cm depth, 5 cm diameter) were collected at random from each plot in October 2020 and pooled together to form one composite sample per plot, for a total of 16 soil samples. After that, the composite samples were passed through a 2 mm sieve to remove any roots or stones. Each soil sample was stored in three parts: one part was air-dried and stored at room temperature for the determination of various chemical properties, e.g., SOC, TN, pH, AP, and AK; the second sample was stored in a refrigerator at

−20 ◦C for the determination of ammonium (NH4 +) and nitrate (NO3 −) content; the third part was stored at −80 ◦C and then subjected to DNA extraction and microbial community analysis. Soil samples were collected from 0–20 cm depths using cutting rings of 100 cm<sup>3</sup> volume to calculate soil bulk density (BD). All maize plants in each plot were collected at the ripening stage, and the crop yield was analyzed as the dry weight of grains.

The SOC and TN contents were determined through combustion using a Vario EL III Elemental Analyzer (Elementar). Soil available phosphorus (AP) and available potassium (AK) were measured according to Hanway and Heidal (1952) and Olsen (1954) [51,52], respectively. Soil pH was measured using a pH electrode (FE20-FiveEasy pH, Mettler Toledo) with a soil-to-deionized-water ratio of 1:2.5 (*w*/*w*). Soil mineral N, including NO3 − and NH4 +, was extracted with 0.01 mol L−<sup>1</sup> KCl solution and analyzed using an auto-analyzer (TRAACS-2000, Bran+Luebbe).

The SOC stock was calculated using the following equation [53]:

$$\text{SOC stocks} \left(\text{Mg} \,\text{C} \,\text{ha}^{-1}\right) = \text{SOC} \left(\text{g} \,\text{kg}^{-1}\right) \times BD \left(\text{g} \,\text{m}^{-3}\right) \times \frac{H(\text{cm})}{100} \tag{1}$$

where *BD* and *H* represent bulk density and soil depth (0–20 cm), respectively.

The soil fertility index was determined using the averaging approach as follows: to prevent the differences in data scale, TN, AP, AK, NO3 −, and NH4 <sup>+</sup> estimated in this investigation were normalized from 0 to 1 using the "max–min" approach, and the soil fertility index was then calculated as follows [54]:

$$
\textit{Fertlikty index} = (-\textit{x}\_{i,min}) / (\textit{x}\_{i,max} - \textit{x}\_{i,min}) \tag{2}
$$

where *xi* is the measured soil properties, *xi*,*min* is the minimum of soil properties *I*, *xi*,*max* is the maximum of soil properties *i*.

#### *2.3. Solid-State 13C NMR Spectroscopy*

The chemical composition of SOM was investigated by determining the relative abundance of functional groups using solid-state 13C cross-polarization and magic-anglespin (CPMAS) NMR. Soil samples were treated with a 10% hydrofluoric acid (HF) solution to concentrate the organic matter and remove paramagnetic minerals [55]. The 13C-CPMAS NMR spectra were acquired using an AVANCE III 400 WB spectrometer (Bruker, Billerica, MA, USA) at 100 MHz for 13C and 400 MHz for 1H with a spinning rate of 5 kHz, an acquisition of 20 ms, a recycle time of 1 s, and a contact time of 1 ms.

Bruker TopSpin (v4.1.1) was used to compare the 13C-CPMAS NMR spectra of different samples; peak areas were calculated and integrated to estimate their relative proportions. The methods outlined by Bonanoni et al. were used to select spectral regions and identify C functional groups (chemical structures) [56]: 0–50 ppm, Alkyl C; 50–60 ppm, Methoxyl C; 60–95 ppm, O-alkyl C; 95–110 ppm, Di-O-alkyl C; 110–145 ppm, Aryl C; 145–160 ppm, Phenolic C; and 160–190 ppm, Carboxyl C. In this study, the Methoxyl C, O-alkyl C, Di-Oalkyl C, and Carboxyl C groups were classified as labile C groups, whereas Alkyl C, Aryl C, and Phenolic C groups were classified as recalcitrant C groups [57,58]. Different indices of SOM stability were calculated as follows [59,60]:

The *A/O-A* (*Alkyl C/O-alkyl C*) represents the degree of decomposition and was calculated as follows:

$$A/\text{O-A ratio} = \frac{Alkyl\text{ C}}{Method\text{ C} + \text{O-alkyl C} + Di\text{-O-allkyl C}}\tag{3}$$

The Aliphatic C/Aromatic C (*Alip/Arom*) indicates the degree of aliphaticity, and was calculated using the following equation:

$$Alip/Arom = \frac{Alkyl\,\,\text{C} + Methoxyl\,\,\text{C} + O\text{-}alkyl\,\,\text{C} + Di\text{-}O\text{-}alkyl\,\,\text{C}}{Ayl\,\,\text{C} + Phnelic\,\,\text{C}}\tag{4}$$

The *HB/HI* (Hydrophobic C/Hydrophilic C) indicates the hydrophobicity degree, and was calculated using the following equation:

$$\text{HB/HI} = \frac{\text{Alkyl } \text{C} + \text{Aryl } \text{C} + \text{Ph} \text{moli} \text{ } \text{C}}{\text{Metbo} \text{xyl } \text{C} + \text{O-Alkyl } \text{C} + \text{Di-O-Alkyl } \text{C} + \text{C} \text{arboxyl } \text{C}} \tag{5}$$

The *Aromaticity* response to soil SOC stability was calculated as follows:

$$Aromaticity = \frac{Aryl \text{ } \text{C} + \text{Phenylic} \text{ C}}{Alkyl \text{ } \text{C} + Methoxyl \text{ } \text{C} + \text{O-alkyl } \text{C} + Di\text{-O-alkyl } \text{C}}\tag{6}$$

#### *2.4. DNA Extraction and Amplicon Sequencing*

DNA was extracted from fresh soil (0.25 g) according to the manufacturer's instructions using a PowerSoil kit (MoBio Laboratories, Carlsbad, CA, USA). The DNA extract was tested on a 1% agarose gel, and the concentration and purity of the DNA were determined using a NanoDrop 2000 UV-vis spectrophotometer (ThermoScientific, Wilmington, USA). The hypervariable region V3-V4 of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5 -ACTCCTACGGGAGGCAGCAG-3 ) and 806R (5- GGACTACHVGGGTWTCTAAT-3 ) using an ABI GeneAmp® 9700 PCR thermocycler (ABI, Carlsbad, CA, USA). The PCR amplification of the 16S rRNA gene was performed as follows: initial denaturation at 95 ◦C for 3 min, followed by 30 cycles of denaturing at 95 ◦C for 30 s, annealing at 55 ◦C for 30 s, and extension at 72 ◦C for 45 s, then single extension at 72 ◦C for 10 min, and a final hold at 10 ◦C. The PCR mixtures contained 4 μL of 5X TransStart FastPfu buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of forward primer (5 μM), 0.8 μL of reverse primer (5 μM), 0.4 μL of TransStart FastPfu DNA Polymerase, template DNA 10 ng, and finally ddH2O up to 20 μL. The PCRs were performed in triplicate. The PCR product was extracted from 2% agarose gel and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer's instructions and quantified using Quantus™ Fluorometer (Promega, Madison, WI, USA).

Raw amplicon sequences were subjected to quality control using the following criteria: (1) the low-quality sequencing reads which had an average quality score of <20 or contained ambiguous nucleotides were filtered using Mohur version 1.31.1 using the UPARSE (version 7.1) software [61]; and (2) only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence. The maximum mismatch ratio of the overlap region was 0.2. Reads that could not be assembled were discarded; sequences with ≥97% similarity were assigned to one operational taxonomic unit (OTU), and chimeric sequences were identified and removed. Taxonomy was assigned to each OTU by the RDP (Ribosomal Database Project) classifier [62]. The final sequencing products contained 43,750 sequences with an average length of 415 bp per sample for downstream analysis.

#### *2.5. Statistical Analyses*

SOC chemical functional groups, bacterial community composition, and other variables were analyzed by SPSS 22.0 (IBM Corp., Armonk, NY, USA) using one-way ANOVA with a randomized group design. Differences were considered significant at *p* < 0.05, and a post hoc least-significant-difference test was carried out to compare the differences among farmland management systems. The normal distribution of data was tested by using the "Shapiro. test" function of the stats package in R v.4.1.2 (R Core Team, 2021).

Changes in bacterial community composition were evaluated by nonmetric multidimensional scaling (NMDS), based on the Bray–Curtis distance calculation method. To further determine significant differences in bacterial community structure for any pair of samples, Adonis analysis was performed using NMDS and statistical analyses were performed using R (version 4.1.2) with the vegan package (R Core Team, 2021).

Linear discriminant analysis (LDA) effect size (LEfSe) was applied to determine if there were significant differences in bacterial taxa among the four farmland management systems [63]. We performed LDA in combination with Kruskal–Wallis (KW) test to identify species with significant differences in abundance between management systems setting log LDA scores > 2.0 (http://huttenhower.sph.harvard.edu/galaxy/, accessed on 16 November 2022). Circus (https://hiplot.com.cn/, accessed on 16 November 2022) was used to aid in the identification and analysis of similarities and differences resulting from phyla relative abundance comparisons [64], which emphasizes the statistical significance and biological relevance [63].

The co-occurrence networks of bacterial communities with SOC chemical functional groups and bacterial co-occurrence networks under each management system were compared by setting the same metrics (dissimilarity threshold for KLD matrix maxima and Spearman correlation threshold of 0.6) and running the networks using all OTU taxa for the four farmland management systems. For each edge and measure, the reciprocal and bootstrap distributions were generated with 100 iterations. The *p*-value of the measure was calculated as the area of the mean of the bootstrap distribution and the standard deviation generated by the mean of the Gaussian curve under the reciprocal distribution. The *p*-values were adjusted using the Benjamini–Hochberg procedure [65,66]. Only the most significant interactions with strong linear connections (r > 0.6) are shown in the network diagram. The nodes in the constructed network represent genus and SOC chemical functional groups, and the edges represent strong and significant correlations between the nodes. Network visualization was conducted using Gephi (version 0.9.2) and Cytoscape (version 3.8.2)(Ideker, 2011). The topological characteristics of the calculated network include positive and negative correlations, nodes, edges, network density, closeness centrality, and betweenness centrality (Tables S4 and S5).

To further elucidate the pathways through which all factors regulate crop yield, partial least-squares path modeling (PLS–PM) and Pearson correlations were performed. Model evaluation was performed based on the goodness of fit (GOF), with a GOF > 0.7 considered an acceptable value. The models were constructed using R (version 4.1.2) (R Core Team, 2021). Ordinary least-squares regression was performed to test the correlation between the factors.

#### **3. Results**

#### *3.1. Soil Organic Matter, Soil Fertility, and Crop Yield*

Compared with that in the FP, the SOC content increased by 15% in IP, 33% in HY, and 87% in ISSM (*p* < 0.05; Figure 1a). Correspondingly, the SOC stock increased by 11%, 31%, and 77% under IP, HY, and ISSM, respectively, compared to that in FP (*p* < 0.05; Figure 1b). Compared to the other three farmland management systems, ISSM increased TN, AP, AK, and NO3 − contents (*p* < 0.05; Table 2). Consequently, ISSM increased the soil fertility (*p* < 0.05; Figure 1c) approximately four times than FP. The highest crop yield was obtained under ISSM, with a 9.4% increase than FP (Figure 1d).

**Table 2.** Effects of management systems on soil physicochemical properties in a 12-year field experiment.


Error bars indicate standard errors of the mean (n = 4). Letters within the same row indicate significant differences among farmland management systems at *p* < 0.05. SOC: soil organic carbon; TN: total nitrogen; AP: available phosphorus; AK: available potassium; NH4 +: ammonium N; NO3 −: nitrate N; BD: soil bulk density. FP: farmer practice; IP: improved farmer practice; HY: high yield system; ISSM: Integrated soil–crop system management.

#### *3.2. Soil Organic Matter Functional Groups Depending on Management Systems*

Across all farmland management systems, the proportion of O-alkyl C (27–30%) was the largest, followed by Alkyl C (24–26%), Aryl C (16–21%), and Carboxyl C groups (7.9–9.2%). These four functional groups dominated SOM, accounting for more than 70% of the total functional group spectrum across all management systems (Table S1). Compared with FP, the ISSM increased the levels of Methoxyl C and O-alkyl C groups by 17% and 8.1%, respectively, and decreased Aryl C group levels by 21% (*p* < 0.05; Figure 2c,d). In general, the ISSM system increased the labile C groups (*p* < 0.05; Figure 2a). The A/O-A was similar in the FP and IP (0.55–0.60) and lower in ISSM (0.52) (Figure 2e). Among all management systems, the highest Alip/Arom was under ISSM (3.4) (*p* < 0.05; Figure 2e).

**Figure 1.** (**a**,**b**) Soil organic carbon content and stock (0–20 cm) under four farmland management systems, (**c**,**d**)soil fertility index, and crop yield after 12 years under four farmland management systems: FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. Letters represent significant differences (*p* < 0.05) among management systems. Error bars indicate standard errors of the mean (n = 4). The fertility index was calculated based on the standardized "max–min" approach of five soil properties.

#### *3.3. Bacterial Community Structure Depending on Management Systems*

The Shannon index of the bacterial community increased more in ISSM than in FP and HY (*p* < 0.05; Figure 3a). Management systems changed the soil bacteria community structure, as indicated by NMDS and Adonis tests (*p* = 0.01; Figure 3b). The pairwise comparison revealed that the soil bacterial community under ISSM was different from that under FP and HY (Table S3). A Mantel test revealed that the major factors shaping bacterial community structure were soil pH (r = 0.74, *p* < 0.05), AP (r = 0.67, *p* < 0.05), AK (r = 0.55, *p* < 0.05), SOC (r = 0.54, *p* < 0.05), TN (r = 0.40, *p* < 0.05), and O-alkyl C (r = 0.34, *p* < 0.05) (Table S4).

**Figure 2.** 13C CPMAS NMR spectra of functional groups of SOC (**a**), and their proportions (**b**–**d**) and ratios (**e**) under four farmland management systems: FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. Letters represent significant differences (*p* < 0.05) among farmer practices. Error bars indicate standard errors of the mean (n = 4). A/O-A, Alkyl C/O-alkyl C; Alip/Arom, Aliphatic C/Aromatic C; HB/HI, Hydrophobic C/Hydrophilic C.

**Figure 3.** Bacterial community profiles in soil under four farmland management systems: FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. (**a**) Bacterial alpha-diversity (Shannon), (**b**) nonmetric multidimensional scaling (NMDS) analysis of bacterial community structure based on Bray−Curtis distance, (**c**) relative abundances of bacterial taxa at the phylum level, (**d**) the linear discriminant analysis effect size (LEfSe) analysis showed the significantly different taxa of bacterial communities. Letters represent significant differences (*p* < 0.05) among four farmland management systems. Error bars indicate standard errors of the mean (n = 4).

Across the farmland management systems, the bacterial community composition was dominated by Actinobacteria (37%), Proteobacteria (26%), and Chloroflexi (10%) (Figure 3c, Table S2). The similarities in the relative abundance of dominant phyla and classes among farmland management systems were further investigated (Figure 3c, Table S2). The relative abundances of *r*-strategists such as Firmicutes increased 4-fold, whereas those of Myxobacteria and Bacteroidetes (phylum level) increased 0.5-fold in ISSM as compared with those of the other three systems (*p* < 0.05; Figure 3c, Table S2). The LEfSe also showed that the relative abundances of Firmicutes and Myxobacteria were the highest in ISSM (*p* < 0.05; Figures 3d and S1) than in other systems. At the class level, the relative abundances of Bacteroidia, Polyangia, Clostridia, Chloroflexi, and Bacilli were higher in ISSM (*p* < 0.05; Figure 3c, Table S2). Further comparisons revealed that the relative abundances of taxa Clostridia (belonging to the phylum Firmicutes) and Paeniclostridium, Romboutsia, and Terrisporobacter (belonging to the class Clostridia) were higher under ISSM than in other farmland management systems (*p* < 0.05; Figures 3c and S1, Table S2).

#### *3.4. Bacterial Network and Linkage with SOM Functional Groups Depend on Management Strategy*

The node and edge numbers of the bacterial co-occurrence network were high in ISSM, indicating a greater complexity of bacterial communities (Figure 4). The positive correlations between nodes were greater in the ISSM network (73%) than in the FP (16%), IP (30%), and HY (29%) networks (Figure 4a,b). ISSM increased the average clustering coefficient (*avg CC*), average closeness centrality (*avg C*), and average betweenness centrality (*avg BC*), but decreased the average path length (GD) and modularity compared to FP (*p* < 0.05; Table S5).

**Figure 4.** Co-occurrence networks of bacterial communities under four farmland management systems: FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. Node colors indicate phyla (**a**) and modularity classes (**b**). The size of each node is proportional to the number of degrees. For visual clarity, co-occurrence networks only show nodes with at least 10 degrees.

Co-occurrence network analysis provided evidence of the correlation between SOC chemical functional groups and bacterial community composition (Figure 5). The network pattern indicated that the correlation between labile C (O-alkyl C, Methoxyl C, and Di-O-alkyl C) and bacterial genera was 1.4 times higher under ISSM than under FP. In

particular, positive correlations increased by 34% and negative correlations decreased by 55% (Figure 5a,d; Tables S6 and S7). Notably, labile C and Aryl C were positively correlated with Firmicutes (Clostridia and Bacilli) in ISSM and FP (*p* < 0.05; Figure 5, Table S7). Additionally, *Paeniclostridium*, *Romboutsia* (Clostridia), *Psychrobacillus*, and *Solibacillus* (Bacilli), positively correlated with labile C (O-alkyl C, Methoxyl C, and Di-O-alkyl C), were stronger in ISSM than in other systems (*p* < 0.05; Figure 5d, Table S2). The relative abundances of *Fluviicola* (Bacteroidetes) and *Anaeromyxobacter* (Myxobecteria) were positively correlated with Methoxyl C, and the relative abundance of *Haliangium* (Myxobecteria) was negatively correlated with Aryl C (*p* < 0.05; Figure 5d, Table S2).

**Figure 5.** Network analysis revealing the associations between bacterial taxa and labile (**a**–**d**) and recalcitrant (**e**–**h**) SOC functional groups. The co-occurrence network was constructed at the genus level and colored according to the phylum classification. Seven functional groups of SOC functional groups are indicated with red triangles. Red and green lines represent strong linear connections (r > 0.6), respectively. FP: farmer practice; IP: improved farmer practice; HY: high-yield system; ISSM: integrated soil–crop system management. Labile C groups: Methoxyl C, O-alkyl C, Di-O-alkyl C, Carboxyl C. Recalcitrant C groups: Alkyl C, Aryl C, Phenolic C.

#### *3.5. Linking SOM Quantity and Quality, Microbial Diversity, and Soil Fertility with Crop Yield*

The relative abundance of Bacteroidetes, Firmicutes, and Myxobacteria were linked to soil fertility and labile C (Methoxyl C, O-alkyl C, and Di-O-alkyl C), which were positively correlated with SOC (*p* < 0.05; Figure 6a). PLS-PM analysis showed that crop yield was directly dependent on soil fertility (path coefficient = 0.58; Figure 6b), which was also supported by the Pearson correlations (Figure 6d). The SOC showed the largest effect on fertility via direct (path coefficient = 0.57) and indirect effects (path coefficient = 0.28) on the relative abundance of *r*-strategists and labile C pools (Figures 6c and S2). Soil fertility and *r*-strategists increased with SOC and labile C content (Figure 6c). There were corresponding strong positive correlations between the relative abundances of *r*-strategist bacteria and the labile C groups (*p* < 0.05; Figures 6d and S3). The effects of labile C and *r*-strategists on crop yield were indirect rather than direct (Figures 6b,c and S2). Specifically, crop yield increases were directly dependent on soil fertility. The relative abundance of *r*-strategists, SOC, and soil fertility positively correlated with crop yield (*p* < 0.05; Figure 6d). The SOC, labile C, and *r*-strategists accounted for 75% of the soil fertility index, and the soil fertility index explained 33% of the variance in total crop yield (Figure 6b).

**Figure 6.** Correlations of crop yield, fertility index, soil properties, SOM chemical structure, and bacterial communities (**a**). Partial least-squares path (PLS−PM) analysis for crop yield, showing the relationships among soil fertility index, SOM chemical structure, and bacterial communities. GOF, goodness of fit (**b**). Black and red solid arrows indicate positive and negative associations, gray dashed arrows indicate insignificant correlations. Standardized total effects of each variable from PLS−PM (**c**). Regression of multiple factors to study the correlation between the factors (**d**). The solid line was fitted using ordinary least-squares regression and the shaded area corresponds to the 95% confidence interval. Significance codes: \*\*\* *p* < 0.001, \*\* *p* < 0.01, \* *p* < 0.05. The fertility index was calculated based on the standardized "max–min" approach of five soil properties. The *r*−strategists are the first axis of the principal component analysis based on the three bacterial taxa (Support Infographic Figure S3). SOC: soil organic carbon; TN: total nitrogen; AP: available phosphorus; AK: available potassium; NH4 +: ammonium N; NO3 −: nitrate N.

#### **4. Discussion**

#### *4.1. Integrated Soil–Crop Management System Increases SOM Quantity and Quality*

The ISSM increased SOC by 15–87% compared to the other three systems (*p* < 0.05; Figure 1a), particularly the labile C functional groups (Methoxy C and O-alkyl C) (Figures 2c and 7). Methoxyl C and O-alkyl C are considered readily degradable components [67]. O-alkyl C groups are usually derived from carbohydrates from fresh plant material; for instance, the anomeric C is derived from cellulose and hemicellulose, whereas Methoxyl C groups are derived from fresh lignin and carbohydrates [68]. The increase in SOM, particularly labile functional C groups, may be caused by the following reasons. First, the high planting density and crop residues in ISSM can provide high amounts of labile C such as O-alkyl C [69–71]. This is consistent with the highest yield being recorded under the ISSM system (Figure 1). Second, the addition of organic fertilizers in this system resulted in a surplus of carbohydrates, including labile C [72]. Third, organic fertilizers stimulate the degradation of Aryl C groups and consequently form new Methoxyl C and O-alkyl C groups [73,74]. The high accumulation of labile C resulted in high SOC stock (Figures 1b and 2c), supporting our first hypothesis. Additionally, a lower A/O-A ratio

in ISSM than in FP (Figure 2e) indicates that the SOC decomposition was slow in ISSM, resulting in more C being accumulated in the soil. These results confirm that the ISSM system reduces CO2 emissions. Consistent with the results estimated by Cui et al. that the CO2 emissions were reduced by 12.9% for the ISSM system-based interventions compared to FP [43]. In addition, the lower A/O-A ratio and higher Alip/Arom (Figure 2e) indicate that the ISSM system not only increases the SOC quantity, but also the quality [75]. In this study, a large amount of organic fertilizer input in ISSM systems has increased the yield and SOC content, but the problem of high nutrient losses to the environment, such as ammonia emissions, nitrate leaching, and denitrification losses, is undeniable. According to previous studies, the comprehensive results showed that ISSM system-based interventions reduced reactive nitrogen losses by 13.3–21.9% and GHG emissions by 4.6–13.2%. The yield-scaled nitrogen footprint averaged 4.6kg reactive nitrogen loss per Mg of maize produced, compared to 6.1 kg Mg−<sup>1</sup> without intervention. Similarly, yield-scaled GHG emissions were 328 kg compared to 422 kg CO2 equivalent per Mg for maize, respectively. Farmer practice applied 300 kg−<sup>1</sup> ha−<sup>1</sup> of N fertilizer, while in ISSM we reduced the N fertilizer application rate to 195 kg ha−<sup>1</sup> [43]. According to unpublished data from our group, the N fertilizer utilization rates of the four farm management systems were 37.3 kg kg<sup>−</sup>1, 59.7 kg kg−1, 39.4 kg kg<sup>−</sup>1, and 62.7 kg kg−1, respectively, showing that ISSM had the highest N fertilizer utilization rate and relatively reduced N losses. Therefore, ISSM systems are beneficial to C sequestration and environmental protection, and thus may serve as a sustainable farmland management practice.

#### *4.2. Integrated Soil–Crop Management System Links Labile SOC with r-Strategists*

Bacterial community composition depends on a set of environmental factors [76]. As shown in Figure 6, the relative abundances of Bacteroidetes, Firmicutes, and Myxobacteria were highly correlated with soil pH, SOC, TN, AP, and AK (Table S5). Soil pH is a vital edaphic factor affecting the diversity and composition of soil bacteria in agricultural and forest soils [77–79]. In addition to SOC and pH, other edaphic factors (e.g., TN and AK contents) strongly influenced bacterial diversity (Figure 6a); this is congruent with previous studies [80–83].

The bacterial co-occurrence network was altered by management systems. The high bacterial network complexity in ISSM (Figure 4, Table S5) is owing to the availability of a diverse range of resources (e.g., increased availability of soil carbon and nutrients, Table 2 and Figure 1a). This is in agreement with the results of studies indicating that increased soil microbial network complexity is tied to resource availability and diversity [84–86]. In general, the increased positive correlation and the reduction in the average path lengths in the ISSM co-occurrence network (Figure 4) were indicative of efficient interactions between microorganisms, promoting a denser co-occurrence network [87].

The phyla Firmicutes and Bacteroidetes prefer substrates rich in available C [31,80]. The relative abundance of Firmicutes, Myxobacteria, and Bacteroidetes increased 0.5–4 times more in ISSM than in FP (Figure 3c, Table S2). The higher relative abundance of *r*-strategists was due to the 25% increase in labile C functional groups (Figure 6d). In particular, *Paeniclostridium*, *Romboutsia* (Firmicutes), and *Bacteroidetes* were positively correlated with labile C (Methoxy C and O-alkyl C) (Figures 5 and 6a), confirming our second hypothesis. The Bacillus (Firmicutes) dominates microbial communities under labile C and available P inputs [88,89]. The relative abundance of *r*- (copiotrophs) strategists (e.g., *Bacteroidetes*) increased when the available substrate content in the soil was high (Figure 6a,d). Therefore, the enrichment of *r*-strategists under ISSM was influenced by substrate efficiency and C quality (Figure 6b). Microorganisms considered "opportunists" (*r*-strategists) preferentially exploit less-complex organic compounds [12,31,32]. This is also supported by the positive correlations between labile C pools and abundances of *r*-strategies microbial groups (Figure 6d).

**Figure 7.** Concept of the effects of farmland management systems on the SOC content and compositions, and bacteria community structure and yield of maize. La: labile C; Re: recalcitrant C; BAC: bacterial community; SOM: soil organic matter; TN: total nitrogen; AP: available phosphorus; AK: available potassium; NO3 −: nitrate N.

#### *4.3. SOM Quantity and Quality Increased Crop Yield by Regulating Bacterial Community Composition and Enhancing Soil Fertility*

Soil quality is the key to achieve high crop productivity, and soil quality is tied to SOC content and composition [90]. ISSM improved soil fertility by increasing the SOC content and abundance of *r*-strategists (Figure 6b,c), supporting the improvement of soil quality, which is more commonly associated with a transformation of bacterial community structure [91]. Therefore, the bacterial community structure is an effective biomarker for assessing the environmental impact of a wide range of agricultural practices on soil quality [92]. Specifically, significant positive correlations between SOC and soil chemistry were observed (e.g., TN, AP, AK, and NO3 −), which is suggestive of the coupling between C sequestration and nutrient cycling, a link that explains the increase in soil fertility under ISSM (Figure 6a). Previous studies indicate that combining chemical fertilizers with organic matter such as straw and manure is an effective way to maintain soil productivity [93,94]. Congruent with these results, mixing NPK fertilizer with manure resulted in yields 25% higher than that when the NPK fertilizer was used alone. This result is directly due to enhanced soil nutrients, which leads to increased bacterial diversity [95].

Soil fertility and crop yield were significantly positively correlated (Figure 6b). Moreover, crop yield directly depended on soil fertility, and SOC and labile C indirectly affected crop yield by mediating *r*-strategists to enhance soil fertility (Figure 6b,c). Our results demonstrate that bacterial communities in carbon-rich soil primarily consist of r-strategist bacteria. These microorganisms are involved in soil nutrient cycling, which can regulate soil fertility and indirectly improve crop yields. The decomposition of plant residues and the soil's own organic matter in the soil is a biochemical process, through which carbon is returned to the atmosphere in the form of CO2; and nitrogen, phosphorus, sulfur and trace elements are released into the soil in inorganic form for use by higher plants; some nutrients are assimilated by soil microorganisms into microbial biomass and participate in the rapid turnover process of soil microorganisms. The long-term substitution of straw and manure not only accelerated nutrient cycling, but also increased soil quality and crop yields by increasing bacterial diversity and changing the composition and chemical functional groups, which equally support our results [3,96]. Thus, reasonable agricultural management practices to improve soil fertility and C quality and further improve microbial diversity can help maintain sustainable crop production in Northeastern China [97–99].

#### **5. Conclusions**

Integrated soil–crop system management changed SOC chemical functional group levels, bacterial community structures, and their relationships with crop yield. ISSM enhanced SOC content and improved its quality by increasing labile Methoxyl C and O-alkyl C groups, resulting in a high Alip/Arom ratio. Optimization of farm management systems promotes the correlations between SOC functional groups and bacterial communities, particularly the labile Methoxyl C and O-alkyl C, with Firmicute. SOC content is one of the determinants of soil fertility, while r-strategists and labile C increase with SOC. Crop yield depends directly on soil fertility, while SOC and labile C indirectly affect crop yield by mediating *r*-strategists to enhance soil fertility. In conclusion, ISSM changed SOC quantity and quality and altered microbial community. Overall, the combined effects of these factors culminated in improved soil quality and productivity. Therefore, it is essential to optimize management strategies to maintain soil quality and long-term sustainability of agricultural production.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/agriculture13010134/s1, Table S1: Effects of farmland management systems on functional groups of SOC. Table S2: Relative abundance of soil bacteria at the phyla and class level. Table S3: Results of Adonis test for bacterial communities under four farmland management systems. Table S4: Spearman's rank correlation (R-value) between soil chemical properties and chemical structure of organic carbon and bacterial communities based on Mantel-test. Table S5: Topological properties of the empirical molecular ecological networks (MENs) of bacteria communities at four farmland management systems. Table S6: Topological properties of network of between bacterial taxa and SOC functional groups parameters of soils. Table S7: Detailed network characteristics of interactions between bacterial taxa and SOC functional groups parameters of soils. Figure S1: Linear discriminant analysis (LDA) to identify the taxa leading to community differences for bacteria. FP: farmer practice; IP: improved farmer practice; HY: high yield system; ISSM: Integrated soil-crop system management. Figure S2: Standardized effects of each variable from the partial least squares path analysis (PLS-PM). (a) and (b) represent the standardized direct effects to crop yield and fertility index, (c) and (d) represent the standardized indirect effects. The fertility index was calculated based on the standardized "Max-Min" approach of five soil properties. SOC: soil organic carbon. Labile C: Methoxyl C, O-alkyl C, Di-O-alkyl C, Carboxyl C. The relative abundances of Firmicutes, Myxobacteria, and Bacteroidetes were chosen to *r*-strategists by principal component analysis. Figure S3: Principal component analysis (PCA) analyzing. The scores of the first PCA axis were used as *r*-strategists bacteria. The relative abundances of Firmicutes, Myxobacteria and Bacteroidetes were chosen to *r*-strategists bacteria principal components for analysis. Shaded area is 95% of the confidence interval.

**Author Contributions:** Conceptualization, J.T. and Q.L.; methodology, validation, and writing original draft preparation, visualization, Q.L.; writing—review and editing, visualization, J.T, A.K. and Y.K.; resources, Z.S. and Q.G.; supervision, funding acquisition, F.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (2022YFD1901300), "Research and Application of Key Technologies for Cultivation of Arable Land Quality and Green Development of Agriculture" program (Z191100004019013), the National Natural Science Foundation of China [grant nos. 32071629], and Agriculture Carbon Neutral Account Establishment Program in Quzhou (202127). We would also like to thank the Beijing Science and Technology Program, Beijing Advanced Disciplines, and the RUDN University Strategic Academic Leadership Program.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare that they have no conflicts of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Guofeng Wang 1, Pu Liu 1, Jinmiao Hu <sup>1</sup> and Fan Zhang 2,\***


**Abstract:** Explaining the methane emission pattern of Chinese agriculture and the influencing factors of its spatiotemporal differentiation is of great theoretical and practical significance for carbon neutrality. This paper uses the IPCC coefficient method to measure and analyze the spatial and temporal differentiation characteristics of agricultural methane emission, clarify the dynamic evolution trend of the kernel density function, and reveal the key influencing factors of agricultural methane emission with geographical detectors. The results show that China's agricultural methane emissions showed a first increasing and then declining trend. Agricultural methane emissions decreased from 21.4587 million tons to 17.6864 million tons, with an upward trend from 2000 to 2005, a significant decline in 2006, a slow change from 2007 to 2015, and a significant decline from 2015 to 2019. In addition, the emissions pattern of the three major grain functional areas is characteristic; in 2019, agricultural methane emissions from main producing area, main sales area, and balance area were 10.8406 million tons, 1.2471 million tons, and 5.599 million tons, respectively. The main grain producing area is the main area of methane emissions, and the emission pattern will not change in the short term. The variability of grain functional areas is the decisive factor for the difference in agricultural methane emissions. The state of industrial structure is the key influencing factor for adjusting the spatial distribution—the explanatory power of the industrial structure to the main producing areas reached 0.549; the level of agricultural development is the most core influencing factor of the spatial pattern of the main grain sales area—the explanatory power reached 0.292; and the level of industrialization and the industrial structure are the core influencing factors of the spatial pattern of the balance area—the explanatory power reached 0.545 and 0.479, respectively.

**Keywords:** agriculture; methane emission; spatiotemporal pattern; kernel density; influencing factors

#### **1. Introduction**

Climate change has become the biggest threat to global sustainable development [1–4]. The Paris Agreement stipulates that the state parties should keep the global average warming to 2 degrees Celsius higher than the pre-industrial revolution level and strive to limit it to 1.5 degrees Celsius [5]. The "climate critical points" of 1.5 and 2 degrees Celsius are the key time points for catastrophic climate events but are also related to the threat of human survival, which may produce a chain reaction once exceeded [6]. In addition, frequent occurrence of global extreme climate caused by climate change has brought significant economic and social losses to human beings [7]. How to mitigate and adapt to climate change has become a major issue for countries all over the world. Carbon peaking and carbon neutrality is a solemn commitment of China to address global climate change, and a major declaration of the transformation of social and economic development mode [8]. To achieve carbon peaking by 2030 and strive to carbon neutrality by 2060 is not only an essential part of China's high-quality growth but also a reform requirement for major adjustments in different industries [9].

**Citation:** Wang, G.; Liu, P.; Hu, J.; Zhang, F. Spatiotemporal Patterns and Influencing Factors of Agriculture Methane Emissions in China. *Agriculture* **2022**, *12*, 1573. https://doi.org/10.3390/ agriculture12101573

Academic Editors: Yinglong Chen, Masanori Saito and Etelvino Henrique Novotny

Received: 29 August 2022 Accepted: 24 September 2022 Published: 29 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Methane is the second largest greenhouse gas responsible for climate change after carbon dioxide [10]. Strengthening methane emissions reduction has become a necessary item in the 21st century [11]. In 2020, global methane emissions were 570 million tons, of which human activities produced 340 million tons. IPCC pointed out in Working Group I of the Sixth Assessment Report that within 20 years after emission, the greenhouse effect of 1 ton of methane is comparable to 84 tons of carbon dioxide and its warming effect is still 28 times that of carbon dioxide even after 100 years [12]. In addition, since it is easier to reduce methane than carbon dioxide, the International Energy Agency (IEA) notes that 75% of methane leaks in the global oil and gas industry can be controlled with existing technologies, and 50% of methane emissions are reduced at net zero cost. [13]. Methane has a shorter lifetime in the atmosphere than carbon dioxide and can be a priority for emissions reduction. Emission reductions of non-carbon dioxide greenhouse gases such as methane are essential for the goal of controlling global warming to 1.5 ◦C by the end of this century [14–16].

As agriculture is the main source of methane emissions, it is of great significance to clarify its spatial distribution and influencing factors [17]. Methane emissions from agriculture account for about 1/5 of the total global emissions, and the main sources of agricultural methane emissions are farming and animal husbandry. The carbon emissions produced by the planting industry are mainly generated by paddy planting, while the carbon emissions produced by the breeding industry are mainly generated by the intestinal fermentation and manure management of livestock and poultry [18].

Scholars from all over the world have studied agricultural methane emissions from different perspectives. Some scholars have estimated the intestinal methane emissions of livestock and poultry in the Southern African Development Community (SADC) [19]. Some scholars have studied methane emissions from wetland rice fields [20]. In addition, some researchers have compared the intestinal methane emissions of livestock and poultry in Germany in the late 19th century with the current situation [21]. As a major agricultural country, China's methane emissions from agricultural activities account for 50.15% of the country's total emissions. Agricultural methane emission patterns and influencing factors have become a hot topic in academic research [22–27]. After estimating the methane emissions of various regions in China, it is indicated that the total methane emissions from agricultural interaction in 2018 were 18.22 million tons, among which the intestinal fermentation emissions from livestock and poultry were the largest, accounting for 50.69%, the emissions from paddy were 35.17%, and the emissions from livestock and poultry management were 14.14% [28]. By measuring methane and nitrous oxide emissions from 1980 to 2018, Li Yang pointed out that the methane missions generated from agriculture increased from 0.56 × 109 CO2-eq to 0.73 × <sup>10</sup><sup>9</sup> CO2-eq, which were mainly affected by efficiency factor, structure factor, and population size factor [29]. In addition, scholars' research focused on the provincial level, and calculated the changes in Jiangxi, Guangdong, Heilongjiang, and other provinces. Some scholars have also measured methane emissions from lake farms and their response to ecological restoration [30], which provided a solid foundation to study agricultural methane emissions. However, at present, there are relatively few estimates based on the national level, and there is a lack of comparative studies on the differences among main grain producing area, main sales area, and balance area; thus, relevant influencing factors need to be further clarified. In this paper, the spatial and temporal patterns of methane emission in China were analyzed by using the IPCC method and the calculation method of the Provincial Greenhouse Gas Inventory Guide. The influencing factors were calculated by using the analysis method of the geographic detector, and the emission reduction countermeasures and suggestions were put forth to provide some references for methane emission reduction in China.

#### **2. Materials and Methods**

#### *2.1. Measurement Method of Agricultural Methane Emission*

The CH4 emissions of the national agricultural system mainly come from paddy fields and livestock and poultry manure management and intestinal fermentation [28]. Among them, the CH4 emissions from paddy fields generally follow the basic methods and requirements determined by the IPCC guidelines, and the formula is calculated as:

$$E\_{CH\_4ric\varepsilon} = \sum EF\_{\hat{1}} \times AD\_{\hat{1}} \tag{1}$$

where *ECH*4*rice* represents CH4 emissions from paddy fields (×104 t), *EF*<sup>i</sup> represents methane emission factors from paddy fields (kg/hm2), and *AD*<sup>i</sup> represents the sown area of this type of methane emission factor (×103 hm2),The methane emission factors of paddy fields are listed in Table 1.

**Table 1.** Methane emission factors from paddy fields.


The calculation formula of CH4 generated by livestock and poultry manure management is calculated as follows:

$$E\_{CH\_4, numure, i} = EF\_{CH\_4, numure, i} \times AP\_i \times 10^{-7} \tag{2}$$

where *ECH*4,*manure*,*<sup>i</sup>* refers to the amount of methane produced by manure management of species *<sup>i</sup>* (×104 t), *EFCH*4,*manure*,*<sup>i</sup>* refers to manure management methane emission factor of species *i*, and *APi* refers to the number of species *i*, the methane emission factors for manure management are listed in Table 2.


**Table 2.** Methane emission factors from manure management.

The calculation formula of CH4 produced by intestinal fermentation of livestock and poultry is calculated as follows:

$$E\_{\text{CH4,metric},i} = \sum EF\_{\text{CH4,carrier},i,j} \times AP\_i \times R\_j \times 10^{-7} \tag{3}$$

where *ECH*4,*enteric*,*<sup>i</sup>* is the amount of methane produced by intestinal fermentation of species *<sup>i</sup>* (×104 t), *EFCH*4,*enteric*,*i*,*<sup>j</sup>* is methane emission factor from intestinal fermentation of livestock and poultry of species I, *APi* is the number of species I, and *Rj* is the breeding proportion of this livestock and poultry, The methane emission factor from enteric fermentation of livestock and poultry are listed in Table 3.


**Table 3.** Methane emission factor from intestinal fermentation of livestock and poultry.

#### *2.2. Kernel Density Function*

The kernel density function, due to its non-parametric estimation power of the probability density, can be used to measure the distribution form of random variables, with the expression of:

$$f(\mathbf{x}) = \frac{1}{nh} \sum\_{i=1}^{n} K(\frac{\mathbf{x} - \mathbf{x}\_i}{h}) \tag{4}$$

where *n* represents the number of the sample, *h* represents bandwidth, *h* = 0.9*SN*4/5 (*N* represents the number of the sample and *S* represents sample standard deviation). *K*( *<sup>x</sup>*−*xi <sup>h</sup>* ) represents kernel density function, and the Epanechnikov kernel density form is used. The dynamic evolution of methane emission zones can be reflected by the distribution interval, morphology, and kurtosis extension of kernel density function. If the function presents a "single peak" as a whole, there are no multiple equilibrium states; if the function appears a "double peak" or "multi-peak" state, there are two or more equilibrium points.

#### *2.3. Geographic Detector*

To deeply explore the causes of regional differences in agricultural methane emissions, geographic detector analysis is used in this study. Geographic detectors can not only test the spatial differentiation of the delay variables but also excavate the influencing factors of spatial differentiation, explain the decisive role of influencing factors, and provide panoramic display for the exploration of the spatial differences. In this study, the factor detection method is adopted to determine the fundamental factor of spatial differentiation of agricultural methane emissions through the determinant force index *q*, in which *Y* represents the agricultural methane emission, *X* = {*Xm*} represents influencing factors,(*m* = 1, 2, . . . *L*;), *L* refers to the number of partitions of factor *X*, and *Xm* refers to the different partitions of factor *X*. The determinant force index *q* of factor *X* on *Y* is as follows:

$$\begin{aligned} q &= 1 - \frac{\frac{L}{m} N\_m \sigma\_m^2}{N \sigma^2} = 1 - \frac{SSW}{SST} \\ SSW &= \sum\_m^L N\_m \sigma\_m^2 \\ SST &= N \sigma^2 \end{aligned} \tag{5}$$

where *N* is the number of provinces in the study, *Nm* is the number of provinces contained in the *m*-th partition of factor *X*, *σ*<sup>2</sup> is the variance of region *Y*, *σ*<sup>2</sup> *<sup>m</sup>* is the variance of the driver *X* in the m subdomain, and *SSW* and *SST* represent the sum of variance in each grain function and the total variance in the whole region, respectively. In general, the larger the *q*-value of factor *X*, the stronger the factor of driving force for the spatial analysis. When the driving factors *X* and *Y* have driving effects on each other, the sum of the variances within regions are usually less than the sum of the variances between regions. In order to compare whether the cumulative variance of different grain production regions is significantly different from the overall variance of the whole region, the F-statistic test is introduced in this paper.

#### *2.4. Data*

In this paper, the paddy fields area, livestock and poultry quantity, and breeding proportion are from the China Statistical Yearbook (2000–2019), State Administration of Grain, China Agriculture Yearbook, and China Animal Husbandry and Veterinary Yearbook (2000–2019). Methane emission factors in different paddy fields, from livestock and poultry manure management in different regions, and from livestock and poultry intestinal methane emission factors under different feeding methods are from the Provincial Guidelines for Greenhouse Gas (National Development and Reform Commission, 2011), and some missing data are supplemented by a linear interpolation method.

#### **3. Results**

#### *3.1. Measurement of Agricultural Methane Emission in China* 3.1.1. Overall Agricultural Methane Sequential Variation

Agricultural methane emissions showed a trend of first increasing and then decreasing. From 2000 to 2019 (Figure 1), agricultural methane emissions decreased from 21.46 million tons to 17.69 million tons, with an average annual drop of 0.88%. Overall, China's green development in agriculture has been significantly improved. From 2000 to 2005, agricultural methane emissions continued to increase, with a significant decline in 2006.

**Figure 1.** Agricultural methane emissions in China, 2000–2019.

From 2006 to 2015, there was a slow change period, with a slight increase in methane emissions, and since 2015, methane emissions have fallen significantly. From the perspective of evolution, the No.1 central document in 2004 "opinions on several policies of promoting farmers' income" was put forward after 18 years, and its main body returned to "the field of agriculture, rural areas and farmers". The document put forward the "three subsidies" of direct subsidies for farmers, subsidies for improved varieties, and subsidies for agricultural machinery. It also proposed some measures to reduce the agricultural tax burden, which effectively enhanced the enthusiasm of farmers to grow grain and prompted significant improvement in the level of green development in 2006. In 2015, in order to cope with the problems such as the "ceiling" of prices, the rise of the "floor" of production costs, and the intensified challenges of the "hard constraints" of resources and the environment, the General Office of the State Council of the People's Republic of China issued "the Opinions on Accelerating the Transformation of Agricultural Development Mode", which provided an important starting point for promoting agricultural development on quantity, quality, and efficiency. Under the guidance of the document, the structural adjustment to the agriculture led to a further decline in methane emissions. In terms of the sources of methane emissions, from 2000 to 2019, emissions from paddy fields fell from 6.51 million tons to 6.29 million tons, emissions from livestock and poultry manure decreased from 2.72 million tons to 1.88 million tons, and emissions from livestock and poultry intestinal fermentation dropped from 12.23 million tons to 9.53 million tons (Figure 2). The decline

of livestock and poultry intestinal emissions played a crucial role in reducing agricultural methane emissions.

**Figure 2.** Agricultural methane emissions percentage in China, 2000–2019.

#### 3.1.2. Spatial Analysis of Different Production Areas

Considering the differences between agricultural development among different regions in China, the analysis between spatial regions is more targeted. Based on the classification of agricultural grain functional areas, the whole country is divided into main producing area, main sales area, and balance area.

Agricultural methane emissions in the main producing area, main sales area, and balance area showed a downward trend, but the changes had their own characteristics (Figure 3). The main producing area was the region with the highest agricultural methane emissions. From 2000 to 2019, the main producing area decreased the most, from 12.79 million tons to 10.84 million tons. From the perspective of the planting industry, the sown area of double-cropping early rice and double-cropping late rice in the main producing areas decreased by 27.61% and 31.83%, respectively, but the sown area of single-crop rice increased by 48.38%, which led to an increase in methane emissions from paddy fields in the main producing areas. From the perspective of the breeding industry, although the methane emissions from paddy fields in the main producing areas increased to a certain extent, an increase of 14.36%, the methane emissions from livestock and poultry manure management and enteric fermentation decreased by 507,500 tons and 2,028,000 tons, respectively. From 2000 to 2019, the number of livestock and poultry breeding dropped significantly, and the number of non-dairy cows, goats, pigs, poultry, horses, camels, and donkeys decreased by 43.83%, 20.43%, 29.04%, 16.78%, 62.97%, 88.65%, and 70.96%, respectively. The number of poultry breeding fell the most. The main production areas bear the heavy responsibility of China's food security, and the decline in emissions shows that the industrial structure has been improved, to a certain extent, and agricultural green production has gradually improved. From the perspective of planting industry, the sown area of double-cropping early rice and double-cropping late rice in the main producing areas decreased by 27.61% and 31.83%, respectively, but the sown area of single-crop rice increased by 48.38%, which led to an increase in methane emissions from paddy fields in the main producing areas. From the perspective of the breeding industry, although the methane emissions from paddy fields in the main producing areas increased to a certain extent, with an increase of 14.36%, the methane emissions from livestock and poultry manure management and enteric fermentation decreased by 507,500 tons and 2,028,000 tons, respectively. From 2000 to 2019, the number of livestock and poultry breeding dropped significantly, and the number of non-dairy cows, goats, pigs, poultry, horses, camels, and

donkeys decreased by 43.83%, 20.43%, 29.04%, 16.78%, 62.97%, 88.65%, and 70.96%, respectively, and the number of poultry breeding fell the most. The main production areas bear the heavy responsibility of China's food security, and the decline in emissions has shown the improved industrial structure and increased agricultural green production. Methane emissions in the balance area fell from 6.29 million tons to 5.60 million tons. In the balance area, methane emissions from rice fields, livestock manure management, and livestock intestinal fermentation all decreased by 23,600 tons, 169,800 tons and 285,400 tons, respectively. The balance area is located mainly in Central and Western China, and the change in emissions was partly due to the implementation of China's western development strategy. In addition, methane emissions in the main sales area dropped from 2.37 million tons to 1.25 million tons, decreasing by 47.48%. The methane emissions from paddy fields, livestock manure management, and livestock enteric fermentation in the main sales area decreased by 41%, 48.14%, and 61.25%, respectively. Except for dairy cows, the number of livestock and poultry has decreased significantly. This was mainly due to the fact that the main sales areas, such as Beijing, Tianjin, and Shanghai, are relatively developed areas in China. With limited land and other planting resources and a large food shortage, it is necessary for these developed areas to rely on other provinces to transport food, forcing the region to adjust the industrial structure and objectively promoting the reduction of agricultural methane emissions. The decline in methane emissions in main producing and sales regions was the main reason for the overall decline. Therefore, the main producing areas have high methane emissions and are key areas for the green transformation and development of agriculture.

From 2000 to 2019, the emissions of the main producing areas, main sales areas, and balanced areas were all trending downward, of which the largest decline rate was the main sales area (47.48%), and the largest decline amount was the main producing area (1.9539 million tons). Although they all showed a downward trend, the changes had their own characteristics. For example, the main grain producing areas were the areas with the highest agricultural methane emissions which have dropped by approximately 15.27% in the past 20 years, with an average annual decline of 0.83%, of which 2000–2005 and 2010–2015 were the rising periods, with increases of 10.78% and 2.34%, respectively; 2005–2010 and 2015–2019 were the decreasing periods, and the decline was relatively large, with 13.14% and 13.96%, respectively. Although the total amount of emissions had a downward trend, the discharge of paddy fields had an upward trend, and the proportion of emissions was also different from the total amount of emissions. 2000–2005 and 2015–2019 were decreasing periods, 2005–2015 was the rising period, where the total amount increased from 4.0507 million tons in 2000 to 4.6324 million tons in 2019, and the proportion of emissions increased from 31.66% to 42.73%. Manure management emissions and intestinal fermentation emissions were on a downward trend, decreasing by approximately 30.02% and 28.75%, respectively. In short, the emission structure of the main producing areas was optimized, where livestock and poultry intestinal emissions were the main source of emissions; the proportion of paddy field emissions and intestinal emissions was close, with 42.73% and 46.36%, respectively in 2019. Paddy field management will continue to be an important measure to reduce emissions in major producing areas. Agricultural methane emissions from the main sales areas were the lowest, and although the agricultural methane emissions from paddy fields were not the highest among the three grain divisions, the ratio of emission sources in the main sales areas was the highest, from 58.67% in 2000 to 65.91% in 2019. In addition, 2000–2005 was a decline period, reducing from 58.67% to 52.99%, followed by an increase period in 2005–2019, with a growth of 8.91%. The trend of changes in the proportion of manure management emissions and intestinal emissions was similar, showing the trend of first rising and then declining, of which the total amount of manure management emissions did not change much in 2000–2010, but the proportion of emissions increased from 14.42% to 19.32%. Intestinal fermentation emissions fell by 23,500 tons between 2000 and 2005, but the proportion of emissions increased from 26.91% to 29.99%, and then in 2019, the proportion of emissions fell to 19.85%. Overall, intestinal fermentation

emissions fell by 391,300 tons, with the proportion decreasing by 7.05%. The total decline (691,200 tons) and the decline rate (10.99%) in the balance area in 20 years were the lowest, similar to the change trajectory of the main producing areas where the increase periods were 2000–2005 and 2010–2015, and the decline periods were 2005–2010 and 2015–2019. The emission sources were similar to those in the main producing areas, with intestinal fermentation emissions as the main emission source, but the change amplitude of the three emission sources was not high. In addition, it should be noted that the proportions of intestinal emissions and paddy field emissions were different from the change trend of the main producing areas and the main sales areas; the change range from 2000 to 2019 was between 3.81% and −2.12%, and the three emission sources were reduced by 236,000 tons, 169,800 tons, and 285,400 tons, respectively, from which can be seen that the emission reduction effect of the balance area was not good. In addition, in 2000, the main producing area and the main sales area were 5.39 times and 2.65 times, respectively, the total emissions of the balanced area and by 2019, they had developed into 8.69 times and 4.49 times, respectively; this, to a certain extent, also reflects the widening of spatial differences.

#### 3.1.3. Spatial Analysis of Different Provinces

In view of the different emphases of industrial structure and economic development in each province, the changes in agricultural methane emissions in each province are analyzed at the provincial level (Figure 4). In 2000, the top provinces in agricultural methane emissions were Hunan, Henan, Sichuan, Guangxi, and Shandong, and by 2010, the top provinces were Hunan, Sichuan, Henan, Inner Mongolia, and Guangxi. The emissions of Shandong decreased from 5.70% to 3.90%, and the emissions of Inner Mongolia increased from 3.27% to 5.68%. In addition, Heilongjiang rose significantly from 3.19% to 5.01%. In 2019, the top emission provinces were Hunan, Sichuan, Inner Mongolia, Heilongjiang, and Jiangxi. The reason for the highest agricultural methane emissions in Hunan was that it is a large paddy production and marketing province in China, with the paddy planting area of 3.8556 million hm2 in 2019. The large-scale breeding rate of cattle and sheep in Sichuan Province was low, and the proportion of large-scale pig breeding in the first half of 2020 (53.8%) exceeded the free-range ratio (46.2%) for the first time, resulting in high intestinal methane emissions from livestock and poultry; however, intensive production is gradually becoming the mainstream of production, and greening has been further improved. In 2018, Henan guided farmers to implement the plan of advocating intensive agricultural production, and at the end of 2018, 980,000 livestock and poultry free-range households had been "withdrawn from the village", which caused agricultural methane emissions in Henan to decrease by 2.31% from 2010 to 2019. The number of sheep raised in Inner Mongolia was much higher than that of other provinces; the large-scale breeding rate needs to be improved, which was the main reason for the high emission of agricultural methane. Therefore, scale is the only way for the transformation and development of animal husbandry. It is impossible for retail farming to completely withdraw, and financial means should be strengthened to support and promote new industrial models. For example, building a comprehensive platform for investment funds and technologies by industry experts and attracting retail investors to join the company can not only reduce the risk and cost of retail investors but also increase the rate of large-scale breeding and effectively reduce methane emissions. The agricultural methane emissions of Tianjin, Beijing, Shanghai, and other Eastern developed provinces were relatively small and remained at the bottom of the emission proportion in 2000, 2010, and 2019. The reason lay in the composition of industrial structure. For example, in 2019, the proportion of the first industry in the GDP of the three provinces was 1.3%, 0.3%, and 0.3%, respectively. The number of livestock and poultry in Ningxia and Hainan was low; in 2019, the number of livestock and poultry was 20.707 million and 63.913 million, respectively, far lower than the national average of 233.2885 million, so the methane emission was low. Liaoning and Jilin, also major grain-producing regions in Northeast China, had relatively low emissions, ranking lower (20th and 21st, respectively, in 2019). Compared with Heilongjiang, the two provinces had

a lower paddy sown area, with only 22% and 13%, respectively, than Heilongjiang in 2019. In addition, among the main producing areas, Hebei accounted for 4.05% of emissions in 2000, and the number of livestock and poultry decreased by 37.26% in 2010 compared with 2000, which was the main reason for its low emissions, ranking 18th in 2019.

**Figure 3.** Agricultural methane emissions by region in China, 2000–2019. (**a**) 2000; (**b**) 2005; (**c**) 2010; (**d**) 2015; and (**e**) 2019.

From the perspective of emission sources, in 2019, the agricultural methane emission in 12 of the 31 provinces was dominated by paddy field emissions, represented by Jiangsu (82.73%), and the remaining 19 provinces were dominated by intestinal fermentation emissions, represented by Qinghai, Gansu, Tibet, Xinjiang, and Inner Mongolia, which accounted for more than 90%, and the intestinal fermentation emissions caused by them reached 18.57% of the national total emissions.

#### 3.1.4. Dynamic Evolution of Spatial Differentiation

To effectively reflect the spatial characteristics of agricultural methane emissions, the dynamic distribution of 2000–2019 and 2010–2019 was analyzed by non-parametric kernel density. As can be seen from Figure 5, the kernel density showed two peaks in the two time periods, and then decreased continuously. The distribution curve showed the right tail, the peak moved significantly to the right, and the gap of provincial agricultural methane emissions showed a widening trend.

#### *3.2. Influencing Factors of Agricultural Methane Emission in China* 3.2.1. Model Building of Influencing Factors

Referring to relevant references, the factors affecting agricultural methane emissions include agricultural development level, residents' living conditions, industrial structure, urbanization level, industrialization level, environmental regulation, agricultural machinery input, and agricultural science and technology R&D investment, etc.

(**a**)

**Figure 4.** *Cont*.

(**b**)

(**c**)

**Figure 4.** Agricultural methane emissions by provinces in China, 2000–2019. (**a**) 2000; (**b**) 2010; and (**c**) 2019.

**Figure 5.** Epanechnikov kernel density function diagram.

Agricultural development level. There is a double-sided effect between agricultural development level and agricultural methane emissions. The improvement of agricultural development level is often accompanied by the expansion of production scale or structural adjustment, which may produce more methane emissions in this process. On the contrary, the improvement of agricultural development level drives the improvement of infrastructure construction and agricultural production capacity and promotes the decline of methane emissions. In this paper, (agr) is used to represent the per capita gross output value of agriculture, forestry, animal husbandry and fishery. Residents' living conditions. There is an "inverted U-shaped" relationship between residents' living conditions and agricultural methane emissions. With the improvement of residents' living conditions, a greener and better ecological environment is needed. Therefore, there is a negative relationship between them, which is represented by (inc), the per capita income of rural residents. Industrial structure. Methane emissions come from planting and breeding. Therefore, the internal structure of the agricultural industry will affect methane emissions, which is represented by (ind), the proportion of planting industry in the total output value of agriculture, forestry, animal husbandry, and fishery. Urbanization level. Urbanization will occupy agricultural land, thus reducing agricultural methane emissions, which is measured by the proportion of urban population in total population and represented by (ubr). Industrialization level. With the development of industrialization, the emergence of green technology will reduce methane emissions in the process of agricultural production, which is represented by (iind), the proportion of non-agricultural total output value. Environmental regulations. Environmental regulations will have an impact on agricultural production and will encourage farmers engaged in agricultural production, especially aquaculture production, to adopt green technology and improve new technologies in the process of agricultural production. In this paper, (enr) is used to characterize the proportion of pollution control projects completed in each region in the year in the regional GDP. Agricultural machinery input. The input of agricultural machinery is mainly used to reduce methane emissions in agricultural planting. The input of great horsepower machinery can effectively improve planting efficiency, which is represented in this paper by (am), the total power of agricultural machinery. Agricultural science and technology R&D investment. Investment in agricultural science and technology R&D is the key to promoting agricultural technology progress. In this paper, (ati) is used to represent the proportion of government investment and financial support to agriculture in the total expenditure of provincial researchers.

3.2.2. Analysis of Influencing Factors of Spatial Differences of Agricultural Methane

On the basis of the geographical detector, this paper takes the three grain production areas as the main hierarchical method, uses the natural breakpoint method to treat variables by stratification, and divides the impact factors into four categories(Table 4). It can be seen that 9 factors play a driving role in the spatial differentiation of methane emissions, and that the influencing factors are divided into key and secondary influence factors according to the size of the driving factor explanatory power (q-statistical value). Among them, spatial factors (spa), industrial structure status (ind), industrialization level (iind), and agricultural science and technology R&D investment (ati) are the key influencing factors; agricultural development level (agr), residents' living conditions (inc), urbanization level (ubr), environmental regulation (enr), and agricultural machinery input (am) are the secondary influencing factors. The explanatory power of spatial factors reaches 0.438, which is the most key influencing factor, indicating that the production differences in different grain production zones played a decisive role in agricultural methane emissions. The explanatory power of industrial structure, industrialization level, and agricultural R&D investment are 0.122, 0.241, and 0.206, respectively. These variables represent the driving effect of agricultural transformation and upgrading on agricultural methane reduction from multiple aspects. In addition, the explanatory power of residents' living conditions, urbanization level, environmental regulation, and agricultural machinery input are all above 0.03, which cannot be ignored even though it is not the most critical factor influencing regional differences during the investigation period.

**Table 4.** Detection results of influencing factors of agricultural methane emissions in China.


3.2.3. Influencing Factors of Agricultural Methane Emission in Different Production Areas

Agricultural science and technology R&D investment has an impact on the spatial distribution of the three regions. It provides important technical support for methane emission reduction, and continuously strengthening the government's investment in agricultural science and technology is an important driver to achieve green development. The key influencing factors of main producing area, main sales area, and balance area are different(Table 5). Specifically, for the main producing area, the industrial structure is the most critical factor, and its explanatory power can reach 0.549. Main producing areas undertake the due obligations to ensure food security. Under the safety bottom line, adjusting the proportion between agriculture, forestry, animal husbandry, and fishery, especially the reasonable industrial structure adjustment process under the guidance of nutrition, will play an important role in methane emission reduction. The explanatory power of industrialization level and agricultural machinery input is also above 0.08, which is the main influencing factor of the spatial distribution of agricultural methane emission in main producing areas.

**Table 5.** Factor detection results of methane emission in different production areas.


For the main sales area, agricultural development level is the most important factor, and its explanatory power reaches 0.292. The level of agricultural development in main sales area will affect the income of residents and their demand for a better life. The better the agricultural development level, the stronger the demand for green and lowcarbon products, which will reduce the methane emissions. The explanatory power of residents' living standards and agricultural mechanization input is also higher than 0.1, indicating that these factors are also important factors influencing the spatial distribution of agricultural methane.

For the balance area, the industrial structure is the main influencing factor with the explanatory power of 0.479, and the industrialization level is also the main influencing factor with the explanatory power of 0.545. The evolution of industrial structure is that the transformation from planting to the integration of primary, secondary, and tertiary industries, which will effectively reduce agricultural methane emissions. The acceleration of industrialization will shift the rural labor force and release a greater demand to replace chemical fertilizers with organic fertilizer. Agricultural development level and investment in agricultural science and technology R&D are also key factors affecting the spatial distribution of balance area.

#### **4. Discussion**

Agriculture is both a source of greenhouse gas emissions and one of the most sensitive sectors to the impact of temperature change. In recent years, agricultural methane emission reduction has achieved some success, but the impact of greenhouse gas emissions caused by agricultural activities on climate change cannot be ignored. The year 2019 was the fifth warmest year since 1951, the national average temperature was 0.79 ◦C higher than usual, the average number of high temperature days was 11.8 days, with 4.1 days more than usual, compared with the second warmest year. Methane is the second largest greenhouse gas, and the recyclability of methane shows that controlling agricultural methane emission is an effective way to achieve the dual carbon goal. From the mechanism of agricultural methane emission, the formation ways of agricultural methane include paddy field emissions from the planting industry, livestock and poultry manure emission from the breeding industry, and livestock and poultry intestinal fermentation emission. Among them, the emission of livestock and poultry intestinal fermentation is the main emission source, accounting for 51.55–60.71%. Therefore, it is important to focus on controlling the emission of livestock and poultry intestinal fermentation to control agricultural methane emission, adjust the livestock breeding structure according to the actual situation of each region, and promote standardized breeding is an effective means. Paddy field emissions from planting industry is the second largest emission source after livestock and poultry intestinal fermentation, accounting for 26.45–35.54%. As the largest paddy producing country, China should choose appropriate paddy varieties and cultivation methods, promote large-scale planting, improve the utilization capacity of straw returning to the field, and other measures to reduce the emission coefficient of paddies. Compared with the previous literature, this paper adds the analysis of grain functional areas, uses nuclear density to investigate the differences and changes, utilizes geographic detectors to analyze the influencing factors, and then puts forward targeted emission reduction measures and countermeasures.

On the basis of the above analysis, the countermeasures to reduce agriculture methane emissions are provided as follows.

Firstly, it is important to carry out spatial cooperation and explore the spatiotemporal absorptive mechanism of agricultural methane emission reduction. The administrative boundary should be broken, the spatial association network should be constructed, and the free flow of elements (manpower, capital, etc.) in space should be realized. In addition, the temporal and spatial benefit protection and incentive policy support system needs to be constructed. Initially, it is necessary to clarify the spatial constraint characteristics of agricultural development, especially the differences between natural resources endowments. From the national level, legal, economic, and social tools should be adopted to

provide institutional guarantee for the coordination mechanism of cross-regional green development of agriculture. Then, a national cross-regional cooperation platform should be established. Consistent with the development of the national carbon market, the platform of methane emission reduction also needs to be built. In terms of agricultural methane emission reduction, attention should be paid to reducing the differences between regions, such as setting up a regional agricultural methane emission reduction fund, to form a platform for cooperation and exchange between regions.

Secondly, it is crucial to strengthen the research and development of green and lowcarbon agricultural technologies. The investment in agricultural science and technology has a significant impact on the spatial distribution of agricultural methane in the three grain functional areas; therefore, increasing investment in agricultural science and technology research and development is the most effective measure to form China's unique development pattern of high-value agriculture and low-carbon agriculture. In the input, the role of the promising government and the effective market shall be highlighted, social capital shall be guided into basic research and development, and an integrated guidance mechanism of breeding, reproduction, and promotion shall be established. In the output, the enthusiasm of various agricultural business entities needs to be stimulated, and a new pattern of green development should be formed by combining key technologies such as big data and smart agriculture. In addition, methane can be further utilized by recycling. The liquid storage process can be reduced for manure with high emission potential, and methane can be recovered through anaerobic fermentation to reduce greenhouse gas emissions. Methane can be recycled through the construction of biogas projects, and emissions can be reduced by changing the wet manure into dry manure and changing the way of manure is stored by covering.

Thirdly, the mechanism of reducing agricultural methane emission by time, land, and classification should be formed. For main producing area and balance area, emphasis should be placed on the optimization of industrial structure, especially on the premise of ensuring food security to promote the green transformation of industrial structure. The adjustment of industrial structure should be carried out from the aspects of spatial distribution, product structure, and industrial chain upgrading. The spatial layout should mainly highlight local characteristics, optimize the spatial layout, take into account variety improvement and quality improvement in product structure, improve regional quality agricultural products, solve the problem of "high quality but not good price", upgrade the industrial chain to achieve vertical extension, and promote advanced planting and breeding, advanced planting and fishing, ecological circular agriculture development, and realize the integrated development of primary, secondary, and tertiary industries. For the main grain sales area, the leading function of agricultural green transformation should be brought into play, and residents' awareness of green agricultural products should be initiated to drive the improvement and promotion of the entire agricultural industry chain.

This paper selects the emission coefficient published by IPCC for calculation, which is a widely used research method; however, due to the different actual conditions in various regions of the country, there is a problem of low accuracy in the face of the change of emission coefficient. At the same time, this paper lacks specific analysis at the provincial level, which is the direction of improvement in the future.

#### **5. Conclusions**

On the basis of the construction of the agricultural methane measurement method, the total amount and structure of the agricultural methane emissions in China is calculated and the spatiotemporal differentiation features of the agricultural methane emission through the grain production areas are analyzed. Then, through the interpretation of the kernel density function curve, distribution location, sample state, and kurtosis extension—the key factors affecting the spatial differentiation of agricultural methane emission—are analyzed by using geographic detectors. The main conclusions are as follows.

Firstly, agricultural methane emissions showed a trend of increasing first and then decreasing from 2000 to 2019; agricultural methane emissions showed a downward trend in 2006 and have been stable since then. During the 13th Five-Year Plan, the greening level was further improved, and methane emissions declined steadily. Agricultural methane emission reduction is closely related to the green development policy.

Secondly, the characteristics of agricultural methane emissions in the three major grain areas are basically the same as those in the whole country, but the changes in the main producing area, main sales area, and balance area have their own characteristics. From 2000 to 2019, the main reason for the reduction of agricultural methane emissions in the main producing area was the decline of intestinal fermentation emissions of livestock and poultry, which exceeded the increase of emissions from paddy fields. The paddy field, livestock and poultry intestinal fermentation, and manure management emissions in the main sales area and balance area fell, which generally reflects that the overall effect of green development in China is gradually emerging. At the provincial level, emphasis should be placed on reducing agricultural methane emissions in major producing provinces such as Hunan, Sichuan, Inner Mongolia, Heilongjiang, and Jiangxi.

Thirdly, the spatial patterns of agricultural methane emissions in China will remain unchanged in the short term, and it is difficult to achieve convergence.

Fourthly, the production heterogeneity of grain production areas has a decisive impact on the difference in agricultural methane emissions. Industrial structure, the level of industrialization, the investment in agricultural science, and technology R&D are the key influencing factors. Agricultural development level, residents' living conditions, urbanization level, environmental regulation, and agricultural machinery input are secondary influencing factors. In particular, the investment in agricultural science and technology R&D is influential in the agricultural methane emission pattern in the three main grain areas. For the main grain sales area, the agricultural development level is the core driving factor for the main sales area, and the industrialization level and industrial structure are the main driving factors for the balance area.

**Author Contributions:** Conceptualization, J.H.; data curation, P.L.; formal analysis, P.L.; investigation, J.H.; methodology, G.W.; writing—original draft, G.W.; writing—review & editing, F.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. 72221002) and the National Natural Science Foundation of China [72003111].

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** Data available from the authors upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

