Hotter and drier climate made the Mediterranean Europe and Northern Africa region a shrubbier landscape

A shift to higher temperatures has left the Mediterranean Europe and Northern Africa (MENA) region more vulnerable to drought and land degradation. We used MODIS LAI (leaf area index) and GPP (gross primary production) deficits, the differences between actual and historical-maximum values, to describe vegetation structural and functional changes and consequential landcover change in response to changing climate conditions during 2001–2019 in the area (20° W–45° E, 20° N–45° N). We found that 1) the vegetation responses varied significantly among eight landcover types with the decreasing importance: forests, savannas, a mosaic of cropland and natural vegetation (CNV), croplands, permanent wetlands, urban land, grasslands, and shrublands, each with distinctive yet overlapping signatures over the ranges of the climate conditions considered. 2) Forests, occupying the coolest and wettest niche, showed the strongest response to severe drought with a lag of 1–3 years and a legacy effect for 10 years. Shrubs, occupying the hottest and driest niche, were the most resilient under a hotter and drier climate. 3) The total areas of savannas and CNV increased by 394,994 and 404,592 km2, respectively, while that of forests decreased by 33,091 km2. Shrublands extended by 287,134 km2 while grasslands and croplands retreated by 490,644 and 225,263 km2. The area of wetlands increased by 49,192 km2, and that of urban land increased by 39,570 km2. A total of 57,649 km2 of barren land became vegetated over the years. Along with higher temperature and more extended period of drought, MENA has evolved towards a shrubbier landscape.


Introduction
Landcover, especially vegetation dynamics over land, plays a critical role in atmospheric processes and the water and carbon cycles globally (Foley et al. 1996;Sitch et al. 2003;Assal et al. 2016). An analysis of FLUXNET data revealed that the exchanges of carbon, water, and energy between terrestrial ecosystems and the atmosphere are limited primarily by water availability when the mean annual temperature is above a threshold of 16 °C (Yi et al. 2010). Climate records from the past 6 decades show that the annual mean temperature of a significant part of the Mediterranean, Europe, and Communicated by Paul Stoy.
Northern Africa (MENA) region have shifted from below to above this temperature threshold (Yi et al. 2014;Allen et al. 2015;Huang et al. 2015). This shift to higher temperatures has left the largely arid and semiarid MENA region more vulnerable to drought and land degradation (Somot et al. 2008;Segui et al. 2010;Friend 2010;IPCC 2013;Zhang et al. 2020;Peñuelas and Sardans 2021). Indeed, some studies have reported drought-induced forest impacts and diebacks in the Mediterranean region (Allen et al. 2010;Carnicer et al. 2011;Peñuelas et al. 2001;Martínez-Vilalta and Piñol 2002;Matusick et al. 2013), as well as shifts in vegetation composition (Jump and Peñuelas 2005;Anderegg et al. 2012).
The drivers of changing vegetation patterns are complex, and human activity as well as climate change is involved (Alexandrian et al. 1998). Indeed, changing climate is clearly linked to changes in human land use in the MENA region (Cramer et al. 2018). However, in this study, we seek to identify statistical correlations between patterns of vegetation and climate in the MENA region. Climate is a strong contributor to variability and change of vegetation and is a matter of both ecological and economic concern, as strong sensitivity to climate can result in rapid land use change (Vanacker et al. 2005;Serra et al. 2008a, b). Over longer time scales, relatively small shifts in background climate can impact the distribution of ecosystems and potentially the viability of agricultural and pastoral systems (Ciais et al. 2005;Vicente-Serrano et al. 2012a, b). Interest in the climate sensitivities of vegetation in semiarid regions is evident in the substantial body of research devoted to characterizing the relationships between precipitation, soil type, land management, and vegetation growth in these water-stressed regions (Zaitchik et al. 2007;Vicente-Serrano et al. 2012a, b). Previous studies have shown that, since 1948, a warmer climate has moved the 16 °C isotherm poleward, leading to a predicted northward shift of vegetation distribution of the Northern hemisphere (Yi et al. 2014). Vicente-Serrano et al. (2012a, b) also showed a decline of the average vegetation cover in the semiarid Mediterranean environments during 1984-2008.
At regional to continental scales, satellite observation offers a feasible and effective approach for monitoring vegetation dynamics. One of the most common observations of vegetation dynamics is based on the normalized difference vegetation index (NDVI). The NDVI, computed as a normalized ratio of reflectance in the near-infrared and red portions of the electromagnetic spectrum, provides a measure of chlorophyll abundance at the Earth's surface and offers an indirect measure of energy absorption and vegetation density (Myneni et al. 1995a, b;Kerr and Ostrovsky 2003).
Several vegetation indices with direct ecological meanings are derived from NDVI using additional algorithms and information. One of them is the leaf area index (LAI). LAI is a fundamental parameter that reflects vegetation structure involved in the processes of fixing atmospheric CO 2 into organic matter. LAI is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as half the total needle surface area per unit ground area in coniferous canopies. Another widely used vegetation index is gross primary production (GPP), which is a measure of vegetation function in the processes of fixing atmospheric CO 2 into organic matter. GPP correlates closely with LAI where LAI is ~ 4 or less, suggesting that leaf area is a critical determinant of GPP in most of the MENA region (Chapin et al. 2011).
This study took advantage of interannual climatic variability during 2001-2019 (Sprintsin et al. 2009;Fischer and Knutti 2015) to characterize the climatic sensitivities of vegetation in the MENA region. We applied the "perfect deficit" approach (Yi et al. 2012) to identify the primary climatic drivers of LAI variability over vegetation across the entire MENA domain aggregated by landcover class.
We used both LAI deficit and GPP deficit as direct measures of climate and non-climate stress on ecosystem structure and function experienced by the vegetation of the semiarid MENA area. We focused on climate stress expressed by some climate indices describing thermal and water conditions and their combined effects. We attempted to establish empirical connections of LAI deficit and GPP deficit with a few climate variables/indices (i.e. annual average temperature (T), annual average precipitation (mm), Dryness Index, Temperature-Precipitation (TP) Index, and Standardized Precipitation-Evapotranspiration Index (SPEI)). Our aim was to understand how landcover changes, especially vegetation structural and functional changes and variabilities, respond to these climate variables: Specifically, we examined which of these variable(s) are the dominant driver(s) for the interannual variability of LAI and GPP deficit. When vegetation types were considered individually, we hypothesized that LAI deficit and GPP deficit of different vegetation types would respond to climate variables differently, and that the distribution of vegetation types may extend or shrink in response to climate change. The direction and intensity of a possible landcover shift were further explored with landcover transition matrices. Such ecological assessment at the temporal scale of decades is much needed for the strategic planning of resource management, for optimizing vegetation productivity and ecosystem services.

Climate and remote sensing data
We collected remote sensing data such as MODIS LAI, GPP and landcover, and various climate measures for the 1 3 MENA region, defined as the bounding box of 20° N-45° N and 20° W-45° E. This is a box-shaped area around the Mediterranean Forests, Woodlands & Scrub biome defined by the World Wide Fund for Nature (WWF) in the MENA region ( Fig. 1).
For the remote sensing datasets of LAI and GPP suitable for this study, we compared the spatial resolution across all available products for the study area and noticed that 500 m was the highest spatial resolution with continuous time series coverage. Using LAI as an example, there are alternative datasets such as the NOAA CDR AVHRR LAI FAPAR: Leaf Area Index and Fraction of Absorbed Photosynthetically Active Radiation, Version 5 and GCOM-C/ SGLI L3 Leaf Area Index (V2). However, the spatial resolution of the NOAA dataset is 0.05 degree, approximately 5.5 km while that of the GCOM-C/SGLI is 2.5 arc minutes, which is about 4.58 km. Both are much coarser than what we employed. Given the large size of our study area, the pixel size of 500 m is assumed to convey sufficient information to address our research questions.
Although the mosaic of plant communities also responded to local geologic, topographic and soil heterogeneity, for which a finer resolution of 1-100 m landcover dataset might be more appropriate, a coarser landcover dataset with 500 m resolution would suffice for our focus on vegetation-climatevegetation interactions and feedback loops in this paper.

Climate indices
We chose two primary instrumental measures (temperature and precipitation) and three compound climate indices which transform basic climate variables such as temperature and precipitation in different ways to assess water stress, one with a linear function, one with an exponential function, and one being standardized, for this study.
We used monthly global surface air temperature (average air temperature at 2 m height) and precipitation (monthly sum) data (01-Jan-2001-31-Dec-2019) at 0.25° resolution provided by ECMWF/Copernicus Climate Change Service within Google Earth Engine. The source of these data ERA5 is the fifth generation ECMWF atmospheric reanalysis of the global climate that combines model data with observations from across the world into a globally complete and physically consistent dataset.
To calculate dryness index, we collected monthly average net radiation (W m −2 ) data (01-Jan-2001-31-Dec-2019) at 0.1 degree resolution from the FLDAS (Famine Early Warning System Network (FEWS NET) Land Data Assimilation System) dataset (Mcnally et al. 2017) provided by the NASA GES DISC at NASA Goddard Space Flight Center within Google Earth Engine. This FLDAS dataset uses Noah version 3.6.1 surface model with CHIRPS-6 hourly precipitation that has been downscaled using the NASA Land Surface Data Toolkit.
We first resampled the radiation data with the same projection system and resolution of the temperature and precipitation data. Net radiation (R n ) is defined as the sum of net shortwave radiation and net longwave radiation in a downward direction.
The dryness index of Budyko (1961) was calculated by: where R n (MJ m −2 year −1 ) and P (mm year −1 ) are, respectively, monthly mean net radiation and precipitation for each grid cell and L is a constant 2.5 MJ kg −1 , the enthalpy of vaporization. Because temperature is often positively correlated with R n, dryness increases in response to increasing T or decreasing P. When there is no precipitation in a month, the dryness index value could be infinitely large. On the global-scale the average net radiation was between 98 and 112.6 W m −2 (Liang 2018) and average precipitation multiplied by L was comparable, therefore, the dryness index values (annual average of the monthly values) were capped under 100 for later analysis. Temperature-Precipitation (TP) Index was derived initially from a model of soil nitrogen variations (Jenny 1984) and further adapted to model vegetation responses to a changing climate (Yi et al. 1996). It is calculated by: where T (°C) and P (mm year −1 ) are annual mean temperature and precipitation for each grid cell, respectively (Yi et al. 1996). TP Index increases in response to decreasing T or increasing P. SPEI (Standardized Precipitation-Evapotranspiration Index) is a measure of water balance, calculated by the standardization of water deficit D, where P (mm year −1 ) and PET (mm year −1 ) are precipitation and potential evapotranspiration for each grid cell, respectively (Vicente-Serrano et al. 2010). The SPEI can measure drought severity according to its intensity and duration. It can identify the onset and end of drought episodes as well. The SPEI is a standardized variable, and it can, therefore, be compared with other SPEI values over time and space. An SPEI of 0 indicates a value corresponding to 50% of the cumulative probability of D, according to a log-logistic distribution (a continuous probability distribution for a nonnegative random variable) of D, over a reference period, while positive values indicate wetter than typical conditions and negative values indicate abnormally dry conditions.
A globally gridded (0.5° resolution) SPEI index calculated from the CRU TS3.23 precipitation and reference evapotranspiration was provided by the SPEI website (https:// digit al. csic. es/ handle/ 10261/ 202305). The SPEI used in the analyses was at the 12-month temporal scale ending December of the previous year.

Perfects, deficits, and relative deficits
Many large-scale vegetation and landcover studies are formulated as statistical analyses of climatological anomalies. This study took a vegetation-centric approach, focusing on the differences in vegetation performance under optimal and suboptimal conditions. Yi et al. (2012) introduced a "perfect-deficit" approach for quantifying links between climate extremes and variations in carbon storage. We applied this approach in developing a climate stress indicator based on our MODIS LAI and GPP data in the MENA area, available for every 8-day period throughout a year, then converted to monthly values. For each pixel, the "perfect" LAI value of the month is defined as the maximum LAI value for this month over the 19 years, and the "deficit" LAI value of the month for each year is defined as the difference between the observed value from the "perfect" value of that month. Therefore, for the entire studied time span (2001-2019) for the ith month, and for the ith month of the jth year, Essentially, we were looking at how the observed LAI for a given month under suboptimal climate and ecological conditions departed from the highest value observed for that time of the year over the 19 years. The perfect value is presumed to correspond to the optimal growing conditions of climate, edaphic features, and other impinging ecological factors for that particular area.
The monthly LAI or GPP deficit values were first averaged within each year to give us an annual mean deficit for each pixel, then were averaged spatially (pixel by pixel) across each of the landcover classes for each year.

Data analyses
There were four sequential steps of data analyses (Supplementary Fig. 1).
Step 1-eight vegetated landcover types with no change over 18 years were extracted to conduct multiple regression analyses; Step 2-regressions of vegetation structural (LAI deficit) and functional (GPP deficit) responses over climate drivers (including temperature, precipitation, TP index, dryness index, and SPEI) were conducted over all vegetated pixels and within each of eight vegetation types; Step 3-results of ten regression analyses and probability density function (PDF) curves of seven variables all pointed to the strong constraints within each vegetation type, which shed new light on landcover changes over the years.
Step 4-using the characteristic order along the spectrums of climate variables and vegetation response variables, direction and intensity of landcover change were quantified using time series and transition matrices.
All LAI, GPP, and landcover data were retrieved within Google Earth Engine from MODIS products (MOD15A2 and MOD17A2, MCD12Q1). All climate datasets (global surface air temperature, precipitation, and net shortwave radiation data with temporal resolution of 1 month) were resolved to a common projection at 0.25° resolution. All annual LAI deficits, annual GPP deficits, and all climate indices were calculated within Google Earth Engine using JavaScript. Because we are looking for long-term largescale patterns, all data were aggregated to annual averages of eight landcover types. All processed spatial data were then exported to Geotiff files with WGS84 and 500 m resolution for further statistical tests and graphing in R-Studio (version 1.2.1335). All maps were made in ArcGIS Desktop 10.7.1.11595.

Climate drivers
The first objective of the study was to find the climate indices that the vegetation structure indicator LAI and function indicator GPP responded to the most, and how those relationships vary across different vegetation types. We found that overall LAI deficit was negatively correlated with annual average temperature (Fig. 2a, adjusted R 2 = 0.529, p < 0.001) and dryness (Fig. 2g, adjusted R 2 = 0.583, p < 0.001), positively correlated with precipitation ( Fig. 2c, adjusted R 2 = 0.613, p < 0.001), and TP Index (Fig. 2g, adjusted R 2 = 0.491, p < 0.001), and not correlated with SPEI ( Fig. 2i, adjusted R 2 = 0.007, p = 0.155) ( Table 1). Among the eight landcover categories, only forest consistently showed a similar pattern to the overall vegetation structural responses while shrubland showed an opposite trend or no trend.
Among all the eight landcover categories: only forest showed an almost consistently similar pattern to the overall vegetation structural responses while shrubland showed an opposite trend or no trend. The only exception was in precipitation. Within each vegetation type, GPP deficit decreased when precipitation increased. In another word, the more precipitation, the less reduction of GPP within each vegetation  (LAI deficit) and functional (GPP deficit) responses over climate drivers including a, b annual average temperature (°C, T), c, d accumulative precipitation (mm, P), e, f Temperature-Precipitation (TP) index, g, h dryness index, i, j Standardized Precipitation Evapotranspiration Index (SPEI) for all vegetated pixels (black lines, solid if p < 0.01, dotted if p > 0.05) and within each of 8 vegetation types (color coded with 95% confidence intervals in gray); for each panel, probability density function (PDF) curves were plotted with the same scheme of color code, horizontal on the top for the climate variables (T, P, TP index, dryness or SPEI), and vertical on the right for the response variables (LAI or GPP deficit) type as we would normally expect under semiarid climates. The overall correlation between GPP deficit and precipitation was positive.

Landcover-specific signatures over the spectrums of climate conditions
Based on the probability density function (PDF) curves for both dependent variables (i.e. vertical PDF graphs for LAI and GPP deficits) and independent variables (i.e. horizontal PDF graphs for climate indices) (Fig. 3), each vegetation cover demonstrated unique and sometimes overlapping distribution signatures over the span of structural and functional responses, and climate conditions. For example, forests had by far the highest LAI deficits (11.3 ± 1.1 m 2 m −2 ), followed by savannas (7.1 ± 0.7 m 2 m −2 ), CNV mosaic (6.8 ± 0.6 m 2 m −2 ), croplands (6.8 ± 0.5 m 2 m −2 ), wetlands (6.1 ± 0.5 m 2 m −2 ), urban land (3.9 ± 0.3 m 2 m −2 ), grasslands (3.3 ± 0.4 m 2 m −2 ), and shrublands (1.5 ± 0.2 m 2 m −2 ) in descending order ( Table 2). The LAI deficit of shrublands was less than half of that of grasslands, which was the second lowest among all. Likewise, forests had the highest GPP deficits, followed by the same descending order (Table 2). This order reflected the differential responses, both structural and functional, to these climate variables among these vegetation types, and their differential responses to climate stress.
When looking into the span of climate conditions, stratification of landcover classes by climate conditions is evident in the data record. Forests occupied the coolest (10.0 ± 0.5 °C) and wettest (1102.3 ± 87.3 mm) niche while shrublands occurred in the hottest (17.5 ± 0.3 °C) and driest (342.2 ± 49.2 mm) area, polar opposite to forests. Each vegetation landcover type is constrained within a limited range of climate conditions but each has a certain level of overlap with others. Within a constrained range of climate conditions, alternative steady states (vegetation types) may exist.

Remotely sensed areal change of landcover
According to the MODIS annual landcover product during 2001-2018, the total areas of savannas and CNV mosaics in MENA increased by 394,994 km 2 (5.3% of savanna area in 2001) and 404,592 km 2 (29.7% of CNV area in 2001), with a significantly steady increasing trend (R 2 = 0.930 and 0.913, respectively), while the area of forests decreased by 33,091 km 2 (1.6%) (Fig. 3a-c) despite the fertilizer effect of elevated ambient CO 2 .
More specifically, the area identified as forest in this analysis decreased by 41,401 km 2 (a 2.9% drop) over 3 years following the drought/heat wave of 2002, and did not recover to the pre-drought level until 2014, 11 years after the drought. It is unclear whether other droughts (e.g. 2005) also impacted vegetation dynamics and prolonged the recovery. A second major drop happened between 2016 and 2017, leaving an overall areal decline between 2001 and 2018.

Matrices of landcover change
To quantify the direction and intensity of change, we used the area of each vegetation type that had changed to make transition matrices. A transition matrix could illustrate the detailed vegetation transition during either a short-term period (i.e. during 2002-2003 right after the 2002 heatwave/drought) or a long-term period (i.e. during 2001-2018) ( Table 3; Fig. 4).
In Table 3, the diagonal cell areas remained in the same vegetation types. The upper right triangle portion indicates transitions to a less productive vegetation type, while the lower-left triangle indicates transitions to a more productive vegetation type, based on what we have learned from PDF curves of LAI and GPP deficits. During 2002During -2003,775 km 2 of the forested area became less productive and remote sensing identified them as savannas, while a total of 22,680 km 2 of barren land became vegetated (turned into grassland and shrubland). As a result, a total of 52,986 km 2 became less productive based on the four categories with over 5000 km 2 (i.e. from forest to savanna, from savanna to grassland, from cropland to grassland, and from grassland to shrubland) while a total of 56,814 km 2 became more productive based the five categories listed in the lower triangle (i.e. from grassland to savanna, from grassland to cropland, from shrubland to grassland, and from barren land to either grassland or shrubland).
During 2001-2018, 55,948 km 2 of the forested area became less productive savannas (Fig. 4, visible in southern France and western Portugal) while 54,475 km 2 of savannas had the reversed transition, leaving a net loss of forests by a mere 1473 km 2 . A total of 85,257 km 2 of barren land became vegetated (mostly turned into grassland and shrubland) and 27,609 km 2 of vegetated pixels became barren, leaving a net gain of 57,649 km 2 of the vegetated area from barren land (0.9%). Such green-up is visible in parts of northern Algeria and Tunisia, eastern Iran, northern Saudi Arabia and northern Egypt (Fig. 4). There was a total area of 413,922 km 2 that became more productive and 260,560 km 2 became less productive (based only on categories of transition with 20,000 km 2 or more).

Discussion
One of the most climate sensitive regions of the world is the Mediterranean basin (MENA) (IPCC 2013). Changes in human population distribution, conflict, land management and land use have occurred concurrent with changes in climate and have altered vegetation patterns across the region, with feedback to the regional climate (Tanrivermis 2003;Serra et al. 2008;Bajocco et al. 2012;Millán 2014;Luyssaert et al. 2014;Perugini et al. 2017;World Bank 2021;Wolpert et al. 2020;Ruiz and Sanz-Sánchez 2020). We recognize that in evaluating causes of regional vegetational changes it may be difficult to separate those induced by climate regionally from other causes of landscape-scale changes such as political and economic changes, particularly in regions of the globe that lack extensive data on land use and management, and considering the interactions of these societal factors with climate change is an important goal for future analyses, although difficult to do at this large spatial scale. Well-known examples of land use and management change not directly driven by climate in the MENA region include afforestation in Israel (Rotenberg and Yakir 2011), marsh drainage in southern Iraq (Hashim et al. 2019), irrigation intensification in several countries (Krakauer et al. 2020), and widespread urbanization (Rozalis et al. 2010). This report, however, is focused particularly on those potential causes that can be evaluated by satellite remote sensing and climate data and operate with some consistency across political boundaries, with an eye toward evaluating sub-continental to global-scale changes elsewhere in the world for which these data are readily available, while more fine-scale ground-based observational data (e.g. land use and management) are more difficult to retrieve. While policy and management can have an impact on local scale and stepwise changes, regional scale change in a long-term, which is our focus, should be dominated by climatic changes.
The aim of the study was to use LAI and GPP deficits to identify and better understand the impacts of the climate on landcover in the MENA region during 2001-2018. Specifically, we would like to 1) determine critical climate drivers for the variability of the landcover as a whole and several ecosystem types individually, 2) explore landcover-specific signatures over the span of climate conditions, and 3) identify landcover changes over time and relate these changes to climate changes (and other non-climatic drivers).

Vegetation-centric approach of LAI and GPP deficits
LAI and GPP deficits provide a different perspective into the climate-vegetation interaction, in contrast to the traditional climate anomaly approach.
Through examining the responses of LAI deficit and GPP deficit of all vegetation types combined and individually to various climate indices including T, P, TP Index, dryness and SPEI, both LAI and GPP deficits proved to be effective indicators of how ecosystem photosynthetic structures and functions respond to climate stress. This study took this as a starting point to further explore the responses of the two deficits to interannual variations of the climatic conditions for the MENA region.
The perfect-deficit approach was used at a monthly scale derived from the original 16-day product, which is much less sensitive to extreme value than a daily or hourly product. Technical documentation of MODIS LAI and GPP products listed detailed uncertainty evaluation rules. The MODIS LAI/FPAR algorithm consists of a main look-up-table (LUT) based procedure (Knyazikhin et al. 1998) that exploits the spectral information content of the MODIS red (648 nm) and near-infrared (NIR 858 nm) surface reflectance, and the backup algorithm that uses empirical relationships between Normalized Difference Vegetation Index (NDVI) and canopy LAI and FPAR. The theoretical estimates of uncertainties (%) in the BRFs used in the C6 LAI/FPAR algorithm are 20-30% for red and 5-15% for NIR. The validation strategy for MODIS GPP was initially established based on Running et al. (1999) and has been widely accepted and applied in numerous publications (e.g. Turner et al. 2003Turner et al. , 2004Turner et al. , 2006. In general, regressions of LAI deficits and GPP deficits behaved very similarly, except to precipitation (Fig. 2c, d). Both Forest and overall LAI deficits (m 2 m −2 ) had a positive correlation with annual precipitation (mm), which was quite counter-intuitive. To understand this, we can start with Fig. 2d.
Within each vegetation type, GPP deficits (kg C m −2 ) decreased when annual precipitation (mm) increased, which meant less GPP loss within the same year with more precipitation. This is consistent with existing literature (Wu et al. 2011;Nielsen and Ball 2015). The overall regression of GPP deficits over precipitation, however, had a counterintuitively positive slope, which was overwhelmingly dominated by intrinsic variation among different vegetation types (i.e. forests generally have the highest GPP base. Therefore, its absolute increase in response to increasing precipitation should be the highest). The overall pattern was more strongly driven by intrinsic productivities of each vegetation type rather than local responses to climate variables alone.
Likely, as in Fig. 2c, the overall regression of LAI deficits over precipitation had a positive slope, with the same reason as for GPP deficits-they are tightly constrained within Table 2 Means and standard deviations (SD) of vegetation responses (i.e. annual Leaf Area Index (LAI) deficit and gross primary productivity (GPP, in g C m −2 year −1 ) deficit), climate conditions (i.e. annual mean temperature (°C), precipitation (mm month −1 ), temperature-precipitation (TP) index, dryness index, Standardized Precipitation Evapotranspiration Index (SPEI)), and area (100,000 km 2 ) for each landcover type that remained as the same landcover type over 19 years (pixels at 500 m resolution) each vegetation type. However, within forests, the regression slope was also positive. To understand this pattern, we should also look at what has happened to the area ( Fig. 3a; Table 3). In 2002, the MENA region experienced a widespread severe drought. The total area of forests decreased drastically from 2,102,011 km 2 in 2002 to 2,040,610 km 2 in 2003, and it continued to decline for two more years that followed, with 2005 having the smallest area of forests. It took another 10 years (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) for the area of forests to restore to the pre-drought level. In other words, there were 1 3 3 years of lag and 10 years of legacy effects in terms of the distribution area of forests. Similarly, forest area was at its highest in 2017, then declined drastically in 2018, possibly due to another extreme event.

Years
If we looked at the time series of LAI deficits, similar lag and legacy effects existed. The decline of LAI happened in the following few years, therefore, LAI was at its highest when precipitation was at its lowest, or LAI deficit was the lowest when climate conditions (such as precipitation) suddenly became less optimal. Therefore, we see this counterintuitive positive slope for regression of LAI deficits over annual precipitation.

Five climate indices
There are many climate indices available. For this study, we used two primary instrumental measures (temperature and precipitation) and three compound climate indices, each presenting a different group of indices that emphasize certain aspects of hydro and/or thermal conditions. Dryness index is expressed in a linear relationship between net radiation and precipitation; TP index is expressed with an exponential formula, containing temperature and precipitation; SPEI is a standardized index allowing flexibility of temporal scales.
Because Dryness, TP index, and SPEI all contain the information of precipitation and are not independent of each other, multiple univariate regressions might be easier to interpret than a single multivariate regression.
When response variables LAI deficit and GPP deficit were plotted against each of the five climate indices, all except SPEI were tightly clustered within each vegetation type ( Fig. 2; Table 2). Because SPEI is a standardized index, site specificity has been removed during the calculation. On the other hand, these tightly distributed clusters may indicate strong resistance to transitions between vegetation types and alternative stable states caused by the changing climate (Hirota et al. 2011;Chapin et al. 2011). Although they currently seem to be within the resilience ranges, further studies with finer spatial and temporal scales, and with a finer categorization of vegetation types and functional groups could tell a different story. Patterns operating at finer scales will be different from what has been looked at (bigger scales), but this may not necessarily lead to a better understanding of the processes operating at the coarser resolution we have studied (Liao et al. 2020;Chen 2021;Chen et al. 2021;Guo et al. 2021).

Eight vegetation covers and their PDFs
Despite their structural and functional responses to climate and ecological conditions, most vegetation types are tightly constrained, representing steady states for specific combinations of climate conditions and their immediate alternative steady states if strong enough forces change them one way or another (Fig. 2).
Opposite to patterns in the forests, shrubland LAI deficit was positively correlated with dryness index. However, that does not necessarily mean if the area experienced prolonged directional change related to temperature, water availability and solar radiation, a transition of vegetation coverage between forests and shrubland could happen. When we looked deeply into the ranges of annual average temperature (10.0 ± 0.5 vs 17.5 ± 0.3 °C), precipitation (91.86 ± 7.28 vs 28.52 ± 4.1 mm month −1 ), TP index (0.51 ± 0.03 vs 0.19 ± 0.01), Dryness (1.00 ± 0.25 vs 11.11 ± 1.83), and SPEI (− 0.34 ± 0.61 vs − 0.64 ± 0.62) of forest and shrubland (Table 1), forest, and shrubland are experiencing very different local climate changes, thus leading to polar opposite responses. The results are consistent with the "16 °C threshold" prediction (Yi et al. 2010), which stated that the exchanges of carbon, water, and energy between terrestrial ecosystems and the atmosphere are limited primarily by water availability when the mean annual temperature is above a threshold of 16 °C and by temperature when below 16 °C. Other studies also have looked into the changes in landcover under the influence of climate change. For example, Vidal-Macua et al. (2017) studied factors affecting forest dynamics in the Iberian Peninsula from 1987 to 2012 based on Landsat scenes, and found that the geographical transition from shrubland to forests is closely related to higher soil moisture (Topographic Wetness Index-TWI) and lower winter solar radiation. Meanwhile, strategies aiming at decreasing the risk of decline and promoting resistance to abrupt stress in the short term may not enhance long-term resilience (Vilà-Cabrera et al. 2018).

Landcover changes and their implications
Combining with transition maps, the transition matrix is a powerful tool to specify the direction and intensity of landscape changes over the region. Comparing short-term and long-term transition matrices can also indicate how vegetation responds to extreme events (such as the 2002 drought) immediately and in a long run.
Even though the percentage of vegetated area of forest (4.7%) was far less than cropland (41.2%), grassland (23.9%), savanna (17.4%), and shrubland (6.0%), the patterns of LAI deficit of forest were always consistent with the LAI and GPP deficit of all vegetation types, which indicates the strong forest LAI response to climate and its dominating influence over the response of the integrated vegetated landcover in our analysis. The MENA area was still benefiting from the increase of ambient temperature and CO 2 during 2001-2018, as indicated by the green-up of barren land and the gradual recovery of forested areas after the 2002 drought. Vegetation response to elevated CO2 is conserved across a broad range of productivity (Norby et al. 2005), but what was initially observed in the temperate free-air CO2 enrichment (FACE) experiments (Hendrey 1992;Hendrey et al. 1999) may not be representative of other regions (Hickler et al. 2008), particularly in degrading landscapes, and belowground response could be more continuous than aboveground response to CO 2 enrichment (Jackson et al. 2009) which cannot be detected by satellite data.
This analysis of available data indicates that overall, forests, cropland, and grassland are more vulnerable to climate stress, thus declining over the past 18 years, while savanna and shrubland are more resilient and their distribution expanded. Expansion of shrubland seems to be counter to the need for improved agriculture and forestry goals of land managers (World Bank 2021), supporting our thesis that climate change is a key driver of the trends described here.
The two dominant vegetated landcovers in this area, cropland (41.2%) and grassland (23.9%) ( Table 2), fluctuated in opposite patterns over the years (Fig. 3). This is likely due to the expansion of cropland that was usually an encroachment into grassland, while fallow croplands generally returned to grassland.
One of the most striking patterns was the areas of increased shrubs that replaced cropland and grasses (upstream of productivity) and barren land (downstream of productivity). According to the areal changes of vegetation covers during 2001-2018 (partly in Table3b), a total area of 38,938 km 2 was converted from grassland to shrubland and 22,897 km 2 from cropland to shrubland. At the same time, 37,853 km 2 of shrubland was converted back to grassland and 14,163 km 2 back to cropland, perhaps as a result of improved land management. Further, a total area of 52,376 km 2 was converted from barren land to shrubland, while 15,232 km 2 of shrubland were converted back to barren land, and neither of these directional changes would seem to be a desirable outcome for current land management practices (World Bank 2021). There have been various reports on local shrubland dieback in south Spain (Lloret et al. 2016) and the overall trend of shrubland cover based on satellite data for the entire MENA region shows increasing dieback.
With higher temperature and more extended period of drought, MENA is becoming less productive (with fewer forests and croplands) with more shrubby vegetation covers, with just a small fraction of barren land (less than 1%) that became vegetated. Although some authors (e.g. Bastin et al. 2019) have recently advocated for tree plantation to combat desertification derived from global change, our study suggests that a more resilient (and short) vegetation might be more suitable for restoration programs in areas like MENA.
The central location and average altitude of forest vegetation have not yet changed despite the changes in temperature and dryness. On the other hand, the mean and median latitude of cropland decreased over the last decade, likely due to human effort in improving irrigation systems for agriculture purposes. The mean latitude and median elevation of grassland also decreased over the last decade. Any increase in grassland probably benefited from the fertilization effect of elevated ambient CO 2 concentration. This also might have a tight association with the escalating incidence of wildfire near arid and semiarid areas throughout the world, some of which have caused devastating losses (Bladon 2018;Bowman et al. 2020;European Commission 2021).
For 2001-2018, the mean latitude (R 2 = 0.864, p < 0.001) and median latitude (R 2 = 0.513, p < 0.05) of cropland both decreased significantly. Such a counter-intuitive change is probably due to the huge effort of human intervention: establishing irrigation systems and building dams for agriculture (United Nations 1999;FAO 2008).
The mean latitude (R 2 = 0.779, p < 0.001) and median elevation (R 2 = 0.473, p < 0.05) of grassland also decreased significantly. Prior field manipulation experiments have shown that elevated ambient CO 2 concentration can stimulate the growth and accumulation of standing biomass of grassland in semiarid areas, acting like a carbon fertilizer and a booster of water use efficiency (Dijkstra et al. 2010). Grasses dieback during the dry seasons, which naturally turn into standing fuel, causing increasing wildfire risk and devastating loss. This has been happening throughout the world, especially in Mediterranean ecoregions such as California, US (Abatzoglou and Williams 2016;Parks et al. 2017;Goss et al. 2020).
The area of permanent wetlands increased by 49,192 km 2 (24.9%) most likely due to sea level rise. At the same time, urban land also increased by 39,570 km 2 (2.9%) (Fig. 3g, h), which reflects the urbanization in the region.

Conclusions
• LAI and GPP deficits, the vegetation-centric approach, provide a useful means to differentiate structural and function responses among different vegetation types under climate stress in the MENA regions, showing the areas and productivities of different vegetation types have experienced significant short-term and long-term changes in response to varying climate and non-climatic (e.g. land management) conditions during 2001-2019. • Over the study period, the areas of savannas and CNV mosaics increased by 394,994 km 2 and 404,592 km 2 , respectively, while that of forests decreased by 33,091 km 2 . Meanwhile, shrublands extended by 287,134 km 2 while grasslands and croplands retreated by 490,644 km 2 and 225,263 km 2 , respectively. The area of wetlands increased by 49,192 km 2 , and that of urban land increased by 39,570 km 2 . Finally, 57,649 km 2 of barren land became vegetated over the years.
• Vegetation responses to climate variations depend on vegetation types with distinctive yet overlapping signatures over the span of climate conditions considered. The climate sensitivity decreases in the following order: forests, savannas, a mosaic of cropland and natural vegetation (CNV), croplands, permanent wetlands, urban land, grasslands, and shrublands. Shrubs were the most resilient under a hotter and drier climate. Forests showed the strongest and most dominating response to severe drought with a lag of 1-3 years and a legacy effect for 10 years.