Significance

Numerous empirical analyses have revealed how net primary production responds to climate change and variability over time. However, linear approaches predominantly used in previous studies often generate a weak relationship between climate factors and productivity, which greatly limits our understanding of the mechanism associated with this fundamental ecological process. We demonstrated dryland sensitivity to climate change and variability by incorporating nonlinear dynamics. Dryland sensitivity patterns revealed in this study are largely unrecognized and partly counterintuitive; however, the underlying mechanisms are inferable and cannot be fully revealed by linear approaches. Our work highlights the importance of nonlinear, state-dependent sensitivity of productivity to climate change and variability, accurately forecasting potential biosphere feedback to the climate system.

Abstract

Primary productivity response to climatic drivers varies temporally, indicating state-dependent interactions between climate and productivity. Previous studies primarily employed equation-based approaches to clarify this relationship, ignoring the state-dependent nature of ecological dynamics. Here, using 40 y of climate and productivity data from 48 grassland sites across Mongolia, we applied an equation-free, nonlinear time-series analysis to reveal sensitivity patterns of productivity to climate change and variability and clarify underlying mechanisms. We showed that productivity responded positively to annual precipitation in mesic regions but negatively in arid regions, with the opposite pattern observed for annual mean temperature. Furthermore, productivity responded negatively to decreasing annual aridity that integrated precipitation and temperature across Mongolia. Productivity responded negatively to interannual variability in precipitation and aridity in mesic regions but positively in arid regions. Overall, interannual temperature variability enhanced productivity. These response patterns are largely unrecognized; however, two mechanisms are inferable. First, time-delayed climate effects modify annual productivity responses to annual climate conditions. Notably, our results suggest that the sensitivity of annual productivity to increasing annual precipitation and decreasing annual aridity can even be negative when the negative time-delayed effects of annual precipitation and aridity on productivity prevail across time. Second, the proportion of plant species resistant to water and temperature stresses at a site determines the sensitivity of productivity to climate variability. Thus, we highlight the importance of nonlinear, state-dependent sensitivity of productivity to climate change and variability, accurately forecasting potential biosphere feedback to the climate system.
Climate warming, precipitation changes, and increasing aridity are major facets of contemporary and future climate change across terrestrial ecosystems (14), particularly in drylands that cover ~45% of Earth’s land surface and support >2 billion people (5). Ongoing climate change has led to increasing climate variability via increased frequency of extreme climatic events, such as droughts and heat waves (1, 2, 6). Empirical studies have primarily used linear or occasionally nonlinear equation–based approaches to long-term (>10 y) observational data (711) to determine how climate change and variability affect plant primary productivity (hereafter, productivity), which is a key ecosystem process regulating the global carbon cycle (12). In recent years, abrupt changes in plant productivity driven by increasing aridity across space (13) and time (14) in drylands have been more common, highlighting dryland vulnerability to future aridification (2, 3). Some of these studies have further quantified the sensitivity of productivity to climate change and variability at large spatial scales to forecast potential biosphere feedback to the climate system (7, 911). Nevertheless, such quantifications have largely relied on a linear approach wherein the linear regression slope of productivity against climate variables was used to represent sensitivity (811). Furthermore, most precipitation–productivity relationships at any given site are weak and account for a relatively small percentage of variation over time (15). However, because natural systems are often complex and dynamic (16, 17), alternative approaches are needed to determine the relationship between productivity and climate at a site over time. Herein, we used an equation-free, data-driven approach incorporating complex nonlinear dynamics (16, 18) to reveal how productivity responds to variability in precipitation, temperature, and aridity over 40 y at 48 sites across Mongolia. Such spatial and temporal information is widely lacking; hence, it is crucial to improve forecasts of dryland sensitivity to future climate change.
The effects of climatic drivers on productivity fluctuate over time (8, 19, 20); therefore, system behavior is state-dependent (i.e., nonlinear). State dependency is a defining hallmark of complex nonlinear systems, where the relationships among interacting variables vary in time according to different states of the dynamical system (16, 17, 21). For example, increasing precipitation can increase productivity during some years but decrease productivity during other years because of varying system states and depending on the state of a third variable, such as temperature (19, 22). Consequently, “mirage correlations” are common in long time-series data, wherein linear correlations between variables can appear and disappear or even change sign across time (16, 17, 23) (SI Appendix, Fig. S1 AD). This may lead to the previous notion that productivity is less sensitive to climate variations at a site over time, compared with across sites along a mean annual precipitation (MAP) gradient (24, 25). In addition, previous studies often reported a considerably weaker or nonsignificant temporal linear relationship between productivity and temperature, compared with that between productivity and precipitation (10, 26), emphasizing the primary role of precipitation in driving annual changes in productivity. Notably, however, even the absence of a linear correlation (SI Appendix, Fig. S1D) does not imply lack of causal effects of climatic drivers other than precipitation on productivity (16). Such nonlinear dynamic properties of a system cannot be described by predefined linear or nonlinear equations (16, 17), that is, nonlinearity in this case is not defined by asking whether the underlying equations are linear or nonlinear (17). Statistical power to test the effects and relative importance of precipitation and temperature as well as changes in their interannual variabilities on productivity over time is thus limited by previous linear or nonlinear parametric approaches.
A critical step to account for the state dependency of nonlinear dynamic systems is to perform more accurate forecasting of dryland sensitivity to climate change and variability. Empirical evidence has shown that annual productivity responds positively to increasing annual precipitation and decreasing aridity and that the magnitude in such responses increases in more arid regions (7, 8, 10). Yet, increasing annual precipitation and decreasing aridity do not necessarily enhance annual productivity according to varying system states. In drylands, the state dependency of the system can also arise from the time-delayed effects of annual variation in precipitation and aridity on primary production, which are referred to as legacy effects (11, 15, 27). Anomalous dry years can have negative effects on productivity following drought (15, 27), and a subsequent wet year after a dry year may not necessarily increase productivity as expected (SI Appendix, Fig. S1E). When negative legacy effects of precipitation on productivity prevail across time, the sensitivity of annual productivity to annual precipitation can sometimes be negative. The effects of increasing annual mean temperature on annual productivity could be positive or negative, depending on the vegetation response to warming-induced changes in water availability and lengthening of the growing season (19, 22, 26, 28). A regional ecosystem model in drylands predicted that an extended growing season because of warming can increase productivity despite increasing aridity (22). In contrast, some studies reported that warming effects on productivity can be negative in grasslands dominated by C3 species, which are generally less resistant to temperature and water stresses than grasslands dominated by C4 species (26, 29). Although empirical evidence is limited (7, 9, 11), the overall effects of interannual climatic variability on mean productivity over a certain period will depend on the relative magnitude of positive and negative ecosystem responses (including immediate and delayed responses) to annual climate conditions (11). Previous studies suggested positive effects of interannual precipitation variability in drylands with <300 mm MAP (9, 11), but this may need to be revisited after incorporating nonlinear dynamics.
We applied empirical dynamic modeling (EDM) (16, 18), an equation-free, data-driven approach to aboveground net primary productivity (ANPP) and climate data from over 40 y (1978 to 2017) at 48 grassland sites (i.e., ~1,920 records of ANPP and climate) that cover extensive gradients in the long-term MAP and mean annual temperature (MAT) across Mongolia (Fig. 1 A and B and SI Appendix, Fig. S2 and Table S1). Such long-term ground-based data enable a direct coupling of climate and productivity and should be suitable for the quantification of nonlinear dynamics in drylands. First, we examined the overall effects and relative importance of annual precipitation, annual mean temperature, summer mean temperature, and annual aridity for the water year (September–August), as well as their interannual variabilities on productivity. Second, we compared an EDM approach with a linear approach in quantifying the sensitivity of productivity to climate change and variability across sites.
Fig. 1.
Geographic patterns of climate and the strength of the causal effects of climate variables on aboveground net primary productivity (ANPP) across 48 meteorological sites in Mongolia. (A) Long-term (1978 to 2017) MAP and (B) MAT across sites (indicated by black circles). Geographic patterns in MAP and MAT were visualized using spatial interpolation via ordinary kriging. (C and D) The cross-map skill ρ (indicating the causal strength) for observed ANPP predicting climate variables with maximum time-series length at each site was shown as joint violin and box plots to present the kernel probability density of the values and the interquartile range and median of the values. (C) The causal strength of annual precipitation, annual temperature, summer temperature, and SPEI on ANPP. (D) The causal strength of interannual precipitation variability, interannual temperature variability, interannual summer temperature variability, and interannual SPEI variability on mean ANPP in 6-y moving windows. Interannual variabilities of precipitation, temperature, summer temperature, and SPEI were calculated as the coefficient of variation of annual precipitation, and SDs of annual mean temperature, summer mean temperature and annual SPEI, respectively. Cross-mapping significance was determined by comparing ρ with the maximum time-series length as well as convergence (a difference between ρ at the maximum and minimum time-series lengths) between original and surrogate time-series data (we used 1,000 surrogate time-series data obtained by randomizing the phases of a Fourier transform of climate variables; refer to Materials and Methods for details). The P-value was estimated for each site as the number of surrogates showing a higher ρ with the maximum time-series length as well as a higher convergence, divided by the total number of surrogates. The metasignificance (meta P-value below each violin plot) was then calculated using a recently proposed method to combine the P-values [harmonic mean P-values (28)]. Different letters above violin plots indicate significant differences in the causal strength among climate variables by paired Wilcoxon tests (P < 0.05).
We hypothesized that negative legacy effects of annual precipitation on productivity would prevail across time under high interannual precipitation variability in the drier southern half of Mongolia (15). We expected that the increasing proportion of annual species that survived the dry period as seeds in drier grasslands (30, 31) would explain the predominance of negative effects of drought legacy. In contrast, positive legacy effects of annual precipitation would be common, and/or the magnitude of positive effects of precipitation would be immediate and relatively large in the wetter northern half of Mongolia (11, 15). This might result in a shift from positive to negative in the sensitivity of productivity to increasing annual precipitation and decreasing annual aridity along climatic gradients (Fig. 1 A and B) in Mongolia. We expected an opposite pattern in the sensitivity of productivity to increasing annual temperature. Indeed, experimental evidence from northern Mongolian grasslands (32, 33) suggested that vegetation in cold environments would generally have low resistance to warming-induced soil water deficiency, and productivity would thus respond negatively to increasing temperature. Given that the proportion of C4 species that have higher resistance to water stress than C3 species increased with aridity in Mongolia (34), the negative effects of increasing temperature might be limited in southern Mongolian grasslands. Productivity sensitivity to interannual climatic variability would also reflect shifts in the dominance of drought resistance (26, 29) and avoidance (30, 31, 35) traits in vegetation along the north–south climatic gradient in Mongolia (34). We expected that the balance between productivity gain caused by positive climate extremes and loss caused by negative climate extremes would change depending on vegetation gradients across Mongolia, which may lead to contrasting patterns in the sensitivity of productivity to interannual climate variability (9, 11). Using the EDM approach, we detected geographic patterns in the sensitivity of productivity to climate across Mongolia and identified potential mechanisms underlying ecosystem sensitivity to year-to-year variability in climate.

Results and Discussion

We used convergent cross-mapping (CCM) (16), an EDM method for detecting causality in nonlinear dynamic systems, to examine the degree to which precipitation, temperature, and aridity, and their variabilities, forced ANPP across the 48 sites. The basic idea of this method is that if a candidate effect, variable Y, can predict the current or previous state of a candidate cause, variable X, then X causally influences Y. The cross-map skill was evaluated using Pearson’s correlation coefficient (ρ) between the cross-map estimates and observations. Note that the values of the selected best θ (the tuning parameter indicating nonlinearity of the dynamical system) in the S-maps were >0 in most cases of our CCMs (SI Appendix, Table S3), indicating that grassland productivity generally displayed nonlinear dynamics across the study area (17, 18).
Annual precipitation, annual mean temperature, summer mean temperature, and annual aridity [quantified as a 12-mo integrated standardized precipitation evapotranspiration index; SPEI (36)] showed significant causal forcing across Mongolia (Fig. 1C and SI Appendix, Table S4). The results were metasignificant for all climate variables, as determined using a recently proposed method to combine the P-values [harmonic mean P-values (37)]. The strength of the causal effect of annual mean temperature on ANPP was less, on average, than that of annual precipitation, summer mean temperature, and SPEI (Fig. 1C). There were no significant differences between the strengths of the causal effects of annual precipitation, summer mean temperature, and SPEI on ANPP. Previous studies reported a considerably weaker effect of temperature on productivity than precipitation in drylands (10, 24), probably because of the state-dependent behavior of drylands across time. However, our analysis incorporated the state dependency highlights the importance of summer mean temperature as well as that of annual precipitation and aridity in driving annual ANPP. At most sites, climate variables affected ANPP with a time delay of 1 to 3 y (SI Appendix, Table S4).
All interannual climate variabilities causally influenced mean ANPP in a 6-y moving window across Mongolia (Fig. 1D). Interannual variabilities in precipitation, temperature, summer temperature, and SPEI were calculated as the coefficient of variation of annual precipitation, and SDs of annual mean temperature, summer mean temperature, and annual SPEI, respectively. The results were metasignificant based on the harmonic mean P-values. In this case, the strength of the causal effects of interannual variability in temperature and summer temperature on ANPP was stronger, on average, than that of interannual variations in precipitation and SPEI on ANPP (Fig. 1D), which were not significantly different. At most sites, climate variability affected mean ANPP with a delay of 1 to 3 moving windows (SI Appendix, Table S5). These results were robust because they were qualitatively similar to the results observed with 5- or 4-y moving window (SI Appendix, Fig. S5 and Tables S6 and S7).
The considerable, but until now underappreciated, importance of temperature and its interannual variability in driving dryland productivity are particularly notable given that the future increase in aridity will be driven by a steady increase in global surface temperature in response to elevated atmospheric CO2 levels (2, 3, 38). Our observational records indicated increasing aridity over time because of continued warming (SI Appendix, Fig. S3 B, C, and E). We also observed increases in interannual variabilities in temperature and aridity (SI Appendix, Fig. S3 G and H). Plant growth in cold climates across Mongolia (SI Appendix, Table S1) can be stimulated in warmer years by an extended growing season (22, 39, 40) (as suggested by increases in summer mean temperature and growing degree days over time; SI Appendix, Fig. S3 C and D) and by ameliorated water availability in wetter years (41). Taken together, our study provided general evidence that grassland productivity was significantly driven by both precipitation and temperature at annual and interannual time scales.
Next, we assessed the sensitivity of ANPP to climate and its variability at each site using scenario exploration analysis with multivariate EDM (17, 42, 43). For each historical time point (i.e., a year or moving window), we predicted hypothetical changes in ANPP (ΔANPP) with small perturbations (ΔZ) in climate variables or their variabilities. A higher positive ΔANPP/ΔZ value suggested a more sensitive positive causal effect of climatic drivers and vice versa. The calculated values at each historical time point were averaged across the time series to determine the system-level sensitivity of ANPP over the observation period.
We detected a clear geographic pattern in the predicted map using kriging (Fig. 2A), showing a positive effect of annual precipitation on ANPP in the northern half and negative effect in the southern half of Mongolia. The sensitivity of ANPP to changes in annual precipitation was significantly positively related to long-term MAP and negatively related to long-term MAT at each site (SI Appendix, Fig. S6 A and B). The negative effect of precipitation on productivity could be attributed to the time-delayed effects of precipitation on productivity (11, 15, 27), which are negative after a dry year or positive after a wet year (11, 15). Indeed, we generally detected time-delayed effects ranging from 1 to 3 y across the study area (SI Appendix, Table S4). Negative time-delayed effects are more intense and dominant over the observation period in the highly arid southern half of Mongolia. We observed an increase in the proportion of annual species that are capable of surviving drought as seeds (30, 31) as long-term MAP decreased and the MAT increased across Mongolia (SI Appendix, Fig. S4 E and F). This explains the predominance of negative time-delayed effects of annual precipitation in the southern half of Mongolia. Conversely, positive time-delayed effects are more common, and/or the magnitude of positive effects is more immediate and greater in the northern half of Mongolia.
Fig. 2.
Maps of the sensitivity of aboveground net primary productivity (ANPP) to changes in annual precipitation, annual mean temperature, summer mean temperature, and annual SPEI. The system-level sensitivity of ANPP was interpolated among 48 sites via ordinary kriging. (AD) Maps of the sensitivity of ANPP to each climate variable evaluated using scenario exploration analysis (i.e., a nonlinear approach). (EH) Maps of the sensitivity of ANPP to each climate variable evaluated using a GLS regression (i.e., a linear approach).
In contrast, ANPP responded negatively to annual mean temperature in the north-eastern part and positively in the rest of Mongolia (Fig. 2B). Although the pattern in the sensitivity of ANPP to summer mean temperature slightly differed from that to annual mean temperature (Fig. 2C), we observed a positive relationship between the sensitivity of ANPP to summer mean temperature and long-term MAT at each site (SI Appendix, Fig. S6F). Both results highlight the contrasting responses of ANPP to temperature between warmer and colder regions in Mongolia. Plant growth in cold climates across Mongolia (SI Appendix, Table S1) can be stimulated by lengthening the growing season in warmer years (22, 39). However, the total response of productivity to changes in temperature depends on the degree to which the extended growing season interacts with warming-induced stress (22, 26, 29). In north-eastern Mongolia, vegetation exposed to colder temperatures and increased precipitation may generally have low resistance to soil water deficiency caused by warming-induced evapotranspiration. This creates drought stress in warmer years (29, 32, 33), resulting in the negative response of productivity to increasing temperature. In contrast, the arid and warmer grasslands in Mongolia are more frequently exposed to drought stress and are dominated by vegetation that is resistant to such harsh climate conditions (26, 29). Therefore, we propose that the negative effects of increasing temperature are limited in these arid and warmer grasslands. We observed that the proportion of C4 species that generally have higher resistance to warming-induced drought stress (29, 44, 45) increased along the climatic gradients in Mongolia (SI Appendix, Fig. S4 A and B), which was consistent with a previous report (34) that supports these arguments. Currently, there are large uncertainties on the effects of increasing temperature on productivity, which can be attributed to limitations associated with the number of experiments (41, 46), long-term observations (47), and models (22). Our study provides long-term, large-scale empirical assessments of productivity sensitivity to temperature changes in drylands.
Generally, an increase in SPEI (i.e., decreasing aridity), an index that integrates precipitation and temperature, negatively affected ANPP across Mongolia (Fig. 2D). Although the overall negative effect of SPEI on productivity is counterintuitive, given the well-known positive relationship between decreased aridity and productivity in drylands (8, 10, 22, 23), we propose that the vegetation across Mongolian grasslands is governed by negative time-delayed effects (i.e., drought legacy) (11, 15) of aridity. In addition, the magnitude of negative effects of SPEI on productivity varied spatially, suggesting that the strength of time-delayed effects differed with areas. We thus suggest that the amelioration of soil water deficiency resulting from the balance between precipitation and temperature in a given year does not necessarily enhance annual productivity.
The sensitivity of ANPP to interannual precipitation variability exhibited a clear geographic pattern, going from negative in the north-eastern part to positive in the other regions of Mongolia (Fig. 3A). We observed a slightly decreasing trend in the sensitivity of ANPP to interannual SPEI variability along the site precipitation gradient (SI Appendix, Fig. S7G). Sensitivity of ANPP to interannual SPEI variability demonstrated a clear geographic pattern similar to that of interannual precipitation variability (Fig. 3D). Previous studies suggested that the response of productivity to interannual precipitation variability switches from negative to positive at a long-term MAP of ~300 mm, with more positive responses at drier sites (9, 11). Even below this threshold (most of our sites received < 300 mm; SI Appendix, Table S1), ecosystem sensitivity to interannual variability in precipitation and aridity exhibited a contrasting geographic pattern reflecting shifts in dominance of drought resistance (26, 29, 34) and avoidance (30, 31) traits in vegetation along the north–south climatic gradient in Mongolia (SI Appendix, Fig. S4). The balance between productivity gain caused by positive extremes and loss caused by negative extremes would change depending on such vegetation shifts, resulting in contrasting patterns of productivity responses to interannual variability in precipitation and aridity (9, 11).
Fig. 3.
Maps of the sensitivity of mean aboveground net primary productivity (ANPP) to changes in climate variability (interannual precipitation variability, interannual temperature variability, interannual summer temperature variability, and interannual SPEI variability) in 6-y moving windows. Interannual variabilities of precipitation, temperature, summer temperature, and SPEI were calculated as the coefficient of variation of annual precipitation, and SDs of annual mean temperature, summer mean temperature and annual SPEI, respectively. The system-level sensitivity of ANPP was interpolated among 48 sites via ordinary kriging. (AD) Maps of the sensitivity of ANPP to climate variabilities evaluated using scenario exploration analysis (i.e., a nonlinear approach). (EH) Maps of the sensitivity of ANPP to climate variabilities evaluated using a GLS regression (i.e., a linear approach).
We observed an increasing trend in the sensitivity of ANPP to interannual temperature variability along the site precipitation gradient (SI Appendix, Fig. S7 C and E; however, the trend in the sensitivity of ANPP to interannual summer temperature variability against long-term MAP was marginal, P = 0.067). The predicted map indicated more positive sensitivity of ANPP to interannual temperature and summer temperature variability in northern than in southern Mongolia (Fig. 3 B and C). Across Mongolia, with increasing interannual temperature variability, the magnitude of the positive effects in warmer and wetter years caused by a longer growing season and improved water availability (3941) was larger than the negative effects in warmer and drier years that constrained plant physiological activity (41). Although the degree to which the sensitivity of ANPP to climate variability differed slightly, depending on the size of the moving window (SI Appendix, Figs. S8 and S9), the geographical pattern of the sensitivity of ANPP to interannual variability in both temperature and SPEI was relatively robust to changes in the size of the moving window. Nonetheless, uncertainties still remain with regard to ecosystem sensitivity to climate variability in very dry grasslands, such as those in Mongolia. Longer time series that include a greater number of extreme climatic events are needed to generate data-driven confirmation of the effects of climate variability on primary production.
Further, we showed large discrepancies (Figs. 2 and 3) between predictions of dryland sensitivity to climatic change and variability based on our approach (i.e., scenario exploration), compared with those based on traditional regression analyses (8, 10, 11). Notably, when we used a linear approach, the observed geographic patterns of productivity sensitivity to annual climate variations are largely in line with well-documented patterns (7, 8, 10). Although the CCM for the causal effect of annual precipitation on ANPP did not perform better than a simple linear model, our CCMs generally performed better than simple linear models based on prediction accuracy in most cases (SI Appendix, Tables S4‒S7). However, because ANPP essentially displayed nonlinear dynamics across Mongolia (SI Appendix, Table S3), a linear approach is in principle ill-posed, and a significant linear correlation does not imply causation. Indeed, the detection of significant causal effects of climate variables on ANPP did not necessarily match with the significance detected by the linear models (SI Appendix, Tables S4‒S7). This highlights the need to fully capture the nonlinear, state-dependent sensitivity of productivity to climate change and variability across time (21, 43, 48) when quantifying system-level sensitivity. We thus demonstrated dryland sensitivity to climate change (Fig. 2 AD) and variability (Fig. 3 AD) using an EDM framework (16, 18) that acknowledges the complex nonlinear dynamics of ecological systems. Generally, dryland productivity often displays state-dependent behavior, as presented in SI Appendix, Table S3, thus making inferences regarding dryland sensitivity to climate change and variability from linear models, in particular, more challenging.
Using an equation-free, nonlinear time-series analyses (16) and the most extensive (48 sites) field-based data of climate and productivity collected over the longest period (over 40 y) to date, our study provided findings on dryland sensitivity to climate change and variability. Although dryland sensitivity patterns that we found are largely unrecognized and partly counterintuitive, at least two underlying mechanisms are inferable. First, a time-delayed climate effect modifies the responses of annual productivity to annual climate conditions. Second, the proportion of plant species resistant to water and temperature stresses at a given site determines productivity sensitivity to variation in climate. Our results also highlighted that the effects of climate on vegetation should be evaluated across multiple time scales (i.e., annual and interannual time scales) by incorporating the potential occurrence of time-delayed effects. To fully understand the sensitivity of global drylands to climate change and variability using an EDM approach, there is a need for prolonged long time-series data (16, 17) spanning across the globe; therefore, more extensive analysis is warranted in future studies. In particular, it is necessary to clarify whether the state-dependent, nonlinear responses of productivity to climate change and variability are specific to cold and water-limited drylands, such as Mongolia. In addition, global-scale manipulative experiments incorporating climatic extremes (4951) and model simulations (9, 22) are required to elucidate the fundamental mechanisms underlying the observed patterns in dryland sensitivity to climate. Another consideration is the need for an explicit statistical test on the synergistic causal effects of temperature and precipitation on productivity under further methodological development, as the CCM permits the identification of causal relationships between two time-series variables (16, 17, 52). Because grazing is the dominant land use in drylands and can interact with climate to impact ecosystem functioning (53), future studies should also evaluate how grazing mediates dryland responses to climate change and variability. Given the free-range livestock management strategy in Mongolia, this approach might be the key to effectively address and adapt to challenges associated with climate change and variability; however, it requires a holistic understanding of the impacts of grazing and climate change on drylands (54, 55). Nonetheless, our spatially explicit regional assessment and visualization of how productivity responds to climate change and variability can be used for validation as well as for important improvements incorporating nonlinear vegetation dynamics in global ecosystem models. Finally, our findings emphasize the importance of considering nonlinear dynamics of drylands to accurately forecast potential biosphere feedback to the climate system.

Materials and Methods

Datasets.

The Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) of Mongolia monitored climate variables (daily mean air temperature and precipitation) and ANPP at 48 meteorological observation sites. These sites are widely distributed across Mongolian grasslands (spanning 2,160 km in linear distance from the western to eastern Mongolia and 1,050 km from the northern to southern Mongolia; SI Appendix, Fig. S2 and Table S1) in 1978 to 2017 (note: at some sites, the time series started after 1979). ANPP was calculated as the mean of aboveground plant biomass harvested from four “1 × 1 m” quadrats randomly placed within each station at the beginning, middle, and end of August, within the peak season of plant growth (56, 57). At each harvest, the locations of the sampling quadrats were moved within each site to avoid the effects of harvesting on the ANPP measures. In addition, we harvested only the current-year shoot of the shrub species to ensure that harvested plant materials can approximate ANPP. Although the number of samples at one time was relatively small, repeated measurements of aboveground plant biomass during the peak growing season would allow the estimation of ANPP at each site. Each station was fenced to avoid damage to observation equipment by livestock. Our time-series analyses require continuous and equidistant data (16, 17); therefore, we imputed missing values for ANPP up to a maximum of five consecutive years using a Kalman smoother algorithm (58). Sites that had more than five consecutive missing years were removed from the analysis, resulting in 48 of 70 candidate meteorological sites. Toward lower latitude sites, the long-term MAP generally decreased whereas the MAT increased (Fig. 1 A and B and SI Appendix, Table S1). Mongolia is characterized by a dryland climate, and ~70% of annual precipitation falls between June and August. In addition, we observed increasing growing degree days with rising annual mean temperature and summer mean temperature over time (SI Appendix, Fig. S3). As such, our climate measures are conducted on an annual basis, including extreme climatic events in summer, such as drought and extreme precipitation over a long period of time. The vegetation type across all sites was grasslands, including forest steppe, steppe, and desert steppe (SI Appendix, Table S1). Across the datasets of all the 48 sites, 92.7% of ANPP values (ranging from 77.5 to 100%) remained fixed (not imputed) over the time series (SI Appendix, Table S1). There were no missing values for the climate variables.

Data Processing.

At each meteorological site, annual precipitation was calculated as the sum of daily precipitation across a year. Annual mean temperature and summer mean temperature were calculated as the average of mean daily air temperature across a year and the average of mean daily air temperature between June and August (i.e., a growing season), respectively. Here, a year was defined from the beginning of September in the preceding year to the end of August in the focal year. We focused on the role of annual mean temperature as well as summer mean temperature in driving ANPP in the following analyses. Previous studies have suggested that both annual and summer temperatures can affect plant productivity, especially in cold regions (26, 33, 59). We also calculated growing degree days (°) as the sum of mean daily temperature above 10 ° in a year.
We calculated the aridity index, known as SPEI, which integrates both precipitation and temperature (36). To achieve this, we used monthly potential evapotranspiration and precipitation data and calculated monthly SPEI with a 12-mo integration period, as previously suggested (10). Negative SPEI values indicated more arid conditions. Next, we calculated the annual SPEI as the mean of the monthly SPEI from the beginning of September in the preceding year to the end of August in the focal year.
To quantify climate variabilities, we aggregated yearly data in 6-y moving windows and calculated the coefficient of variation of annual precipitation and SDs of annual mean temperature, summer mean temperature, and annual SPEI. We used the SDs of annual mean temperature, summer mean temperature, and annual SPEI to avoid interference caused by divisions by their mean values closer to zero. The mean ANPP for each window was also calculated. This approach (11) allowed us to directly test the causal effects of climate variability on plant primary productivity. We decided to use overlapping window analyses because nonoverlapping windows reduce the time-series length and may undermine the efficacy of time-series analysis (16, 17). The potential effects of the overlapping nature of the data on our results were addressed by developing a null test with a surrogate time-series (refer to the next section). To examine the robustness of our results to different time windows, we repeated the time-series analysis using data with 4- and 5-y moving windows.

Long-Term Trends in Regional Climate.

Long-term trends in climate variables and their variabilities across Mongolia were visualized using a generalized least-squares (GLS) regression with the site as a random effect and the first-order autoregressive processes (AR1). Annual precipitation did not change significantly across the 48 sites from 1978 to 2017 (SI Appendix, Fig. S3A). Over the same period, annual mean temperature and summer mean temperature increased (SI Appendix, Fig. S3 B and C). Growing degree days also increased steadily over time (SI Appendix, Fig. S3D). Because of rising temperature despite no changes in precipitation amount, annual SPEI has declined, indicating that aridity has increased over time (SI Appendix, Fig. S3E). Furthermore, interannual precipitation variability in 6-y moving windows showed no directional change (SI Appendix, Fig. S3F). Interannual variabilities of annual mean temperature and SPEI have increased significantly over time (SI Appendix, Fig. S3 G and H).

Vegetation Shifts along Climatic Gradients.

Vegetation data are not available for the same time-series of ANPP and climate because vegetation surveys were not conducted simultaneously with ANPP measurements. Instead, we used recent observational data of vegetation between 2012 and 2019 at or near the meteorological sites (SI Appendix, Table S2). The sampling locations near the meteorological sites were located on flat areas or gentle slopes where no livestock dung or plants damaged by grazing was evident (i.e., a similar condition to the fenced meteorological sites). At each location, plant species composition was sampled annually at the beginning of August between 2012 and 2019 (note: observation years slightly varied among the locations). Plant species composition was determined along two 50-m permanent lines. Along each line, a 50-cm pole was placed at 25-cm intervals (i.e., 201 points in total), and all plant individuals whose leaves and/or stems touched the pole were identified, and the total number of touches were counted. The abundance of each plant species was then determined by dividing the total number of touches for each species by the total number of points along a line. Data on plant species composition along the two lines were pooled across observation years. To confirm the potential vegetation shifts along climatic gradients across Mongolia (34), we used a generalized linear model with a quasibinomial error structure and a logit link function to analyze the changes in the relative abundance of C3/C4 species, annual/perennial species, and grass/forb/shrub species (defined based on the existing literature) (34, 60) according to long-term MAP and MAT at each site (SI Appendix, Fig. S4). In accordance with the previous report across Mongolia (34), the relative abundance of C4 (SI Appendix, Fig. S4 A and B) and annual species (SI Appendix, Fig. S4 E and F) increased with decreasing precipitation and increasing temperature.

Data Analysis.

To determine the potential causal effects of climate variables and their variabilities on ANPP and mean ANPP in moving windows, we used the EDM method to detect causality [i.e., CCM (16)]. The basic concept of CCM is to use the prediction between variables as a test for causality. If variable X (e.g., precipitation) had a causal effect on variable Y (e.g., productivity), then the causal information of variable X should be present in Y. The attractor recovered for variable Y should be able to predict the state of variable X. The prediction skill of cross-mapping (cross-map skill) was measured using Pearson’s correlation coefficient (ρ) between the predicted and observed X values. This procedure was repeated using a subset of the time series X with different lengths. When the cross-map skill from variable Y to X improves with the time-series length (i.e., convergence), variable X has a causal effect on Y (16). Theoretical details of the CCM algorithm can be found in previous studies (16, 17).
Prior to CCM analysis, all time-series data were standardized to have zero mean and unit variance. Throughout the analysis, the prediction skill was evaluated based on Pearson’s correlation coefficient (ρ) using leave-one-out cross-validation. First, we used the S-map (sequential locally weighted global linear maps) (18) to determine the appropriate embedding dimension (E) and to check the nonlinearity (indicated by a nonlinear localization parameter θ) of the time series of ANPP and mean ANPP in moving windows. According to Takens’ Theorem, the true dynamics of the system can be reconstructed by a set of time-lagged coordinates of a single time series (16, 17). The embedding dimension is the number of time-lagged coordinates used for reconstructing the state space (16, 17). In the majority of real-world scenarios, the embedding dimension is not known beforehand and requires estimation. The appropriate embedding dimension (E), that is the embedding dimension showing the best performance of EDM predictions, depends on several factors, including system complexity, time-series length, and noise (16). The S-map is a model-free method that predicts the near future using local linear regression under a certain embedding dimension, and is a tool for evaluating whether a system displays linear or nonlinear dynamics over time (18). S-map forecasts the trajectory of a target system state by a locally weighted linear regression using the all data points in the state space (17, 18).
We evaluated the predictability of the S-map for all possible combinations of E = 2–8 (with an increment of 1) and θ = 0–10 (with an increment of 0.1). We selected the optimal values of E and θ, under which the best prediction skill (ρ) was obtained. The nonlinear localization parameter θ determines the degree to which the points are weighted when fitting a local linear map. When θ = 0, all points are equally weighted such that the local linear map is identical for different points in the reconstructed state space (17, 18). In this case, the S-map is identical to the global linear map. When θ > 0, nearby points in the state space receive a larger weight, and the local linear map can vary in the state space to accommodate state-dependent, nonlinear behavior (17, 18). The weighting function in the S-map is defined as the form of an exponential decay kernel (18), w(d) = exp (−θd/dm). Here, d is the Euclidean distance between the predictee and each time point, and dm is the mean Euclidean distance of all paired time points. We confirmed that ANPP and mean ANPP in all (4-, 5-, and 6-y) moving windows generally displayed nonlinear dynamics (SI Appendix, Table S3).
The significance of CCM was judged by comparing the cross-map skill ρ at the maximum time-series length as well as convergence (the difference between ρ at the maximum and minimum time-series lengths) between the original and surrogate time-series data (17, 61, 62). We randomly generated 1,000 surrogate time-series of causal climate variables by randomizing the phases of a Fourier transform, which preserved the power spectra or autocorrelation of the cross-mapping target time series (63). To apply CCM with the original and surrogate time-series of climate variables to predict the time series of productivity, we calculated the cross-map skill ρ at the minimum (optimal embedding dimension E+1) and maximum time lengths for 1,000 surrogate time-series and the original time series. The P-value was then estimated as the number of surrogate time-series data showing a higher ρ with the maximum time-series length as well as a higher convergence than those for the original time-series data, divided by the total number of surrogate data (61). Here, we added 1 to the numerator and denominator to correct for finite sampling. In addition, because climate would exhibit time-delayed effects on productivity (often termed legacy effects) (11, 15, 27, 64), we considered time lags (time to prediction tp of a potential causal time-series) in cross-mapping, i.e., lagged CCM (52). We examined the time-delayed effect from 0 to −3 time lags (years or moving windows). We report the best result (highest ρ with the maximum time length) among the lagged CCM analyses (SI Appendix, Tables S4‒S7).
To determine the metasignificance of CCM tests across Mongolian grasslands, we applied a recently proposed method to combine P-values and harmonic mean P-values (37). Note that the metasignificance may be relatively conservative compared to other methods to combine P-values because this method controls for family-wise error rate and, importantly, is robust to interdependent samples of P-values. In addition, we compared the cross-map skill with the maximum time length among the causal climate variables using paired Wilcoxon tests.
We further tested the directionality of the responses of plant productivity to climate variables using scenario exploration with multivariate EDM (17, 42, 43, 65). In nonlinear systems, drivers generally have an effect that is state-dependent; the strength and direction of the effect depend on the current state of the system (17, 21, 42, 65).
For each historical time point t (a year or a moving window), we used S-maps (18) to predict ANPP at time t + 1 with a small increase (+ΔZ/2) and a small decrease (−ΔZ/2) in the climate variables at time t, Z(t). We used ΔANPP/ΔZ to approximate the effects of climate variables at time t. A higher positive ΔANPP/ΔZ value suggested a more sensitive positive causal effect of the climate variables, and vice versa. We then determined the average of ΔANPP/ΔZ at each historical time point across the time series to represent the system-level sensitivity of ANPP to climate at each site, as proposed in a previous study (48). We used 50% of the SD of the observed climate variables, Z(t), as ΔZ. This value corresponded to ΔZ = 15.6–48.4 mm (with a mean of 27.7 mm) and 0.37 to 0.73 ° (with a mean of 0.51 °) across sites for annual precipitation and annual mean temperature, respectively. These values were generally comparable to the prediction range of IPCC AR6 (66) that MAP and MAT between 2080 and 2100 increase by 10 to 30% and 1.5 to 5 °, respectively, over middle latitudes compared to MAP and MAT between 1850 and 1900. The magnitude of change, ΔZ, in scenario exploration can thus be determined depending on the research contexts (17).
Because we treated time-series across all sites equivalently, we a priori selected parameters in the S-map procedure that worked across all the time series, as previously suggested (43): E = 3 and θ = 4 to predict ANPP changes against small perturbations of climate in a given year, and E = 3 and θ = 0.3, to predict mean ANPP changes against small perturbations of climate variabilities at a given moving window.
To compare the results based on scenario exploration with those using a linear-based approach, we used a GLS regression to quantify the responses of plant productivity to climate variables. We estimated the sensitivity of productivity to climate variations as the linear regression slope of productivity against climate variables to be comparable with the methods used in previous studies (8, 10, 11, 24). Considering the autocorrelations among observations across time, we included a term that corrected for first-order AR1 in the models. We used the “gls” function in the “nlme” package in R to perform these analyses. In addition, we compared the prediction accuracy based on CCM and GLS, and we examined the difference between terminal ρ from CCM and the prediction accuracy (a correlation coefficient between observed and predicted values) from GLS using paired Wilcoxon test.
Finally, we visualized the geographic patterns in the sensitivity of ANPP to climate and its variability across Mongolia by interpolating the sensitivity values across sites. All spatial interpolations were performed with the “autoKrige” function in the automap package in R. All other data analyses were performed with R software (67) using the “rEDM,” “harmonicmeanp,” “SPEI,” “nlme,” and “automap” packages.

Data, Materials, and Software Availability

The datasets generated during and/or analyzed during the current study are available through figshare (68).

Acknowledgments

We thank the staff at the IRIMHE for their long-term efforts to maintain the meteorological stations and collect field data. This work was financially supported by: Fostering Joint International Research A (no. 19KK0393) from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The Joint Research Program of Arid Land Research Center, Tottori University to T.S. (no. 30F2002 and 02F2002). U.S. NSF, LTER (DEB # 1655499) and LTREB (DEB-1856383) programs to S.L.C. and J.A.R.

Author contributions

T.S., S.L.C., J.A.R., G.B., E.B., and T.K. designed research; T.S., S.L.C., J.A.R., G.B., E.B., and T.K. performed research; T.S. analyzed data; and T.S., S.L.C., J.A.R., G.B., E.B., and T.K. wrote the paper.

Competing interests

The authors declare no competing interest.

Supporting Information

Appendix 01 (PDF)

References

2
J. Huang et al., Dryland climate change: Recent progress and challenges. Rev. Geophys. 55, 719–778 (2017).
3
J. Huang, H. Yu, X. Guan, G. Wang, R. Guo, Accelerated dryland expansion under climate change. Nat. Clim. Change 6, 166–171 (2016).
1800
1801
1802
1803
1804