Diversifying Crop Rotation Increases Food Production, Reduces Net Greenhouse Gas Emissions and Improves Soil Health
Diversifying Crop Rotation Increases Food Production, Reduces Net Greenhouse Gas Emissions and Improves Soil Health
Diversifying Crop Rotation Increases Food Production, Reduces Net Greenhouse Gas Emissions and Improves Soil Health
1038/s41467-023-44464-9
Global food production faces challenges in balancing the need for increased
yields with environmental sustainability. This study presents a six-year field
experiment in the North China Plain, demonstrating the benefits of diversify-
ing traditional cereal monoculture (wheat–maize) with cash crops (sweet
potato) and legumes (peanut and soybean). The diversified rotations increase
equivalent yield by up to 38%, reduce N2O emissions by 39%, and improve the
system’s greenhouse gas balance by 88%. Furthermore, including legumes in
crop rotations stimulates soil microbial activities, increases soil organic carbon
stocks by 8%, and enhances soil health (indexed with the selected soil phy-
siochemical and biological properties) by 45%. The large-scale adoption of
diversified cropping systems in the North China Plain could increase cereal
production by 32% when wheat–maize follows alternative crops in rotation
and farmer income by 20% while benefiting the environment. This study
provides an example of sustainable food production practices, emphasizing
the significance of crop diversification for long-term agricultural resilience and
soil health.
Producing more nutritious food to alleviate world hunger1 while safe- systems have also emitted great amounts of greenhouse gases (GHG)
guarding the environment2–4 is a significant challenge for humanity. and caused environmental degradation4–6. For instance, China’s crop
The challenge is much pronounced in highly-populated countries and production has increased by 74% in the past 30 years (1986–2016), but
regions where agricultural resources are limited. Conventional inten- this came with a >300% increase in synthetic fertilizer use7. In 2019,
sified food production systems’ reliance on major inputs of synthetic China’s GHG emissions associated with food production reached
agrochemicals, such as fertilizers and pesticides, have boosted food 2.4 Gt CO2-equivalent (CO2-eq)8, with about 50% emitted during the
production significantly since the ‘Green Revolution’, but these crop production stage8,9. The loss of soil fertility which may go along
1
State Key Laboratory of Efficient Utilization of Agricultural Water Resources, Beijing 100083, China. 2College of Water Resources & Civil Engineering, China
Agricultural University, Beijing 100083, China. 3School of Tropical Agriculture and Forestry, Hainan University, Haikou 570228, China. 4College of Life and
Environmental Sciences, Wenzhou University, Wenzhou, Zhejiang 325035, China. 5The μBC-Soil Group, Tallus Heights, Kelowna, BC, Canada. 6State Key
Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China. 7Key Laboratory of Agri-
cultural Water Resources, Centre for Agricultural Resources Research, Institute of Genetics and Developmental Biology, Chinese Academy of Sciences,
Shijiazhuang 050022, China. 8Department of Biological and Environmental Engineering, Riley-Robb Hall, Cornell University, Ithaca, NY 14853, USA. 9The UWA
Institute of Agriculture, The University of Western Australia, Perth, WA 6001, Australia. 10Land-CRAFT, Department of Agroecology, Aarhus University,
Aarhus, Denmark. 11Institute of Meteorology and Climate Research, Atmospheric Environmental Research (IMK-IFU), Karlsruhe Institute of Technology (KIT),
Garmisch Partenkirchen, Germany. e-mail: [email protected]; [email protected]; [email protected]; [email protected]
with the intensification of crop production further complicates food available phosphorus (AP), microbial biomass carbon (MBC), microbial
production10–12 and exposes it to climate risks13,14 and environmental biomass nitrogen (MBN), and microbial community composition and
health concerns15,16. diversity. These measurements enabled us to test the hypotheses that
Innovative concepts like integrated farming systems with diver- (a) cash-crop diversified system increases farmers’ net income without
sified crop rotations17,18 have emerged to address these food produc- jeopardizing crop yields, (b) legume diversified systems reduce field
tion and environmental sustainability challenges. These systems offer scale GHG emissions, and (c) integrating diversified rotations increases
a range of food crops to meet consumer demand for plant-based, food production, reduces GHG emissions, and benefits soil health
healthy food19—an increasing dietary trend in high and upper-middle- (Fig. 1). Our results demonstrate that (a) instructive findings from
income countries20—while providing other agricultural products such newly designed, tested, and validated diversified systems could guide
as animal feed, industrial fiber, or multi-purpose biofuels21,22. More- the North China Plain in establishing a more sustainable system to
over, integrated food production systems help increase farmers’ maintain or increase grain and protein production with reducing the
income and deliver socio-economic benefits20,23. However, little is damage to the environment and soil ecosystems, and (b) the results
known about how cash-crop and legume-diversified cropping systems from such a representatively intensive food producing region may
can achieve the triple goals—increasing food yields, reducing envir- provide a guide for the countries/regions with similar agricultural
onmental footprint, and benefiting soil health. environments to follow on an expanded scale.
In this context, we conducted a 6-year (2016–2022) field study in
the North China Plain—the food basket of China—one of the most Results
intensively cultivated regions in the world, where crop production is Diversified crop rotations increase ecosystem productivity
dominated by simple winter wheat (Triticum aestivum L.) and summer The food crops evaluated in the study are from different genera
maize (Zea mays L.) double cropping (wheat–maize, WM; two crops 1 families; thus, the yield of each crop was converted to the ‘equivalent’
year) which occupies 70% of the area’s arable land, delivering about product yield to wheat for a compatible comparison (detailed in
23% of China’s total cereal food24. The specific objective of the case Methods). The sweet potato (Ipomoea batatas L.) →winter
study was to asses comprehensively several diversified cropping sys- wheat–summer maize rotation (SpWM) increased the annual equiva-
tems in terms of food production, GHG balance, soil health benefits, lent yield by 38% compared to the conventional winter wheat–summer
and farmers’ income. We tested various diversified cropping systems, maize double-cropping which yielded 13,185 kg ha–1 annually (Fig. 2a).
incorporating cash-crops, legumes, other cereals, or forages into the Rotations diversified with sweet potato (SpWM), peanut (Arachis
conventional wheat–maize system. Major performance metrics inves- hypogea L.) (PWM), or soybean (Glycine max L.) (SWM) significantly
tigated were: (i) plant biomass, grain and protein yields, and farmers’ increased the annual economic benefit (net income) compared to the
net incomes; (ii) soil GHG fluxes, soil carbon (C) sequestration, and net winter wheat–summer maize control (Fig. 2b), with sweet potato→-
GHG emissions; and (iii) soil health parameters including soil pH, bulk winter wheat–summer maize rotation increasing 60%, and pea-
density, soil water content, total nitrogen (TN), dissolved organic nut→winter wheat–summer maize rotation and soybean→winter
carbon (DOC), soil nitrate-N (NO3–-N), ammonium-N (NH4+-N), wheat–summer maize rotation increasing by 13–22% (P < 0.05). Protein
Respiration
idue
idue
idue
Topsoil
Carbon
Fig. 1 | Schematic illustration of system integration from issues to outcomes. In legume crops can maintain crop yields, increase farmers’ income, and reduce GHG
the North China Plain—the case study area, traditional cereal monoculture (such as emissions due to the biological N2 fixation by legumes partly substituting for syn-
wheat–maize double-cropping, i.e., two cereal crops per year) requires inputs of thetic N inputs. Legume-included rotations can also enhance soil health by stimu-
synthetic agrichemicals and irrigation in food production, causing large green- lating soil microbial activities, increasing carbon sequestration, and enhancing
house gas (GHG) emissions; in contrast, rotation systems diversified with cash and nutrient cycles.
a 30 b 60000 c 3500
yr 1)
yr 1) 25 50000 3000
yr 1)
a
1
a a b bc c bc b
bc b
1
b
Equivalent yield (t ha
20 40000 c 2500
15 30000 2000 c
c
10 20000 1500
5 10000 1000
0 0 500
WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM
d 30
e f 3500
Winter wheat Winter wheat Winter wheat
yr 1)
50000
Summer maize Summer maize Summer maize
3000
yr 1)
1
yr 1)
25
1
a a a 2500 ab
20 ab ab b
Grain yield (t ha
b
b b 30000 b b 2000
15
1500
20000
10
1000
10000
5 500
0 0 0
WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM
Fig. 2 | Agroecosystem productivity. The entire rotation system productivity: lowercase letters denote significant differences between the rotations at P < 0.05.
a equivalent yield to wheat, b economic benefit – net income, and c protein yield. The exact P values: P < 0.001 in a–c, P = 0.001 in d, P = 0.032 in e and P = 0.004 in f.
The productivity of WM preceded by alternative crops in rotation: d annual grain In b and e, $1 US = 6.95 Chinese Yuan (as of May, 2023). For a–c, n = 9; d–f, n = 18.
yield, e economic benefit, and f protein yield. In a–c, the black bar within violin Treatment abbreviations: WM winter wheat–summer maize (control case), SpWM
plots shows the 25 and 75 percentiles, the whiskers beyond the bar represent 95 and sweet potato→winter wheat–summer maize rotation, PWM peanut→winter
5 percentiles, and the inner white dot indicates the median. The violin shaped area wheat–summer maize rotation, SWM soybean→winter wheat–summer maize
reveals the data distribution. The dashed lines in d–f are the baseline values of the rotation, SmWM spring maize→winter wheat–summer maize rotation, RSWM
WM rotation, and the error bars are the standard deviation of annual value summed ryegrass–sorghum→winter wheat–summer maize rotation. Source data are pro-
over winter wheat and summer maize. In a–f, one-way ANOVA with two-sided and vided as a Source Data file.
post-hoc test was conducted to determine significant differences. Different
yield varied among the rotation systems, with ryegrass (Lolium multi- control which was at 8.9 ± 1.0 kg N ha–1 (Table S1, Fig. 3a), decreasing
florum L.)–sorghum (Sorghum bicolor L.) rotation (RSWM) and soy- by 30%, 42%, and 49% in the peanut→wheat–maize rotation,
bean→winter wheat–summer maize rotation having the highest soybean→wheat–maize rotation, and sweet potato→wheat–maize
protein yields at 2274 and 1883 kg ha–1 yr–1, respectively, being 8–31% rotation, respectively. These rotations with peanut, soybean, and
higher than winter wheat–summer maize rotation (Fig. 2c). The sweet potato received a reduced dosage of N fertilizer as compared to
ryegrass-sweet sorghum→winter wheat–summer maize rotation, with the wheat–maize control (Table S1). The soils of all crop rotations,
four crops in 2 years, produced the most aboveground biomass for including wheat–maize, acted as net sinks for atmospheric CH4.
livestock, while the soybean→winter wheat–summer maize rotation However, the diversified crop rotations increased the sink strength by
produced the most protein for human nutrition. In contrast, an 33–76% above wheat–maize (Fig. 3b). The wheat–maize control had
entirely cereal-based system (i.e., spring maize→ winter the highest annual global warming potential (GWP) for soil N2O and
wheat–summer maize rotation, SmWM) decreased equivalent yield by CH4 fluxes (3764 kg CO2-eq ha–1 yr–1, Fig. 3c). The cereal-including
16% and protein yield by 23% compared to the winter wheat–summer rotations (ryegrass-sweet sorghum→wheat–maize rotation and spring
maize control. maize→wheat–maize rotation) had 22% and 19%, respectively, lower
The preceding crops in rotation significantly affected grain yield GWPs than wheat–maize rotation, while the rotations including
(Fig. 2d), net income (Fig. 2e), and protein yield (Fig. 2f) of the suc- legumes (peanut, soybean) or sweet potato had 32%, 43%, and 51%,
ceeding winter wheat–summer maize crops as the rotation continued. respectively, lower GWPs than wheat–maize rotation. Besides field-
The diversified crop rotations with sweet potato, peanut, and soybean scale direct GHG emissions, indirect GHG emissions associated with
had a positive carryover effect on the productivity of the succeeding the use of agrochemicals and irrigation were highest for wheat–maize
winter wheat and summer maize. On average, winter wheat–summer rotation (8802 kg CO2-eq ha–1 yr–1) and 34–41% lower for rotations with
maize preceded by non-cereal crops in rotation showed increased sweet potato, peanut, soybean and spring maize (Fig. 3d).
grain yield, economic benefit, and protein yield by 26–32%, 39–46%, The soil carbon stocks (SOC) in the 0–90 cm layer significantly
and 25–29%, respectively, as compared to the winter wheat–summer increased from 2016 to 2022 for all treatments, but the magnitude of
maize cereal monoculture. the change differed among rotation systems. Six years after the
initiation of the experiment (in 2022), the peanut→wheat–maize
Diversified crop rotations reduce net GHG emissions rotation had the highest soil C sequestration (2.03 t ha–1 yr–1),
The amounts of supplementary irrigation and fertilizers applied to followed by soybean→wheat–maize (1.91 t ha–1 yr–1) and sweet
each rotation varied among them, as did the N2O and CH4 emissions potato→wheat–maize (1.44 t ha–1 yr–1), while ryegrass-sweet
measured weekly from October 2016 to October 2022 (Supplementary sorghum→wheat–maize, spring maize→wheat–maize rotation, and
Figs. S1–2). The diversified crop rotations significantly reduced wheat–maize rotations had significantly lower C sequestration rates
annual cumulative N2O emissions compared to the wheat–maize (0.21–0.69 t C ha–1 yr–1) (Fig. 3e). Soil C sequestration offset total GHG
a 12 b 0 c 7000
yr 1)
a
Annual mean cumulative N2O
yr 1)
1
1
emission (kg N2O-N ha
bc
1
8
(kg CO2-eq ha
cd b b
b 4000
d bc b
6 c bc
bc cd
-2 3000
d
4
2000
2 -3 1000
WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM
d N P2O5 K2 O
e 100
0-30cm
f 20
70-90cm
a b
yr 1)
(t CO2-eq ha 1 yr 1)
10 cd cd
d
GHG emissions
e c
7500
1
60
(kg CO2-eq ha
5
d e e
5000 0
40 f
d e
-5
c
2500 20 Direct GHGs a b
-10 Indirect GHGs
Soil carbon sequestration
Net GHG emissions
0 0 2016 2022 2016 2022 2016 2022 2016 2022 2016 2022 2016 2022
-15
WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM WM SpWM PWM SWM SmWM RSWM
Fig. 3 | Direct and indirect greenhouse gas (GHG) emissions in the diversified between the rotations at P < 0.05. The exact P values: P < 0.001 in a–c and f. In e,
crop rotations. a Annual mean N2O emissions; b CH4 emissions; c GWPN2O+CH4; ** denotes significant changes in soil carbon stock in the 0–90 cm soil layer
d indirect GHG emissions; e changes in soil carbon stock in the 0–90 cm depth; between 2016 and 2022 at P < 0.01 (P = 0.005 for SpWM, P = 0.002 for PWM and
f net GHG emissions. In a–c, diamond box plots show triangles spanning 25 and 75 SWM). Treatment abbreviations: WM winter wheat–summer maize (control case),
percentiles meeting at the median. Small squares indicate means. Whisker lines SpWM sweet potato→winter wheat–summer maize rotation, PWM peanut→winter
extend to the minimum and maximum. In e and f, error bars in each column wheat–summer maize rotation, SWM soybean→winter wheat–summer maize
represent the standard deviation of whole-profile sums of three replicates. In rotation, SmWM spring maize→winter wheat–summer maize rotation, RSWM
a–c, e, and f, one-way ANOVA with two-sided and post-hoc test was conducted for ryegrass–sorghum→winter wheat–summer maize rotation. Source data are pro-
statistical significance; different lowercase letters denote significant differences vided as a Source Data file.
emissions by 75–89% in the rotations with sweet potato, peanut, and sorghum→wheat–maize rotation had the lowest SOC and DOC due to
soybean and 7–21% in the spring maize→wheat–maize, ryegrass-sweet the removal of forage for silage during the growing season, which
sorghum→wheat–maize and wheat–maize rotations (Fig. 3f). Conse- stimulated soil microbial activity, leading to the highest MBC (Fig. S3).
quently, wheat–maize rotation had the largest net GHG emissions Soil microbial community composition and diversity—an indi-
totaling 10,025 kg CO2-eq ha–1 yr–1, and the diversified rotations with cator of the biological mechanism of soil health—were determined in
sweet potato, soybean and peanut lowered net GHG emissions by 83%, 2016 and 2022 to evaluate the cumulative changes after 6 years’
90%, and 92%, respectively, as compared to wheat–maize (P < 0.05). rotational experimentation. It showed that the diversification of crop
rotations had positive effects (Fig. 5). For example, the Shannon index
Diversified crop rotations enhance soil health and microbial —a popular index for quantifying microorganism diversity taking into
diversity account the number of species and relative abundance of each species
Soil health is often disregarded in the assessment of food crop pro- in a sample—significantly increased by 7–10% in the rotations
duction systems. In the present study, we used the Cornell Soil Health with sweet potato, peanut, and soybean rotations, but did not
Assessment (CSHA) scoring method and applied Principal Component change in wheat–maize or the graminaceous crop-based spring
Analysis (PCA) to evaluate how crop diversification influenced soil maize→wheat–maize and ryegrass-sweet sorghum→wheat–maize sys-
health. A soil health score was calculated for each rotation based on tems from 2016 to 2022 (Fig. 5a). Inclusion of the legumes peanut or
key indicators (detailed in Methods). The PCA revealed that the first soybean in rotations increased bacterial community operational
two principal components accounted for 70.5% of the cumulative taxonomic units (OTU) richness (Fig. S4a–f) and fungal community
percent variability (PC1 = 51.1% and PC2 = 19.4%) in soil health (Fig. 4a), Shannon index, Chao 1, and OTU richness (P < 0.05) compared to
with the measured soil health indicators clustered in distinct groups. wheat–maize (Fig. S4g–l). Many of the soil health indicators are
At the end of the 6-year study, the peanut→wheat–maize rotation had dynamic and may change weekly, seasonally, or annually. However,
the highest soil health score (59.5), followed by soybean→wheat–maize after 6 years of ‘rotation nourishment,’ soils in the rotations with sweet
rotation (56.8) and sweet potato→wheat–maize rotation (52.9) potato, peanut, and soybean had significantly higher values for Chao 1,
(Fig. 4b), being 41–59% higher than wheat–maize (33.2). abundance-based coverage estimator (ACE), and community richness,
A closer examination of the soil indicators that contributed to soil compared to the cereal monoculture spring maize→wheat–maize sys-
health scores revealed that the soil in the peanut→wheat–maize rota- tem (Fig. 5b–d).
tion had a 6.5% higher SOC concentration and 29.7% higher dissolved Redundancy analyses, depicting the association of soil physio-
organic carbon (DOC) but 5.7% lower total N and 15.4% lower available chemical properties and crop rotations in the bacterial (Fig. 5e) and
P (AP) concentration than wheat–maize (Fig. S3). These changes were fungal (Fig. 5f) communities, found that the first two sets of compo-
due to reduced N fertilizer (−41%) and associated higher nutrient use nents explained 54.6% (bacterial) and 43.5% (fungal) of the total var-
efficiencies by crops in the diversified rotations. The ryegrass-sweet iation in the two communities, respectively. The relationship was
complex, but it was evident that bacterial and fungal community Discussion
composition changed from 2016 to 2022. The rotations with sweet The empirical evidence from the 6-year field experiment demon-
potato, peanut, and soybean had distinct bacterial community com- strated significant positive impacts of diversified crop rotations on
positions from the other three systems (Fig. 5e), whereas the five stu- various agroecosystem functions and services. Based on the study’s
died diversified rotations had fungal community compositions findings and double-cropping planted area25, we estimated that
distinctly different from the wheat–maize control (Fig. 5f). The 6 years diversified cropping systems could reduce net CO2-eq emissions by
of rotational nourishment affected soil physiochemical properties 106.8 ± 31.7 million tonnes annually, offsetting 5.6% of the annual GHG
interactively with soil microbial communities. The rotations with sweet emissions associated with China’s food system (1.9 billion tonnes in
potato, peanut, and soybean bacterial communities were closely 2020 reported by FAO26). The primary reason behind this reduction is
associated with SOC and AP, whereas those of wheat–maize, spring the decrease in synthetic N fertilizer use by 3.6 million tonnes. Farmers
maize→wheat–maize, and ryegrass-sweet sorghum→wheat–maize in the North China Plain region could benefit from 20% increases in
rotations were closely associated with TN and microbial biomass car- annual net income, equivalent to 84 billion Yuan ($11.6 US billion) in
bon (MBC). In contrast, the rotations with sweet potato, peanut, and total. Moreover, these diversified rotations will increase Cornell Soil
soybean fungal communities were associated with SOC, pH, and DOC. Health Assessment scores by 45%, due to increased C sequestration
rates compared to the winter wheat–summer maize rotation. Fur-
Multiple functions assessments of diversified crop rotations thermore, wheat and maize yields may increase by 32%, equivalent to
We used the comprehensive evaluation index concept (CEI) to assess 73.5 million tonnes per year if planted following alternative crops such
the synergies and trade-offs of the different crop rotations related to as sweet potato, soybean or peanuts; this would make about 36.1
yield, nutritional value, soil-related indicators (health, C sequestration, million tonnes of additional straw biomass available annually for
microbial biodiversity), net GHG emissions, and economic benefit alternative uses, such as feed, bioenergy, or enhancing soil carbon
(Eqs. 11–20 in Methods) (Fig. 6a, d). The results showed that the stocks.
wheat–maize rotation had the lowest CEI value, averaging 0.19 The three-dimensional benefits (grain and protein yields, GHG
(Fig. 6b), while the sweet potato→wheat–maize rotation had the emissions, and soil health) of crop diversification in the winter
highest CEI value of 0.81, followed by peanut→wheat–maize rotation of wheat–summer maize rotation system on the North China Plain will
0.75, and soybean→wheat–maize rotation of 0.69. The diversified add significant ‘social novelty’ to the challenge of bringing essential
rotations with sweet potato, peanut, and soybean had CEI values more food nutrients to dining tables without adversely affecting soil health.
than triple relative to wheat–maize and were significantly higher than The social novelty is far beyond the agricultural benefits of crop
the cereal-based spring maize→wheat–maize and ryegrass-sweet diversification that are broadly recognized worldwide. For example,
sorghum→wheat–maize rotations (0.25–0.29). Significant relation- long-term studies in North America show that rotational diversification
ships occurred among the key indicators across the rotation systems increases crop yields, reduces the impact of adverse weather on eco-
(Fig. 6c). Soil health was positively corelated with crop yield (r = 0.79), system productivity27, and enhances system robustness28. Studies in
economic benefits (r = 0.85), and nutrition score (r = 0.90), but was Europe and Africa demonstrate that field-scale diversity has a sub-
negatively correlated with net GHG emissions (r = –0.83). Soil health stitutive interaction with fertilization leading to increased yields at low
was also positively correlated with soil carbon sequestration (r = 0.79) N fertilizer doses29–32. Several Chinese studies implementing crop
and soil biodiversity (r = 0.67), contributing substantially to the trade- diversification in various forms or scales show that farmers can adopt
off against net GHG emissions (r = –0.83). multiple strategies to increase farming resilience and economic
100
a b
1.0
WM SWM
SpWM SmWM
PWM RSWM NO3-
80 +59%
0.5
+52%
Soil health score
60 +54%
MBC c +43%
BD c
TN d
0.0
DOC +13%
40 f e
MBN
-0.5
20
AP
-1.0
0
-1.0 -0.5 0.0 0.5 1.0 2016 WM SpWM PWM SWM SmWM RSWM
PC1 (51.09%)
Fig. 4 | Soil health assessment for different rotation systems. a Principal com- Different lowercase letters denote significantly different groups of crop rotations at
ponent analysis (PCA) for eigenvalues of the ten key soil properties, including soil P < 0.05. The red dashed line with red percentages is the baseline of the WM
organic carbon (SOC), dissolved organic carbon (DOC), total nitrogen (TN), NO3–-N, rotation. The black dashed line with black percentages is the baseline of 2016.
available P (AP), pH, microbial biomass N (MBN), microbial biomass carbon (MBC), Treatment abbreviations: WM winter wheat–summer maize (control case), SpWM
bulk density (BD), and soil water content (SWC). Each soil property was normalized sweet potato→winter wheat–summer maize rotation, PWM peanut→winter
as an individual CSHA (Cornell Soil Health Assessment) contributor to the overall wheat–summer maize rotation, SWM soybean→winter wheat–summer maize
soil health score. b CSHA score for each rotation. In b, the whiskers in each column rotation, SmWM spring maize→winter wheat–summer maize rotation, RSWM
represent the standard deviation of three replications; circles are the corre- ryegrass–sorghum→winter wheat–summer maize rotation. Source data are pro-
sponding data points for the bar chart, n = 3. One-way ANOVA with two-sided and vided as a Source Data file.
post-hoc test was conducted for the significance test among different treatments.
a Shannon b Chao1
7.0 8000
+22%
6.5 +10%
+18% +20%
7000
+8% +18%
+7% +12% +14%
d p = 0.002 p = 0.001
c
5.0 5000
2016 WM SpWM PWM SWM SmWM RSWM 2016 WM SpWM PWM SWM SmWM RSWM
2022 2022
c ACE d Richness
8000 7000
+22%
+18% +20%
+24%
+18% +22%
7000 +13% +14% 6000 +20%
+20% +16% +16%
ab ab a
ab ab ab
a
ab ab
b ab
6000 5000 b
c p = 0.001 p = 0.0003
c
5000 4000
2016 WM SpWM PWM SWM SmWM RSWM 2016 WM SpWM PWM SWM SmWM RSWM
2022 2022
e f
Bacteria Fungi
0.8
0.8
DOC**
pH**
DOC NO3--N**
AP
0.4
0.4
pH NH4+-N
SOC*
RDA2 (17.70%)
RDA2 (25.15%)
TN NO3--N
SOC*
0.0
0.0
WM
SpWM TN* WM
SpWM
PWM MBC/MBN* MBC* PWM
SWM SWM
SmWM SmWM
-0.8
-0.8
RSWM RSWM
-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0
RDA1 (29.40%) RDA1 (25.81%)
Fig. 5 | Soil microbe community diversity and soil properties in response to indicating soil properties including soil organic carbon (SOC), dissolved organic
rotation. a Shannon Index; b Chao 1; c ACE index; d OTU richness measured in the carbon (DOC), total nitrogen (TN), NO3–-N, NH4+-N, available P (AP), pH, microbial
six crop rotations from 2016 to 2022; e Redundancy analysis (RDA) of bacterial OTU biomass carbon (MBC), microbial biomass N (MBN), and the MBC to MBN ratio.
and soil properties; f RDA of fungal OTU and soil properties. The percentage values Treatment abbreviations: WM winter wheat–summer maize (control case), SpWM
in a–d represent the increment from 2016 to 2022; dot is the mean value, whisper sweet potato→winter wheat–summer maize rotation, PWM peanut→winter
line bars on each dot are the standard deviations for three replications. In a–d, one- wheat–summer maize rotation, SWM soybean→winter wheat–summer maize
way ANOVA with two-sided and post-hoc test was conducted to determine sig- rotation, SmWM spring maize→winter wheat–summer maize rotation, RSWM
nificant differences among treatments. Different lowercase letters below pillars ryegrass–sorghum→winter wheat–summer maize rotation. Source data are pro-
indicate significant differences between rotations at P < 0.05. Charts e and f show vided as a Source Data file.
redundancy analysis with colored dots and triangles indicating rotation, vectors
a b
1.0
Equivalent yield (+)
0.8
Economic
benifit (+) Soil health (+)
CEI
0.6
0.2
0.0
WM SpWM PWM SWM SmWM RSWM
Nutrition Soil c
yield (+) biodiversity 1
Equivalent yield Equivalent
(+) yied *** * * ** 0.8
emissions (-) sequestration (+) Soil biodiversity 0.40 0.52 0.65 -0.95 0.93 Soil
biodiversity
-0.6
-0.8
Soil
Soil health 0.79 0.85 0.90 -0.83 0.79 0.67
health
WM SpWM PWM SWM SmWM RSWM -1
Equivalent Economic Nutrition Net GHG Soil Soil Soil
yield benefit yield emissions carbon biodiversity health
sequestration
d
30 Equivalent Equivalent Equivalent
yield yield yield
20
Economic Soil health Economic Soil health Economic Soil health
10 benefit benefit benefit
Net GHG Soil carbon Net GHG Soil carbon Net GHG Soil carbon
SWM emssions sequestration SmWM emssions sequestration RSWM emssions sequestration
Fig. 6 | Comprehensive assessments of the different rotation systems. a radar soil carbon stock in the 0–90 cm layer from October 2016 to October 2022; Net
map of multiple objective analysis to assess the various functions of the rotations; GHG emissions denotes the net GHG budgets including direct and indirect GHG
b comprehensive evaluation index (CEI, Eqs. 15–24 in Method, supplementary emissions minus soil carbon sequestration; Soil biodiversity is represented by the
material); c synergy and trade-off relationship among the key variables across the Shannon Index of bacterial and fungal communities; Soil health score is calculated
different crop rotations; d, Nightingale Rose Charts of multiple objectives analysis from 10 indicators using the overall Cornell Soil Health Assessment (CSHA) matrix
to assess the detailed functions 1for each crop rotation. In a, values are normalized. combined with PCA81–83. Treatment abbreviations: WM winter wheat–summer
(+), positive variables, (–), negative variables. Equivalent yield re-expresses crop maize (control case), SpWM sweet potato→winter wheat–summer maize rotation,
yield as equivalent to the grain yield of wheat; Economic benefit denotes the crop PWM peanut→winter wheat–summer maize rotation, SWM soybean→winter
yield multiplied by its market price minus consumable cost; Nutrition yield includes wheat–summer maize rotation, SmWM spring maize→winter wheat–summer maize
17 nutritional components of each crop such as crude protein, fat content, macro- rotation, RSWM ryegrass–sorghum→winter wheat–summer maize rotation. Source
and micro-elements, and vitamins68. Soil carbon sequestration denotes changes in data are provided as a Source Data file.
benefits28,33. Overall, cropping diversification offers a comprehensive potato, peanut and soybean significantly increased soil C after 6 years
systems’ approach to enhance agroecosystem productivity and of rotation nourishment compared to the traditional winter
adaptability to changing climates worldwide. wheat–summer maize double-cropping system. Introducing sweet
A major factor in designing diversified cropping systems is to potatoes and legumes into these rotations has been shown to stimu-
ensure an increase in soil C sequestration—as this plays a crucial role late microbial growth, increased C use efficiency36, and promoted the
for soi health, in climate change mitigation and achieving C formation of highly stable mineral-associated soil C complexes37.
neutrality34,35. We found that the diversified rotations with sweet Labile litter compounds from these rotational crops are the dominant
source of microbial products that affect stable SOC with a high capa- increased fungal community diversity compared to winter
city for C stabilization38. Furthermore, our findings highlight a close wheat–summer maize, shaping soil nutrient profiles with different
correlation between soil health and nutrition yield, underscoring the carbon substrates that impacted microbial diversity. Moreover,
interconnectedness of both dimensions and their reliance on agri- diversified crop rotations potentially reduce soil agrochemical con-
cultural practices39. tamination due to decreased application of synthetic N fertilizer and
A noteworthy aspect of the study is that we measured SOC con- herbicides, favoring soil microenvironments54,55. Other studies have
centrations to 90 cm depth. Researchers often focus on topsoil layers demonstrated the positive influence of diversified rotations on soil
(0–20 cm) when studying SOC status, potentially masking some health. Maize yields and SOC content in diversified rotations were
treatment effects, especially in crops with different root masses at more resilient than those found in monoculture56, potentially
various soil depths. We found that SOC in the 0–20 cm soil layer increasing adaptability to climate change57. Diversified agroecosys-
accounted for 51% of the total SOC, 0–30 cm layer accounted for 67%, tems are considered an economically viable alternative to business-as-
and 0–50 cm layer accounted for 80%. Our results highlight the usual maize–soybean rotations in the Midwest US58, with cover crops,
importance of subsoil C stocks when assessing soil C storage for a perennials, and small grain cereals enhancing soil health by 32–49%
more precise estimation of the contributions of soil C sequestration in and crop productivity by 16–29%59.
mitigating climate change. The novelty of this study is its highly diversified crop rotations,
Lowering GHG emissions without reducing crop yields is a major with cash crops and legumes replacing cereal monocultures, together
challenge for global agriculture. We found that the legume-based with rotating shallow-rooted crops with deep-rooted crops. It posi-
(peanut→wheat–maize, soybean→wheat–maize) rotations significantly tively affects ecosystem services and functions, partly compensating
reduced N2O emissions by 30–42% compared to the conventional for the lower services in less diversified cereal-based rotations. Crop
cereal-based wheat–maize rotation while significantly increasing sys- diversification offers significant benefits in reducing GHG emissions
tem productivity (grain and protein yields). Legume crops improve soil and has a synergistic effect on plant biomass and protein production,
N cycling40, but may also increase N2O emissions, specifically following soil health, and microbial community biodiversity. Diversification
residue incorporation41,42. However, legume cultivation requires sig- reduced crop inputs, particularly N fertilizer and irrigation amounts,
nificantly less fertilizer than cereal crops such as maize or wheat. minimizing environmental costs while increasing farmers’ incomes.
Therefore, the legume-based rotations with peanut and soybean Thus, diversified crop rotations can serve as an effective strategy for
required 37% less fertilizer than conventional wheat–maize, explaining sustainable food production in areas worldwide with environments
the overall reduction in N2O emissions for the legume-based rotations similar to the North China Plain.
as N fertilizers were the main source of N2O emissions in our study40. The scientific evidence presented in this study holds significant
The manufacture of inorganic N fertilizer remained the most sig- potential for informing and reinforcing current policies and incentive
nificant contributor to total GHG emissions (about 36–39%). Similarly, schemes to improve the green transition of farming systems in China.
we found that the diversified rotations with sweet potato, soybean, and Specifically, China set objectives to promote soybean production
peanut significantly reduced net GHG emissions by 75–92% compared (initiated in 2019), implement an action plan for zero increase in fer-
to the wheat–maize, spring maize→wheat–maize, and ryegrass-sweet tilizer and pesticide use (launched in 2015), and develop strategies to
sorghum→wheat–maize systems. In a study in Western Canada, cereals reduce agricultural GHG emissions with the ultimate goal of trans-
rotated with legumes decreased net GHG emissions by 17–35% com- forming agricultural systems into zero ‘emitters’ by 2060. When
pared to cereal monocultures43 due to improved N uptake associated designing effective diversified cropping systems, it is important to
with organic N mineralization released from legume residues44. In consider local environmental conditions to balance land use for food
Europe, increasing legume production by introducing legumes into production with the other functions of an ecosystem. Optimizing crop
cereal rotations effectively improved plant-based protein while redu- configuration is essential—through social, economic, and environ-
cing environmental impacts45,46. In the US corn belt, diversified maize mental dimensional planning within a complex framework of various
and soybean systems maintained or increased maize yields while indicators to ensure long-term sustainability. We recommend that
reducing N losses to the environment47. developing and adopting diversified cropping systems should be a key
Introducing shallow-rooted crops like sweet potato and annual consideration in agricultural policy setting and a top priority for on-
legumes in rotation with deep-rooted crops such as winter wheat farm decision-making, as they are critical for achieving long-term food
enables the use of soil nutrients across the entire rooting zone, production resilience, maintaining soil health, and ensuring environ-
improving nutrient use efficiency and reducing soil nutrient losses mental sustainability. Given that climate variability will increase, multi-
through leaching. Alternating deep- and shallow-rooted crops also year assessments of the sustainability of agricultural production sys-
enhances root distribution, increases soil porosity and permeability, tems remain tentative, and funding should be sought for additional
reduces soil bulk density48, and improves soil aggregate stability49 to observation sites for long-term monitoring.
anchor SOC50. Furthermore, diversified crop mixtures alter the quality
and quantity of organic matter entering the soil as plant litter, Methods
enhancing the stabilization of microbial-derived compounds that The field experiment was established in October 2016 at the Luan-
provide positive feedback for plant growth. The distinct microbial cheng Agro-Ecosystem Station of the Chinese Academy of Sciences
communities associated with diverse crop mixtures can promote soil C (37°50′ N, 114°40′ E; altitude, 50.1 m) in Luancheng County, Hebei
cycling51. In our study, diversified legume–cereal and sweet Province (Fig. S5a), which represents the agricultural production and
potato–cereal rotations significantly increased the α-diversity of soil climate conditions of the North China Plain. The experimental area
bacterial and fungal communities after 6 years of rotation. The experiences a warm, temperate zone, semi-humid, monsoon climate
increased fungal diversity resulted from a higher number of species, with a frost-free period of 200 days. The annual mean air temperature
while the increased Shannon’s diversity index in the bacterial com- was 14.7 ± 1.0 °C over the last 20 years, and the annual average pre-
munity resulted from changes in the relevant abundance of certain cipitation was 472 ± 161 mm over the last 62 years (Fig. S5b). Fig. S5c
species. Diverse crops stimulate bacterial community activities in the presents the average monthly precipitation and air temperatures from
rhizosphere52, as the residue root exudates and plant litter from rota- 2016 to 2020. The experimental site had loam soil with sandy loam in
tional crops provide a greater diversity of residue carbon substrates, the surface layers, light/medium loam at 40–80 cm depth, and light
supporting the growth of diverse soil microorganisms53. We found that clay below 80 cm. The 0–10 cm soil layer before the experiment start
the rotations with sweet potato, peanut, and soybean significantly has pH 7.6 ± 0.2, 240 ± 42 µS cm–1 EC, 1.49 ± 0.07 g cm–3 bulk density,
0.42 ± 0.03 cm3 cm–3 field capacity, 0.10 ± 0.02 cm3 cm–3 wilting point,
CPI j
1.06 ± 0.14 g kg–1 total N, 11.5 ± 0.4 g kg–1 SOC, 9.3 ± 3.1 mg kg–1 available EBij = Y ield ij × Priceij × Cos t ij ð2Þ
phosphorous, and 109.6 ± 24.1 mg kg–1 available potassium. Prior to CPI 2008
the experiment, the core food crops on the same fields were wheat and
maize grown in a 12-month double-cropping system that is popular PC i = Y ield i × βi ð3Þ
regionally. From the start of the experiment, the sites were managed
with optimized irrigation and seeding techniques, with all residues where EYij is equivalent yield of crop i at year j (kg ha–1), wheat pricej is
returned to the soil during the winter wheat–summer maize year. price of winter wheat in year j adjusted for China National Consumer
The experiment comprised eight crops: cereals winter wheat (cv. Price Index (CPI, relative Yuan), economic benefit (EB) is within-year
Kenong 2009), summer maize (cv. Jundan 20), and spring maize (cv. crop yield multiplied by its corresponding price in Chinese Yuan minus
Jundan 20); legumes soybean (cv. Shidou 12) and peanut (cv. Jihua 4); crop production costs (Yuan ha–1), j is a given year, PCi is total protein
cash/forage: sweet potato (cv. Shangshu 19), ryegrass (cv. Dongmu 70) content of crop i (kg ha–1), and βi is protein concentration per kg crop
and sorghum (cv. Jintianza 3). There were five newly designed diversified product (%) referred from previous studies67,68. (For comparison, the
crop rotations: [sweet potato→winter wheat–summer maize rotation following exchange rate was used: $1 US = 6.95 Chinese Yuan and 1
(SpWM, 2-year cycle); peanut→winter wheat–summer maize rotation € = 7.60 Yuan in May 2023.) Each indicator for one crop rotation was
(PWM, 2-year cycle); soybean→winter wheat–summer maize rotation summed across the crops within one crop rotation cycle. The CPI was
(SWM, 2-year cycle); ryegrass-sweet sorghum→winter wheat–summer used to adjust crop prices across years to eliminate inflation or
maize rotation (RSWM, 2-year cycle); spring maize→winter deflation as a source of variability between cases in different years.
wheat–summer maize rotation (SpWM, 2-year cycle)]. The regionally Table S2 lists the crop prices and CPI values for individual years and
dominant winter wheat–summer maize double-cropping system (WM, crops, and Table S3 lists each crop’s protein content.
1-year cycle) was the control. ‘→’ denotes crop rotation across years and
‘–’ denotes crop rotation within a year. Unlike pure winter Gas flux and net GHG emissions
wheat–summer maize, the alternative rotations included fallow seasons Soil GHG (N2O and CH4) fluxes were measured using the static cham-
for parts of their non-winter wheat–summer maize years (Fig. S6). ber method69. Plastic frames with water-sealing grooves were inserted
Each rotation was replicated three times, resulting in 18 plots 5 cm into the soil between crop rows, and chambers (20 × 20 × 30 cm)
established using a randomized complete block design (Fig. S6). All were placed on top of the frames during measurements. For flux
cycles were repeated on their respective plots from October 2016 to measurements, gas samples (>30 mL per sample) were taken between
October 2022. The winter wheat–summer maize double-cropping 9:00 and 10:00 am at specific time intervals (0, 12, 24, and 36 min after
control rotation completed six cycles, while the other five rotations chamber closure) using 100 mL plastic syringes attached to a three-
completed three rotation cycles. Supplementary Table S1 lists the total way stopcock. Gas samples were collected about once every 7 days
amounts of inputs, including N, P, K, irrigation, diesel for tillage, from sowing to harvest, with more frequent measurements (every
electricity for irrigation, pesticide, seeds, and labor for each crop in 3 days) after N fertilizer application and precipitation events. The
each rotation. Supplementary Data 1 details fertilizer and irrigation headspace air temperature in the chamber was recorded with a ther-
amounts and timings for each crop in each rotation in each year. Basal mometer. Gas samples were analyzed using a gas chromatograph
N applications (as urea) were broadcast and plowed (20 cm depth) into (Agilent 7890 A, Agilent Inc., USA) equipped with a flame-ionization
the soil before seeding. Topdressing of N was performed at the key detector for measuring CH4 at 200 °C and an electron capture detector
growth stages with specific amounts summarized in Supplemen- for measuring N2O at 330 °C. The column temperature was set to 55 °C.
tary Data 1. CH4 was separated with a Porapak Q column and detected by the
flame-ionization detector. N2O was detected by a micro electron cap-
Equivalent yield and economic benefit ture detector (μECD). The carrier gas was high-purity N2 (99.999%)
Crops were harvested at maturity with a conventional combine har- with a flow velocity of 20 mL min–1. A purge gas (5% CO2 in N2) was used
vester with an automated weighing system, which measured the total to avoid interference with CO2. The gas flux was calculated using the
crop yield of each plot. After air-drying to a constant moisture content linear approach. Soil N2O and CH4 fluxes for each rotation across the 6
(13%), wheat and maize grain yields were recorded60,61. Sorghum bio- years are illustrated in Figs S1 and S2.
mass yield (fresh weight) per plot was weighed to represent silage For the net GHG emissions (kg CO2-eq ha–1) estimation, we con-
moisture of 30% (70% dry matter)61. Ryegrass dry matter was deter- sidered indirect GHG emissions of crop inputs from a life cycle
mined in 1 × 1 m sub-sample after drying at 85 °C until constant assessment, direct N2O and CH4 emissions from the soil, and changes
weight62. Sweet potato yield (fresh weight) was obtained by weighing in soil carbon stock from 0–90 cm during the 6-year experiment. The
the fresh root tubers harvested from each plot using a single row system boundary in the life cycle assessment was set from manu-
harvester63. All mature pods in each plot were collected to determine facture, storage, and delivery of external inputs such as fertilizer to
peanut yield accurately. Ten uniform peanut plants were used for crop harvest, including soil carbon storage change (Fig. S7). Potential
investigating the pod numbers per plant. All harvested peanut pods emissions related to the logistics of transporting, exporting, or mar-
were air-dried to a water content of approximately 14% and weighed to keting grain products beyond the farm gate were excluded as they
calculate the pod yield64. After harvesting each plot, soybean seeds were considered outside the system boundary.
were cleaned and weighed, with the grain yields adjusted to 13.0% Indirect GHG emissions (CE) result from manufacturing, storing,
moisture content65. transporting, and delivering agricultural inputs, such as fertilizer (N, P,
For yield comparisons among different crop products, the K), pesticides, diesel, and electricity to the farm gate, as follows:
equivalent yields of the six crop rotations were calculated by multi-
plying the crop yield by the market price relative to that of winter X
n
CE = ðDk × C k Þ ð4Þ
wheat in the same year66,67. Each cropping system’s equivalent yield,
k =1
economic benefit, and protein yield were determined as follows:
where CE (kg CO2-eq ha−1 yr–1) is the indirect GHG emissions for one
crop within a given rotation, Dk is the amount of the kth input (Table S1),
Priceij
EY ij = Y ield ij × ð1Þ and Ck is the corresponding CO2-eq emission coefficients (see Table S4
W heat pricej for the manufacture, storage, and delivery of agricultural inputs).
Direct N2O and CH4 fluxes from soil were calculated from October samples were determined using a Multi 3100 N/C TOC analyzer
2016 to October 2022, as follows: (Analytik Jena, Germany). A modified Bremner’s standard protocol
was used to extract soil NO3–-N and NH4+-N using 2 M KCl75, mea-
M T 0 P dc sured with a UV-1800 spectrophotometer (Mapada Instruments,
f= H × 60 ð5Þ
V 0 T P0 dt Shanghai, China). Soil available P concentrations were determined
using 0.5 mol L–1 NaHCO3 (pH 8.5) and the molybdenum-blue col-
where f (ug m–2 h–1) is GHG emission fluxes of N2O and CH4 emissions, orimetric method76.
M is the molecular mass of the specified gas (44 g mol–1 for N2O and Soil organic carbon (SOC) concentrations in the 0–90 cm soil
16 g mol–1 for CH4), V0 (m3mol−1) is the molar volume of an ideal gas, T0 layer (0–10, 10–20, 20–30, 30–50, 50–70, and 70–90 cm increments)
(K) is the ideal gas temperature, T (K) and P (kPa) are the air tem- were measured after each crop harvest from October 2016 to October
perature and pressure in the chamber at the sampling time, P0 (kPa) is 2022, in samples collected by auger at five random locations per plot.
the standard air pressure, H (m) is the height of chamber, dc/dt (ppm Soil samples from each sampling depth within a plot were combined,
min–1) is the slope of the linear regression curve for the change of dried at 105 °C for 24 h, sieved (2 mm), pretreated with 0.5 M HCl to
headspace gas concentrations during times of chamber closure, and remove carbonates77, and ball-milled. SOC concentrations were
60 is min h−1. Total GHG emissions during each crop growth and fallow determined using a vario Macro CNS Analyzer (Elementar, Germany).
period were calculated using linear interpolation70. Soil bulk density in all soil layers down to 0.9 m in all plots was mea-
Cumulative N2O and CH4 emissions during each crop’s growing sured using the gravimetric method78.
season or fallow period were calculated as: Soil samples were also collected before the start of the
experiment (October 2016) and at the final summer maize harvest
n
X
ðf i+1+ f iÞ (October 2022) for DNA sequencing and microbial community
F= × ðt i + 1 t i Þ × 24 ð6Þ
2 composition and diversity measurements. Subsamples for DNA
i=1
sequencing were immediately placed on dry ice and stored at –80 °C
where F is the cumulative emission of one of the two gases from one until DNA extraction, while the other subsamples were taken to the
experimental plot (kg ha–1) over a specified period, f is the gas emission laboratory for chemical analysis. Soil DNA was extracted from 0.5 g
flux (kg ha–1 h–1), t is the sampling time in days, i is the ith sampling from freeze-dried soil using the Powersoil DNA Isolation Kit (Mo Bio
the plot, and 24 is the conversion coefficient between hours and days Laboratories, Carlsbad, CA, USA) according to the manufacturer’s
(h day–1). protocol. DNA quality was assessed according to the 260/280 nm
The global warming potential (GWP) of seasonal soil N2O and CH4 and 260/230 nm absorbance ratios using a NanoDrop ND2000
emissions, combining Eq. (7) results and re-expressed in CO2-eq71, was spectrophotometer (NanoDrop, ND2000, Thermo Scientific, 111
calculated as: Wilmington, DE) and then stored at –80 °C for further molecular
biology analysis.
GW P N 2 O + CH 4 = N 2 O × 273 + CH 4 × 27 ð7Þ Primer pairs—barcode-515F/806 R (5′-GTGCCAGCMGCCGCGGTAA-
3′/5′-GCACTACHVGGGTWTCTAAT-3′)—were used to amplify the V3 + V4
The annual change in SOC stock (ΔC; kg CO2-eq ha–1 yr–1) from region of the bacterial 16 S rRNA gene79, yielding accurate taxonomic
October 2016 to October 2022 was estimated for each crop rotation, as information with few biases among various bacterial taxa. Primers—ITS1
follows: (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′-TGCGTTCTTCATC-
GATGC-3′)—were used to amplify the ITS1 region of the fungal rRNA
ðT 1 T 2 Þ × 44 gene. The PCR was carried out in a mixture (final volume, 50 µL) com-
ΔC = 12 ð8Þ
n prising 2 µL DNA template (1–10 ng), 5 µL of each 2 μM primer, 25 μL
Premix Ex Taq (Takara Biotechnology), and 13 μL sterilized water.
where T1 and T2 (t ha–1) are the total soil carbon storage values in the Thermal-cycling conditions were as follows: initial denaturation of 3 min
0–90 cm soil layer in October 2016 and October 2022, respectively, n is at 94 °C, six touchdown cycles of 45 s at 94 °C, 60 s from 65 to 58 °C,
the number of years of the study period, and 44/12 is the conversion 70 s at 72 °C, followed by 22 cycles of 45 s at 94 °C, 60 s at 58 °C, 60 s at
coefficient from C into CO2. 72 °C, with a final elongation of 72 °C for 10 min. Tag-encoded high-
Net GHG emissions were equal to the sum of indirect (CE) and throughput sequencing of the 16 S and ITS genes from purified and
direct (GWP of N2O and CH4) GHG emissions minus soil carbon quantified PCR products was performed by the Magigene Company
sequestration (ΔC) when soil carbon was the sink, calculated as follows: (Beijing, China) using the Illumina MiSeq platform.
Raw sequences were quality processed using the Quantitative
Net GHGs = CE + GW P N 2 O + CH 4 ΔC ð9Þ Insights into Microbial Ecology (QIIME) pipeline (version 2.1)79.
Resampling operational taxonomic units (OTUs) was based on the
minimum sequence numbers to correct the sampling effort before
Soil physiochemical properties, microbial diversity and soil further analysis. The OTUs were clustered with a 97% similarity cutoff
health scoring using the UPARSE pipeline (version 7.1)80. Taxonomic assignment was
Soil samples (0–20 cm depth) were collected at the beginning of the carried out using OTUs with SILVA (16 S) and Unite (ITS).
experiment (October 2016) and at the end of the summer maize har- Alpha diversity was characterized by six indices (Shannon-Wea-
vest in October 2022 to measure soil physicochemical properties: ver, Chao1, ACE, Richness, Simpson, and Pielou), calculated using the
including soil pH, bulk density, soil water content, total nitrogen (TN), OTU numbers for bacteria and fungi. (see Supplementary Information
dissolved organic carbon (DOC), soil nitrate-N (NO3–-N), ammonium-N for detailed calculations).
(NH4+-N), available phosphorus (AP), and microbial biomass carbon The Cornell Soil Health Assessment (CSHA)81–83 combined with
and nitrogen (MBC and MBN, respectively). principal component analysis (PCA)84,85 was used to determine health
Soil pH was measured in a 1:5 (w/v) soil: water suspension as scores for soils sampled in 2016 and 2022, which allowed comparing
described by McLean (1982)72. TN was determined using the mod- each plot’s changes in soil health after 6 years of ‘rotation nourish-
ified Kjeldahl method73. MBC and MBN were determined using ment’. The ten soil indicators (physical attributes: bulk density, soil
the chloroform fumigation and extraction method74. DOC con- water content; chemical attributes: soil pH, TN, SOC, DOC, NO3–-N, AP;
centrations extracted from fumigated and unfumigated soil biological attributes: MBC and MBN) were normalized as individual
CSHA scores81. The weights of individual CSHA scores were based on (3) Calculate the entropy weight (Wj) of the evaluation index:
PCA of all soil indicators (Table S5), representing the sum of the
eigenvectors derived from the first three principal components, which 1 Hj
Wj = Pn ð16Þ
were selected based on the inflection point from a Scree plot and n Hi
i=1
Kaiser’s cut-off (eigenvalues > 1) (Fig. S8). These first three principal
components accounted for 51.09%, 19.44%, and 15.04% of the data’s (4) Construct a weighted normalized decision matrix (Z) using the
total variation, with a cumulative variance of 85.57%, capturing most of normalized matrix fij and the weights of each index Wj:
the variation among the soil indicators. The soil health overall score (%) 2 3
was computed as a weighted average of all individual CSHA scores, z 11 z 12 z 1m
6z z 22 z 2m 7
calculated as follows: 6 21 7
Z = ðW j × f ij Þn × m = 6 7 ð17Þ
4 5
ðA1 × w1 Þ + ðA2 × w2 Þ + + ðAn × wn Þ z n1 z n2 z nm
Soil health score = ð10Þ
w1 + w2 + + wn
(5) Determine the positive ideal solution vector Z+ and the negative
where A is the CSHA score for each individual soil indicator, and w is ideal solution vector Z− using the positive ideal solution (Zj+) and
the weighting factor for each soil indicator (Table S5). the negative ideal solution (Zj−):
most positive and furthest from the negative ideal solution86–88. The j=1
uncertainty associated with extreme weather events and climate 20. Lucas, E., Guo, M. & Guillén-Gosálbez, G. Low-carbon diets can
change. reduce global ecological and health costs. Nat. Food 4,
394–406 (2023).
Reporting summary 21. Ray, D. K., Ramankutty, N., Mueller, N. D., West, P. C. & Foley, J. A.
Further information on research design is available in the Nature Recent patterns of crop yield growth and stagnation. Nat. Commun.
Portfolio Reporting Summary linked to this article. 3, 1293 (2012).
22. Ranganathan, J., Daniels, R. J. R., Chandran, M. D. S., Ehrlich, P. R. &
Data availability Daily, G. C. Sustaining biodiversity in ancient tropical countryside.
The raw reads from Illumina sequencing described in this study, are Proc. Natl Acad. Sci. USA 105, 17852–17854 (2008).
available at NCBI under the accession no. SRR26083954-26083974 23. Crippa, M. et al. Food systems are responsible for a third of global
(16 S rRNA) and SRR26083997-26084017 (ITS). The dataset generated anthropogenic GHG emissions. Nat. Food 2, 198–209 (2021).
in this study has been deposited in the database under accession code 24. Yang, X. L. et al. Diversified crop rotations enhance groundwater
on Figshare at https://doi.org/10.6084/m9.figshare.24793563. Source and economic sustainability of food production. Food Energy Secur.
data are provided with this paper. 10, e311 (2021).
25. Liu, Z. J. et al. Patterns and causes of winter wheat and summer
Code availability maize rotation area change over the North China Plain. Environ. Res.
Code scripts used in this study are available at: https://figshare.com/s/ Lett. 17, 044056 (2022).
5c07d2ac7c8dc4a5b678. 26. FAO. FAOSTAT. (2023). https://www.fao.org/faostat/en/#data/EM/
visualize, last access: 19 Jul 2023.
References 27. Bowles, T. M. et al. Long-term evidence shows that crop-rotation
1. Webb, P. et al. The urgency of food system transformation is now diversification increases agricultural resilience to adverse growing
irrefutable. Nat. Food 1, 584–585 (2020). conditions in North America. One Earth 2, 284–293 (2020).
2. Gerten, D. et al. Feeding ten billion people is possible within four 28. Li, J. X. et al. Diversifying crop rotation improves system robustness.
terrestrial planetary boundaries. Nat. Sustain. 3, 200–208 (2020). Agron. Sustain. Dev. 39, 38 (2019).
3. Beerling, D. J. et al. Potential for large-scale CO2 removal via 29. Giller, K. E. The food security conundrum of sub-Saharan Africa.
enhanced rock weathering with croplands. Nature 583, Glob. Food Secur. 26, 100431 (2020).
242–248 (2020). 30. Hijbeek, R. et al. Do organic inputs matter – a meta-analysis of
4. Cui, Z. L. et al. Pursuing sustainable productivity with millions of additional yield effects for arable crops in Europe. Plant Soil 411,
smallholder farmers. Nature 555, 363–366 (2018). 293–303 (2017).
5. Garnett, T. et al. Sustainable intensification in agriculture: Premises 31. Mhlanga, B., Ercoli, L., Pellegrino, E., Onofri, A. & Thierfelder, C. The
and policies. Science 341, 33–34 (2013). crucial role of mulch to enhance the stability and resilience of
6. Gu, B. J. et al. Cost-effective mitigation of nitrogen pollution from cropping systems in southern Africa. Agron. Sustain. Dev. 41,
global croplands. Nature 613, 77–84 (2023). 29 (2021).
7. National Bureau of Statistics of China. Agriculture. China Statistical 32. MacLaren, C. et al. Long-term evidence for ecological intensifica-
Yearbook − 2021. (China Statistics Press, Beijing, China, 2021). tion as a pathway to sustainable agriculture. Nat. Sustain. 5,
8. Liu, G., Zhang, F. & Deng, X. Z. Half of the greenhouse gas emissions 770–779 (2022).
from China’s food system occur during food production. Commun. 33. Lin, B. B. Resilience in agriculture through crop diversification:
Earth Environ. 4, 161 (2023). adaptive management for environmental change. BioScience 61,
9. Friedlingstein, P. et al. Global carbon budget 2021. Earth Syst. Sci. 183–193 (2011).
Data 14, 1917–2005 (2022). 34. Robertson, G. P., Paul, E. A. & Harwood, R. R. Greenhouse gases in
10. Cumming, G. S. et al. Implications of agricultural transitions and intensive agriculture: contributions of individual gases to the
urbanization for ecosystem services. Nature 515, 50–57 (2014). radiative forcing of the atmosphere. Science 289,
11. Marris, E. A world without soil: The past, present, and precarious 1922–1925 (2000).
future of the earth beneath our feet. Nature 601, 503–504 (2022). 35. Powlson, D. S. et al. Limited potential of no-till agriculture for cli-
12. Xu, H. G. et al. Ensuring effective implementation of the post-2020 mate change mitigation. Nat. Clim. Change 4, 678–683 (2014).
global biodiversity targets. Nat. Ecol. Evol. 5, 411–418 (2021). 36. Liu, K., Bandara, M., Hamel, C., Knight, J. D. & Gan, Y. T. Intensifying
13. Wang, B. et al. Sources of uncertainty for wheat yield projections crop rotations with pulse crops enhances system productivity and
under future climate are site-specific. Nat. Food 1, 720–728 (2020). soil organic carbon in semi-arid environments. Field Crops Res.
14. Shi, Y. et al. Building social resilience in North Korea can mitigate 248, 107657 (2020).
the impacts of climate change on food security. Nat. Food 3, 37. Chenu, C. et al. Increasing organic stocks in agricultural soils:
499–511 (2022). knowledge gaps and potential innovations. Soil Tillage Res. 188,
15. Suweis, S., Carr, J. A., Maritan, A., Rinaldo, A. & D’Odorico, P. Resi- 41–52 (2019).
lience and reactivity of global food security. Proc. Natl Acad. Sci. 38. Drinkwater, L. E., Wagoner, P. & Sarrantonio, M. Legume-based
USA 112, 6902–6907 (2015). cropping systems have reduced carbon and nitrogen losses. Nature
16. Tilman, D. et al. Future threats to biodiversity and pathways to their 396, 262–265 (1998).
prevention. Nature 546, 73–81 (2017). 39. Montgomery, D. R. & Biklé, A. Soil health and nutrient density:
17. Chai, Q. et al. Integrated farming with intercropping increases food beyond organic vs. conventional farming. Front. Sustain. Food Syst.
production while reducing environmental footprint. Proc. Natl 5, 417 (2021).
Acad. Sci. USA 118, e2106382118 (2021). 40. Guinet, M., Nicolardot, B. & Voisin, A.-S. Nitrogen benefits of ten
18. Xie, W. et al. Crop switching can enhance environmental sustain- legume pre-crops for wheat assessed by field measurements and
ability and farmer incomes in China. Nature 616, 300–305 (2023). modelling. Eur. J. Agron. 120, 126151 (2020).
19. Guo, Y. X. et al. Air quality, nitrogen use efficiency and food security 41. Abdalla, M. et al. A critical review of the impacts of cover crops on
in China are improved by cost-effective agricultural nitrogen man- nitrogen leaching, net greenhouse gas balance and crop pro-
agement. Nat. Food 1, 648–658 (2020). ductivity. Glob. Change Biol. 25, 2530–2543 (2019).
42. Olesen, J. E. et al. Challenges of accounting nitrous oxide emissions 62. Hu, Z. Y. et al. Italian ryegrass–rice rotation system for biomass
from agricultural crop residues. Glob. Change Biol. 00, 1–10 (2023). production and cadmium removal from contaminated paddy fields.
43. Gan, Y. T. et al. Improving farming practices reduces the carbon J. Soils Sediment. 20, 874–882 (2020).
footprint of spring wheat production. Nat. Commun. 5, 5012 (2014). 63. Song, D. G. et al. Influence of planting distance and density on the
44. Ma, B. L., Liang, B. C., Biswas, D. K., Morrison, M. J. & McLaughlin, N. yield and photosynthetic traits of sweet potato (Ipomoea balatas L.)
B. The carbon footprint of maize production as affected by nitrogen under an intercropping system with walnut (Juglans regia) saplings.
fertilizer and maize-legume rotations. Nutr. Cycl. Agroecosyst. 94, Soil Tillage Res. 196, 104484 (2020).
15–31 (2012). 64. Zhao, J. H. et al. The impacts of soil tillage combined with plastic
45. Karlsson, H., Ahlgren, S., Strid, I. & Hansson, P.-A. Faba beans for film management practices on soil quality, carbon footprint, and
biorefinery feedstock or feed? Greenhouse gas and energy bal- peanut yield. Eur. J. Agron. 148, 126881 (2023).
ances of different applications. Agric. Syst. 141, 138–148 (2015). 65. Liu, X. B. et al. Yield response of continuous soybean to one-season
46. Plaza-Bonilla, D., Nolot, J.-M., Raffaillac, D. & Justes, E. Innovative crop disturbance in a previous continuous soybean field in North-
cropping systems to reduce N inputs and maintain wheat yields by east China. Field Crops Res. 138, 52–56 (2012).
inserting grain legumes and cover crops in southwestern France. 66. Alam, M. K., Bell, R. W. & Biswas, W. K. Decreasing the carbon
Eur. J. Agron. 82, 331–341 (2017). footprint of an intensive rice-based cropping system using con-
47. Roesch-McNally, G. E., Arbuckle, J. G. & Tyndall, J. C. Barriers to servation agriculture on the Eastern Gangetic Plains. J. Clean. Prod.
implementing climate resilient agricultural strategies: The case of 218, 259–272 (2019).
crop diversification in the U.S. Corn Belt. Glob. Environ. Change 48, 67. Wang, S. Q. et al. Diversified crop rotations reduce groundwater use
206–215 (2018). and enhance system resilience. Agric. Water Manag. 276,
48. Shah, K. K., Modi, B., Pandey, H. P., Subedi, A. & Shrestha, J. 108067 (2023).
Diversified crop rotation: an approach for sustainable agriculture 68. Groot, J. C. J. & Yang, X. L. Trade-offs in the design of sustainable
production. Adv. Agric. 2021, 1–9 (2021). cropping systems at a regional level: a case study on the North
49. Liang, J. P. & Shi, W. J. Cotton/halophytes intercropping decreases China Plain. FASE 9, 14 (2022).
salt accumulation and improves soil physicochemical properties 69. Zheng, X. H. et al. Quantification of N2O fluxes from soil–plant
and crop productivity in saline-alkali soils under mulched drip irri- systems may be biased by the applied gas chromatograph meth-
gation: a three-year field experiment. Field Crops Res. 262, odology. Plant Soil 311, 211–234 (2008).
108027 (2021). 70. Huang, J. X., Chen, Y. Q., Sui, P. & Gao, W. S. Estimation of net
50. Hazra, K. K. et al. Diversification of maize-wheat cropping system greenhouse gas balance using crop- and soil-based approaches:
with legumes and integrated nutrient management increases soil Two case studies. Sci. Total Environ. 456–457, 299–306 (2013).
aggregation and carbon sequestration. Geoderma 353, 71. IPCC. Climate change 2021: The physical science basis. Contribu-
308–319 (2019). tion of working group i to the sixth assessment report of the inter-
51. Mooshammer, M. et al. Microbial feedbacks on soil organic matter governmental panel on climate change. (ed. Masson-Delmotte, V.
dynamics underlying the legacy effect of diversified cropping sys- P. Z. A.) (Cambridge University Press, Cambridge, United Kingdom
tems. Soil Biol. Biochem. 167, 108584 (2022). and New York, NY, USA, 2021). https://www.ipcc.ch/report/sixth-
52. Dias, T., Dukes, A. & Antunes, P. M. Accounting for soil biotic effects assessment-report-working-group-i/.
on soil health and crop productivity in the design of crop rotations. 72. Mclean, E. O. Soil pH and Lime Requirement. in Methods of Soil
J. Sci. Food Agric. 95, 447–454 (2015). Analysis: Part 2 Chemical and Microbiological Properties, 9.2.2,
53. Garbeva, P., Veen, J. A. V. & Elsas, J. D. V. Microbial diversity in soil: Second Edition (ed. Page, A. L.) 199-224 (American Society of Agr-
Selection of microbial populations by plant and soil type and onomy, Inc. and Soil Science Society of America, Inc., Madison,
implications for disease suppressiveness. Annu. Rev. Phytopathol. Wisconsin, USA, 1982).
42, 243–270 (2004). 73. Bremner, J. M. & Shaw, K. Denitrification in soil. I. Methods. Investig.
54. Selim, M. M. A review of advantages, disadvantages and challenges J. Agric. Sci. 51, 22–39 (1958).
of crop rotations. Egypt. J. Agron. 41, 1–10 (2019). 74. Ladd, J. N. & Amato, M. Relationship between microbial biomass
55. Garland, G. et al. Crop cover is more important than rotational carbon in soils and absorbance (260 nm) of extracts of fumigated
diversity for soil multifunctionality and cereal yields in European soils. Soil Biol. Biochem. 21, 457–459 (1989).
cropping systems. Nat. Food 2, 28–37 (2021). 75. Bremner, J. M. Inorganic forms of nitrogen. in Methods of Soil Ana-
56. Yang, Y. H. et al. Optimizing crop rotation increases soil carbon and lysis: Part 2 Chemical and Microbiological Properties (ed. Norman,
reduces GHG emissions without sacrificing yields. Agric. Ecosyst. A. G.) 1179–1237 (American Society of Agronomy Madison, 1965).
Environ. 342, 108220 (2023). 76. Olsen, S. R. & Sommers, L. E. Phosphorus. in Methods of Soil Ana-
57. Jarecki, M. et al. Long-term trends in corn yields and soil carbon lysis: Part 2 Chemical and Microbiological Properties, 9.2.2, 2nd
under diversified crop rotations. J. Environ. Qual. 47, Edition (ed. Page, A. L.) 403–430 (American Society of Agronomy,
635–643 (2018). Inc. and Soil Science Society of America, Inc., 1982).
58. Baldwin-Kordick, R. et al. Comprehensive impacts of diversified 77. Harris, D., Horwáth, W. R. & van Kessel, C. Acid fumigation of soils to
cropping on soil health and sustainability. Agroecol. Sustain. Food remove carbonates prior to total organic carbon or CARBON-13
Syst. 46, 331–363 (2022). isotopic analysis. Soil Sci. Soc. Am. J. 65, 1853–1856 (2001).
59. Chahal, I., Hooker, D. C., Deen, B., Janovicek, K. & Van Eerd, L. L. 78. Grossman, R. B. & Reinsch, T. G. 2.1 Bulk density and linear exten-
Long-term effects of crop rotation, tillage, and fertilizer nitrogen on sibility. in Methods of Soil Analysis: Part 4 Physical Methods (eds.
soil health indicators and crop productivity in a temperate climate. Dane, J. H. & Topp, G. C.) 201–228 (Soil Science Society of America,
Soil Tillage Res. 213, 105121 (2021). Inc., 2002).
60. Fang, Q., Zhang, X. Y., Chen, S. Y., Shao, L. & Sun, H. Selecting traits 79. Caporaso, J. G. et al. QIIME allows analysis of high-throughput
to increase winter wheat yield under climate change in the North community sequencing data. Nat. Methods 7, 335–336 (2010).
China Plain. Field Crops Res. 207, 30–41 (2017). 80. Edgar, R. C. UPARSE: Highly accurate OTU sequences from micro-
61. da Silva, P. C. G. et al. No-tillage and crop rotation increase crop bial amplicon reads. Nat. Methods 10, 996–998 (2013).
yields and nitrogen stocks in sandy soils under agroclimatic risk. 81. Gugino, B. K. et al. Cornell Soil Health Assessment Training Manual,
Field Crops Res. 258, 107947 (2020). Edition 1.2, (Cornell University Press, Cornell University, 2007).
82. Idowu, O. J. et al. Farmer-oriented assessment of soil quality using Author contributions
field, laboratory, and VNIR spectroscopy methods. Plant Soil 307, X.Y., S.K., T.D., and X.J. designed the study. X.Y. designed the field
243–253 (2008). experiment. X.Y., J.X., S.L., L.X. and Y.S. facilitated data collection. X.Y.,
83. Moebius-Clune, B. N. et al. Comprehensive Assessment of Soil J.X., X.J. and K.B.B. do the data analysis. X.Y., X.J., K.B.B., Y.G. and
Health – The Cornell Framework, Edition 3.2, (Cornell University K.H.M.S. wrote the manuscript. Y.G., K.H.M.S., S.P. and T.S.S. improved
Press, Cornell University, 2016). and revised the writings. All authors discussed the results and con-
84. Chang, T. T., Feng, G., Paul, V., Adeli, A. & Brooks, J. P. Chapter Three tributed to the final manuscript.
- Soil health assessment methods: Progress, applications and
comparison. in Adv. Agron., Vol. 172 (ed. Sparks, D. L.) 129–210 Competing interests
(Academic Press, 2022). The authors declare no competing interests.
85. Congreves, K. A., Hayes, A., Verhallen, E. A. & Van Eerd, L. L. Long-
term impact of tillage and crop rotation on soil health at four tem- Additional information
perate agroecosystems. Soil Tillage Res. 152, 17–28 (2015). Supplementary information The online version contains
86. Wang, H. D. et al. Multi-objective optimization of water and fertilizer supplementary material available at
management for potato production in sandy areas of northern https://doi.org/10.1038/s41467-023-44464-9.
China based on TOPSIS. Field Crops Res. 240, 55–68 (2019).
87. Chen, Y., Zhu, M. K., Lu, J. L., Zhou, Q. & Ma, W. B. Evaluation of Correspondence and requests for materials should be addressed to
ecological city and analysis of obstacle factors under the back- Xiaolin Yang, Taisheng Du, Xiaotang Ju or Yantai Gan.
ground of high-quality development: Taking cities in the Yellow
River Basin as examples. Ecol. Indic. 118, 106771 (2020). Peer review information Nature Communications thanks Baojing Gu,
88. Liu, X. G., Peng, Y. L., Yang, Q. L., Wang, X. K. & Cui, N. B. Deter- Robert Rees and Feng Zhou for their contribution to the peer review of
mining optimal deficit irrigation and fertilization to increase mango this work. A peer review file is available.
yield, quality, and WUE in a dry hot environment based on TOPSIS.
Agric. Water Manag. 245, 106650 (2021). Reprints and permissions information is available at
89. Delgado, A. & Romero, I. Environmental conflict analysis using an http://www.nature.com/reprints
integrated grey clustering and entropy-weight method: A case
study of a mining project in Peru. Environ. Modell. Softw. 77, Publisher’s note Springer Nature remains neutral with regard to jur-
108–121 (2016). isdictional claims in published maps and institutional affiliations.
90. Chen, C. W., Wang, J. H., Wang, J. C. & Shen, Z. H. Developing
indicators for sustainable campuses in Taiwan using fuzzy Delphi Open Access This article is licensed under a Creative Commons
method and analytic hierarchy process. J. Clean. Prod. 193, Attribution 4.0 International License, which permits use, sharing,
661–671 (2018). adaptation, distribution and reproduction in any medium or format, as
91. R Core Team. R: A language and environment for statistical com- long as you give appropriate credit to the original author(s) and the
puting. in Computing (R foundation for statistical computing, source, provide a link to the Creative Commons licence, and indicate if
Vienna, Austria, 2014). changes were made. The images or other third party material in this
article are included in the article’s Creative Commons licence, unless
Acknowledgements indicated otherwise in a credit line to the material. If material is not
We thank J.Z., N.Q., P.H., F.L., Y.C. and R.C. from China Agricultural included in the article’s Creative Commons licence and your intended
University for field work contributions, and Dr. Y.S. for assisting the use is not permitted by statutory regulation or exceeds the permitted
analysis of soil microbial diversity. This work was financed jointly by the use, you will need to obtain permission directly from the copyright
National Natural Science Foundation of China (Grant No.32071975 to holder. To view a copy of this licence, visit http://creativecommons.org/
X.Y., 52239002 to T.D. and 41830751 to X.J.), the Hebei Province Key licenses/by/4.0/.
Research and Development Program of China (20326411D-1 to X.Y.) and
Pioneer Center for Research in Sustainable Agricultural Futures (Land- © The Author(s) 2024
CRAFT) (DNRF grant number P2 to K.B.B.).