Climate-driven variations in productivity reveal adaptive strategies in Iberian cork oak agroforestry systems

José Crlos Pérez-Girón ,Emilio Rfel Díz-Vrel ,Pedro Álvrez-Álvrez

aDepartment of Organisms and Systems Biology,Polytechnic School of Mieres,University of Oviedo,E-33600,Mieres,Asturias,Spain

bResearch Group on Planning and Management in Complex Adaptive Socio-Ecological Systems (COMPASSES).School of Engineering,University of Santiago de Compostela,E-27002,Lugo,Spain

Keywords:Net primary production Carbon use efficiency Climate Quercus suber Agroforestry system

ABSTRACT Background: Cork oak agroforestry systems (AFS) have been managed for centuries by humans to produce cork and other goods and services and have recently been recognised as an important reservoir for biodiversity improvement and conservation.However,despite having recently been included as a natural habitat of community-wide interest within the EU Habitats Directive,these systems are in a critical situation of decline.Among other factors,they are strongly threatened by climate change,the effects of which are also expected to be particularly severe in the Mediterranean region.In this study,we aimed to evaluate the influence of climate variability by examining primary production indicators and also to analyse whether the geographical location may have a role in the incidence of the adverse effects of climate.Methods: Cork oak AFS were identified in the Forest Map of Spain and the Land use map of Portugal and categorized on the basis of canopy cover.Seasonal climate data from 2001 to 2020 were used to model relationships with climate predictors and proximity to the coast.Hotspot analysis was conducted to identify significant spatial clusters of high-and low-efficiency areas.Results:The responses to the influence of climatic conditions differed among the various cork oak AFS categories,particularly in the forest category,which was less dependent on climate variations.Relative humidity and water availability were the main drivers of net primary production (NPP).Carbon use efficiency (CUE) was limited by relative humidity and spring temperature in open ecosystems.Proximity to the coast proved beneficial,especially in years with adverse weather conditions,but was not a limiting factor for survival of the ecosystem.Finally,the results of the hotspot analysis supported the other findings,highlighting high-efficiency areas close to the coast and cold spots grouped in specific areas or dispersed inland.Conclusions: Canopy plays a key role in the influence of climatic conditions,particularly in forest categories in which a high density seems to generate microclimate conditions.Water availability,both via the soil and air moisture,is the main driver of primary production,reflecting different adaptive strategies.The oceanic atmosphere may act as a buffer in years of extreme drought.

1.Introduction

Climate is considered one of the main drivers of biodiversity and ecosystem change (Millennium Ecosystem Assessment,2005;IPBES,2019) and is expected to become a major stressor (Bellard et al.,2012;Urban,2015;Willeit et al.,2019).The Mediterranean region has been considered one of the hotspots in future climate change (Giorgi,2006),and a pronounced increase in temperature(4°C-5°C)and a considerable decrease in mean precipitation(of around 25%-30%),mainly in summer,are expected (Giorgi and Lionello,2008),thus seriously increasing the effects of summer drought in the region.Droughts will occur more frequently,and be of greater duration and intensity,with the added aggravation in the Iberian Peninsula that they might not be followed by wet winters (Böhnisch et al.,2021).Furthermore,this trend has been confirmed by the latest IPCC report,which indicates that the Mediterranean region has been getting warmer and drier in the last few years(IPCC,2021).This may constitute a threat to the future of cork oak ecosystems,which are already in a critical situation of decline (Costa et al.,2009;Pinto-Correia et al.,2011) in the Mediterranean region,where water availability is the main driver of primary production(Reichstein et al.,2007;Garbulsky et al.,2010)and is a major ecosystem function sensitive to changes in climate(Huang et al.,2019;Tang et al.,2019;Stocker et al.,2019).

Cork oak (Quercus suberL.) is a typically evergreen Mediterranean tree species,with a range of distribution expanding around the western Mediterranean basin,where the largest populations of this species are found,specifically in the southwest of the Iberian Peninsula (Díaz-Fernández et al.,1995;Pereira,2007;Navarro Cerrillo et al.,2013).The species is well adapted to mild Mediterranean climates with Atlantic influence,i.e.mild winters and hot and dry summers with high relative humidity(Pereira,2007;Quero et al.,2008).The presence of the species in continental areas depends on there being some oceanic influence(Navarro Cerrillo et al.,2013).The optimum mean annual temperature ranges between 13°C and 19°C,and the species tolerates cold poorly and does not withstand periods of frost,especially below -5°C (Pereira,2007;Gil and Varela,2008;Navarro Cerrillo et al.,2013).The rainfall regime tolerated is within a wide range,mainly greater than 500 mm and reaching up to 2,400 mm (Gil and Varela,2008),and the species can withstand up to 4 months of summer drought due to its powerful root system.However,cork oak is very sensitive to waterlogging (Pereira,2007).To deal with summer drought,cork trees have developed deep root systems that enable groundwater uptake(David et al.,2007;Piayda et al.,2014)and physiological mechanisms that prevent water loss,such as efficient stomatal control of transpiration (Nardini et al.,1999;Mediavilla and Escudero,2004;Pérez et al.,2005;David et al.,2007;Besson et al.,2014).

Centuries-long management ofQ.suber,frequently associated with holm oak(Q.ilexL.andQ.rotundifoliaLam.)and to a lesser extent with other oaks (Q.fagineaLam.andQ.pyrenaicaWilld.),has given rise to multifunctional agroforestry systems (AFS) (Joffre et al.,1999;Costa et al.,2009;Pinto-Correia et al.,2011).In these systems,known asmontadosin Portugal anddehesasin Spain,cork production prevails together with the production of other goods and services,such as crops and other non-timber productions(mainly grazing and hunting)(Pereira,2007).For simplicity,we will refer to these systems hereafter asdehesas.The land cover patterns ofdehesasare similar to those of savanna,characterized by the presence of scattered trees in varying densities(Aronson et al.,2009;Fonseca and Pinto-Correia,2009;Pinto-Correia et al.,2011).However,the tree density is generally low,with the presence of herbaceous or shrub vegetation in the understory(Pereira et al.,2007;Correia et al.,2014,2016;Piayda et al.,2014).These systems have been key elements of the landscape since time immemorial (Eichhorn et al.,2006;Fonseca and Pinto-Correia,2009),making valuable contributions to the landscape and the environment,fulfilling fundamental functions and processes such as primary production,soil formation and regulation of nutrient cycles or hydrological flows,in addition to being an important reservoir for biodiversity improvement and conservation(Plieninger and Wlbrand,2001;Torralba et al.,2016).Dehesashave recently been included as a natural habitat type of community-wide interest within the EU Habitats Directive;in view of their conservation status,they have been categorized as in serious danger of disappearance.

Although the problem is not new,already having been reported in the mid-twentieth century,cork oak AFS are in a critical situation of decline(Costa et al.,2009),and their viability has been seriously threatened and affected by various factors(both natural and human-induced)in addition to climate change(Aguilera et al.,2020),such as the proliferation of pests and diseases(Brasier et al.,1993;González et al.,2020),fire recurrence(Silva and Catry,2006;Guiomar et al.,2015),lack of regeneration,change in land use,and land abandonment (Pinto-Correia and Mascarenhas,1999;Bugalho et al.,2011;Godinho et al.,2016).Thus,following a cascade reaction(Haines-Young and Potschin,2010),the loss of ecological functions and processes due to the impact of these factors(Schröter et al.,2005,2019;Mooney et al.,2009),may trigger a reduction in the capacity to provide ecosystem services,and consequently a risk to human well-being.

The current state of cork oak AFS and the effects of climate on the development of these ecosystems can be assessed using indicators.Primary production indicators,such as gross primary production(GPP)and net primary production (NPP),have been widely modelled in carbon cycle-related studies at the global scale as they are climate-sensitive(Huang et al.,2019;Tang et al.,2019;Stocker et al.,2019).GPP is the total amount of carbon stored by plants,which takes into account autotrophic respiration (Collalti and Prentice,2019;Collalti et al.,2020).Subtraction of autotrophic respiration gives us the net carbon transformed into biomass (NPP).The NPP/GPP ratio is a measure of carbon use efficiency (CUE),which represents the efficiency of plants to sequester carbon from the atmosphere through photosynthesis.In this type of assessment,the use of open source remote sensing data is very useful as it provides continuous,valuable information on ecosystem productivity over large areas.

The aims of the present study were (i) to evaluate the influence of climate variability on different cork oak AFS in the Iberian Peninsula(categorized by canopy cover) using production indicators (NPP and CUE)and(ii)to analyse whether geographical location may play a role in the incidence of the adverse effects of climate on these ecosystems.

2.Materials and methods

2.1.Study area and climate data

The study focused on mainland Portugal and Spain (Iberian Peninsula) and did not include the islands (Fig.1).The study region is surrounded by water,mainly by the Atlantic Ocean and the Mediterranean Sea,including 1,793 and 4,964 km of coastline of Portugal and Spain,respectively.The region also comprises 65% of the total distribution range ofQuercus suber(Caudullo et al.,2019).

The area is very heterogeneous in terms of climate and broadly speaking can be divided into three zones:dry climate zones(widespread in the south and southeast);temperate zones with dry,hot summers(most of the Iberian peninsula,i.e.approximately 40%of its surface);and temperate zones with dry,temperate summer climates (most of the northeast of the Peninsula,as well as almost all of the west coast of mainland Portugal)(AEMET,2011).

Climate data were obtained by combining ERA5-Land monthly averaged data from 1981 to present (Mu～noz Sabater,2019) and ERA5 monthly averaged data on pressure levels from 1979 to present(Hersbach et al.,2019).We download monthly mean variables from 2001 to 2020,for temperature,precipitation,relative humidity,total evaporation,potential evaporation and volume of soil water at different layers(layer 1:0-7 cm;layer 2:7-28 cm;layer 3:28-100 cm;layer 4:100-289 cm).All variables were provided at 0.1°× 0.1°spatial resolution (ca.9 km × 9 km pixel size),except relative humidity,which was distributed at 0.25°× 0.25°spatial resolution (ca.31 km pixel size).These data were summarized by season,i.e.for winter,spring,summer and autumn months.

Although proximity to the coast is not a climatic variableper se,it is a factor that affects climate conditions,with temperatures being higher or lower and coastal areas being wetter than inland areas.Thus,the shortest distance between eachQ.suberplot and the coast was computed.

2.2.Vegetation productivity data

The global MODIS data collection was obtained from the Land Processes Distributed Active Archive Center (LP DAAC) data pool.We used the MOD17A2HGF.006 and MOD17A3HGF.006 products(Running and Zhao,2019a,2019b),which provide gross primary production(GPP)and net primary production(NPP)data(in kg carbon per m2)respectively,from 2001 to 2020,at 500-m resolution.GPP and NPP values for non-vegetated or artificial areas were excluded from the analysis(Zhang et al.,2014),and the land pixel values were multiplied by a scale factor of 0.0001(Running and Zhao,2015),as ordered in the metadata file,to return the original value at those pixels.

Fig.1.Range of distribution of Q.suber in the Iberian Peninsula and the different types of study plots.

The GPP data set (originally one for every 8 days) was used to calculate annual mean values and carbon use efficiency (CUE)(Pérez-Girón et al.,2020).Annual CUE values were calculated as the NPP/GPP ratio,representing the efficiency of plants to sequester carbon from the atmosphere through photosynthesis.

2.3.Plot selection

The digital maps available for Portugal and Spain were different,due to the different purposes of each.For Portugal,we used theCarta de Uso e Ocupaç～ao do Solo de Portugal Continental2018(COS),developed with the aim of characterizing land cover in the country in 2018.In Spain,we used the Forest Map of Spain (MFE),at scale 1:25000 (MFE25) or 1:50000(MFE50),depending on availability,as the map was created with the information captured in national forest inventories carried out at different times:the MFE50 project was completed between 1997 and 2006,and the MF25 was developed after 2007.

We selected areas designated as“Florestas de sobreiro”y“Superfícies agroflorestais (SAF) de sobreiro”from the COS maps as those withQ.suberpresence.To harmonize the selection criteria,we reproduced the classification criteria used by COS in MFE,selecting patches classified as“Alcornocales”or“Dehesa”,predominated byQ.suberand with a canopy cover greater than 10%,the threshold established in the COS maps.

Finally,to ensure that GPP and NPP values corresponded toQ.suberAFS,we only selected the areas if more than 80% of the MODIS pixels were fully occupied by areas in whichQ.suberAFS was present.

2.4.Categorization

The social,economic and ecological importance of the cork oak has led to its management over centuries,giving rise to different tree densities depending on the aim of the AFS(Aronson et al.,2009;Fonseca and Pinto-Correia,2009).Tree density seems to control edaphic (Gallardo,2003) and climatic conditions in this type of ecosystem (Joffre et al.,1999).It has recently been demonstrated that forest canopy has a buffering effect on climatic conditions and their variations,even generating microclimates (De Frenne et al.,2021;Haesen et al.,2021).

Thus,in order to address the different densities in the AFS,we used the Tree Cover Density (TCD) 2018 from the high-resolution products provided by the Copernicus Land Monitoring Service (European Environment Agency) to assign to each previously selected MODIS pixel the average value of the TCD pixels on which the MODIS pixel overlaps.The selected pixels were then classified according to TCD into the following categories (Costa et al.,2006):low-density AFS (<10%),moderate-density AFS (10%-25%),high-density AFS (25%-50%) and forest stands (>50%).According to Costa et al.(2006) the average expected stand densities (in trees per hectare) are 25-35 (±18) for low-density AFS,36-42 (±17) for moderate-density AFS,48-55 (±18)for high-density AFS,while the stand density for forest category will be much higher.To check the relationship between TCD and stand density,data from the 3rdNational Forest Inventory (IFN3;from its acronym in Spanish)were used to assess whether the stand density(number of trees per hectare)increased with the TCD.Thus,based on TCD,the first three types (low-,moderate-and high-density AFS) can be considered open forest systems,with low tree densities,while the forest category more closely resembles a forest stand structure,with canopy closure expected to create a microclimate.

2.5.Statistical analysis

We used multiple linear regression to determine how cork oak AFS primary production was related to climate predictors and proximity to the coast.Hotspot analysis was used to identify significant spatial clusters of high-and low-efficiency areas.

Multiple linear regression was used to model the relationships between primary production indicators (NPP and CUE) and the climatic variables and proximity to the coast.An exhaustive search for the best predictor subsets was performed using the branch-and-bound algorithm,to test all possible combinations of predictors for final selection of the best model(Narendra and Fukunaga,1977).As this algorithm returns the best model of each size,we limited the number of predictors to 5.Candidate models were compared using the adjusted coefficient of determinationand the root mean square error (RMSE).To check multicollinearity,we computed the variance inflation factor (VIF) and excluded predictors yielding VIF>10(Mandeville,2008).

For proximity to the coast,we used a univariate linear regression approach and applied linear-log models.Spearman correlation coefficient was used to determine the strength of the linear relationship between proximity to the coast and climatic variables.A highRvalue indicates a stronger relationship,while a lowRvalue indicates the opposite.PositiveRvalues reveal the same trend,while lowRvalues reveal the opposite trend(Spearman,1904).

Hotspot analysis was conducted by applying the local Getis-Ord Gi＊statistic(Getis and Ord,1992)to identify high and low CUE areas(hotspots and cold spot areas)within the selected pixels.The analysis only included the CUE as this value establishes the threshold to the physiological activity of plants(Amthor,2000;Van Iersel,2003;Keith et al.,2010).Local spatial autocorrelation was measured to assess how each data point is surrounded by other data points (neighbourhood) with similarly high or low values.The method returns aZ-score andp-value for each data point assuming a normal distribution.Z-scores greater than 1.96 and less than -1.96 are statistically significant atp<0.05,while larger positive and negative values indicate greater clustering.

The leaps library (Lumley,2013) implemented in the R software environment (R Core Team,2020) was used to fit multiple linear regression based on branch-and-bound algorithms.Graphical analyses were conducted with the ggplot2 package(Wickham,2009).Significant differences were determined using the Wilcoxon-Mann-Whitney test (at α=0.01).

3.Results

3.1.Primary production indicator trends

The 870420 AFS plots were distributed as follows:11.6% lowdensity,28.9% moderate-density,44.4% high-density and 15.1% forest plots.Tree density differed significantly between the different groups based on forest inventories(Figure S1 in supplementary material),except between the low-density and moderate-density groups.A large difference in tree density between high-density plots(TCD between 25%and 50%)and forest plots(TCD>50%)was noted,and therefore large differences in the influence of climatic conditions were expected.

The trends in GPP and NPP in relation to TCD showed that the indicator value and dispersion increased with canopy cover (Fig.2).The mean NPP ranged from 0.655 kg C•m-2•yr-1for low-density plots to 0.985 kg C•m-2•yr-1for the forest plots,with maximum values of up to 2.21 kg C•m-2•yr-1reached.Inversely,the CUE value decreased as the canopy cover increased,with mean values ranging from 0.64 in lowdensity plots to 0.58 in forest plots and maximum values reaching up to 0.72 were reached in each of these types of plots;the minimum value was 0.30 in high-density and forest plots.

Regarding annual trends,the mean values were maintained(Fig.3a).In general,the dispersion range and average NPP increased with canopy cover,while CUE decreased as canopy cover increased.However,generalized slight decreases in NPP were observed,particularly in 2005,but were not related to CUE.The minimum mean NPP values of the historical series were recorded in 2005,reaching values of 0.474 kg C•m-2•yr-1in the low-density plots,0.544 kg C•m-2•yr-1in moderatedensity plots,0.611 kg C•m-2•yr-1in the high-density plots and 0.898 kg C•m-2•yr-1in the forest plots.This represents decreases of more than 20%for low-,moderate-and high-density plots,while the decrease in the forest plots was only 11%.The decrease in NPP values was immediately recovered in the following year with increases of the same magnitude as the decreases.The mean CUE values for this period suffered decreases of less than 2% that lasted 2 years,i.e.recovery of the values before the large decrease in NPP occurred in 2007 (not the largest decrease in the historical range).In 2012 and without any precedent,CUE values below 0.3 were observed for a total of 97 samples in the forest category.Similar findings were observed in 2015,2017 and 2020,in respectively 172,87 and 201 samples in the forest category.Rather than being an isolated event,this seems to have become a more frequent trend in recent years,first affecting the high-density plots in 2015.

3.2.Influence of climate factors

The parameter estimates and goodness of fit statistics for the multiple linear regression models based on canopy cover are summarized in Table 1.NPP showed a good fit in all categories,with the variance explained increasing with canopy cover(ranging from 0.64 to 0.72),except for forest category,which yielded the lowest variance explainedOn the other hand,CUE also showed a good fit in low-,moderate-,and high-density plots,with the variance explained increasing with canopy cover (from 0.46 to 0.57) as observed for NPP.The goodness of fit for CUE in forest category was very low,which may indicate that the CUE is influenced by factors other than climate factors.

Relative humidity in the summer months(rhm_Summer)and volume of soil water in the first two layers of soil(swl1 and swl2)in the summer months were the most important predictors for NPP models for low-,moderate-and high-density plots.The influence was similar in models:NPP benefited from an increase in relative humidity and volume of soil water of layer 2,but was negatively influenced by an increase in volume of soil water in layer 1.The most important predictors in the forest plots changed slightly.Relative humidity in the summer months was maintained,but the volume of soil water was replaced by total precipitation in summer (tp_Summer).In addition to changing the predictors,the influence also changed,with the lower rainfall during the summer months leading to higher NPP values.

Fig.2.Box-and-whisker plot showing the mean trends in GPP,NPP and CUE.The dots represent outliers.

Fig.3.Box-and-whisker plot showing the annual trends in GPP,NPP and CUE.The dots represent outliers.Dashed black lines indicate the threshold of ecosystems at potential risk (CUE=0.3).

Regarding CUE,both relative humidity in the summer months(rhm_Summer) and mean temperature of spring months (t2m_Spring)were constant throughout the models,and relative humidity had a positive influence,as occurred with NPP.The mean temperature of spring months negatively affected the CUE.The variables that changed were total evaporation in summer (e_Summer) in the low-density model and volume of soil water in layer 3 in spring(swl3_Spring)in the high-density model.

3.3.Proximity to the coast

The relationship between the primary production indicator(NPP and CUE) and the proximity to the coast decreased sharply for short distances,slowing down or becoming asymptotic as the distance increased.Therefore,a linear-log model was constructed:the goodness of fit statistics for NPP and CUE are shown in Fig.4.In general,theRa2djvalues for each indicator and group are similar to theRa2djvalues obtained for climatic factors.In the case of NPP,the outcomes were very clear,with the indicator values increasing as the distance to the coast decreased.A weaker relationship was obtained for CUE,although the trend continued to be observed.Finally,forest CUE was not related to the proximity to coast,apparently due to a large number of points near the coast with low CUE values,which may also explain the lack of any relationship between CUE and climatic factors.Annual models (Figures S2,S3,S4 and S5 in Supplementary material) followed the same trends shown here.However,dependency was higher in some years,i.e.2005 and 2017.

Proximity to the coast was highly and negatively correlated with spring and summer relative humidity (r=-0.67 and -0.81,respectively),negatively correlated with average temperature of autumn and winter (r=-0.52 and -0.54,respectively) and positively correlated with summer average temperature (r=0.48).The variable was not correlated or only slightly correlated with precipitation and soil water content.

3.4.Hotspot analysis

The Gi＊statistics revealed clusters ofQ.suberAFS with high and low CUE within the study area,respectively corresponding to hotspot and cold spot areas(Fig.5).In general,the largest hotspot areas occurred in the southwest of Portugal and south and northeast of Spain,always close to the coast,with large clusters in low-,moderate-and high-density plots,mainly in the Alentejo and Algarve regions (Portugal).The cold spot areas were generally located in inland areas,usually sparsely distributed,but forming large clusters in high-density and forest plots.Here we can distinguish three large clusters mainly located in Los Alcornocales natural park (south of Spain),Sierra de Hornachuelos and Sierra Norte de Sevilla natural parks (north of the Guadalquivir valley) and south of Caceres (central Spain).Forest cold spots were almost absent from AFS plots in Portugal.

Comparative box-and-whisker plots indicate significant differences between hotspot and cold spot CUE areas regarding the main independent variables selected in CUE multiple linear regression models(Fig.6).The trend shown is the same for the different types of AFS,and the differences between hotspots and cold spots were slightly smaller only in the forest category.Hotspots were identified in coastal zones while cold spots were found in inland areas (Fig.6a).The average distance to the coast was about 10 km for hotspots and more than 120 km for cold spots,except for the forest category,for which the distances were shorter,around 30 km.The dispersion was also remarkable as hotspots varied very little and the values were grouped very close to the mean,with the inverse observed for in cold spots.Related to the above,relative humidity during the summer months was higher in hotspots(Fig.6b),with average values around 68%-70%,while for cold spots it was around 47%-48%,except for the forest category,for which the values were slightly higher(61%).Finally,the trends in average spring temperature (Fig.6c) were very similar,regardless of the canopy cover,with values between 16.5°C and 18°C in the hotspots and between 1°C and 2°C in the cold spots.

4.Discussion

4.1.Changes with stand density

The multifunctional character of traditionally managed cork oak agroforestry systems has conditioned their structure and composition(Bugalho et al.,2011),mainly characterized by a low tree density with the presence of herbaceous or shrub vegetation in the understory.Consequently,the contribution of each part of the structure to primary production is different.In the dehesa ecosystem,understory species contribute between one third and one half of the total GPP,and consequently,about half or two-thirds of the GPP is contributed by trees(Dubbert et al.,2014;Correia et al.,2016).Therefore,in low-,moderateand high-density categories the contribution to biomass should be considered in relation to the ecosystem as a whole,rather than in relation to isolated species as e.g.in monospecificQ.suberstands.However,in the forest category,trees contribute more than half the value of the primary production indicator and thus establish a relationship with the species.

Fig.4.Linear-log models for NPP and CUE grouped by canopy cover categories.Note that“x”in the model formula indicates log(x).

Fig.5.Location of hotspot (red points) and cold spot (blue points) areas of Q.suber AFS based on CUE and grouped by canopy cover categories:a) low-density,b)moderate-density,c) high-density and d) forest.(For interpretation of the references to colour in this figure legend,the reader is referred to the Web version of this article.)

Fig.6.Box-and-whisker plot comparing a) proximity to the coast,b) summer relative humidity and c) spring temperature in hotspots and cold spot CUE areas.Statistical significance:ns:p>0.05;＊:p ≤0.05;＊＊:p ≤0.01;＊＊＊:p ≤0.001;＊＊＊＊:p ≤0.0001.The black dots represent outliers.The red dots represent mean values.(For interpretation of the references to colour in this figure legend,the reader is referred to the Web version of this article.)

The results obtained are sensitive to changes in canopy cover,highlighting different trends between the different categories as the canopy cover increases and especially with the forest category,where the effects may be attributed to the species.Both GPP and NPP values increased with canopy cover.The highest GPP and NPP values were associated with the forest category and related to a large number of trees and a greater contribution to the total amount of plant biomass of trees to the indicators than herbaceous or shrub vegetation.By contrast,the opposite effect was observed for influence of the canopy cover on the CUE,with average CUE value decreasing and the total number of low-value outliers increasing as canopy cover increased.These slightly low CUE values,relative to the categories with lower tree canopy density,are consistent with the stand development,as halts in growth are not followed by halted respiration,thus leading to a decrease in CUE (Collalti and Prentice,2019).

However,the increasingly common very low extreme CUE values in 2012,2015 and 2020 in high-density AFS and forest categories are of particular concern.These values are very close to the established threshold of 0.2,which may limit the physiological activity of plants(Amthor,2000;Van Iersel,2003;Keith et al.,2010) by increasing the additional respiratory cost,reducing growth and new tissue formation,or even leading to plant collapse,thus indicating that ecosystems will potentially be at risk (Pérez-Girón et al.,2020).However,such low values were not observed in the categories with lower canopy cover(lowand moderate-density categories),which may be related to the scarcity of trees.The understorey vegetation in cork oak AFS may mainly comprise annual species or crops,the growth period of which depends on the amount of rainfall and its seasonal distribution,or shrub vegetation well adapted to the Mediterranean summer drought stress (Correia et al.,2016).However,this is not possible with trees,as the weather conditions do not determine whether the trees are alive or dead,as trees can tolerate a range of conditions,with modified maintenance costs affecting carbon assimilation.This difference would therefore explain the absence of low CUE values in the low-and moderate-density categories.

The study findings also highlight the higher CUE of managed ecosystems than of unmanaged ecosystems (Fernández-Martínez et al.,2014;Campioli et al.,2015),suggesting that the conservation and future of these systems -as well as the biodiversity they maintain and the ecosystem services they provide -are strongly dependent on human management (Pereira and da Fonseca,2003;Bugalho et al.,2011).However,some changes appear to be taking place in these ecosystems as the same pattern was not observed in forest plots and high-density AFS within the historical period considered,indicating that a change in some factor has been accentuated in recent years.Here,different drivers such as economic(e.g.declining profitability of traditional dehesa products or agricultural intensity),socio-cultural (e.g.rural exodus),political (e.g.availability of access to direct payments of EU Common Agricultural Policy(CAP)),technological(e.g.new technics)and natural(e.g.climate change) factors may be affecting these systems in a complex and simultaneous way(Pinto-Correia et al.,2011;Plieninger et al.,2021).

4.2.Responses of NPP and CUE to climate influence

Previous studies have claimed that water is the main driver of primary production in cork oak AFS,mainly linking productivity to the limitations caused by summer droughts (Pereira et al.,2007;Piayda et al.,2014;Ramos et al.,2015;Correia et al.,2016).Our findings also indicate water limitations in summer as the main driver for NPP models in alldehesacategories,but to a greater extent in open ecosystems,i.e.categories with low canopy cover(<50%),than in systems with a forest stand type structure.This suggests that summer is a critical season for cork oak AFS,regardless of the canopy cover.Furthermore,although water availability was also a limiting factor in the CUE models in our study,fewer limitations were observed than in the NPP models,as other factors were involved.

The productivity of any plant depends on its ability to maintain photosynthetic tissues with an adequate water level.When water is limiting,water loss in plants is minimized via transpiration,with closure of stomata restricting the entry of CO2and thus also limiting photosynthesis (Mediavilla and Escudero,2004;Rzigui et al.,2018;David-Schwartz et al.,2019;Grossiord et al.,2020).Although this is a typical adaptive mechanism in Mediterranean species,such asQ.suberandQ.ilex,protecting against summer drought,trees generally display stomatal control over transpiration (Nardini et al.,1999;Pérez et al.,2005;David et al.,2007;Besson et al.,2014);drought also induces leaf senescence in the understorey vegetation,and all of these factors lead to a decrease in primary production levels in the ecosystem(Pereira et al.,2007).Relative humidity plays a key role in this process(closely related to vapor pressure deficit:Grossiord et al.,2020),because the transpiration rate falls as relative humidity increases.This explains why relative humidity in the summer months is the most important variable in all of our models and why it is positively related to NPP.The relationship with CUE is clear,as the maintenance cost increases with the transpiration rate(Amthor,2000;Van Iersel,2003),and the maintenance cost will therefore be lower at high relative humidity,and the CUE will increase.

When soil water is available to plants,transpiration rates are mainly controlled by climatic factors.However,when the soil water becomes limiting,the transpiration rate falls (Gardner and Ehlig,1963).Thus,in the NPP models,in addition to being affected by relative humidity,the low-,moderate-and high-density categories were also positively affected by water volume in the second layer of soil(swl2;between 7 and 28 cm)and negatively by soil water volume in the first 7 cm (swl1).In open ecosystems such as dehesas,the first layer of soil is particularly sensitive to water loss through evaporation,which will form part of the relative humidity of the air;thus,favouring an increase in evaporation would benefit the ecosystem to a greater extent.From the point of view of water absorption by plants,water is more beneficial in deeper layers,where the roots of most trees and of the understory vegetation occur(David et al.,2007;Baldocchi et al.,2010;Correia et al.,2016).This explains the positive relationship with the presence of water in the second layer of soil.Although greater limitation of the amount of water in the third layer of soil was expected due to the location of most roots(from 28 cm to 1 m),we believe that these findings may be explained by daily fluctuations in soil water content due to hydraulic lift processes(David et al.,2007).The water rises to the upper soil layers where it becomes available both to the oak tree roots and understorey vegetation(Mara～nón et al.,2009).

Temperature only appears to influence CUE,especially the average temperature in spring.Temperature plays a fundamental role because it affects both photosynthesis and autotrophic respiration(Ra),the rates of which increase exponentially with temperature,thus increasing the maintenance cost(Ryan,1991;Ryan et al.,1994).As this only occurs in low-,moderate-and high-density plots,and asQ.suberis adapted to Mediterranean climates characterized by high temperatures,this pattern is associated with the vegetative period of the understorey vegetation(Dubbert et al.,2014;Correia et al.,2016).Higher temperatures in early spring may accelerate the activation and germination processes in annual plants,favouring greater photosynthetic activity and therefore greater CUE.However,higher temperatures in late spring are also likely to have the opposite effect,shortening the duration of photosynthetic activity by advancing senescence.Furthermore,temperature affects both the plant part of the ecosystem and also other climatic factors such as relative humidity and soil water content.

Finally,the findings suggest structural and compositional differences between low-,moderate-and high-density AFS and the forest plots.Under forest canopy cover,the amount of photosynthetically active radiation and wind speed are reduced,directly reducing the variations in humidity and temperature and extreme events (De Frenne et al.,2021;Haesen et al.,2021).The effects of canopy throughfall on the soil characteristics are also modified by the contribution of organic matter(Mara～nón et al.,2009).All of this is translated into differences in primary production and the associated factors.Thus,it was only possible to model the NPP response,and the goodness of fit was lower than that of the previous model;it was not possible to model the CUE,possibly because this variable does not depend on the predictors used.Similarly,the unexpected and negative effect of summer rainfall on NPP may be related to the fact that the basic needs of humidity and water are covered in these stands.Heavy storms occur at the end of the summer,which may cause waterlogging that is damaging to the ecosystem (Pereira,2007).The storms can also provoke fluctuations ranging from flooding to water deficiency,the latter of which is favoured in periods of high temperatures and can potentiate the spread and infectivity ofPhytophthoraspecies(González et al.,2020).

4.3.Oceanic influence as a buffer for extreme droughts

The study findings highlighted the influence of oceanic climatic conditions on the NPP of cork oak AFS,with proximity to the coast providing a clear benefit.Similarly,the decrease in NPP with increasing distance from the coast does not seem to be limiting for these ecosystems but rather is asymptotic,either because the tree species is within its distribution range (Gil and Varela,2008;Caudullo et al.,2019) or because it is adapted to the driest climatic conditions(David et al.,2007;Besson et al.,2014).

The study findings reflect a gradient of change along the dehesa categories.In the low-density category,the highest NPP values,mainly contributed by the understory(Dubbert et al.,2014;Correia et al.,2016),are lower than for the category immediately above the same distance from the coast,up to the highest values found in the forest category.The forest category,coinciding with the ecological response ofQ.suber,is benefited by a mild oceanic climate with high relative humidity(between 65% and 80%) (Quero et al.,2008),which is also consistent with the observed differences between hotspots and cold spots(Fig.6).Therefore,under optimal climate conditions a greater number of trees will be more productive.In addition,the high and negative correlation between the distance to the coast and the summer relative humidity again suggests that relative humidity in the summer months is a limiting factor for all cork oak AFS as it is closely linked to stomatal closure (David-Schwartz et al.,2019).A similar but slightly weaker pattern was observed for CUE,the levels of which varied with the proximity to the coast within ranges of distance that generally do not affect the physiological functioning of the ecosystems.In this last case,the microclimate generated by the forest(De Frenne et al.,2021;Haesen et al.,2021) seems to be sufficient to maintain the conditions necessary for adequate carbon assimilation,without depending on the coastal influence.

The findings of the hotspot analysis were consistent with these previous results,highlighting high-efficiency areas near the coastline and cold spots grouped in specific areas or dispersed inland.However,considering the CUE value,some clusters of points suggest some degree of risk to the ecosystem(Pérez-Girón et al.,2020).These CUE values are very low(generally below 0.3)both in high-density and forest categories(Fig.4)and are particularly notable in the latter at short distances to the coastline,in the south of Spain(Fig.5),and considering the annual distribution was more pronounced in 2012,2015,2016,2017 and 2020(Fig.3).According to the historic drought database for Spain (Vicente-Serrano et al.,2017),droughts were recorded in 2005,2012,2015,2017 and 2019,with those occurring in 2005 and 2017 being the most extreme.For these two years in particular,the goodness-of-fit of our linear-log models increased significantly (Figures S2,S3,S4 and S5 in Supplementary material),which suggests that in extreme drought years,proximity to the coast may buffer cork oak AFS from extreme climate conditions.This further confirms the suggestions of Piayda et al.(2014),who argued that the high vapor pressure deficit(related to temperature and relative humidity (Grossiord et al.,2020) found in cork oak ecosystems in Portugal may be a consequence of the lack of entry of oceanic air masses.

5.Conclusions

The responses of the different TCD-based categories of cork oak AFS to primary production indicators suggest that canopy cover and hence tree density play a key role in the influence of climatic conditions.Forest plots can maintain microclimatic conditions that make them less dependent on environmental conditions,while AFS plots with an open ecosystem (lower densities) depend on macroclimate conditions.Therefore,within the same ecosystem,the response to climate change may vary depending on tree density.

Regarding the influence of climate variability,our findings showed that the responses of the ecosystems reflect the ecological traits and the different adaptive strategies used by the component trees and understory plants to survive drought seasons,where water(soil or air moisture)is the main driver of primary production.Relative humidity is associated with transpiration and water loss through closure of stomata,which will vary depending on the severity of water deficiency,limiting photosynthesis to a greater or lesser extent.At the same time,both the hydraulic lifting processes and the deep roots allow the trees to take advantage of the groundwater from the deeper layers and make it available to understory vegetation with shallow,but not superficial,roots.Temperature only seems to influence CUE in open ecosystems (low-,moderate-,and highdensity) in which the understory layer makes a greater contribution to primary production.In particular,an increase in spring temperature could advance the growing season but could also shorten the growth period of annual plants and increase the maintenance cost.

Several factors affect these ecosystems and do so in a complex way,and it is therefore difficult to isolate the individual effects.For example,relative humidity and the proximity to the coast are closely related,thus influencing the carbon balance in cork oak AFS.Our findings show that proximity to the coast improves productivity levels and may also buffer climate conditions in extreme drought years,reducing the associated adverse effects and the risk to the ecosystem.Therefore,in future climate change scenarios,in which the Mediterranean region is expected to be one of the most severely affected and droughts are expected to be more frequent,prolonged,and intense,an important risk of loss of ecosystems and their associated functions will appear.This will affect all ecosystems,although inland ecosystems -where the first disturbances have already been detected -may not be buffered against the oceanic influence and would therefore be particularly affected,while coastal areas,such as southwestern Portugal,may cope better.

Ethics approval and consent to participate

We have no ethical concerns to declare.All authors consented to participate on this manuscript.

Consent for publication

All authors agree with the content of this manuscript and its publication inForest Ecosystems.

Availability of data and material

The datasets analysed during the current study are freely available from the following hosts.ERA5-Land monthly averaged data from 1981 to present and ERA5 monthly averaged data on pressure levels from 1979 to present can be accessed from the Copernicus Climate Change Service(C3S) Climate Data Store (CDS) (https://doi.org/10.24381/cds.68d2bb30 and https://doi.org/10.24381/cds.6860a573).MODIS MOD17A2HGF.006 and MOD17A3HGF.006 products can be accessed from the Land Processes Distributed Active Archive Center (LP DAAC)data pool (https://lpdaac.usgs.gov/products/mod17a2hgfv006/and https://lpdaac.usgs.gov/products/mod17a3hgfv006/).Forest Map of Spain can be accessed from the Ministerio para la Transición Ecológica y el Reto Demográfico webpage (https://www.miteco.gob.es/es/carto grafia-y-sig/ide/descargas/biodiversidad/mfe.aspx).Carta de Uso e Ocupaç～ao do Solo de Portugal Continental 2018 can be accessed from the Sistema Nacional de Informaç～ao Geográfica Of Direç～ao-Geral do Território(https://snig.dgterritorio.gov.pt/rndg/srv/por/catalog.search).

Funding

JCPG is in receipt of a“Severo Ochoa”PhD Grant provided by the Government of Principado de Asturias(PA-18-PF-BP17-026).

Authors’contributions

JCPG,PAA,and EDV together conceived the research idea;JCPG collected and analysed data and wrote the initial draft of the manuscript;all authors contributed to the manuscript review and editing.All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Acknowledgements

Not applicable.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://do i.org/10.1016/j.fecs.2022.100008.