ARTICLE

Conservation agriculture increases the soil
resilience and cotton yield stability in climate
extremes of the southeast US

Amin Nouri

1,7

✉

, Daniel C. Yoder

1

, Mohammad Raji

2,8

, Safak Ceylan

1

, Sindhu Jagadamma

1

, Jaehoon Lee

1

,

Forbes R. Walker

1

, Xinhua Yin

3

, Judith Fitzpatrick

4

, Brady Trexler

4

, Prakash Arelli

5

& Arnold M. Saxton

6

Climate extremes pose a global threat to crop security. Conservation agriculture is expected

to offer substantial climate adaptation bene

fi

ts. However, synergistic effects of conservation

practices on yield during normal versus extreme climates and underlying regulatory

mechanisms remain elusive. Here, we analyze 29-years of climate data, cotton (Gossypium

hirsutum L.) yield, and soil data under 32 management practices in Tennessee, USA. We

fi

nd

that long-term no-tillage enhanced agroecosystem resilience and yield stability under climate

extremes and maximized yield under favorable climate. We demonstrate that no-tillage

bene

fi

ts are tied with enhanced soil structural stability and organic carbon. No-tillage

enhanced the effectiveness of legume cover crop in stabilizing cotton yield during relatively

dry or wet, and dry years, while nitrogen fertilizer rate and precipitation timing, controlled

yield stability in wetter years. Our

fi

ndings provide evidence-based insights into how man-

agement strategies can enhance agroecosystem resilience and production stability in climate

extremes.

https://doi.org/10.1038/s43247-021-00223-6

OPEN

1

Department of Biosystems Engineering and Soil Science, The University of Tennessee, Knoxville, TN, USA.

2

Department of Electrical Engineering and

Computer Science, The University of Tennessee, Knoxville, TN, USA.

3

Department of Plant Sciences, The University of Tennessee, Jackson, TN, USA.

4

Proli

fi

c Earth Sciences, Inc., Montgomery, NY, USA.

5

USDA-ARS, Crop Genetics Research Unit, Jackson, TN, USA.

6

Department of Animal Science, The

University of Tennessee, Knoxville, TN, USA.

7

Present address: Hermiston Agricultural Research and Extension Center, Oregon State University, Hermiston,

OR, USA.

8

Present address: Google, Sunnyvale, CA, USA.

✉

email:

amin.nouri.g@gmail.com

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

1

123456789

0():,;

C

limate extremes pose a global threat to agricultural produc-
tion and future food, feed, and

fi

ber security. The economic

cost of crop yield losses in

fl

icted by extreme climatic events

between 2003 and 2014 in developing countries around the world
was about $80 billion

1

. In the U.S., the number of extreme weather

incidents costing more than $1 billion has steadily increased in recent
decades

2

. Projections of climate change suggest that with rising

average global temperature and the impact of warming weather on
the water cycle, the frequency of dry spells will increase in subtropical
and lower mid-latitude terrestrial regions. While the average pre-
cipitation is expected to decline in most regions, light and moderate
rainfall will shift to less frequent higher-intensity events

3

. This is

expected to increase the precipitation contrast between dry and wet
regions and seasons, leading to further increases and decreases in
storm erosivity in wet and dry regions, respectively

4

,

5

. It is thus likely

that the adverse impact of large seasonal and annual weather varia-
bility on natural resources and

fi

eld crop production will be more

pronounced, even in the near term

6

.

The U.S. southeastern and lower southern plains

—

a region

known as the cotton belt

—

comprise almost 17% of global upland

cotton production, representing over $6 billion in total revenue
and approximately 38% of the world

’

s exports of raw cotton

fi

ber

7

. The cotton belt also outranks the rest of the country in

billion-dollar extreme weather disasters

2

. Rising summer tem-

peratures (especially during the night) are also evident in the vast
majority of the cotton belt

8

. Although cotton is considered a

drought-tolerant plant, warm summer nights have a documented
impact on yield losses by stimulating respiration more than
photosynthesis

9

. Altered rainfall distribution has also led to more

frequent soil moisture de

fi

cits during the growing seasons

3

, while

higher off-season heavy rainfalls have raised concerns about
potential erosion and nonpoint source pollution.

Far-reaching agroecological and socioeconomic consequences

of climate extremes demand a systems approach to climate-
resilient agroecosystems

10

. As an integral component of the

complex multidimensional ecological resilience concept increased
agroecosystem resilience to climate perturbations requires greater
emphasis on improving the system

’

s self-regulatory capacity

11

.

Soil, as a central linkage of interconnected atmospheric and ter-
restrial domains, is a primary element of agroecosystem resilience
to adverse climate conditions. The soil system does climate-
buffering through regulating the hydrologic cycle, energy
exchange, and temperature

12

,

13

, thereby alleviating stresses on

plants and biota. The climate adaptability of the soil system is
strongly related to intrinsic attributes while careful land use and
judicious soil management have a documented impact on soil
structural stability, yielding greater permeability during heavy
rainfall, enhanced water storage in drought, and improved gas
exchange for biological respiration and heat acclimation

14

,

15

.

Such management practices offer a great potential to enhance the
soil system

’

s capacity to resist and recover from climate pertur-

bations, thereby maintaining vital soil functions and services.

Among the numerous bene

fi

ts of conservation agricultural

practices such as minimized soil disturbance, cover cropping, and
crop diversi

fi

cation,

16

has suggested that improved soil organic

carbon (SOC), soil aggregation, nutrient cycling, and biodiversity
are the most important factors enhancing soil resilience and
recovery potential, either directly or by invoking other processes.
However, the extent to which those qualities can uphold soil
productivity under various types, severities, duration, and timing
of climate extremes is poorly understood. The historical and
prospective impact of increasing climate extremes on agroeco-
system productivity and environmental services have been
extensively reported in the literature

17

–

19

. There is still a lack of

advanced knowledge, however, about the potential of agricultural
practices to enhance soil system adaptability to climate extremes

6

.

Soil resilience comprises a key component of both

“

soil quality

”

(or

“

soil health

”

) and

“

soil degradability

”

, although certain dis-

crepancies exist. For example, the common soil quality assess-
ment frameworks barely consider the prolonged consequences of
soil degradative processes (e.g., soil erosion) and changes in
short- and long-term soil productive, environmental, and capital
functions when exposed to extreme climate events

6

. In these

models, therefore, anthropogenic soil disturbances such as tillage,
machinery traf

fi

c, or biomass removal are the primary stressors

that agroecosystems and their managers must deal with to
maintain

“

soil quality

”

10

,

20

. Likewise, indicators of soil degrad-

ability rarely address the soil

’

s capacity to recover from envir-

onmental stresses. Therefore, under average climate conditions
with a low proportion of extreme to near-normal years, soil
quality indicators may well represent the sustainability of soil
management systems. However, with the manifest increase in
extreme weather conditions, integrating the soil degradation
process into soil quality can provide operational insight into
whole-system attributes of soil resilience and sustainability.

Several methods have been recommended for assessing soil

resilience, the majority of which are based on soil

“

malleability

”

,

de

fi

ned as changes from the initial state in individual soil quality

indicators after exposure to stresses

16

,

21

,

22

. This approach can be

useful when a static snapshot of a soil-speci

fi

c situation at a given

time and space is desired. However, it neglects many spatio-
temporally dynamic characteristics of the soil system and stres-
sors, as well as the external factors in

fl

uencing soil resilience.

Among those factors is the interaction between multiple soil
properties (soil function), recovery rate (elasticity), terrain char-
acteristics, stressor attributes (e.g., severity, magnitude, timing,
duration, and target), and soil surface cover. A quite different
approach to characterizing soil resilience is evaluating soil capa-
city to sustain its productivity

23

, environmental

24

, and capital

25

functions under adverse climate conditions. The soil productivity
function (e.g., biomass) is a surrogate for the multi-faceted
interactions among soil

–

landscape

–

climate properties and pro-

cesses that determine the capacity of soil to uphold production
stability in unfavorable climates. In addition, soil productivity
functions can represent the effect of physical damage due to cli-
mate extremes on yield. Finally, external inputs such as fertilizers,
lime, and irrigation can enhance production stability in the face
of climate stresses without considerably contributing to the soil
self-regulatory capacity. As such, a measure of soil resilience
based on productivity functions should be informed by detailed
knowledge of the type and level of external inputs.

Climate extremes are complex natural phenomena with a wide

range of de

fi

nitions based on the exposed resource and spatio-

temporal extent. The literature displays large inconsistencies
regarding the conceptual and operational de

fi

nition of drought,

making it dif

fi

cult to compare and synthesize results

26

. De

fi

ning

drought based on the local divergence of climate parameters from
long-term normals can considerably reduce the applicability of
study outcomes across scales. The use of standardized drought
indices has been suggested by

26

as an effective approach to

resolve this issue. Drought indices assimilate data on precipita-
tion, evapotranspiration, stream

fl

ow, and other weather inputs

into a single numerical value, allowing for multi-scalar evaluation
of weather anomalies. These indices are currently used by pro-
grams such as the U.S. Drought Monitor and the European
Drought Observatory for drought monitoring and early warning
at regional to global scales. However, these climate-based early
warning systems may not be well-suited to predict the risk of crop
failure across plants and management systems. Plants vary largely
in drought/wetness tolerance and heat and moisture demand
across phenological stages

27

. The USDM outcomes are utilized by

the United States Department of Agriculture (USDA), the Farm

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

2

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

Service Agency (FSA), and the Internal Revenue Service (IRS) to
support decisions on eligibility for low-interest loans, livestock
forage programs, and tax deferrals

28

. Assessment and prediction

of production vulnerability, and of responsive environmental,

fi

nancial, and social programs require the drought assessment

systems to be calibrated for crop- and management-speci

fi

c cli-

mate adaptation potentials. However, it should be acknowledged
that developing such systems requires information input from
multiple rather than a single long-term experiment as
presented here.

Developing climate-resilient agricultural practices is the anchor

of many national and global climate adaptation frameworks

29

,

30

,

but the recommended practices are often not based on experi-
mental evidence and do not provide operational details, so remain
generic and impractical

6

. The major limitation in developing

evidence-based strategies is the lack of long-term soil and crop
data under diverse management systems across a representative
range of different climate conditions. Long-term

fi

eld experi-

ments with supervised management systems and controlled biotic
stresses

—

e.g., pest, diseases, weed pressure

—

are essential in

understanding the dynamic effect of management systems on soil
adaptive capacity under favorable and unfavorable environmental
conditions.

In this study we address the long-term impacts of various crop-

management systems on the stability of cotton yields in the face

of climate extremes. We do this by utilizing twenty-nine years of
cotton yield data across thirty-two management systems, along
with shorter-term measured soil biophysical properties, and a
suite of the most common drought indices to (1) identify the
most relevant climatic variables explaining long-term

fl

uctuation

in cotton yield; (2) identify the management systems resulting in
the greatest yield stability across extreme climate conditions; and
(3) better understand underlying soil mechanisms controlling soil
resilience and consequently upholding yield stability in extreme
climate conditions. Here we show that among climate indices,
cumulative standardized precipitation index (SPI) and standar-
dized precipitation evapotranspiration index (SPEI) during the
growing season best explains the long-term mean cotton yield
variabilities. Among the raw climate variables, mean soil tem-
perature (MST) during the entire growing season, and especially
low values of MST in August was the most critical yield-limiting
factor. Here we show that throughout the May

–

October growing

seasons, drought severity in August has the most detrimental
effect on cotton yield. Our results suggest that long-term no-
tillage (NT) is a practice that maintains yield stability in unfa-
vorable climate conditions while enhancing crop yield in favor-
able years. Across N rates, legume cover crops

—

particularly hairy

vetch (HV)

—

maintained higher yield than cereal and no cover

crop treatments in near-normal and relatively dry (RD)/wet years.
However, they were not effective in maintaining yield stability in

May

June

Jul

Aug

Sep

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

May

June

Jul

Aug

Sep

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

 

PHDI

 PMDI

 SPEI

 PNP

 SPI

 PET

 Z-index

 DD

 sc-PDSI

 

PDSI

Coefficient of determinaon-R

2

a

)

Planng

Emergence

Fir

st squares

Fir

st Bloom

Open Boll

Peak bloom

Peak bloom

Harv

e

st

Harv

e

st

 

PRCP

 Tmax

 Tmin

 Tmean

 STmax

 STmin

 STmean

Coefficient of determinaon-R

2

b

)

Planng

Emergence

Fir

st squares

Fir

st Bloom

Open Boll

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

 STmin

 STmax

 PRCP

 Tmean

 Tmax

 Tmin

 STmean

Growing season

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

 PHDI

 PDSI

 PNP

 sc-PDSI

 PET

 Z-index

 SPEI

 SPI

 DD

 PMDI

Growing season

Fig. 1 The coef

fi

cients of determination for quadratic regressions, explaining yield variability based on the monthly and seasonal variability of drought

indices and raw climate variables.

The

R

2

values represent the strength of the relationship between annual mean crop yield and monthly and cumulative

growing season (GS) values of

a

drought indices, including standardized precipitation index (SPI), standardized precipitation evapotranspiration index

(SPEI), palmer drought severity index (PDSI), self-calibrated palmer drought severity index (sc-PDSI), palmer hydrological drought index (PHDI), palmer
moisture anomaly index (Z-index), palmer modi

fi

ed drought index (PMDI), and percentage of normal precipitation (PNP) and

b

raw climate variables

including cumulative precipitation (PRCP), maximum temperature (Tmax), minimum temperature (Tmin), mean temperature (Tmean), maximum soil
temperature (STmax), minimum soil temperature (STmin), and mean soil temperature (STmean).

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

3

dry and wet years. High SOC concentration is a common char-
acteristic of high-yielding practices in near-normal and RD/wet
years. However, yield stability in dry and wet years is more
strongly linked to higher wet aggregate stability (WAS) and N
content.

Results

Relationship between climatic factors and long-term yield
variability

. The relationship between the time series of overall

cotton yields and climate variables varied largely across drought
indices/climate indicators and cotton developmental stages
(Fig.

1

). Drought indices with a varying degree of association

showed a quadratic relationship with the long-term mean cotton
yield. In general, the coef

fi

cient of determination was higher for

management systems with higher yields. Among all indices,
cumulative SPI (

R

2

=

0.68) and SPEI (

R

2

=

0.63) during the

growing season exhibited the strongest relationship with the
mean cotton yield across management systems (Fig.

2

). By

accounting for the effect of potential evapotranspiration (PET) in
water balance, SPEI led to one unit greater drought severity
compared with SPI, resulting in the

x

-axis difference in Fig.

2

.

Throughout the May

–

October growing seasons, drought severity

(represented by SPI, SPEI, and Palmer Z-index) in August had
the most detrimental effect on cotton yield (Fig.

1

a). Long-term

cotton yield showed an identical responses to PET and degree
days (DD

60

’

s). The relationship between long-term cotton yield

and raw climate variables including monthly air temperature and
cumulative precipitation was similar to those of drought indices,
but with lower

R

2

values (Fig.

1

.b). However, yield responded

differently to air and soil temperatures. MST during the entire
growing season (

R

2

=

0.45), and especially low MST in August

(

R

2

=

0.38) was the most critical yield-limiting factor among the

raw climate variables. Air temperature in either individual
months or during the entire growing season was poorly related to
the cotton yield. Among the possible 37 input variables, the
MARS procedure selected SPI-GS (SPI in the growing season),
SPEI-GS, and MST-GS to explain 78% of the variability in long-
term cotton yield. Diagnostic plots and model equations are given
in Fig. S1. The selected variables were utilized for the hierarchical
clustering of the experimental years into different climate con-
ditions. The resulting dendrogram grouped experimental years
into six clusters, including dry (D; cluster A; 2 years), RD (cluster
B; 3 years), near-normal (NN; cluster C; 9 years), high-yielding
near-normal (NNhy; cluster D; 7 years), relatively wet (RW;
cluster E; 5 years), and wet (W; cluster F; 3 years) years (Fig. S2).

The impact of management system on yield stability across the
climate conditions

. Long-term mean cotton yield in the current

experiment showed a strong correlation with the published mean
yield across Tennessee (

r

=

0.80) and to a smaller extent with the

mean yield across the U.S. Cotton Belt (

r

=

0.60)

31

(Fig. S3).

Considering the yield outcome in low- and high-yielding years
and disregarding the source of environmental stress, management
systems with NT showed the most consistent effect on yield
stability in poor environmental conditions (Fig.

3

e, f, i, j). NT also

increased yield potential for management systems with legume
cover crops (Fig.

3

a, b). Legume cover crops, in turn, decreased

the yield gap among the management systems with a range of N
application rates (Fig.

3

g, h). Across drought severities, NT

maintained a higher cotton yield than CT. The cover crop effect
was signi

fi

cant in RD, near normal (NN), and near-normal high

yield (NNhy) years, but not in dry (D) and wet (W) years. The N
rate and N rate by cover crop interaction were signi

fi

cant in all

climates except for the dry years (Table

1

and Table

2

). In dry

years, where the drought severities exceeded

−

3.4 (SPI) and

−

5.1

(SPEI) during the growing seasons, cotton yields under all
management systems fell below 50% of the global mean yield
(1 Mg ha

−

1

) (Fig.

4

). Within the lower 10 percentile of global

mean yield (up to 0.57 Mg ha

−

1

), NT signi

fi

cantly lowered the

risk of yield losses compared with CT, whereas the effect of cover
crop, N rate, and their interaction was not signi

fi

cant (Table

1

). In

general, the lowest probability of crop failure was observed within
RD, NN, and NNhy climates with drought/wetness severities
ranging between SPI scores of

−

1.7 to 0.46 and SPEI scores of

−

2.2 to 0.78 (Fig.

4

). Within these climates, higher N rates and

legume cover crops predominantly controlled higher yield out-
comes. However, in extreme climates, both in dry (D) and wet
(W) years, NT was the management that maintained higher
cotton yield, as shown by the

F

-values in Table

2

. The most

signi

fi

cant contribution of NT to higher yield stability occurred in

wet climates, where NT systems yielded on average 22% more
cotton than CT systems. In W (SPI averaging 2.92 and SPEI
averaging 3.61) and RW (SPI averaging 1.27 and SPEI averaging
1.90) years, the cover crop effect was not signi

fi

cant, but a

minimum N rate of 34 kg ha

−

1

remained critical to maintaining

yield stability. Figure

5

demonstrates the yield bene

fi

ts of tillage,

cover crop, and N rate compared with the corresponding control
treatments in different climates. In all climates except for wet
years, HV planted on either CT or NT resulted in greater yield

-5

-4

-3

-2

-1

0

1

2

3

4

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

-6

-5

-4

-3

-2

-1

0

1

2

3

4

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

b)

R

2

 

= 0.67

Yield = 1.14 - 0.02x - 0.04x

2

Coon Yield (Mg ha

-1

)

Standardized Precipitaon Index

a)

Coon Yield (Mg ha

-1

)

Standardized Precipitaon Evapotranspiraon Index

R

2

 

= 0.62

Yield = 1.12 - 0.02x - 0.02x

2

Fig. 2 Quadratic relationships between drought severity and detrended
cotton yield during the growing seasons for the 32 management
systems.

The drought severity is represented by

a

standardized precipitation

index (SPI) and

b

standardized precipitation evapotranspiration index (SPEI).

The arched dashed line represents the relationship between mean cotton
yield and drought severity. The shaded band around curves represent the
95% con

fi

dence limits.

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

4

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

bene

fi

ts than NC (control no cover crop treatment). CC, the other

legume cover crop, outperformed NC in D, RD, NNhy, NN, and
RW climates, but only when it was planted in an NT system.
However, management systems with both HV and CC cover
crops exhibited lower yields than those with NC in wet years.

Relationship between soil resilience indicators and yield sta-
bility across drought and wetness severities

. Among soil resi-

lience

metrics,

SOC provided

the strongest

explanation

(

R

2

=

0.87) for yield stability under a range of drought and

wetness severities (Fig.

6

). ANOVA tests showed that intercepts,

0.0

0.5

1.0

1.5

2.0

2.5

 HV

 CC

 WW

 NC

Convenonal llage

a)

No-llage

b)

0.0

0.2

0.4

0.6

0.8

1.0

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

0.0

0.5

1.0

1.5

2.0

 0

 34

 67

 101

c)

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

d)

0.0

0.2

0.4

0.6

0.8

1.0

  C

u

mulave

 distribuon

0.0

0.5

1.0

1.5

2.0

2.5

 NT

 CT

Probability density

e)

Non-Leguminous

Leguminous

0.0

0.2

0.4

0.6

0.8

1.0

f)

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

0.0

0.5

1.0

1.5

2.0

2.5

 0

 34

 67

 101

g)

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

h)

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.5

1.0

1.5

2.0

2.5

High N rate

 NT

 CT

i)

Low N rate

0.0

0.2

0.4

0.6

0.8

1.0

j)

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

0.0

0.5

1.0

1.5

2.0

2.5

 HV

 CC

 WW

 NC

k)

0.

2

0.

4

0.

6

0.

8

1.

0

1.

2

1.

4

1.

6

1.

8

2.

0

Yield (Mg ha-1)

0.0

0.2

0.4

0.6

0.8

1.0

l)

Fig. 3 Cumulative distribution (lines) and probability density (shaded areas) of cotton yield across low- and high-yielding environmental conditions.

The graphs represent the interactive effect of the tillage system (

a

–

d

), cover crop type (

e

–

h

), and N rate (

i

–

l

) with the other two factors. CT conventional

tillage, NT no-tillage, NC no cover crop, WW winter wheat, CC crimson clover, HV hairy vetch.

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

5

linear and quadratic slopes, and their interactions with climates
differed (

P

< 0.001) among climate conditions (Table

2

). The

linear slope indicated 17 and 16% decrease in the probability of
yield losses per unit increase in SOC from 7 up to 10 g kg

−

1

in RD

and NNhy climates, respectively. In RW and NN climates a unit
increase in SOC lowered the probability of yield losses by 12 and
10%, respectively. In W climate the corresponding decrease in the
risk of yield losses was 5%, while SOC had a negligible effect on
increasing yield stability in D climate. Total bound N (TN

b

)

explained 78% of the variability in cotton yield stability across
climate conditions (Fig.

6

), with the intercept, common linear

slope, and linear by climate interaction differing (

P

< 0.05) among

climate conditions. Although the common quadratic slope was
signi

fi

cant, it did not considerably relate to yield stability in the

observed range of TN

b

. The effect of TN

b

on yield stability in

different climates revealed an almost similar effect to that seen in
SOC. During NNhy and RD years, greater TN

b

levels effectively

reduced the risk of yield losses by 6.0 and 5.9% per unit increase
in TN

b

, respectively. In contrast to SOC, increases in TN

b

con-

centration up to 200 mg kg

−

1

resulted in 5.0 and 5.3% reductions

in the probability of yield losses in D and W years, respectively.
Small differences in linear slopes among climate categories indi-
cated that of the soil resilience metrics, WAS was the most critical
factor in decreasing the risk of crop failure in the most extreme
(W and D) climates. WAS explained 71% of the variability in
cotton yields stability across climate conditions. Within the WAS
values from 25 to 39%, a unit increase in WAS lowered the

probability of yield loss by 2.4% in NNhy and RD climates. The
corresponding values for RW, NN, W, and D climates were 2.2,
2.2, 2.1, and 2.0%, respectively (Table S1). The POxC and MB
exhibited similar relationships with the yield stability; with
increasing POxC and MB values a decreasing common linear
yield loss slope was evident (

P

< 0.05), but the slopes did not

differ among different climates. DOC did not signi

fi

cantly relate

to yield stability.

Discussion

The rainfed water supply and controlled biotic stresses in this long-
term experiment provided a unique opportunity to observe the
direct impact of true climate perturbations on yield. In general, soil
temperature anomalies in August (optimum MST

27

–

31 °C) and

rainfall de

fi

ciency in August (optimum

100

–

140 mm) and to

some extent in June (optimum

50

–

100 mm) were the most cri-

tical yield-limiting climate factors. Considering the cotton growth
life cycle in the study region, August is the period of peak bloom to
open boll formation, which coincides with peak water consumption
by cotton plants due to high transpiration and evaporative water
demand. Along with the suf

fi

cient water supply, adequate accu-

mulation of heat units upon the appearance of

fi

rst bloom is

essential for achieving boll maturity and

fi

ber quality

32

. Likewise,

the

fi

rst square to

fi

rst bloom development occurs mostly within

June and is a critical time for avoiding water de

fi

cit stress. Cotton

vegetative growth is rapid within the three weeks following the

Table 1 Means for cotton yield (Mg ha

−

1

) under tillage, cover crop, N rate (nitrogen rate), and N rate × cover crop interaction in

dry (D), relatively dry (RD), near-normal (NN), high-yielding near-normal (NNhy), relatively wet (RW), and wet (W) years and
across all climate conditions.

Treatments

D

RD

NN

NNhy

RW

W

Across climates

Number of years

2

3

9

7

5

3

Tillage means

CT

0.48b

1.10b

1.08b

1.07b

0.91b

0.68b

1.00b

NT

0.52a

1.14a

1.16a

1.11a

0.97a

0.83a

1.04a

Cover crop means

NC

0.50a

1.05c

1.10c

1.04c

0.92a

0.78a

0.98c

WW

0.48a

1.04c

1.05c

1.03c

0.97a

0.78a

0.97c

CC

0.50a

1.15b

1.14b

1.10b

0.92a

0.73a

1.03b

HV

0.53a

1.24a

1.20a

1.18a

0.96a

0.75a

1.10a

N rate means

0

0.48a

0.96c

0.95c

0.94b

0.87b

0.66b

0.88c

34

0.51a

1.10b

1.13b

1.12a

0.96a

0.79a

1.03b

67

0.52a

1.20a

1.19a

1.13a

0.97a

0.78a

1.08a

101

0.49a

1.20a

1.22a

1.16a

0.98a

0.80a

1.09a

Cover × N rate means

0 × NC

0.43a

0.82f

0.84f

0.85f

0.75e

0.61e

0.77f

0 × WW

0.46a

0.80f

0.76f

0.78f

0.81de

0.63e

0.74f

0 × CC

0.48a

1.03de

1.10de

1.02e

0.91bcd

0.70cde

0.96e

0 × HV

0.53a

1.18abc

1.12cde

1.11cde

1.01ab

0.71bcde

1.05cd

34 × NC

0.50a

1.01e

1.06e

1.01e

0.96b

0.77abcd

0.96e

34 × WW

0.49a

1.02de

1.04e

1.06de

0.94bc

0.76abcde

0.97e

34 × CC

0.53a

1.12cd

1.17bcd

1.12bcd

0.93bc

0.78abcd

1.04cd

34 × HV

0.55a

1.27a

1.27a

1.27a

1.01ab

0.84abc

1.16a

67 × NC

0.55a

1.19abc

1.24ab

1.12bcd

0.96bc

0.86ab

1.09bc

67 × W

0.50a

1.19abc

1.20ab

1.13bcd

1.02ab

0.85abc

1.08bc

67 × CC

0.50a

1.20abc

1.10de

1.09cde

0.92bcd

0.65de

1.02de

67 × HV

0.52a

1.24ab

1.22ab

1.20ab

0.99ab

0.76abcde

1.11ab

101 × NC

0.51a

1.15bc

1.24ab

1.16bcd

1.01ab

0.86ab

1.10abc

101 × WW

0.47a

1.16bc

1.21ab

1.16bc

1.10a

0.87a

1.10abc

101 × CC

0.47a

1.25ab

1.22ab

1.17bc

0.94bc

0.79abcd

1.09abc

101 × HV

0.51a

1.25ab

1.19abc

1.14bcd

0.85cde

0.69cde

1.07bcd

CT

conventional tillage,

NT

no-tillage, NC no cover crop (control),

WW

winter wheat,

CC

crimson clover,

HV

hairy vetch. Means with no common letter in a column within the means of tillage, cover crop,

N rate, or cover × N rate are signi

fi

cantly different based on the protected LSD at

P

< 0.05.

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

6

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

emergence of the

fi

rst square, requiring rapid uptake of phos-

phorous (P) and potassium (K) through the aqueous medium

33

.

Water de

fi

ciency within this period can adversely affect the biomass

production and the number of fruiting sites as well as the earliness,
which is particularly important in short-season environments. The
effect of sub-seasonal and seasonal climate

fl

uctuations on the

variability of cotton yield was well-addressed by SPI, SPEI, and
MST indices/indicators. Management systems with higher global
yields generally showed a stronger relationship with climate indices.
This indicates that high-yielding management systems diminish the
adverse effect of other environmental stresses on yield and boost
greater yield outcomes with improving climate conditions. SPEI
index includes the PET, but this did not improve the predictability
of yield variability over that from SPI. This can perhaps be attrib-
uted to the indeterminate relationship between PET and DD

60

’

s.

Though rising air temperature during critical cotton growth stages
may enhance the evaporative loss of water, suf

fi

cient heat unit

accumulation (DD

60

’

s

2300) is essential to supporting many

physiological processes that eventually determine cotton yield. As
shown in Fig.

1

, PET, DD

60

’

s, and Tmean across cotton growth

stages demonstrated weak but very similar relationships with yield
variability. The quadratic relationship between PET and yield
variability implies that, for example, in August, there is an optimum
range of PET that maintains a high yield, while lower and higher
values lead to yield losses. This is unlikely however that low PET in
an RD month like August in such a rainfed system causes yield
decline. This inconsistency is likely since low PET means low degree
days due to the similar effect of air temperature on PET and DD

60

’

s,

but their contrasting effects on yield outcome. Alternately, both
seasonal and sub-seasonal variations in MST were more strongly
related to the yield variability (Fig.

1

). In addition to the incidental

solar radiation, the soil temperature is strongly affected by soil
biophysical properties such as texture, compaction, organic matter
content, and surface cover. Management can modify the majority of
these factors, thereby altering the soil water regime, so soil tem-
perature variations may not be well represented by air temperature
variations

34

. Maximum yield outcome was obtained when MST

within the time from

fi

rst bloom to open boll development varied

between 27 and 31 °C. Soil temperatures outside this range caused
signi

fi

cant yield reductions. Residue cover left by conservation til-

lage systems such as NT can effectively buffer extreme summer heat
by insulating the soil surface and enhancing albedo

35

, though this

effect may gradually diminish throughout the growing season

36

,

37

.

Accelerated residue retention by cover crops can further strengthen
the heat buffering capacity of NT soil. Moreover, reduced soil
manipulation stimulates organic matter accumulation in the topsoil
and increases micropore volume, which collectively contributes to a
higher soil water retention capacity

38

,

39

. Except for the two

extremes of near-saturated and near-dry conditions, within a wide
range of soil moisture contents water and vapor movement occurs
simultaneously with temperature gradient and vice versa. A wetter
soil by virtue of its greater thermal conductivity and heat capacity
can decrease the diurnal temperature cycle amplitude. Therefore,
soil temperature may relatively well represent the management
effect on heat regulation. However, further long-term temperature
monitoring at different soil depths is needed to fully comprehend
the role of conservation practices on soil heat regulations and their
impact on yield variability.

NT was the primary management practice boosting yield in

favorable climates and reducing yield losses in unfavorable cli-
mate conditions. Even though the yield bene

fi

t of NT compared

with CT was smaller in generally favorable climate conditions
than the effect of higher N rates (67 and 101 kg ha

−

1

) or legume

vs non-legume cover crops, NT exhibited a consistent yield
bene

fi

t in all climate conditions (Table

1

and Fig.

5

). In addition,

NT increased the general effectiveness of cover crops

—

in

Table

2

Analysis

of

Variance

for

the

effect

of

tillage,

cover

crop,

N

rate

(nitrogen

rate),

tillage

×

cover

crop,

tillage

×

N

rate,

N

rate

×

cover

crop,

and

tillage

×

N

rate

×

cover

crop

interactions

on

cotton

yield

production

in

dry

(D

),

relatively

dry

(RD

),

near-normal

(NN),

high-yielding

near-normal

(NNhy

),

relatively

wet

(RW),

and

wet

(W)

years

and

across

all

climate

conditions.

Number

of

years

D

R

D

N

N

NNhy

RW

W

All

climates

23

9

7

5

3

F

-value

P

-value

F

-value

P

-value

F

-value

P

-value

F

-value

P

-value

F

-value

P

-value

F

-value

P

-value

F

-value

P

-value

ANOVA

P-values

Tillage

9.79

0.0023

4.84

0.0303

31.08

<0.0001

8.21

0.005

8.77

0.003

24.96

<0.0001

12.37

7.00E

−

04

Cover

crop

1.59

0.1967

27.96

<0.0001

19.37

0.0003

12.7

0.001

1.4

0.243

0.74

0.5309

17.79

4.00E

−

04

N

rate

1.56

0.2039

43.26

<0.0001

66.92

<0.0001

40.85

<0.0001

5.95

6.00E

−

04

5.14

0.002

69.17

<0.0001

T

×

C

0.62

0.6027

0.7

0.5567

2.08

0.1088

1.4

0.2479

0.15

0.9296

0.51

0.6734

1.79

0.1553

T

×

N

0.11

0.9536

0.26

0.8544

1.97

0.1247

2.04

0.1144

2.82

0.04

0.96

0.4138

2.04

0.1145

C

×

N

0.72

0.6859

4.38

<0.0001

12.63

<0.0001

6.2

<0.0001

4.97

<0.0001

1.87

0.056

12.27

<0.0001

T

×

N

×

C

0.42

0.9229

0.46

0.8971

0.95

0.4866

0.93

0.5077

0.35

0.9584

0.43

0.921

0.9

0.5272

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

7

particular legumes

—

in enhancing the yield outcome across cli-

mate conditions

40

(Fig.

5

). Moreover, NT considerably increased

N use ef

fi

ciency compared with CT. For example, management

systems with 34 kg N ha

−

1

on NT delivered on average greater

yield bene

fi

ts than systems with 67 kg N ha

−

1

on CT (Fig.

5

).

Reducing N fertilization by 50% while maintaining a high yield in
all climate conditions (except for RD) has profound environ-
mental and economic implications. For instance, it has been
observed that net agronomic nitrous oxide emissions were gov-
erned by N fertilizer rate

41

. However, in the present study, the

major N bene

fi

t of legume cover crops was associated with RD

and near-normal years. In wet (W) and relatively wet (RW) years
WW with the maximum N rate (101 kg ha

−

1

) was associated with

higher yields (Table

1

). This is likely due to higher N loss through

leaching, denitri

fi

cation, and runoff in wet conditions. All wet

(W) years in this study were characterized by heavy early growing
season rainfall, particularly in May and June. Excessive early
season rainfall

—

especially under cover crops with high residue

mulch

—

is likely to adversely affect crop stand and reduce yield

(Fig.

5

). It is tied with low soil temperature due to lower heat

exchange and high moisture content, which may hinder timely
seed germination or plant establishment. Also, among wet years,
late summer and fall storms in 2018 starting in September
delayed harvest and caused major boll weathering and yield
decline. Generally, the increasing global atmospheric temperature
has been assumed to bene

fi

t crop production by extending

growing seasons in some regions

42

. However, increasing surface

temperatures evaporate more water into the atmosphere, which
may cause yield losses by spawning more intense late summer
and fall storms

43

–

45

. Combining WW cover crop either with NT

or CT did not prompt signi

fi

cant yield bene

fi

ts greater than NC

in any climate condition. In general, the yield on CT management
in most climate conditions remained highly dependent on
enhanced N fertilization. The only exception was dry years, in

which enhancing N rate beyond 67 kg ha

−

1

decreased yield in

both CT and NT systems. A strong relationship between long-
term mean cotton yield variability in the current rainfed experi-
ment and the broader regional yield data indicates that yield at a
regional scale is still dictated by climate conditions, especially in
wet and dry years (Fig. S3). Although this may not be fully
addressed through a single long-term experiment, this relation-
ship signi

fi

es the importance of climate early warning systems

that account for the crop- and management-speci

fi

c vulner-

abilities to climate extremes across the phenological stages. Thus,
the precise auxiliary practices (e.g., irrigation, drainage, fertiliza-
tion, cover crop planting and termination, and harvest schedul-
ing) could be implemented on time to avoid considerable
economic losses.

Our results suggest that greater SOC concentration is key for

improved yield stability in NN and RD and wet climate condi-
tions. The minimum threshold of 10

–

15 g kg

−

1

SOC suggested

by

46

to maintain soil resilience, restore soil biodiversity, and

improve structural stability is in perfect accord with our obser-
vations (Fig.

6

). However, there was a minimal effect of higher

SOC on yield stability in dry and wet years (Table S1). Improved
soil structural stability as represented by WAS proved to be the
primary soil attribute linked to maintained yield stability under
extreme climates, particularly in wet years. Soil aggregates phy-
sically protect terrestrial carbon pool and the level of aggregation
and their stability determine the severity of

fl

uvial, aeolian, and

gaseous soil, carbon, and nutrient losses

47

. Climate forces such as

kinetic rainfall energy or aeolian abrasion are known to disrupt
soil aggregates. The resultant particles are prone to transport,
carrying nutrients and C, and forming surface crusts. This process
can promptly and chronically decrease productivity by reducing
available water and nutrients, depleting SOC, and impeding air
and water circulation. The positive relationship between SOC and
WAS in this study (

r

=

0.65) corroborates previous

fi

ndings on

HV

CC

WW

NC

HV

CC

WW

NC

101

67

34

0

101

67

34

0

101

67

34

0

101

67

34

0

101

67

34

0

101

67

34

0

101

67

34

0

101

67

34

0

CT

NT

0.00 0.25 0.50 0.75

0.00 0.25 0.50 0.75

NN

0.00 0.25 0.50 0.75

< 90

th

 percenle

< 50

th

 percenle

RD

< 10

th

 percenle

0.00 0.25 0.50 0.75 1.00

Cumulave probability

0.00 0.25 0.50 0.75

W

NNhy

D

RW

0.00 0.25 0.50 0.75

Fig. 4 Cumulative probability of cotton yield losses based on Gaussian distribution below 10, 50, and 90th percentile of the global mean yield across
climate conditions.

The plots represent the associated values in dry (D), relatively dry (RD), near-normal (NN), high-yielding near-normal (NNhy),

relatively wet (RW), and wet (W) years. CT conventional tillage, NT no-tillage, NC no cover crop (control), WW winter wheat, CC crimson clover, HV
hairy vetch. SPI and SPEI values during the growing seasons. The blue, yellow, and gray bars represent the cumulative probability of yield losses below 10,
50, and 90th percentile of the global mean yield, respectively.

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

8

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

the linkage between SOC and WAS

38

,

47

,

48

. Soil total bound

nitrogen, TN

b

was strongly related to yield stability in extreme

climate conditions (Fig.

6

and Table S1). Besides being an

essential nutrient for plant growth, nitrogen along with phos-
phorus and sulfur are important elements in the process of bio-
mass humi

fi

cation

49

, which in turn contributes to the longevity of

soil structural stability and resilience. Retention and supply of N
in acute climate conditions may increase soil resilience through
multiple processes, such as the regulation of water and nitrate
in

fl

ux and ef

fl

ux rate, maintaining air

–

water balance for con-

tinued N mineralization and nitri

fi

cation, and heat regulation in

dry periods to avoid N volatilization. The maintenance of the N
balance is thus an important indicator of climate-resilient
agroecosystems.

Previous studies have demonstrated the regulatory in

fl

uence of

conservation practices on key soil functions

—

e.g., water, heat,

and nutrient

fl

uxes, structural hierarchy and stability, carbon

sequestration, and biological vitality and diversity

15

,

37

,

38

,

50

. Such

a soil system is expected to provide improved resilience to
extreme climate events such as torrent rainfalls, heat waves, and
droughts, thereby reducing the climate stresses on plants and
stabilizing crop yield. The effect of improved overall soil health on
crop yield has been shown to become more evident over time

51

.

However, there is a lack of suf

fi

cient long-term crop yield data

from different conservation practices exposed to a representative
range of climate conditions. Consequently, there is little evidence-
based information regarding the role of conservation practices in
the climate adaptability of agroecosystems. Our study provides

insight into the long-term interaction between climate extremes
and yield stability under 32 tillage, cover crop, and nitrogen rate
practices. We demonstrate how management-induced changes in
soil condition, control soil resilience which in turn supports yield
stability during unfavorable climate conditions. Our

fi

ndings

indicate that long-term NT provides a consistent basis for greater
yield stability across climate conditions. In addition, NT improves
the effectiveness of cover crops, particularly legumes, in drought
mitigation and crop yields maintenance. Moreover, NT provides
those bene

fi

ts with relatively lower dependence on nitrogen fer-

tilizer which has important environmental and economic impli-
cations. We de

fi

ned the drought and wetness severities based on

standardized drought indices to enable cross-scale evaluation,
comparison, and synthesis of results. Demonstrated crop-and
management-speci

fi

c vulnerability of crop yield to climate

extremes provides insight into climate-resilient agroecosystem
management strategies. However, it should be acknowledged that
these results represent a limited geographic scale and soil resi-
lience has been quanti

fi

ed based on a limited number of variables.

Therefore, further long-term studies are required for improved
quanti

fi

cation of the soil resilience concept and its linkage with

crop yield stability under climate extremes.

Methods

Long-term experiment

. This study was conducted using data from a long-term

rain-fed plot-scale cotton experiment at the University of Tennessee

’

s West Ten-

nessee AgResearch and Education Center in Jackson, TN, USA (35°37

′

N: 88°51

′

W,

altitude 113 m). The humid subtropical (Köppen classi

fi

cation) climate of the study

region is characterized by a mean annual temperature of 15.5 °C and mean annual
precipitation of 1375 mm. The experiment was originally initiated in 1981 on a 2-
hectare area. The complete set of current treatments has been in place since 1984.
The experiment is located on Lexington silt loam (

fi

ne-silty, mixed, thermic Ultic

Hapludalfs) with a

fi

eld slope ranging from 0 to 2%. The experiment layout

comprises a cross-classi

fi

ed arrangement of four nitrogen (N) fertilizer rates: 0, 34,

67, and 101 kg Na ha

−

1

; four cover crop treatments: winter wheat (WW,

Triticum

aestivum L

.), HV (

Vicia villosa

), crimson clover (CC,

Trifolium incarnatum L

.), and

no cover crop (NC,

“

winter weeds

”

); and two tillage systems: conventional tillage

(CT) and NT; thus totaling 32 treatment combination with four replications
established in 128 plots in a randomized complete block design with split

–

split

plot. The experiment was initiated in 1981, and the complete set of current
treatments have been in place since 1984. Cover crops are terminated in late April
to early May with 3.51 L ha

−

1

of Gramoxone

®

SL 2.0 (N, N-dimethyl-4,4-bipyr-

idinium dichloride) (Syngenta Global). Within 10 days, plots are broadcast with
ammonium nitrate (NH

4

NO

3

) at the assigned rates. The CT treatments are double

disked to a depth of 10 cm and harrow leveled to prepare the seedbed. Cotton is
seeded uniformly in all plots. Cotton cultivars planted through the course of the
experiment included Deltapine 50, Stoneville 825, Stoneville 474, Deltapine 425,
Deltapine 451, and Phytogen 375. In October of each year the cotton was har-
vested, ginned, and lint yield was recorded.

Characterization of climate extremes

. Historical daily precipitation and air

temperature data were recorded by a meteorological station at the experimental
site. Weather anomalies were characterized using both raw climatic data and
drought indices. To calculate the drought indices, we used version 1.0.4 of the

“

climate-indices

”

software package for the Python language

52

. The drought indices

used, included SPI, SPEI, palmer drought severity index (PDSI), self-calibrated
PDSI (sc-PDSI), palmer hydrological drought index (PHDI), palmer moisture
anomaly index (Z-index), palmer modi

fi

ed drought index (PMDI), and percentage

of normal precipitation (PNP) (

fi

g. S5). The SPI is the number of standard

deviations of the cumulative precipitation from the long-term median. To calculate
SPI, a three-parameter gamma probability density function is

fi

tted to the time

series of precipitation data. The associated probability density distribution is then
transformed to a normal distribution

53

. The SPI can be calculated at time scales

from ¼- to 48-months to re

fl

ect the availability of different water resources. We

calculated SPI for ¼- and 1-month scales since plants respond relatively quickly to
precipitation anomalies. The SPEI follows the same computational procedure as
SPI but additionally accounts for the effect of temperature variability in drought
assessment. The PET required by SPEI for the assessment of water balance was
calculated using the Thornthwaite

54

model. The PDSI is another well-known index

estimating the duration and severity of drought based on monthly precipitation
and temperature data as well as the soil

’

s plant available water content (PAWC). In

this study, the mean PAWC value (0.18 cm

3

cm

−

3

) across management systems

was obtained from

15

. The sc-PDSI, PHDI, and Z-index

55

,

56

are variants of the

PDSI index with an improved analytical capacity to capture rapidly emerging
drought and moisture anomalies, long-term drought effect on water storage and

101

67

34

WW

CC

HV

-15 -10

-5

0

5

10

15

20

25

30

N rate effect

 D

 RD

 NN

 NNhy

 RW

 W

Convenonal llage

Cover crop effect

g

g

a

)

101

67

34

WW

CC

HV

-15 -10

-5

0

5

10

15

20

25

30

Effect of treament levels (% difference with control)

No-llage

b)

Fig. 5 Yield bene

fi

ts of N rates and cover crops compared with

associated control treatments (0 N rate and no cover crop).

Values are

reported across dry (D), relatively dry (RD), near-normal (NN), high-
yielding near-normal (NNhy), relatively wet (RW), and wet (W) climate
conditions under

a

conventional tillage and

b

no-tillage.

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

9

stream

fl

ow, and extreme weather events in heterogeneous environments including

multiple micro-climatic sub-regions. Soil temperature measured continuously at
10 cm depth was also considered as a variable not associated with any speci

fi

c

index. Finally, cotton growth stages are often de

fi

ned as the number of days after

planting

—

but achieving proper growth at each milestone also strongly depends on

the number of days with optimum air temperature. Degree days (DD

60

’

s) is a

measure of heat units that plants have accumulated before they undergo defoliation
and is often used as an indicator of plant maturity. Therefore, DD

60

’

s was also

calculated as the cumulative average daily air temperature minus 15.6 °C (mini-
mum air temperature in which cotton growth occurs) over monthly periods.

Quanti

fi

cation of soil resilience

. The multi-disciplinary nature of resilience in

ecosystem studies has led to a broad spectrum of de

fi

nitions and assessment

metrics, but there is a lack of proper quanti

fi

cation of agroecosystem resilience,

particularity from the soil perspective. We used a suite of soil physical, chemical,
and biological properties commonly used for soil health assessment to foster
insight into the underlying mechanisms governing soil resilience and consequently
yield stability in extreme climatic conditions. The soil properties included SOC,
permanganate oxidizable carbon (POxC), dissolved organic carbon (DOC), WAS,
totally bound nitrogen (TN

b

), and microbial biomass (MB). During the growing

season of 2019, a composite soil sample including three random sub-samples was
collected from each plot (

n

=

128) at 0

–

10 cm soil depth. Sampling was conducted

by a hand auger except for the aggregate stability analysis, for which clods were
collected by hand and fractured gently into pieces to obtain the natural soil
aggregates. The SOC was measured by dry combustion method at 950 °C using a
CN analyzer (Elementar vario SOC cube in solid mode, Hanau, Germany). Fol-
lowing the method of

57

, POxC was measured by adding 20 mL of 0.02 M potassium

permanganate (KMnO

4

) to 2.5 g of soil sample. The solution was shaken at 120

rpm for 2 min and allowed to settle for 8 min. The supernatant was extracted and
diluted to 1/100. Finally, the absorbency rate of the dilute solution was measured at
550 nm by a microplate reader (Biotek Instruments, Inc., Vermont, USA). To
measure DOC, 10 g soil sample was mixed with 40 mL de-ionized water. The
solution was shaken for 30 min, then centrifuged at 3500 rpm for 20 min and
passed through a 0.45 µm membrane. The DOC concentration was measured by a
SOC analyzer following the method of

58

(Shimadzu Instruments, Inc., North

America). The WAS was measured using a wet sieving apparatus (08.13, Eijkelk-
amp Agrisearch Equipment, Giesbeek, the Netherlands). Natural air-dried soil
aggregates sized between 1 and 2 mm were placed on 0.26 mm sieves and wetted
through capillary water motion from a wet cloth placed on the sieve bottom.
Samples were raised and lowered in distilled water for 3 min with a frequency of 35
times per min

59

. Thereafter, the unstable fraction was oven-dried and weighed and

stable fractions were corrected for the sand content using 2 g L

−

1

of sodium

hydroxide as a dispersant. The fraction of water-stable aggregates was calculated as

WAS

¼

on sieve

sand

ð

Þ

total sample

sand

´

100

ð

1

Þ

Fig. 6 Indicator variable regression (IVR) shows the probability of yield losses below the 50th percentile of the global yield.

Relationships are shown

across dry (D), relatively dry (RD), near-normal (NN), high-yielding near-normal (NNhy), relatively wet (RW), and wet (W) years vs.

a

soil organic carbon

(SOC g kg

−

1

),

b

microbial biomass (MB

µ

g g

−

1

),

c

wet aggregate stability (WAS %),

d

total bound N (TN

b

mg kg

−

1

), and

e

permanganate oxidizable carbon

(POxC mg kg

−

1

).

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

10

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

For MB, soil sampling was conducted in a relatively wet soil condition. Samples
were not sifted unless they displayed rootlets or debris, in which case these were
removed by sifting through a 3 mm soil sieve. Samples were analyzed in triplicate
and the average CV was less than 10%. The assay was by the microBIOMETER

©

(Proli

fi

c Earth Sciences, Montgomery, NY, USA) kit protocol using kit reagents

and supplies. 0.5 mL of soil was packed into the microBIOMETER soil sampler and
added to the 10 mL of the high salt extraction solution which loosens the attach-
ment of the microbes to soil particles. Prior to whisking for 30 seconds, the clay
plug of soil was broken up by mechanically stabbing the base of the tube with a
stainless-steel spatula. Samples were allowed to settle for 20 min, during which time
the soil precipitates, and the

fl

uid above the soil is colored with the microbes which

contain soil pigments. Six drops from the top 13 mm of the extract in the tube were
applied to the test card containing a membrane that traps the microbes on the
surface while allowing the

fl

uid to wick away. The intensity of the color on the

membrane correlates with MB. MB was estimated using the microBIOMETER app
on a cell phone. Brie

fl

y, the phone takes a picture, and the MB is estimated from

the intensity of the color of the test card. At the same time as the microBIOMETER
test was performed, each extraction (three samples) was analyzed by digital
microscopy by making three slides and taking ten pictures of each slide, and
analyzing by particle size. Every soil sample was thus analyzed by microscope and
test card a total of nine times. The microscopic data showed an excellent corre-
lation with the test card data.

Statistical analyses

Yield stability under extreme climate conditions

. The yield data for the

fi

rst 5 years

(1984

–

1988) of the experiment was removed to only consider the more stable

condition of management systems following establishment. So, we considered yield
data from 1989 to 2018 for

n

=

3712

51

. The only exception was 1993, for which

yield information is not available (Supplementary Data 1). To assess seasonal and
sub-seasonal drought in

fl

uences on linearly de-trended cotton yield, all drought

indices, and soil- and plant-speci

fi

c drought indicators were calculated at a

monthly time scale and also accumulated over the growing season (May 1 to
September 30) to obtain the

“

drought severity

”

60

. The resulting 53 variables were

subjected to principal component analysis to reduce multicollinearity. Due to the
non-linear relationships between yield and the drought indicators/indices, we used
multivariate adaptive regression splines (MARS) models

61

to identify the most

relevant variables explaining yield

fl

uctuations. The MARS model is a non-

parametric approach to capture the nonlinearity aspect of singular and interactive
features in multiple regression by assessing cut points (knots). We used the

“

earth

”

function in

R

62

to apply a threefold generalized cross-validation procedure for

“

forward

”

pruning of features to avoid over

fi

tting. The variable selection procedure

was continued until the nature of non-linearity and an optimum number of knots
were identi

fi

ed. The model selected four variables to build a model with

fi

ve

“

basis

functions

”

(including the intercept) up to the second degree of interactions (Fig.

S1). We then applied a hierarchical clustering with the

“

Ward

”

partitioning algo-

rithm to cluster the experimental years based on the drought variables selected by
the MARS procedure. The optimum number of clusters was determined using the

“

elbow

”

method, which restrains the number of clusters as a function of the total

within-cluster sum of squares (WSS), such that adding more clusters does not
improve WSS. The resulting six clusters were marked as two clusters of near-
normal years, and one cluster each of RD, dry, relatively wet, and wet years. The
yield stability was assessed separately within each of these climate clusters. We
de

fi

ned yield stability as the probability that yield of the management system

i

falls

below the speci

fi

ed level of

λ

(10, 50, or 90th percentile of the global mean yield) in

a given environment

j

. This probability was calculated using the equation

Pr

y

ij

≤

λ

¼

Z

λ

1

=

σ

ij

ffiffiffiffiffiffi

2

π

p

exp

λ

μ

i

2

=

2

σ

2

ij

h

i

d

λ

ð

2

Þ

where

μ

i

is the mean yield and

σ

ij

is the variance of the

i

th system. The normality of

yield distribution for all management systems was satis

fi

ed based on the

Shapiro

–

Wilk test (

W

> 0.9).

Soil resilience control on yield stability across climate conditions

. Using

indicator variable regression in SAS

63

, the effect of a suite of soil quality indicators

(Supplementary Data 2) on the probability of yield losses below the

fi

ftieth per-

centile of global mean yield was determined under each climate condition. Indi-
cator variable regression combines analysis of variance and regression. In this
process (also known as dummy regression), a regression line for each treatment
level is modeled, allowing slopes and s to be compared. Here, the climate condi-
tions (e.g., dry and wet) were considered as categorical variables. A full model with
linear, quadratic, and all interactions with treatments was executed, giving

Pr

y

ij

≤

λ

¼

μ

þ

CL

i

þ

β

1

´

SP

jj

þ

β

2

´

SP

jj

´

CL

i

þ

γ

1

´

SP

ij

´

SP

ij

þ

γ

2

´

CL

i

´

SP

ij

´

SP

ij

þ

ϵ

ij

ð

3

Þ

where the probability of yield losses was predicted at

i

th climate condition (CL) in

j

th soil attribute (SP),

μ

is the common intercept,

β

and

γ

are slope parameters, and

ϵ

is the error term. Using a stepwise procedure, non-signi

fi

cant terms (

P

< 0.05)

were dropped. To obtain the parameter estimates, signi

fi

cant terms in the

fi

nal

models were reduced to only one term per component (e.g., intercept, and linear
and quadratic slope) and common intercept was removed using

“

point

”

option

embedded in the

“

DANDA.SAS

”

macro

64

.

Data availability

The datasets relevant to the current study are available at

https://zenodo.org/record/

4959654#.YMkAv6hKhPY

.

Code availability

The R code for statistical analyses including the coding for MARS procedure and
clustering analysis are included in the Supplementary Material of this article.

Received: 9 March 2021; Accepted: 25 June 2021;

References

1.

Livelihoods, To Build Resilient the Impact of Natural Hazards and Disasters
on Agriculture. Food Agric. Organ. UN (2015).

2.

Smith, A., et al. US Billion-Dollar Weather & Climate Disasters 1980-2019.
NOAA National Centers for Environmental Information: Asheville, NC, USA
15 (2019).

3.

Stocker, Thomas F., et al. Climate Change 2013: The physical science basis.
contribution of working group I to the

fi

fth assessment report of IPCC the

intergovernmental panel on climate change. (2014).

4.

Biasutti, M. & Seager, R. Projected changes in US rainfall erosivity.

Hydrology

Earth Syst. Sci.

19

, 2945

–

2961 (2015).

5.

Hoomehr, S., Schwartz, J. S. & Yoder, D. C. Potential changes in rainfall
erosivity under GCM climate change scenarios for the southern Appalachian
region, USA.

Catena

136

, 141

–

151 (2016).

6.

Janowiak, M. et al. Adaptation resources for agriculture: responding to climate
variability and change in the midwest and northeast.

USDA

(2018).

7.

Cooke, B. & Jiang, H. Outlook for US agricultural trade. Electronic Outlook
Report from the Economic Research Service (2020).

8.

Alizadeh, M. R. et al. A century of observations reveals increasing likelihood of
continental-scale compound dry-hot extremes.

Sci. Adv.

6

, eaaz4571 (2020).

9.

Gipson, J. & Joham, H. In

fl

uence of night temperature on growth and

development of cotton (

Gossypium hirsutum

L.). I. Fruiting and Boll

Development 1.

Agron. J.

60

, 292

–

295 (1968).

10. Ludwig, M., Wilmes, P. & Schrader, S. Measuring soil sustainability via soil

resilience.

Sci. Total Environ.

626

, 1484

–

1493 (2018).

11. Lal, R. Rights-of-Soil.

J. Soil Water Conserv.

74

, 81A

–

86A (2019).

12. Hillel, D.

Introduction to Environmental Soil Physics

. (Elsevier, 2003).

13. Yoder, D. C. et al. Soil health: meaning, measurement, and value through a

critical zone lens.

J. Soil Water Conserv.

76

, 103

–

108 (2021).

14. Li, J. et al. Asymmetric responses of soil heterotrophic respiration to rising and

decreasing temperatures.

Soil Biol. Biochem.

106

, 18

–

27 (2017).

15. Nouri, A., Lee, J., Yin, X., Tyler, D. D. & Saxton, A. M. Thirty-four years of no-

tillage and cover crops improve soil quality and increase cotton yield in
Al

fi

sols, Southeastern USA.

Geoderma

337

, 998

–

1008 (2019).

16. Lal, R. Degradation and resilience of soils. Philosophical Transactions of the

Royal Society of London.

Philos. Trans. R. Soc. Lond. Ser. B

352

, 997

–

1010

(1997).

17. Li, Y., Guan, K., Schnitkey, G. D., DeLucia, E. & Peng, B. Excessive rainfall

leads to maize yield loss of a comparable magnitude to extreme drought in the
United States.

Glob. Change Biol.

25

, 2325

–

2337 (2019).

18. McLaughlin, B. C. et al. Weather underground: Subsurface hydrologic

processes mediate tree vulnerability to extreme climatic drought.

Glob. Change

Biol.

26

, 3091

–

3107 (2020).

19. Van Kessel, C. et al. Climate, duration, and N placement determine N2O

emissions in reduced tillage systems: a meta

‐

analysis.

Glob. Change Biol.

19

,

33

–

44 (2013).

20. Döring, T. F. et al. Resilience as a universal criterion of health.

J. Sci. Food

Agric.

95

, 455

–

465 (2015).

21. Corstanje, R., Deeks, L., Whitmore, A., Gregory, A. & Ritz, K. Probing the

basis of soil resilience.

Soil Use Manag.

31

, 72

–

81 (2015).

22. Herrick, J. E. & Wander, M. M. Relationships between soil organic carbon and

soil quality in cropped and rangeland soils: the importance of distribution,
composition, and soil biological activity. in

Soil Processes and the Carbon

Cycle

. 405

–

425 (CRC Press, 2018).

23. Gaudin, A. C. et al. Increasing crop diversity mitigates weather variations and

improves yield stability.

PloS ONE

10

, e0113261 (2015).

24. Biggs, R. et al. Toward principles for enhancing the resilience of ecosystem

services.

Annu. Rev. Environ. Resour.

37

, 421

–

448 (2012).

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv

11

25. Holt-Giménez, E. Measuring farmers

’

agroecological resistance after

Hurricane Mitch in Nicaragua: a case study in participatory, sustainable land
management impact monitoring.

Agric. Ecosyst. Environ.

93

, 87

–

105 (2002).

26. Slette, I. J. et al. How ecologists de

fi

ne drought, and why we should do better.

Glob. Change Biol.

25

, 3193

–

3200 (2019).

27. Tirivarombo, S., Osupile, D. & Eliasson, P. Drought monitoring and analysis:

standardised precipitation evapotranspiration index (SPEI) and standardised
precipitation index (SPI).

Phys. Chem. Earth Parts A/B/C

106

, 1

–

10 (2018).

28. Svoboda, M. et al. The drought monitor.

Bull. Am. Meteorol. Soc.

83

,

1181

–

1190 (2002).

29. Gitz, V., Meybeck, A., Lipper, L., Young, C. D. & Braatz, S. Climate change

and food security: risks and responses. Food and Agriculture Organization of
the United Nations (FAO) Report 110 (2016).

30. Walthall, C. L., Anderson, C. J., Baumgard, L. H., Takle, E. & Wright-Morton,

L. Climate change and agriculture in the United States: Effects and adaptation.
(2013).

31. USDA National Agricultural Statistics Service Cropland Data Layer. Published

crop-speci

fi

c data layer. (2019).

32. Hake, K. & Grimes, D. in Physiology of cotton 255-264 (Springer, 2010).
33. Main, C. L. W288 Cotton Production in Tennessee. (2012).
34. Lal, R. & Shukla, M. K. Principles of soil physics. (CRC Press, 2004).
35. Shen, Y., McLaughlin, N., Zhang, X., Xu, M. & Liang, A. Effect of tillage and

crop residue on soil temperature following planting for a Black soil in
Northeast China.

Sci. Rep.

8

, 1

–

9 (2018).

36. Seneviratne, S. I. et al. Land radiative management as contributor to regional-

scale climate adaptation and mitigation.

Nat. Geosci.

11

, 88

–

96 (2018).

37. Nouri, A. et al. Crop species in no-tillage summer crop rotations affect soil

quality and yield in an Al

fi

sol.

Geoderma

345

, 51

–

62 (2019).

38. Singh, S. et al. Soil organic carbon and aggregation in response to thirty-nine

years of tillage management in the southeastern US.

Soil Tillage Res.

197

,

104523 (2020).

39. Werner, W., Sanderman, J. & Melillo, J. Decreased soil organic matter in a

long

‐

term soil warming experiment lowers soil water holding capacity and

affects soil thermal and hydrological buffering.

J. Geophys. Res.

125

,

e2019JG005158 (2020).

40. Blanco-Canqui, H., Mikha, M. M., Presley, D. R. & Claassen, M. M. Addition

of cover crops enhances no

‐

till potential for improving soil physical

properties.

Soil Sci. Soc. Am. J.

75

, 1471

–

1482 (2011).

41. Jin, V. L. et al. Management controls the net greenhouse gas outcomes of

growing bioenergy feedstocks on marginally productive croplands.

Sci. Adv.

5

,

eaav9318 (2019).

42. Vincent, L., Zhang, X., Mekis, É., Wan, H. & Bush, E. Changes in Canada

’

s

climate: Trends in indices based on daily temperature and precipitation data.

Atmos.-Ocean

56

, 332

–

349 (2018).

43. Houghton, J. T. et al. Climate change 2001: the scienti

fi

c basis. (The Press

Syndicate of the University of Cambridge, 2001).

44. Liu, Q. et al. Extension of the growing season increases vegetation exposure to

frost.

Nat. Commun.

9

, 1

–

8 (2018).

45. Mueller, N. D. et al. Cooling of US Midwest summer temperature extremes

from cropland intensi

fi

cation.

Nat. Clim. Change

6

, 317

–

322 (2016).

46. Lal, R. Restoring soil quality to mitigate soil degradation.

Sustainability

7

,

5875

–

5895 (2015).

47. Qin, S. et al. Temperature sensitivity of SOM decomposition governed by

aggregate protection and microbial communities.

Sci. Adv.

5

, eaau1218 (2019).

48. Blanco-Canqui, H. & Ruis, S. J. No-tillage and soil physical environment.

Geoderma

326

, 164

–

200 (2018).

49. Himes, F. Nitrogen, sulfur, and phosphorus and the sequestering of carbon. in

Soil processes and the carbon cycle

. 315

–

319 (CRC Press, 2018).

50. Li, L. & Schaeffer, S. M. Stabilization mechanisms of isotope-labeled carbon

substrates in soil under moisture pulses and conservation agricultural
management.

Geoderma

380

, 114677 (2020).

51. Nouri, A. et al. Management duration controls the synergistic effect of tillage,

cover crop, and nitrogen rate on cotton yield and yield stability.

Agric. Ecosyst.

Environ.

301

, 107007 (2020).

52. Adams, J. Climate indices, an open source Python library providing reference

implementations of commonly used climate indices. Climate indices in
Python. (2017).

53. McKee, T.B., Doesken, N.J. & Kleist, J. The Relation of Drought Frequency

and Duration to Time Scales.

Proceedings of the 8th Conference on Applied

Climatology.

179-184 (1993).

54. Thornthwaite, C. W. An approach toward a rational classi

fi

cation of climate.

Geogr. Rev.

38

, 55

–

94 (1948).

55. Palmer, W. C.

Meteorological Drought

. Vol. 30 (US Department of Commerce,

Weather Bureau, 1965).

56. Wells, N., Goddard, S. & Hayes, M. J. A self-calibrating Palmer drought

severity index.

J. Clim.

17

, 2335

–

2351 (2004).

57. Culman, S. W. et al. Permanganate oxidizable carbon re

fl

ects a processed soil

fraction that is sensitive to management.

Soil Sci. Soc. Am. J.

76

, 494

–

504

(2012).

58. Jones, D. & Willett, V. Experimental evaluation of methods to quantify

dissolved organic nitrogen (DON) and dissolved organic carbon (DOC) in
soil.

Soil Biol. Biochem.

38

, 991

–

999 (2006).

59. Kemper, W. & Rosenau, R. Aggregate stability and size distribution. Methods

of soil analysis: Part 1 physical and mineralogical.

Methods

5

, 425

–

442 (1986).

60. Fuchs, B. Using the Standardized Precipitation Index (SPI) and the

Standardized Precipitation Evapotranspiration Index (SPEI).

in NDM Center

.

Nebraska-Lincoln, Caribbean Drought Workshop

. 1

–

24 (2012).

61. Friedman, J. H. Multivariate adaptive regression splines.

Ann. Stat.

19

, 1

–

67

(1991).

62. Racine, J. S. RStudio: a platform-independent IDE for R and Sweave. (JSTOR,

2012).

63. Institute, S.

SAS 9.4 Output delivery system: User

’

s Guide

. (SAS Institute, 2014).

64. Saxton, A. A macro for converting mean separation output to letter groupings

in Proc Mixed.

Proc. 23rd SAS Users Group International

, 22

–

25 March 1998,

Nashville, 1243

–

1246 (1998).

Acknowledgements

We fondly remember and honor Dr. Donald Tyler, the founder of this long-term
experiment. We are also grateful to all staff in the West Tennessee AgResearch and
Education Center of the University of Tennessee for many years of effort in collecting
yield data. This project was supported by USDA-ARS, Contract No. 58-6066-8-043. Also,
the

fi

nancial support was in part provided by Agriculture and Food Research Initiative

Competitive Award No. 2015-68007-23212 from the USDA National Institute of Food
and Agriculture. All data needed to evaluate the conclusions in the paper are present in
the paper and/or the Appendix.

Author contributions

A.N. and D.C.Y. conceived the idea. S.J., S.C., J.F., B.T., and X.Y. conducted the
experiment and collected data. A.N., D.C.Y., M.R., and A.M.S. conducted the analysis. All
authors including A.N., D.C.Y., M.R., F.R.W., S.J., S.C., J.F., B.T., P.A., J. L., A.M.S., and
X.Y. contributed equally to the interpretation of results and

fi

nal draft of the paper.

Competing interests

The authors declare no competing interests.

Additional information

Supplementary information

The online version contains supplementary material

available at

https://doi.org/10.1038/s43247-021-00223-6

.

Correspondence

and requests for materials should be addressed to A.N.

Peer review information

Communications Earth

&

Environment

thanks Nele Verhulst

and the other, anonymous, reviewers for their contribution to the peer review of this
work. Primary Handling Editors: Leiyi Chen and Clare Davis.

Reprints and permission information

is available at

http://www.nature.com/reprints

Publisher

’

s note

Springer Nature remains neutral with regard to jurisdictional claims in

published maps and institutional af

fi

liations.

Open Access

This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made. The images or other third party
material in this article are included in the article

’

s Creative Commons license, unless

indicated otherwise in a credit line to the material. If material is not included in the
article

’

s Creative Commons license and your intended use is not permitted by statutory

regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this license, visit

http://creativecommons.org/

licenses/by/4.0/

.

© The Author(s) 2021

ARTICLE

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00223-6

12

COMMUNICATIONS EARTH & ENVIRONMENT

|           (2021) 2:155  | https://doi.org/10.1038/s43247-021-00223-6 | www.nature.com/commsenv