Stem Mortality and Forest Dieback in a 20-Years Experimental Drought in a Mediterranean Holm Oak Forest

Climatic warming is predicted to increase the intensity and duration of extreme weather such as droughts and heat waves. Climate change could therefore increase stem mortality and forest dieback in many ecosystems around the world, especially in semi-arid forests. We investigated the influence of climatic conditions and a throughfall displacement experimental treatment (15% decrease in the amount of soil moisture) on the intensity of stem mortality in a Mediterranean forest. We also investigated the use of remotely sensed data from MODIS as a tool to estimate the intensity of stem mortality and a possible strategy of forest management to mitigate this forest dieback induced by climate change. Stem mortality was higher when mean annual temperature was higher and rainfall was lower, especially during spring and summer. Mortality was higher and more affected by the drought treatment in the holm oak, Quercus ilex, than in co-occurring species of tall shrubs better adapted to drier conditions. Two spectral MODIS indices (NDVI and EVI) were good predictors of stem mortality, but others (GPP, PsnNet, and NPP specifically calculated to predict forest health and development) were not clearly correlated with stem mortality. Selective stem thinning strongly reduced stem mortality (especially in Q. ilex) by buffering the effects of climate change on forest structure. Finally, the possible future substitution of the current dominant species of this forest, Q. ilex, by species of tall shrubs better adapted to drought is also discussed.


INTRODUCTION
Higher air temperatures are projected globally for the coming decades, and climatic patterns around the world are changing: the intensity and duration of extreme climate conditions, such heat waves, floods, and droughts, are predicted to increase (Wang et al., 2012;IPCC, 2018). Higher rates of evapotranspiration in some semi-arid and Mediterranean areas would be induced by these higher temperatures, with no any increase of precipitation (IPCC, 2018). General circulation models thus project an average decrease in soil moisture of 15% in the next 50 years and a return period of extreme droughts 10-fold shorter than in the twentieth century in Mediterranean regions (Bates et al., 2008). Forest ecosystems in these Mediterranean areas, seasonally exposed to water stress, may be particularly vulnerable to even slight increases in water deficits, which can induce reductions in tree growth (Barbeta et al., 2013;Liu et al., 2015Liu et al., , 2018 and crown conditions (Carnicer et al., 2011;Galiano et al., 2012) and increases in tree mortality (Breshears et al., 2005;Allen et al., 2010;Peng et al., 2011;Williams et al., 2012). Precipitation anomalies have also been strongly correlated with the subsequent occurrence of hot extremes in most areas of the world (Mueller and Seneviratne, 2012), and an increased frequency of heat waves is expected in Mediterranean regions, coinciding with summer drought and higher evapotranspiration and thus a lower availability of soil water (Fischer and Schar, 2010).
Several species of tall shrubs associated with holm oak (Quercus ilex L.) forests have lower growth rates but are more resistant to drought conditions than the holm oak, the current dominant species in these forests Ogaya and Peñuelas, 2003). Higher mortality rates and lower seed production and seedling survival of Q. ilex under future drier conditions could drive a progressive substitution of dominant oak species by the co-occurring species more resistant to low water availability (Lloret et al., 2004;Ogaya and Peñuelas, 2007;Liu et al., 2015Liu et al., , 2018. The selective thinning of stems is an ancient practice of forest management widely used to enhance the growth and development of the remaining stems, increasing the availability of resources per stem when stem density decreases. Selective thinning under a scenario of increasing drought could therefore partly offset this water deficit by reducing competition for water resources and thus contribute to the stabilization of forest structure and function. Spectral indices from remotely sensed data are widely used to evaluate the structure and functioning of terrestrial ecosystems . The Normalized Difference Vegetation Index (NDVI) is the most widely used index, mainly as an estimator of the fraction of the photosynthetically active radiation absorbed by the canopy (Tucker et al., 1985). The Enhanced Vegetation Index (EVI) was developed to reduce the noise produced by soil background and atmospheric aerosols and to reduce the saturation of the reflectance signal at increasing levels of green biomass, which is commonly associated with NDVI (Huete et al., 2002). The MODIS system, though, also provides several indices specifically calculated to estimate forest productivity [Gross Primary Production (GPP), Net Primary Production (NPP), and Net Photosynthesis (PsnNet)].
The main objectives of this study were: (1) to identify the relationships between stem mortality and climatic conditions, (2) to determine the effect of an experimental reduction of soilwater availability on stem mortality in the dominant species of the Mediterranean forest, (3) to test the utility of selectively thinning stems to mitigate the effects of the decrease in water availability, and (4) to evaluate the use of remotely sensed data from satellites as a tool for detecting and monitoring forest dieback and stem mortality.

Site Description
The study was conducted in the Prades Mountains, Catalonia, northeastern Spain (41 • 21 ′ N, 1 • 2 ′ E), on a 25 • south-facing slope at an altitude of 950 m a.s.l. The soil is a Dystric Cambisol over Paleozoic schist, ranging in depth from 35 to 100 cm. Annual mean temperature and total rainfall from 1975 to 2018 averaged 11.9 • C and 654 mm, respectively, with hot and dry summers and most of the rainfall concentrated in spring and autumn. The vegetation is a very dense multi-stem forest (16,616 stems ha −1 ), due to an ancient management to obtain charcoal, dominated by Q. ilex, Phillyrea latifolia L., and Arbutus unedo L., with an abundance of other evergreen species well-adapted to dry conditions such as Erica arborea L., Juniperus oxycedrus L., and Cistus albidus L. This forest was not perturbed during the last 70 years, and the maximum heights of the dominant species are about 6-10 m.

Meteorological Data
An automatic meteorological station installed at the study site has monitored temperature, photosynthetically active radiation, air humidity, and precipitation since late 1998. We had previously estimated temperature and rainfall from another meteorological station located in Poblet Monastery, 5.6 km northeast of our study area at 510 m a.s.l. The Poblet station collected data from late 1974 to July 2002, so we had climatic data from both meteorological stations from August 1998 to July 2002. We calculated linear correlations between temperature (R 2 = 0.97) and rainfall (R 2 = 0.75) for the data from these two stations during this period to estimate climatic data for the study site since 1975.

Throughfall Displacement and Stem Mortality
Eleven 15 × 10 m plots were delimited at the study site at the same altitude along the slope; four were control plots, other four were throughfall displacement exclusion (TDE) plots, where the input of rainwater decreased the amount of soil moisture by 15% using permanent plastic strips (Figure 1), and the other three  plots were selectively thinned, where ca. 25% of the basal area of each species was removed ( Table 1). In drought plots, six plastic strips were installed throughout plot slope, about 1 m above soil surface and down of tree canopies. These plastic strips covered also 2 m above plot margin, and 2 m down of plot margin, to avoid edge effect. Each plastic strip consisted of four rigid PVC tubes joined together by transparent PVC between them (Figure 2), these strips covered about 30% of plot surface and partially intercepted throughfall and removed it downslope the plots. Soil moisture was measured each annual season using a time domain reflectometer (Delta-T Devices, Cambridge, UK) in four randomly selected sites per plot. The reflectometer measured soil moisture in the first 15 cm depth. All living stems >2 cm in diameter at a height of 50 cm in the plots were labeled and classified by species at the start of the experiment. All stems that had died during the year were counted in each of the following winters to obtain the rates of mortality per year, plot, and species. The stem-mortality rate (M) was calculated as (Sheil et al., 1995): where N o is the number of stems at the start of the study period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), N t is the number of stems at the end of the study period, and t is the number of years of the study period. The control and drought treatments began in 1999 and ended at the end of 2018, whereas M in the thinned plots was determined from 2012 (immediately after thinning) to 2015.

Remote Sensing
The remotely sensed data for 2000-2012 were obtained from MODIS sensors onboard the Terra and Aqua satellites. The Terra satellite orbits Earth from north to south in the morning, and the Aqua satellite orbits from south to north in the afternoon. We used NDVI and EVI data obtained from the 16-days Terra MOD13Q1 and 16-days Aqua MYD13Q1 products, combined in eight-day series based on a selected 250 × 250 m pixel corresponding to our study site, and including the studied plots and the surrounding area.
We also used GPP and PsnNet data obtained from the 8days Terra MOD17A2 product and NPP data obtained from the annual Terra MOD17A3 product. All these products provided an 1 × 1 km pixel for 2000-2010 centered in the studied area, and including the studied plots and the surrounding area. Mean annual values of the remote-sensing indices were calculated as an annual integration for each growing season covered by NDVI, EVI, GPP, and PsnNet, whereas NPP was calculated directly from annual data.

Statistical Analyses
We used the meteorological data for 1975-2018 to estimate the trends of changes in air temperature and rainfall throughout this period. The difference of soil moisture in control and TDE plots was tested using an ANOVA, where soil moisture was the dependent variable, and treatment the factor. General linear models were also conducted to determine the relationships between the rate of stem mortality and the changes in air temperature and rainfall during the study period. An analysis of variance (ANOVA) was performed, with mortality rate as a dependent variable, and species and treatment (control and TDE plots) as independent factors. A similar ANOVA was then performed to compare stem mortalities in the control and TDE plots. We next calculated the difference in annual M for each year and species and selected 5 extremely dry years and 5 wet years. Another ANOVA was performed, with the difference in mortality rate as a dependent variable and species and treatment as independent factors. Other general linear models were used to determine the relationships of the rate of stem-mortality with mean annual NDVI, EVI, GPP, PsnNet, and NPP. We used linear regressions because they identified the relationships that fitted better with all variables studied. M was arcsin m 0.5 transformed to satisfy the assumptions of a normal distribution. All linear models were performed using StatView (SAS Institute Inc., Cary, USA).

Meteorological Data
Mean annual temperature ranged from 10.7 • C in 1984 to 13.1 • C in 2011, and total annual precipitation ranged from 355 mm in 2015 to 1021 mm in 2018 (Figure 3). Mean annual air temperature during 1975-2018 increased by about 1.7 • C (R 2 = 0.50). This increase was due to the strong increase of 2.7 • C in spring and summer (R 2 = 0.60), but autumn and winter temperatures did not change. Annual precipitation, however, was highly variable, with only a non-significant trend toward a slight decrease from 1975 to 2018, but precipitation during spring and summer decreased ca. 200 mm (R 2 = 0.26).

Stem Mortality
Annual M for all species together was positively correlated with mean annual temperature (R 2 = 0.17) and negatively with spring and summer precipitation (R 2 = 0.18) (Figure 4). Stem mortality was higher in Q. ilex than P. latifolia and A. unedo throughout the study period, and M was significantly higher in the TDE than in the control treatment for all species together  Frontiers in Forests and Global Change | www.frontiersin.org FIGURE 5 | Average stem-mortality rates in the control and TDE plots for each of the three species studied. Error bars correspond to the standard error of the mean (n = 4). Different letters indicate significant differences (P < 0.05).
( Figure 5). A comparison between M and the meteorological data indicated larger increases in stem mortality induced by the TDE treatment in hotter and drier years (Figure 6). The difference in M between the control and TDE plots was larger for Q. ilex than P. latifolia and A. unedo, and this difference for Q. ilex was larger in extremely dry than wet years (Figures 7, 8), even when M for this species was lower in the TDE than in the control plots in wet years. The difference in M for P. latifolia and A. unedo between the control and TDE plots, however, was not significant in both extremely dry and wet years. Finally, selective thinning strongly decreased stem mortality only in Q. ilex, but not in P. latifolia or A. unedo (Figure 9). Stem mortality in the control plots was higher for Q. ilex than P. latifolia and A. unedo, and M in the thinned plots was lowest for Q. ilex.

Remote Sensing
Mean annual NDVI (0.80) was highest in 2003, and mean annual EVI (0.38) was highest in 2008, when annual precipitation for 2000-2012 was highest (926 mm in 2003 and 837 mm in 2008). Mean annual NDVI and EVI (0.74 and 0.33, respectively) were lowest in late autumn 2011 and winter 2011-2012, coinciding with the end of an extremely dry period. Mean annual NDVI and EVI were positively correlated with annual precipitation, and negatively correlated with mean annual temperature and with M ( Table 2).
GPP and PnsNet had another annual pattern than NDVI and EVI, with low values in winter, peaking in spring, decreasing in summer, increasing again in autumn, and a final decrease in late autumn and the following winter. GPP was correlated only with FIGURE 6 | Annual stem-mortality rates in the control and TDE plots for all species studied (A) and annual meteorological data (mean air temperatures and total rainfall) (B). Black dots represent temperature values, and gray bars represent rainfall values. Error bars correspond to the standard error of the mean (n = 4). *P < 0.05. annual precipitation, but PsnNet and NPP were not significantly correlated with M or any climatic variable ( Table 2).

Stem Drought-Induced Mortality
Drought-induced stem mortality is likely to increase in the future. Climatic models forecast an increased frequency of drought together with higher air temperatures (Dai, 2013;Allen et al., 2015). The progressive increase in air temperatures accompanied by a slight decrease in rainfall, especially in spring (when plants flower and develop new leaves) and summer (coinciding with the most critical period for plant development and survival), is exacerbating dry conditions and driving stem mortality and a decline in tree growth in this holm oak forest (Ogaya and Peñuelas, 2007;Barbeta et al., 2013;Ogaya et al., 2016;Liu et al., 2018). The demand for atmospheric moisture FIGURE 7 | Annual stem-mortality rates in the control and TDE plots for Q. ilex (A), P. latifolia (B), and A. unedo (C). Error bars correspond to the standard error of the mean (n = 4). *P < 0.05.
FIGURE 8 | Average differences in stem-mortality rates between the control and TDE plots in dry and wet years for each of the three species studied. Error bars correspond to the standard error of the mean (n = 4). *P < 0.05. has been described as the primary mechanism by which high temperature influences drought-sensitive tree populations (Liu et al., 2013). This mechanism is illustrated by our finding that mean annual temperature was positively correlated with stem-mortality rate, whereas spring and summer precipitation were negatively correlated with M. The deterioration of crown condition has been associated with a reduction in levels of nonstructural carbohydrates (Rosas et al., 2013), that could lead to a carbon starvation with a combination of hydraulic failure (Sevanto et al., 2014;McDowell and Allen, 2015), which may FIGURE 9 | Average stem-mortality rates in the control and thinned plots for each of the three species studied. Error bars correspond to the standard error of the mean (n = 4). *P < 0.05.
contribute to the forest dieback as it was observed in the studied site (Ogaya et al., 2016). As expected, mortality rates were higher for Q. ilex than the co-occurring species of tall shrubs (Figure 3) that are better adapted to drought Ogaya and Peñuelas, 2003). Stem mortality for Q. ilex also strongly increased under low water availability in the TDE treatment or extremely dry years, which did not occur in the other species (Figures 4, 7). Most tree species globally are routinely subject to levels of physiological water stress near their limits of tolerance  (Choat et al., 2012), and experimental evidence suggests high tree mortality induced by drought with relatively modest increases in temperature (Adams et al., 2009;Duan et al., 2014). Higher stem mortality (Ogaya and Peñuelas, 2007;Barbeta et al., 2013;Liu et al., 2015) and lower seedling survival (Lloret et al., 2004) in Q. ilex are leading to a progressive substitution of the current dominant species in this forest by species of tall shrubs better adapted to low water availability (Liu et al., 2018), which could be followed by rapid alterations in ecosystem composition and structure and to the feedbacks between the biosphere and climate, as in many other forests around the world (Wang et al., 2012). Some ecosystem services such as the mitigation of climate change by CO 2 uptake will also be seriously affected.

Selective Thinning
Several studies have reported an increase in stem growth induced by selective thinning, increasing the availability of resources for the remaining stems (Roberts and Harrington, 2008;Keyser, 2010;Ogaya et al., 2019). Higher rates of net photosynthesis and stem growth in semi-arid habitats in response to thinning may be mediated by higher water availability and stomatal conductance (Moreno-Gutiérrez et al., 2011), whereas these responses were observed in temperate forests only in unusually dry years (Herbst et al., 2015). On the other hand, a decrease of tree mortality induced by a selective thinning is also described in several forests around the world, especially for those that are waterlimited (Hood et al., 2016;Neumann et al., 2017;Restaino et al., 2019). An increase in the rate of net photosynthesis can also contribute to mitigating climate change due to an increase in the uptake of atmospheric CO 2 by forests. In addition to the benefits of selective thinning for forest conservation and functioning (Wang et al., 2012), selective thinning has indirect benefits such as decreasing the risk of forest fires and increasing human water supply, because a high stem density can favor the spread of forest fire due to the abundance of dead stems and branches (Ogaya and Peñuelas, 2007), and highly dense vegetation consumes a lot of rainwater (Birot et al., 2011).

Remote Sensing
The NDVI and EVI indices were good predictors of forest dieback and stem mortality. Good relationships between these indices and increased biomass have been previously observed (Garbulsky et al., 2013). Both indices, however, have also strongly decreased under extremely dry conditions (Ogaya et al., 2016). These indices are widely used to estimate vegetation activity, but they can also be used to estimate forest decay or stem mortality. In contrast, annual GPP, PsnNet, and NPP in our study were poorly correlated with increased biomass or stem-mortality rates, and seasonal GPP, and PsnNet did not significantly decrease when many trees died in 2011 (Ogaya et al., 2016). This result was unexpected, because these indices are specifically calculated to directly estimate ecosystem productivity, but GPP, PsnNet, and NPP were not good predictors of either tree growth, or stem mortality in our forest.

CONCLUSIONS
Hotter and drier conditions are driving an increase in forest dieback and stem mortality in many water-limited ecosystems, such as Mediterranean forests, with a progressive substitution of the current dominant species by species of tall shrubs better adapted to future climate. These effects were followed by changes in forest structure and ecosystem services. Selective stem thinning was also identified as an excellent tool for managing this forest in a changing climate; selective thinning not only buffered forest structure and function against increasing drought, but also reduced the risk of forest fire.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
RO collected the data, performed the statistical analyses, and drafted the manuscript. DL and AB collected the data and contributed to the statistical analyses. JP designed the experiment and contributed to the statistical analyses. All authors contributed to the interpretation of the results and the writing of the manuscript.

ACKNOWLEDGMENTS
We are grateful to DARP (Generalitat de Catalunya), Ester Trullols, and Xavier Buqueras for their permission and help to conduct this research in the Poblet Holm Oak Forest.