Introduction 

Worldwide forest ecosystems are characterized by their high biodiversity, harbouring a large number of plant and animal species, which represents an important percentage of the terrestrial biomass ([20]). In Europe, coniferous (Picea abies (L.) H. Karst, Pinus spp., etc.) and broadleaved forests (e.g., Fagus sylvatica L., Quercus robur L., Q. petraea (Matt.) Liebl., Q. canariensis Willd.) are being affected by forest decline episodes since late 20^th century ([8], [42]). Although the first cases were related to biotic (fungi, insects and other pathogens) and abiotic disturbances (bushfires, storms and long drought periods) that caused severe impacts into ecosystem dynamics ([44]) the discussion about forest dieback cannot be nowadays separated from the risks associated with climate change and past management ([45]). Since the Mediterranean region is located in a transition zone between the mild and humid climate of central Europe and the semiarid climate of northern Africa, it has been highlighted as a potentially vulnerable area to climatic shifts ([17]). Specifically, different scenarios related to the climate change impacts on the Iberian Peninsula suggest an important drop in annual rainfall, that oscillates between 18% and 42% depending on the climate scenario ([3]).

The frequency and intensity of these drought episodes, together with other abiotic stresses such as heat, are expected to increase under climate change. In addition to specific differences in the molecular mechanisms underlying plant responses to drought and heat, it is known that the combination of both stresses causes unique physiological and anatomical plant responses, basically relying on changes in the reactive oxygen species (ROS) metabolism and stomatal behavior ([49]). Water scarcity effects can cause delayed forest physiological recovery (“legacy effects” sensu [2]) and escalate when combined with heat stress, causing extensive tree mortality episodes ([1]). It has also been highlighted the interactive effect of drought with other disturbances like bushfire ([4]) or insects and pathogens ([27], [21]). These interactions could cause important alterations in the ecology of plants at different scales, from populations’ demography and species distribution ([30]) or community diversity ([6]) to intraspecific phenotypic responses ([35]). Therefore, the complex relationship between the different stress factors involved in forest dieback ([11], [42]), intensified in a context of climate change, introduces an enormous instability in forest dynamics ([25]) and finally results in a complex scenario not completely understood to date.

In the southern Iberian Peninsula, the forest is shaped mostly with species of the Quercus genus (Fagaceae), being the most widespread the evergreen Holm oak (Q. ilex L.) and cork oak (Q. suber L.) as well as the semi-deciduous Q. faginea Lam. and Q. canariensis Willd., which constitute the typical open woodland agroforestry systems handled by human and cattle for thousands of years ([37]). Climate-change driven distribution shifts are predicted for these forests, but other factors such as intra- and inter-specific competition or tree size are shaping their dynamics ([18]) along with soil fertility and the availability of nutrients. In fact, Q. suber forests with high economic and ecological relevance are already suffering a worrying woodland decrease ([16], [31]) according to similar processes in different global biomes.

During the last decades, several factors other than climate have been identified to threaten the persistence of these oak forests. Among them, the presence of the exotic pathogen Phytophthora cinnamomi Rands., a soil-borne oomycete inducing the root-rot disease ([5], [19]), causes the mortality of a large number of adult trees. Other important factor hindering natural regeneration is the intensive grazing pressure produced by wild and domestic ungulates ([41]). This overgrazing is particularly intense in private lands, where cattle and hunting provides important economic benefits and are actively managed. However, it is still not clear how these management practices will interact with climate and which could be the final consequences for future forest dynamics. In the case of Q. suber, the decline process is specially alarming due to the elevated economic benefits that it provides in the region. Previous research has indicated that climate factors (especially low rainfall and dry spell duration) and past management strategies affect cork oak growth and its demographic dynamics ([16]). The synergies between climate change, lack of regeneration, forest ageing and exotic pathogens have caused a decrease in cork production during last decades, affecting this socioeconomic resource at the southwestern Europe ([26], [50]) and the ecosystem services provided by cork oak woodlands ([29]).

Here, we present the results of a field study on the cork oak forests from the southern Iberian Peninsula aimed to determine how site conditions and management factors affect cork oak forest structure. The objective of this study was twofold: (i) to assess the relationships between site conditions (latitude, elevation and aspect), management (public or private) and local climate with the mortality ratios of cork oak woodlands; and (ii) to disentangle the possible changes in Q. suber demography through studying tree size, density, regeneration and mortality ratios. We hypothesized that private intensively-managed forests show lower tree density, regeneration and mortality than public forests and that southern forests have higher tree size than those of the north due to better climate conditions and past management legacies. These results would allow us to delve into the dynamics of these communities and to shed some light on establishing preventive measures against climate change scenarios as well as to settle optimal management models that allow the conservation of cork oak forest ecosystems.

  Material and methods 

Study area description

The study was carried out at Los Alcornocales Natural Park, located in southern Spain (36° 21′ N, 05° 32′ W - Fig. S1 in Supplementary material). The elevation of this area ranges from the sea level to 1092 m a.s.l. The soil is constituted by sandstone from Oligo-Miocene, with alternate layers of clay materials that constitute steep reliefs in spite of the low altitude. The inceptisols form the most significant edaphic unit of the territory, presenting sandy and permeable soils. In addition to their good drainage, the acidity and the calcium carbonate deficiency is notable. The climate is sub-humid Mediterranean with dry summers and wet autumns and winters. Mean annual temperature is 15.7 °C, being July the warmest month (mean of 24.0 °C). Dominant warm winds arrive from northern Africa (Sahara Desert), taking humidity from the Mediterranean Sea and creating a characteristic fog (crypto precipitation) mainly in the southern mountainous areas, which receive an extra amount of rainfall not recorded by the weather stations, thus generating a humid microclimate that softens the intense summer drought. For the period 1971-2000, the rainfall in the study area oscillated between 600-1400 mm year^-1, with a maximum close to 2000 mm year^-1 in certain high elevations. On average, summer rainfall was 35 mm at the reference weather station for the northern location and 55 mm for the southern one ([26]).

Los Alcornocales Natural Park is dominated by cork oak forests (Q. suber), one of the largest and the best preserved woodlands in Spain ([33]), with 62% of the surface under private ownership. The study area is characterized by a heterogeneous landscape with abundant dense stands (> 300 trees ha^-1) and open woodlands (Dehesa or Montado). At the highest and more humid areas, Q. suber is replaced by Q. canariensis. In the bare zones, wild olives (Olea europea L.) appear above marlstone soils, forming grass areas for cattle and goats. The shrub stratum is dense and characterized by several woody species such as Erica arborea L., E. australis L., E. scoparia L., Phillyrea latifolia L., Genista tridens Cav., G. triacanthos Brot., Calicotome villosa (Poir.) Link, Crataegus monogyna Jacq., Ulex borgiae Rivas Mart., Stauracanthus boivinii (Webb) Samp. and Cistus salviifolius L., among others ([33]). The southern location of this study was the starting area for the exploitation of cork in Andalusia and therefore an area of intense cork activity, which affected the very survival of the trees. In 1865, one of the world’s most important cork factories was installed there, the Armstrong Cork Company, with headquarters in Pittsburgh, PA, USA ([26]).

Sampling design

Field sampling was carried out from January 2019 to February 2020. We established 20 plots of 100 × 100 m (1 ha) across forests dominated by Q. suber, following a north-south dipole and a health status gradient within the Natural Park (Fig. S1a in Supplementary material). The plots were separated by at least 1 km and they showed certain common characteristics: lack of recent bushfires (> 10 years) and reforestation. Thus, 10 plots were sampled in the northern part of the Natural Park limit and 10 in the southern distribution (Fig. S1a). We compared 5 private and 5 public plots within each region, the latter owned by the Forest Administration and local governments.

Each plot was geo-referenced by GPS (GPSMAP 66st^®, Garmin Ltd., USA), and tree density was estimated by the total counting of Q. suber individuals and other tree species. Dead cork oak individuals were also counted, as well as the number of trees with severe defoliation (61-99% of crown surface defoliated), according to the semiquantitative standard scale proposed by ICP guidelines ([14]). Mortality ratio was estimated as the number of dead Q. suber individuals with respect to the total number of living Q. suber trees within the plot. Diameter at breast height (DBH) was measured in at least 30 trees of different diametric classes per plot. DBH of each individual was determined using girthing tape and basal area was calculated based on the formula for the area of a circle [basal area = π × (DBH/2)²]. Moreover, we counted the total number of saplings (height < 50 cm) in order to estimate the regeneration rate in each plot. Tree and shrub covers were visually estimated by two independent observers for each plot using a semi-quantitative scale (0-25%, level 1; 25-50%, level 2; 50-75%, level 3; and >75%, level 4 - [12]).

Local climate data and statistical analyses

In order to evaluate the influence of the latitudinal limits (N vs. S) on the cork oak population dynamics, we obtained total annual rainfall and mean temperature data from nearby local meteorological stations (< 10 km to the study sites - Fig. S1a) for the period 1985-2019: Alcalá de los Gazules in the north (36° 27′ N, 05° 43′ W) and Los Barrios in the south (36° 11′ N, 05° 29′ W). These data were compared by Student’s t statistic (significant differences at P < 0.05) to check if there have been climatic differences between the two groups of plots over the last 3 decades.

We analyzed the relationships of geographic location (N vs. S) and property type (public, hereafter Pu, or private, hereafter Pr) as explanatory variables, with mortality, DBH, basal area, regeneration, tree density as well as tree and shrub covers as response variables. For that purpose, we used a Gaussian function in Generalized Linear Models (GLM), except for the mortality, where we employed a logistic regression since this response variable has a binomial nature (i.e., alive or dead). Likelihood-ratio chi-squared tests were run for every model and the statistic was calculated (L-R χ²). In case of significant (P < 0.05) interactive effects between plot location and property, pairwise differences were detected by post-hoc Tukey’s HSD tests (α = 0.05). If only a significant effect of the geographic location was detected (P < 0.05), we transformed this discrete variable into a continuous one in order to obtain a more detailed view. Thus, we represented the specific latitude of each plot versus the response variable and calculated the coefficient of determination (R²). All statistical analyses were performed in R software ver. 4.0.0 ([38]).

  Results and discussion 

Mortality

Q. suber mortality ratio was significantly lower (L-R χ² = 50.77, df = 1, P < 0.05) in northern plots (8.95%) than in southern ones (15.6%). Moreover, the property influenced the mortality ratios (13.89% in the public vs. 6.86% in private ones; L-R χ² = 19.61, df = 1, P < 0.05 - Fig. 1). Lastly, both factors showed interaction (L-R χ² = 44.94, df = 1, P < 0.05), being the mortality similar in public plots (13.66% in the N vs. 14.39% in the S) but different to the private ones (4.28% in the N vs. 18.08% in the S - Fig. 1). In addition, we detected the presence of dead or very damaged Q. canariensis trees in the 15% of the sampled plots which, as already reported by other studies ([42]), suggests that this semi-deciduous species inhabiting humid habitats and shallow and fertile soils is also susceptible to tree dieback.

Fig. 1 - Mortality ratios (%) of Q. suber individuals in the study plots regarding property type (private or public) and geographic location (North: red; South: blue). Horizontal lines represent the median, and boxes and whiskers represent the interquartile range and the nonoutlier ranges of the box-plots, respectively. Significant effects of geographic location (GL), property type (PT) and/or their interaction are indicated in the corner of the panel (GLM: *P < 0.05). Different letters indicate significant differences between groups (Tukey’s HSD test, α= 0.05).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

It is known that competition between individuals of a population for environmental resources might cause changes in their subsistence under limited resources. Different studies confirmed that competition, tree size and climate are the most important factors implied on oak mortality processes ([31]) and, in some cases, also management intensity and local characteristics ([46]). Our results suggest that mortality processes are strongly influenced by the intensity of the centennial management of Q. suber forests in private properties from southern locations, resulting in higher disturbances due to resource extraction (cork, tannins, coal and wood). In addition, the public sites tended to present higher tree densities (see below), a factor that has been previously related to tree mortality in this species ([31]).

The cork and tannins extraction causes a strong stress on the tree, occasionally leading to death after several harvesting rotations ([34]). Obtaining tannins for the tanning of skins sometimes caused the death of many cork oaks ([9]). In the Forestry Memories of the end of the 19^th century, it was not unusual to read that some of the cork removals carried out in the forest caused serious damage to the trees. Although both northern and southern areas have had a long historical tradition in cork exploitation, a plausible explanation for the higher mortality of Q. suber in the S would be the more intense cork extraction during the last century ([50]). The southern region was the starting point of the cork exploitation in Andalusia and an intense area of cork activity, which affected the survival of remaining trees ([23]).

Regarding the influence of climatic factors on mortality, the increase in mean annual temperature might induce physiological stress of trees due to increasing evapotranspiration rates ([32], [42]). In addition, higher temperatures might enhance water stress and modify the response of trees against certain disturbances ([7]). Our data confirm that the mean annual temperature was significantly higher at the southern location along 1985-2019 period (18.53 ± 0.07 °C in the S vs. 18.01 ± 0.08 °C in the N, P < 0.05 - Fig. S2 in Supplementary material), which would had contributed to the presence of physiological imbalances and would partly explain the higher mortality in the S location (Fig. 1). Annual rainfall was significantly higher (P < 0.05) in the S between 1985 and 2019 (885 ± 64 mm) than in the N location (705 ± 45 mm - Fig. S2). These contrasting conditions in southern locations indicate a more humid and favourable conditions for the propagation of soil pathogens ([11], [21]), which has probably been undermining the health of cork oak trees ([19]), thus causing an increase in the mortality at the S plots. Lastly, the widespread theory according to which plots facing south (more sun-exposed) would show higher mortality rates could not be confirmed. In our study, plot aspect had not a significant effect over mortality ratio (R² = 0.02, P > 0.05 - Tab. 1).

Tree density

Tree density was significantly related to geographic location (L-R χ² = 41.11, df = 1, P < 0.05), with higher mean values in the N (307.2 ± 16.5 trees ha^-1) than in the S (106.3 ± 27.7 trees ha^-1). This result was supported by the high adjustment of the lineal regression of tree density in relation to the latitude of plots (R² = 0.69, P < 0.05 - Fig. 2).

Fig. 2 - Lineal regression showing tree density variations across the sampled latitudinal gradient in Quercus suber study plots (R² = 0.69, P < 0.05). Red dots show sampled plots in the northern location and blue dots in the southern one; crossed and emptied show private and public properties, respectively.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

On one hand, the long tradition of exploitation of resources from Q. suber (i.e., cork and tannin extraction) has modified the natural composition of the forests and has boosted the Q. suber growing at the expense of other species ([26], [47]). The management carried out for the last two centuries has promoted the increase of Q. suber individuals because their presence has been facilitated at the expense of other species in order to obtain greater amounts of cork, leading to forests with higher cork oak density. However, the high intensity of cork harvesting from the end of the 19^th century to the first half of the 20^th century led to numerous mortality episodes ([9], [26]), which would partially explain the lower tree density in the southern location after the more intense harvesting practices in this area.

This legacy of forest structure determines current competition levels and constitutes a crucial factor in demographic processes, like mortality, of many Iberian tree species ([18]). This legacy effect has been also proven to be one of the main factors affecting mortality across the distribution range of Q. suber since intra- and inter-specific competition has modified the ecological conditions of these populations ([31]).

Size structure

Trees differed in size between the two studied locations (DBH = 37.6 ± 0.8 cm in the N vs. 48.2 ± 1.0 cm in the S; L-R χ² = 75.59, df = 1, P < 0.05) and property types (DBH = 37.8 ± 0.7 cm vs. 47.3 ± 1.1 cm in public and private properties, respectively; L-R χ² = 61.76, df = 1, P < 0.05 - Tab. 1), showing a significant interaction between both factors (L-R χ² = 9.84, df = 1, P < 0.05 - Fig. 3a). Thus, differences in DBH within southern plots were higher (41.4 ± 1.0 cm vs. 55.1 ± 1.7 cm in public and private plots, respectively) than within northern plots (34.6 ± 0.8 cm in public plots vs. 40.6 ± 1.3 cm in private plots - Tab. 1).

Fig. 3 - Size structure parameters in Quercus suber study plots. (a) Tree size (DBH) regarding property (public or private) and geographic location (N or S). Values represent means ± standard error. Significant effects of geographic location (GL), property type (PT) and/or their interaction are indicated in the corner of the panel (GLM: *P < 0.05). Different letters indicate significant differences between groups (Tukey’s HSD test: α = 0.05). (b) Lineal regression of basal area across the sampled latitudinal gradient (R² = 0.40; P < 0.05). Red dots show sampled plots in the northern location and blue dots in the southern one; crossed and emptied show private and public properties, respectively.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Forest management has traditionally favoured the presence of larger and thicker Q. suber trees (more productive for the cork industry). Therefore, strong selective pressures for resource extraction have been exerted on cork oak individuals during centuries (see the previous “Tree density” subsection), especially in private plots. In these areas, owners aim to maximize benefits and shorten time between cork extraction events and would have indirectly promoted the increase in trunk thickness. In addition, the higher DBH in southern individuals could be related to microclimate patterns of the Strait of Gibraltar area (mild temperatures due to the proximity to the ocean and high air humidity), that would buffer summer aridity and reduce plant water stress, therefore allowing higher growth rates ([42]). The interaction between microclimate and management would explain the boost of Q. suber DBH under the favourable southern conditions in the more intensively managed private plots, implying a greater DBH variation between property types in the S location.

Basal area was higher in the N than in the S (38.4 ± 4.7 m² ha^-1 in the N vs. 18.5 ± 3.7 m² ha^-1 in the S; L-R χ² = 12.58, df = 1, P < 0.05), in accordance with the higher tree density recorded and higher mortality rates (Fig. 2). Moreover, our results showed a significant linear relationship between the basal area and latitude (R² = 0.40, P < 0.05 - Fig. 3b), with a gradient towards larger basal areas in northern plots. The higher basal area has been also related to an increasing acorn input, hence higher regeneration ([36]).

Regeneration

Regeneration (number of Q. suber saplings per hectare) was significantly affected by the geographic location of the plots (L-R χ² = 9.21, df = 1, P < 0.05), with higher values in the N than in the S (94.9 ± 25.2 vs. 26.0 ± 6.1 saplings ha^-1, respectively - Fig. 4a). In addition, we found a significant interaction between geographic location and property type (L-R χ² = 4.41, df = 1, P < 0.05 - Fig. 4a), being this pattern opposite to that detected for mortality rate and DBH. In this case, significant differences between N and S were found in public plots (Tukey’s HSD tests, α = 0.05; 138.4 ± 38.5 saplings ha^-1 in the N vs. 21.8 ± 11.7 saplings ha^-1 in the S), whereas private ones showed similar values (51.4 ± 20.5 saplings ha^-1 in the N vs. 30.2 ± 4.8 saplings ha^-1 in the S).

Fig. 4 - Regeneration (panel a; Q. suber saplings ha^-1) and shrub cover (panel b) differences between property types (private or public) and locations (N or S). Values represent mean ± standard error. Significant effects of geographic location (GL), property type (PT) and/or their interaction are indicated in the corner of the panel (GLM: *P < 0.05). Different letters indicate significant differences between groups (Tukey’s HSD test: α = 0.05).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The higher regeneration rate at the denser northern forests, compared to the S location, would be caused by the higher proportion of healthy trees per hectare, that produces more fruits (acorns). Although forests in the N location have been exposed to drier conditions during the last 35 years (9 years with a total annual rainfall below 489 mm - Fig. S2b in Supplementary material), which is considered the minimum range under which Q. suber is distributed in Andalusia according to [47]), the detrimental effects of waterlogging for Q. suber sapling survival would explain this difference. Cork oak seedlings have shallower root systems than mature trees and may be affected by excess water in upper soil layers ([28]), as reflected in our results. The loss of trees and regeneration problems of Q. suber forests under a changing climate are altering the species’ population dynamics across its complete distribution range ([31]), a problem that is especially evident at the southern Iberian Peninsula.

The high total annual precipitation (with intense rainfall amount in narrow time windows) is causing soil erosion processes and the loss of soil fertility ([24]), which indirectly contributes to the loss of seeds of different understory species ([15]) and to the decrease in the growth rate of Q. suber saplings in particular. Even though the regeneration rate between public and private plots is similar in the S (and the lowest sapling values registered), this was different in the N, where regeneration is lower in private areas. We attribute it to a slightly higher herbivory pressure, since there is usually an excess of herbivorous populations on private farms due to the high interest of maintaining this economic resource. In the Natural Park, deer (Cervus elaphus L.), fallow deer (Dama dama L.), wild boar (Sus scrofa L.) and other wild ungulates are causing an intense grazing of many shrub species including cork oak saplings, altering the composition of the shrub understory and endangering the regenerative capacity of Q. suber ([41]).

Moreover, it has been verified that Q. suber latitudinal distribution entails demographic changes in the structure of populations ([31]). Those of southern Spain present a smaller density of saplings, oppositely to the populations of the northern distribution edge of the species ([31]). Our results partially confirmed the importance of the geographic gradient (latitudinal) at regional scale in the sapling recruitment. However, we assume that there are no genetic differences between the two provenance regions of Q. suber used in this work, since previous studies of Spanish Q. suber populations pointed to the homogenization of the genetic structure of the species due to the existence of gene flow across populations ([22], [39]). However, differences in certain traits such as survival ratio could be due to local adaptation processes. It has been shown for this species that summer drought and temperature are important selective agents and that it may exhibit some phenotypic plasticity to inter-annual variations in rainfall ([40]).

Shrub cover

Our data indicated an effect of the latitude on shrub cover (L-R χ² = 8.0, df = 1, P < 0.05), with significant lower values in the N (2.3 ± 0.2) than in the S (3.1 ± 0.3). Besides, there was a significant interaction between location and the type of property (L-R χ² = 4.5, df = 1, P < 0.05 - Fig. 4b). Similarly to the regeneration, significant differences were found in the public plots (Tukey’s HSD tests, α = 0.05; 2.2 ± 0.2 in the N vs. 3.6 ± 0.2 in the S) but not in the private ones (2.4 ± 0.2 in the N vs. 2.6 ± 0.4 in the S).

Over the past 30 years, there has been greater and more continuous industrial cork activity in northern plots (both public and private properties) given the high value of this material, implying regular shrub clearing to facilitate cork extraction in this area. Nevertheless, we suggest that the economic crisis of the last 12-13 years has led to a significant decrease in forest clearing works by the forest administration and local governments. That would explain the differences of current shrub cover in public properties, given that the clearings have been focused only on the plots with more cork production (N). In fact, for a few decades now, cork harvesting and clearings are no longer done in some properties of the S (mostly public ones - [26], [10]). Thus, the shrubland could be evolving without the perturbations associated with such activities, remaining undisturbed and increasing in abundance.

Furthermore, we are witnessing the abandonment of grazing in some areas, which allows the shrub cover to slowly increase, due to favourable rainfall condition and mild annual temperatures in the S location. Plant cover generally increases with rainfall but only to a certain limit because of the interactions with other factors such as soil fertility ([48]). In addition, there is less competition in the S between adult trees and shrubs (lower tree density) in terms of capturing light, water and soil nutrients. There is a clear evidence that competition for resources in Mediterranean forests affects plant size, growth, germination and flowering ([13]). Less competition could improve the metabolic processes and allows a higher growth of different shrub species, although we are facing an ecological and edaphic process with long temporal pattern in which other factors such herbivory pressure also interacts.

  Conclusions 

We found an overall influence of property type and geographic location on cork oak population dynamics (including mortality, forest density, regeneration and tree size) likely attributable to differential climatic conditions and historical management. Total rainfall amount showed a latitudinal gradient, with higher values in the south. This extreme of the gradient exhibited a clear negative effect on the two main variables describing the population dynamics of Quercus suber: mortality and regeneration. Likewise, it is proposed that herbivory pressure is an important biotic factor causing a collapse in the population recruitment dynamics. Despite the difficulty of dealing with forest decline in Q. suber forests given the large number of interacting factors and the lack of ongoing monitoring actions, we believe that regeneration can be facilitated by controlling herbivory pressure. Moreover, the introduction of silvicultural practices would improve the stability of these unique Mediterranean forest systems in a scenario of increased aridity. Emphasis should also be placed on the application of preventive treatments to eradicate pathogens in forest soils as well as anti-oomycete fungicide chemicals for controlling diseases in Quercus species. However, protocols of soil improvement and research on mycorrhizal processes should be prioritized as sustainable management strategies. It is therefore critical to implement conservation and management plans adapted to local needs ([43]) in order to improve the stress resilience of trees ([16]) and hence make possible to reduce the vulnerability of Q. suber forests to dieback processes.

  Acknowledgements 

This study has been partially funded by the Jaime González Gordon Foundation, which granted AGR. The Management of Los Alcornocales Natural Park allowed the access to the public plots and facilitated the research. We would also like to thank Héctor Carrasco for his help in the statistical analyses comparing temperature and rainfall data, as well as Dr. Alexandro Leverkus (University of Granada, Spain) and Dr. Pablo González (University of Córdoba, Spain), who made comments on the manuscript draft and provided valuable suggestions. Francisco Conesa, Gabriel Agüera and Virginia Luque helped during fieldworks. Helena Valverde also helped during fieldwork and did the complete translation of the article into English. Raquel Cordero, a technician from the Natural Park, provided comments on insect plagues and tree conservation. We are also thankful to the anonymous reviewer that provided constructive comments on a previous version of the manuscript. JL-J was supported by a predoctoral grant from V-PPIT - U. Sevilla (Fifth Research Plan from the University of Seville), LM by the fellowships II5B and IV2 from VI-PPIT - U. Sevilla and RS-S was supported by VULBOS project (UPO-1263216, FEDER Funds, Regional Government of Andalusia, Consejería de Economía, Conocimiento, Empresas y Universidad 2014-2020) and LESENS (RTI2018-0968 84-B-C33) project from the Spanish Ministry of Science, Innovation and Universities.

VJD conceived the study; VJD and AGR carried out the field sampling; VJD and JL-J wrote the first draft of the manuscript; JL-J performed the data analysis and representation; RS-S and LM helped to interpret the results and provided substantial corrections to several manuscript drafts; FDO provided feedback and discussed ideas in early versions of the manuscript.

  References

  Supplementary Material

Fig. S1 - Map of Los Alcornocales Natural Park and specific location of study plots.
Fig. S2 - Variation in climatic parameters through time at the two different latitudinal study areas.