Climatic limiting factors of North American ecosystems: a remote-sensing based vulnerability analysis

Remote-sensing based vulnerability assessments to climate change are a research priority of critical importance for landscape-scale efforts to prioritize conservation and management of ecosystems. Limiting climatic factors can serve as a proxy for quantifying ecosystem vulnerability, since theory predicts that ecosystems close to critical climate thresholds will be more sensitive to interannual variation in limiting climate factors. Here, we analyze time series of enhanced vegetation index data for continental-scale vulnerability assessments. The analytical approach is a lagged monthly correlation analysis that accounts for memory effects from the previous growing season. Mapping multivariate correlation coefficients reveals that drought vulnerabilities can be found across the continent, including a distinct geographic band across the western boreal forest. The analytical approach reveals climate dependencies at high spatial and temporal resolution, with the direction and strength of correlation coefficients indicating the risk of threshold transgressions at the edge of species and ecosystem tolerance limits. The approach is further useful for hypothesis testing of contributing non-climatic factors to climatic vulnerability, allowing locally targeted management interventions to address climate change.


Introduction
As a tool for assessing the vulnerability of ecosystems to climate change, the strength of correlations between interannual climate variability and vegetation response can be a useful proxy for quantifying ecosystem vulnerability. Ecosystems or species close to critical climate thresholds will be more sensitive to anomalies (Dakos et al 2012), which is an essential element of edge detection theory (Bathiany et al 2020). Critical climate thresholds can apply relatively uniformly to members of natural species assemblages, allowing for ecosystem-level detection of vulnerabilities to climate change (Trisos et al 2020). While natural resource managers and conservationists are primarily concerned with abrupt ecological disruption due to threshold transgressions at the edge of species and ecosystem tolerance limits (Millar and Stephenson 2015), the importance of gradual transitions has recently been pointed out in a metaanalysis of 4600 empirical climate change impact reports (Hillebrand et al 2020). While threshold transgressions threaten profound changes to ecosystems, they are less frequent than gradual ecosystem and species responses to climate change. Comprehensive vulnerability assessments should therefore detect both response types, threshold transgressions and gradual changes, to infer vulnerability of species or ecosystems.
Remote-sensing based vulnerability assessments to climate change is a research priority of critical importance, as few other methodological approaches could guide conservation and management priorities at continental or global scales. Previous works on developing remote-sensing based metrics to assess vulnerability include a number of valuable methods. Seddon et al (2016) developed a vegetation sensitivity index based on seasonally detrended multivariate climate time series, where coefficients of variables were weighted by their significance in a principal component regression. Another widely applied method estimates resistance (response to climate anomaly) and resilience (speed of recovery) metrics from a normalized difference vegetation index (NDVI) time series (de Keersmaecker et al 2015). Li et al (2018) added interpretative values by evaluating resistance and resilience of vegetation in the context of climatic exposure. Li et al (2020) added a stability component to resistance and resilience metrics to arrive at a more comprehensive vulnerability score. Another comparable approach relies on a linear mixed-effects model to quantify the strength of the association between precipitation/aridity and an estimation of primary productivity based on the enhanced vegetation index (EVI) (Maurer et al 2020). All these methods yield a single metric of vulnerability for decision support. While these indices deliver valuable information, they generally do not allow inferring the nature of the vulnerability because the indices are not time-and variable-specific.
A conceptually different approach to vulnerability assessments is to evaluate long-term historical trends in ecosystem health and productivity. In the context of remote sensing, these trends are referred to as greening and browning. For example, Sulla-Menashe et al (2018) observed greening in the eastern North American boreal forests, which is more humid, and browning in the western boreal, where forests are more prone to moisture stress. Pan et al (2018) documented the risk of reversal from long-term greening to browning in the warmer future due to evapotranspiration demand. Both of these remote sensing based observations have also been confirmed through historical tree ring records suggest observed global warming shifted climatic drivers of tree growth from temperature-limited to moisture-constrained (Babst et al 2019), and showing that western North America has seen declines in growth, while the northeast has been described as a climate change refugium (e.g. D'Orangeville et al 2018). The studies generally point to the importance of a multivariate climate gradient from cold/wet to warm/dry conditions, which is a key criterion for assessing ecosystem vulnerability.
In summary, previous remote sensing indices that were developed have focused on evaluating stability based on short term responses to interannual climate variation, or by tracking trends over time. Both approaches have identified evapotranspiration demand as one major climatic drivers of tree growth in a warming world. Furthermore, while the potential of threshold transgression of species and ecosystems is well grounded in ecological theory, gradual responses (positive and negative) as well as their causes should also be quantified. This justifies the need for different approach to assessing ecosystem sensitivity and vulnerability. Our objectives are therefore to: (a) develop a multivariate vegetation response descriptor that conveys where and when forests are likely becoming vulnerable to climate-induced physiological stress associated with drought and high temperatures; (b) apply this multivariate index to the North American continent, mapping the behavior of growth and the underlying climatic driving factors; and (c) to associate these patterns of vegetation response with regional long-term climatic conditions, as well as non-climatic factors, such as forest stand age, elevation, and topo-edaphic factors, with the purpose of testing hypotheses of interactions among climatic and non-climatic factors. We discuss applications for natural resource management in the context of climate change.

Remote sensing data
As a proxy for annual vegetation productivity, we use a remotely sensed vegetation index. Among available vegetation indices with long time series and high resolution, the EVI was selected because it is less sensitive to noise from soil and atmospheric conditions and is less saturated in high-biomass areas than the NDVI (Huete et al 2002). Also, as a measure of vegetation greenness, EVI has been used to derive metrics related to primary productivity and to reflect vegetation responses to drought (e.g. Anderson et al 2010, Seddon et al 2016, Vicca et al 2016. Here, we use 17 years (January 2003-December 2019) of 16 days 500 m EVI records accessed from Moderate Resolution Imaging Spectroradiometer Vegetation Indices (MOD13A1, Collection 6) as described by Huete et al (2002) and Didan et al (2015). The dataset was obtained through NASA Land Processes Distributed Active Archive Centre, U.S. Geological Survey/Earth Resources Observation and Science Centre (https:// lpdaac.usgs.gov), comprising 391 layers (layer for 1 August 2016 was missing from the database, but imputed for this study as explained below).
The quality of EVI records varied among seasons and regions, usually showing a high percentage of missing values across high latitudes areas and during winter months. Records flagged as poor quality (classified as 'Lowest quality' , 'Quality so low that it is not useful' , 'L1B data faulty' , 'Not useful for any other reason/not processed') were excluded from further analysis. Also excluded from this analysis were human-managed ecosystems and grasslands, according to the vegetation classification according to the MCD12Q1 International Geosphere-Biosphere Programme classification (Sulla-Menashe and Friedl 2018). Thus, the analysis includes needle leaf, broadleaf and mixed forests with tree cover above 60%, as well as shrub lands and savannas with at least 10% tree coverage, as defined in classes 1 through 9 in Sulla-Menashe and Friedl (2018). To be included in the analysis, a pixel required at least ten eligible annual classifications for the study period (2003-2019).
After filtering out low-quality observations, agriculture areas and rangeland land cover, we rescaled the spatial resolution from original 500 m to 2 km by averaging the available measurements per temporal interval. This step resulted in more complete time series for a given 2 km pixel. Incomplete time series with less than 300 (out of 391) observations were removed. The remaining missing values of grid cells were estimated with the seasonally decomposed missing value imputation, implemented with the na_seadec() function of the imputeTS package for the R programming environment (R Core Team 2020). This data cleaning approach resulted in gap filling of 1.5% of all data over the entire study area and 17 year time period. To quantify the annual productivity for correlation analysis, we summed the EVI value above 0 per 2 km pixel per year, for each year from 2003 to 2019.

Climate data
We generated spatially interpolated monthly climate grids from 2002 (a year prior to the start of EVI data coverage) to the end of 2019 using the software ClimateNA v7.01 (Wang et al 2016), which is available at http://climatena.ca. The ClimateNA software extracts monthly historical time series of interpolated climate data for particular locations and variables of interest. The interpolated climate grids are based on the Parameter Regression of Independent Slopes Model interpolation method for weather station data (Daly et al 2008). In this study, we used six monthly variables, including three temperature variables: maximum, minimum and average monthly temperature (T max , T min , T avg ), and three sets of precipitation-related variables: monthly precipitation (PPT), monthly relative humidity and a monthly climate moisture index (CMI) developed by Hogg (1997). In total, the correlation analysis relies on 1224 monthly climate grids (6 variables × 12 months × 17 years) for the North American continent.

Analysis of climatic limiting factors
Plant response to climate has lagged effects that vary with vegetation type, soil conditions, variable type and general geographic and macroclimatic region. One widely used approach to flexibly allow for different delayed effects for different variables is a lagged correlation analysis, where Pearson's correlation coefficients between an annual response variable and multiple monthly predictor variables are evaluated. The temporal window of this analysis is usually specified from the end of the last growing season (e.g. last year's September) to the end of current year's growing season (e.g. current August). We implemented this approach in a per pixel analysis of EVI values across 17 year study period. Note that similar types of analysis exist for dendrochronology data, referred to as response function analysis in this field, using a multivariate regression approach, which also allows for statistical inference from individual tree response to populations (Fritts et al 1971, Zang andBiondi 2015). However, in this case we work directly with the statistical population of EVI values and applications of inferential statistics would not be appropriate.
For the purpose of concise reporting, we clustered EVI grid cells with similar response coefficients using recursive partitioning for 14 groups. The cluster analysis was carried out through multivariate recursive partitioning implemented with the mvpart() function of the mvpart package for the R programming environment (R Core Team 2020). The number of groups was pre-set with the rpart.control() option to obtain 14 clusters that explained over 40% of the total variance in response coefficients. To confirm the robustness of this approach, another divisive clustering technique (k-means) was also used with the same number of clusters, yielding similar results. However, we used recursive partitioning for its speed and ability to handle very large datasets.

Climatic response on a continental scale
To visualize continental patterns of limiting climate factors, we mapped remotely sensed vegetation response summarized into 14 clusters of grid cells that show similar vegetation response to interannual climate anomalies (figure 1). To allow a regional discussion of results, we further impose an arbitrary delineation into regions with similar vegetation response to climate anomalies (figure 1(a) and table 1). This delineation is simply meant to communicate a geographic location on the map and is not used in analysis. For a first visualization of broad continental patterns, cluster means from recursive partitioning were colored based on a scale generated by subtracting temperature from precipitation coefficients. Grid cells limited by high temperature and low precipitation are located on the red end of the color ramp and grid cells limited by cold temperatures and high precipitation values are indicated in blue (figures 1(b) and (c)).
Although other combinations are possible, most clusters (or geographic regions) can be ordered along a gradient from being limited by warm temperature and low precipitation (figure 1(c), upper left) to being primarily limited by cold temperature, with precipitation having smaller additional negative effect (figure 1(c), lower right). Two regions deviate to some extent from this diagonal positioning: the Taiga West (TW) region showed positive correlation of both temperature and precipitation with EVI, and Boreal East (BE) has stronger negative association with moisture condition. An important consideration in reading this map correctly is that response coefficients (colors) only indicate how vegetation responds to a local annual climate anomaly, relative to average Figure 1. Vegetation response to climate, summarized into 14 clusters across North America. General geographic areas with similar cluster composition (a) were used to analyze regional variation in cluster composition (table 1). Clusters are colored from most to least limited by high temperatures and low precipitations (b), and subsequently mapped (c). Gray areas represent no data or excluded ecosystems. climate conditions. The coefficients (colors) do not allow an inference on absolute values of vegetation response, nor do they represent the effect absolute climate values at different locations. To give a simplified, univariate example, a positive (blue) temperature coefficient indicates a relative increase of EVI-inferred productivity at this location in response to a year with higher temperature than average for this location.
Broad geographic trends in limiting climate factors are readily apparent ( figure 1(b)), with expected patters of northern and high elevation ecosystems cold-limited and southern interior regions of the continent most drought limited. However, there is also a remarkable regional diversity of clusters, which is summarized in table 1. Generally, there are few landscapes with high uniformity of vegetation response, such as arctic ecoregions that are dominated by only a few clusters with similar thermal limitations, resulting in the lowest within-region variance in vegetation response across the North America study area (table 1, TW, Hudson Plain (HP), and Pacific North (PN)). In contrast, southern and mountainous regions cover almost all clusters from 1 to 14, resulting higher variances (table 1, Sierra US (SU), Temperate South (TS), and Tropical (T)). Given similar latitude, eastern boreal forest ecosystems have a more homogenous response than the western and interior ones, with about half the variance in vegetation response (table 1, Boreal West (BW), Boreal Southwest (BSW), Boreal Central (BC) versus BE, Temperate North (TN), Taiga East (TE)). This is also visible in the map, with discrete moisture-constrained patches in the central and western portions of the boreal forest apparent ( figure 1(b), red patches).

Monthly response function analysis
Monthly response functions for specific climate variables reveal in greater detail at what time of the year specific variables have a positive or negative impact on vegetation productivity (figure 2). While continentwise responses function clusters are provided as supplement figure S1, we focus here on a number of regional examples to illustrate the analytical approach. For instance, the southern savannas and dry forests (region DM: Dry forest, Mexico) are predominantly drought limited. In this region, only drought-limited clusters 1, 2 and 5 occur with significant frequency (table 1, last row). A specific example is shown in figure 2(a) for cluster #1, which has strong negative correlations with temperature and high positive correlations for moisture variables, especially Table 1. for the CMI, which indicates water deficits between March and August.
At the opposite end of the spectrum, the regions Pacific Central (PC), PN, and taiga ecosystems of the Hudson Bay (TH) are dominated by cluster #14 (table 1, upper right), with year-round limitations by cold temperatures, and with precipitation having slightly negative effects in June and July (figures 2(b) and S1 for cluster #14 across all regions). In the more southern regions influenced by Pacific climate, cluster #14 still occurs at lower frequencies, but at lower elevation, there are strong drought limitations with the most frequent cluster being #5 for the sierras of the southern United States (Region SU, figure 2(c)). In these oak or pine dominated forests, drought is a prevalent limiting factor throughout the year, especially in July.
For eastern temperate forest ecosystems from south to north, TS, Temperate Central (TC), TN, a gradient of limiting factors emerges from the vegetation response to climate (figures 2(d)-(f)). Southern pine forests of the United States have diverse response with the dominant clusters being #1 and #3 (other clusters found in this region are discussed below). These clusters are moisture limited from April to August ( figure 2(d)). This moisture limitation is reduced for TC region, where cluster #3 is dominant (figure 2(e)), followed by temperate mixed forests of the Great Lakes, with further reduced moisture limitations (figure 2(f)). For northeastern regions we see increasing cold limitations and negative effects of moisture variables during the growing season ( figure 2(g)).
Another notable continental pattern is the latitudinal sequence of BSW, to BW to TW, with both the northern and southern adjoining regions less limited by evapotranspirative demand (figure 1(b)). The red patches in the western boreal represent cluster #3, indicating growth limitations driven by warm summer temperatures (figure 2(h)).

Inferring vulnerability to climate change
Both our analysis and other research approaches suggest that drought may play an important role in limiting growth of forested ecosystems in the future for a variety of regions of North America. Climate change impacts for high latitudes of North America, where the warming signal is strongest, have already been relatively well documented. While highlatitude boreal ecosystems are generally expected to benefit from warming trends, widespread greening trends that were initially observed in remote sensing in response to climate change have in many cases reversed to browning trends in the most recent decade (Phoenix and Bjerke 2016). For the North American boreal, browning, or trend reversals from greening to browning are prevalent in western Canada and Alaska (de Jong et al 2011, Ju and Masek 2016, Pan et al 2018. This corresponds to diagonal band (dominated by response cluster #3) across the western boreal forest that we also see in this study, where warm summer temperatures and water deficits in spring and fall are limiting factors (figures 1(b) and 2(h)). Our interpretation is that initial greening observed in the 1980s-1990s reversed to browning because of drought limitations by increasing evapotranspiration demands, as inferred by response cluster #3.
Tree ring studies have also pointed to droughtlimitations in the western boreal of Alaska (Barber et al 2000, Trugman et al 2018 and western Canada (Girardin et al 2016, Hogg et al 2017. In contrast to western boreal forest ecosystems, eastern boreal forest suggests the forest productivity could benefit from up to 2 • C warming, and the Northeast might serve as a climate refugia for boreal forest species (D'Orangeville et al 2016(D'Orangeville et al , 2018. This is generally also supported by our analysis that suggests a favorable growth response to higher temperatures in the east (figures 1(b) and 2(g)). Seddon et al (2016) pointed out that a remote sensing index can serve as the first step towards addressing why some regions and locations appear to be more sensitive than others, and with the detailed information from monthly vegetation response to climate, we can interpret the results in the context of candidate variables and test hypotheses regarding the causes of local sensitivity of forest ecosystems and forest stands. A good example is the central boreal forest, where we find a patchwork of forest stands that are cold limited versus summer-drought limited in close proximity. A comparison with stand age (figure 3(a)) reveals a striking 'lock and key' pattern, where drought limitations are most pronounced in medium-aged stands (20-50 years) that are growing fast and have high transpiration rates during the growing season. Both young forest stands (0-20 years) and older forests (>50 years) appear less vulnerable to drought, which could be due to young stands having not yet developed a significant transpiration capacity and older stands having a lower density of larger trees, or possibly better access to groundwater. The response is region-specific, and age class structure has for example no effect on vegetation response in most regions of the eastern boreal forest (figure 3(b)), also apparent at larger scales ( figure 1(b)). Age class structure has also been documented as an important driver of how forest stands respond to climate change in western Canada in research that relies on forest inventory plots (McMillan and Goulden 2008, Peng et al 2011, Luo and Chen 2013. Responses can vary significantly due to region, topo-edaphic factors, and species composition.

Inferring causes of local response
Another example for high local variability in vegetation response can also be found throughout the western montane regions from Alaska to Mexico, which generally contain the same response clusters but at different frequencies ( figure 4). Here, the type of vegetation responses is driven by steep gradients of climatic and topo-edaphic factors associated with elevation gradients. High elevation positions generally fall into thermal-limited clusters (in blue toned colors), but low elevations show stronger drought sensitivity (in orange colors). Such response gradients have, for example been documented for dry

Implications and management applications
Our analysis revealed high spatial heterogeneity of climatic response coefficients across most ecosystems. As a consequence, response to anthropogenic climate change will almost always be locally variable, and mitigation strategies need to be carefully developed considering local circumstances. For example, to mitigate drought sensitivity of fast growing or dense forest stands (e.g. Young et al 2017), silvicultural interventions may including thinning to reduce overall evapotranspiration rates of fast growing forest stands and thereby make the trees overall more resistant to drought episodes (Sohn et al 2016, Wang et al 2019. Such a prescription would, however, only be effective where response coefficients have identified this type of vulnerability during the growing season. The high local and regional diversity in growth climate responses also offers an explanation for contradictory findings in studies from the same general region using similar methods. One example is research across the interior boreal forest of Alaska, where some studies report declines of spruce growth (Barber et al 2000), while others documented overall increases in productivity (e.g. Sullivan et al 2017). We find this particular region very diverse in local response coefficients, so that a diversity of vegetation response from individual research studies would not be unexpected. Climate change impacts and corresponding management interventions to mitigate those impacts would similarly need to be tailored to local circumstances in this region.
This study indicates that time series records of remotely sensed vegetation response are now of sufficient length and quality to infer limiting climatic factors, and map them at high resolution. We further provided examples how vegetation response types can be linked non-climatic factors, such as forest stand age, drainage, nutrient regimes and other topoedaphic factors, which may contribute to resilience or exacerbate vulnerabilities. We believe this analytical approach could prove useful for mapping climate vulnerabilities at spatial scales that can support local management interventions to mitigate climate change impacts, especially when validated by ground observations such as long-term forest inventory plots that is usually available to local forest managers (e.g. Das et al 2021).
As suggested by edge detection theory (Dakos et al 2012), the strength and direction the monthly correlation coefficients may also be interpreted as an indicator of the risk of threshold transgressions at the level of species assemblages (Trisos et al 2020), which are of particular interest to resource managers and conservationists (Millar and Stephenson 2015). A formal analysis of threats due to threshold transgressions would require inclusion of climate change projections, but given average ensembles of CMIP6 projections that project warmer and drier conditions for the Southwest (Mahony et al 2022), climate limitations from this analysis suggest that tree populations in this region are highly vulnerable to continued future declines. This has also already been documented based on past climate change trends over the last several decades (e.g. Floyd et al 2009, Ganey and Vojta 2011, Worrall et al 2013, Byer and Jin 2017.

Conclusions
This study contributes a vegetation response descriptor in the form of multivariate clusters of lagged correlations between cumulative annual EVI values and monthly climate variables. The clusters provide more information than a simple drought vulnerability index, conveying at what times of the year plants are vulnerable to climate-induced physiological stress associated with drought and high temperatures. Both our analysis and prior research suggest that drought may play an important role in limiting growth of forested ecosystems for a variety of regions of North America. We show that the multivariate cluster values can be generated and mapped at continental scales at very high spatial resolution, which reveals that vulnerable forest stands can be found across all ecosystems, including a distinct geographic band across the western boreal forest. The high-resolution analysis highlights the heterogeneity of vegetation response to climatic factors, and is also useful to determine which non-climatic factors, such as stand age or topoedaphic factors, contribute to drought vulnerability. This allows improved targeting of vulnerable forest stands for management interventions to address climate change. The analytical approach can further indicate the risk of threshold transgressions at the edge of species and ecosystem tolerance limits, which are of particular interest to resource managers and conservationists in setting appropriate reforestation and ecosystem restoration objectives under climate change.

Data availability statement
All data needed to evaluate the conclusions in the paper are present in the paper and/or the supplementary materials.
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.6084/m9.figshare.20497992.v2.