Disentangling the climatic and biotic factors driving changes in the dynamics of Quercus suber populations across the species‘ latitudinal range

Impacts of different global change drivers are altering the performance of plant species worldwide. However, these pressures usually differ across the species' distribution range. To properly assess the combined effect of global change at species level, we need to evaluate its consequences across their complete distribution. We focused on recent decline in Cork oak (Quercus suber L.) populations given its high ecological and economic relevance.


| INTRODUC TI ON
The climate alteration recorded during the last decades has the capacity to modify performance and biotic interactions of plant species worldwide (Parmesan, 2006). Although changes in climate might impact species across their distributions, the most dramatic effects are expected at species range edges, where environmental conditions are currently at species' tolerance limit (Hampe & Petit, 2005). It is commonly assumed that populations would expand at the leading edge of the species distribution, whereas population decline or range contraction might occur at the rear edge (Matías & Jump, 2015). However, expanding or declining trends are not always easily detectable under field conditions. At the leading edge of species distributions, populations need enough time to produce and disperse seeds and to effectively establish new individuals beyond current limits once environmental conditions become less limiting, a process known as colonization credit (Jackson & Sax, 2010). On the other extreme of the distribution, by contrast, the longevity and resilience of old trees might assure population maintenance for some time even under suboptimal environmental conditions, also known as extinction debt (Hanski & Ovaskainen, 2002;Jackson & Sax, 2010). Consequently, it is paramount to explicitly consider the differential trends of the different life stages in population dynamics occurring across the distribution of plant species to properly assess their response to changing climatic conditions. However, in spite of its relevance, studies comprising the complete range of a species are really scarce due to the difficulties for compiling reliable field data at large geographical scales (Parmesan, 2006).
Predictions of species population dynamics along geographical gradients are often difficult to obtain for several reasons. It is acknowledged that combined global change drivers can interact resulting in complex effects on population dynamics (Avila, Linares, García-Nogales, Sánchez, & Gómez-Aparicio, 2017;Wong & Daniels, 2017), but the nature of such effects are often unpredictable, limiting the inclusion of these interactions into species distribution models (Leach, Montgomery, & Reid, 2016). In addition, it is commonly observed that seedlings, small trees and large trees respond differently to global change drivers, suggesting an important role of ontogeny in mediating stress responses (Lloret, Peñuelas, Prieto, Llorens, & Estiarte, 2009). Obtaining accurate predictions is even harder because global change drivers might influence population dynamics through its indirect effects on the strength of biotic interactions, such as those coming from intraspecific competition in dense stands (Ruiz-Benito, Lines, Gómez-Aparicio, Zavala, & Coomes, 2013) or pressures from herbivores or pathogens (Gómez-Aparicio et al., 2012;Matías & Jump, 2015). Therefore, a critical step is to obtain field information of the ecological consequences on species population dynamics of interacting climatic and nonclimatic global change drivers across the complete distribution of tree species.
The Mediterranean region is expected to be one of the most vulnerable to global change over the coming decades due to a combination of different pressures (Giorgi & Lionello, 2008).
However, despite the rising number of modelling studies aiming to predict consequences of global change for tree species distributions (Correia, Bugalho, Franco, & Palmeirim, 2018;López-Tirado, Vessella, Schirone, & Hidalgo, 2018), potential trajectories for the future dynamics of typical Mediterranean forests such as oak forests remain largely uncertain (Acácio, Dias, Catry, Rocha, & Moreira, 2017;Lindner et al., 2014). Among oaks, this uncertainty is particularly true for cork oak (Quercus suber L.) (Vessella, López-Tirado, Simeone, Schirone, & Hidalgo, 2017). Q. suber is one of the most important species in the western Mediterranean basin due to its wide distribution (about 2.5 millions of ha) and ecological and economical relevance. Q. suber forests provide multiple important ecosystem services such as enhancement of plant and animal biodiversity, carbon sequestration, protection from soil erosion, provision of wood, timber resources and cultural services (Olea & San Miguel-Ayanz, 2006). In addition, Q. suber forests produce cork, a renewable resource of high economic importance that generates over 30,000 employments related to cork industry (Varela, 2013). However, the persistence of these important forests over the coming decades is uncertain due to the increasing adult mortality rates, overgrazing pressures and lack of regeneration detected in many Q. suber populations (Camilo-Alves, Clara, & Ribeiro, 2013;Ibáñez, Gómez-Aparicio, Ávila, Pérez-Ramos, & Marañón, 2017). Multiple global change drivers have been proposed to cause Q. suber decline, such as increasing drought (Camilo-Alves, Vaz, Clara, & Ribeiro, 2017) and the root-rot disease caused by the exotic soil-borne pathogen Phytophthora cinnamomi Rands. (Brasier, Robredo, & Ferraz, 1993;Sánchez, Caetano, Romero, Navarro, & Trapero, 2006), which interact with the management practices of these forests related to cork harvest and livestock grazing. To date, most studies on the dynamics of Q. suber populations have been conducted at small (i.e., local or regional) spatial scales, which severely limits our understanding of the magnitude and spread of the main threats for the species.
Large-scale studies are urgently needed to identify the most important stressors and susceptible areas, a fundamental step towards the implementation of adequate management strategies that assure the sustainability of these important forest systems.
In this study, we aim to evaluate the main risks for the conservation of Q. suber populations across the complete latitudinal distribution of the species, analysing the relative importance of aridity and exotic soil pathogens (P. cinnamomi) on population structure. For this, we selected 10 different sites across the core distribution of   Table 1). The different populations were selected to maintain aspect, slope and management history as constant as possible (e.g., all populations showed evidence of cork harvesting, but not clearings, reforestations, understorey management or recent fire). Due to geographical constraints, there was only one population replicate at the southernmost site (Toufliht, Morocco). Monthly climatic data (mean, maximum and minimum monthly temperature and total precipitation) were obtained for the different sites from the Climate Research Unit database (Harris, Jones, Osborn, & Lister, 2014) for the 1901-2014 period (CRU TS v.3.23). A xerothermic index (Xi) was calculated from the climatic data following Grossmann, Romane, and Grandjanny (2002), using the formula: Xi = ∑(2TM -P) if 2 TM >P or Xi = 0 if 2TM ≤ P, where TM is the monthly mean of the maximum and minimum temperatures (in °C) and P is the monthly precipitation (in mm).

| Field sampling
In each population, we randomly placed ten transects of 25 m × 10 m separated from each other by at least 20 m, noting the species, number, health status (alive or dead) and size (based on diameter at breast height, DBH) of every tree (young trees: DBH between 5 and 15 cm; old trees: DBH > 15 cm) and the number and demographic stage (seedling: <2 years old; or sapling: >2 years old and DBH < 5 cm) of all recruiting individuals. In addition, we randomly selected 30 adult Q. suber trees per population and quantified the level of defoliation following the standard semiquantitative scale proposed by the ICP Forest Network (Fischer & Lorenz, 2011): no defoliation (0%), light defoliation (1%-25%), moderate defoliation (26%-60%), severe defoliation (61%-99%) and dead (100%) and the crop production, which was visually estimated using the method proposed by Koenig, Mumme, Carmen, and Stanback (1994) for Quercus species, with two independent observers assigning categories from 0 to 4 (0, no acorns; 1, a few seen after close scrutiny; 2, a fair number; 3, a good crop; and 4, a bumper crop, many acorns seen almost everywhere on the tree). Twenty-five acorns per tree were collected for seed weight determination. Field sampling was carried out from September to December 2015, from northern to southern sites following natural acorn maturation. To quantify P. cinnamomi abundance, three soil samples were collected per population beneath the canopy of adult trees showing defoliation signals. Soil samples were air-dried and sieved at 2 mm and analysed by a pathology laboratory specialized F I G U R E 1 Map of the sampled sites along a latitudinal gradient (red dots) over the natural distribution of Q. suber (blue area, source Euforgen; http://www. euforgen.org/). Population codes are available in Table 1 in the studied pathogen (Dpt. Agronomy, University of Córdoba, Spain). Aliquots of 10 g from each soil sample were processed as described in Serrano, Ríos, González, and Sánchez (2015), preparing soil suspensions in 100 ml Water-Agar 0.2%. Aliquots of 1 ml taken from the soil suspensions were plated on NARPH agar Petri dishes (20 dishes per sample). Colonies growing on the plates were morphologically identified and counted. Colonies of Phytophthora cinnamomi are clearly characterized by clustered hyphal swellings and smooth cell-walled chlamydospores (Romero et al., 2007). Results were expressed as colony-forming units per gram of dry soil (cfu g −1 ).
For an overall view of population dynamics across the complete latitudinal range of Q. suber, we calculated the following key ecological indicators of population performance from the collected data:

| Data analysis
Differences in the slope of climatic trends over time across sites were tested using ANCOVA tests. The causes of variation in the different response variables (reproduction, seed weight, demographic structure, sapling dominance, defoliation and mortality) were explored by means of generalized linear mixed models (GLMM). The effects of latitude, xerothermic index (average for the last 50 years) and P. cinnamomi abundance were included as fixed factors in all cases. For those variables calculated at the tree level (i.e., reproduction investment, seed weight and defoliation), tree size and identity were included in the models as additional fixed factors. Finally, stand density was included as fixed factor in the mortality models given its major importance as a driver of tree mortality in Mediterranean forests (Ruiz-Benito et al., 2013). Population was included in all models as random factor, which accounted for the observed interpopulation variability across the latitudinal range. GLMMs were performed using a Poisson error distribution for demography, a Gaussian distribution for reproduction investment, a Gamma distribution for seed weight and a Binomial distribution for defoliation and mortality. The significance of each source of variation included in these models was tested by comparing the values of their Akaike information criterion TA B L E 1 Population sites and codes (within brackets), location, elevation, mean annual precipitation (annual P), mean annual temperature (Temp.), tree density, basal area, regeneration (seedlings per hectare, data are not available for TOU site), mortality percentage and accompanying species of the studied populations across the latitudinal distribution of Q. suber  corrected for small sample size (AIC c ), as well as chi-square tests for the significance of the model (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). We determined the significance of fixed factors by comparing the AICs of a model with only the intercept against the models built, including each fixed factor one at a time and in all appropriate combinations (Bates, 2011). We further calculated Cohen's d to estimate the direction and effect size of the predictive variables included in the best-supported models (Cohen, 1988). Due to the lack of population replication and low tree density, data from the TOU site (southernmost edge, Table 1) were only used for analyses of seed weight and reproduction investment, but were not included in demography (demographic structure, sapling:adult ratio, dominance and mortality) analyses. All analyses were performed using the packages stats, lme4, MuMIn and effsize in r (R Development Core Team, 2016).

| Climate and pathogen abundance
Recent climate (1965-2014 period) differed across sites in mean annual temperature and total annual precipitation, with higher precipitation at the northern edge of the distribution and higher temperatures at central sites, which resulted in higher Xi values in central and southern latitudes (Supporting Information Figure   S1). Climatic conditions have significantly changed during the past century along the latitudinal distribution of Q. suber. Mean annual temperature rose across the studied populations during the century (Supporting Information Figure S2a). This positive trend in temperature was consistent for all the studied populations and TOU) distributions, we detected differential trends in precipitation since the beginning of the 20th Century (Supporting Information Figure S2b). Northern populations are increasing the amount of precipitation at a rate of 16.8 ± 6.7 mm per decade, while there is no change in precipitation at central populations and a negative trend appeared at southern populations (mean decadal reduction of 8.6 ± 3.5 mm).
The soil-borne pathogen P. cinnamomi was detected in the soil of 7 of the 10 studied sites. P. cinnamomi abundance showed a hump-shaped relationship with latitude because it was not detected in either the north (ARN) or south (CHA, TOU) extremes of Q. suber's distribution, but showed higher abundance at central sites ( Figure 4c), with values ranging between 1 and 32 cfu/g.

| Reproductive investment and seed weight
Reproduction investment was influenced by aridity and tree size ( However, seed weight showed a consistent geographical pattern, with heavier seeds as latitude decreased (Figure 2a). In addition, tree identity instead of tree size also explained variance in seed weight (Table 2), pointing to high intrapopulation variability in this trait.

| Demographic structure and sapling dominance
Demographic structure of Q. suber populations strongly varied across the studied sites (Supporting Information Table S1). The best models explaining changes in demographic structure included latitude, Xi and pathogen abundance as predictive variables (Table 2).
We detected a shift in the demographic structure across the latitudinal distribution of the species, as denoted by the significant interaction between age and latitude (χ 2 = 759.3, p < 0.0001). The effect of latitude was not detected for seedling abundance, but a clearer geographical signal appeared as individuals became older ( Figure 3). We detected a marginally significant positive effect of latitude on sapling density, whereas the opposite relationship was found for adult trees (Figure 3). Accordingly, the sapling:adult ratio increased significantly with latitude ( Figure 2b). Seedling density was negatively influenced by Xi (Table 2, Supporting Information   Table S2). Finally, soil pathogens also affected Q. suber demographic structure, as shown by the significant but modest negative relationship between the density of old Q. suber trees in the populations and the abundance of P. cinnamomi in the soil (R 2 = 0.2, p = 0.05).
Q. suber seedlings dominated the seedling bank across the complete latitudinal distribution of the species, with dominance values over 80% in most populations (Figure 2c). However, sapling dominance was influenced by latitude (Table 2, Supporting Information   Table S2). Q. suber sapling dominance was high close to the northern limit of the distribution, but other tree species gained importance as latitude decreased (Table 1).

| Defoliation and mortality
Defoliation was low across the sampled populations, with 90% of the sampled trees showing no or light defoliation. The best-supported F I G U R E 2 Univariate regression models indicating changes in seed weight (a), sapling: adult ratio (b) and regeneration dominance (c) of seedlings (grey dots, grey line) and saplings (open dots, black line) with respect to seedlings and saplings of other tree species across the latitudinal distribution of Q. suber. Significance values of correlations are noted as *p < 0.05, **p < 0.001 and ***p < 0.0001 model (i.e., with lower AIC c ) showed an effect of tree size and pathogen abundance on tree defoliation across the species distribution (Table 2). Defoliation increased with tree size and P. cinnamomi abundance, but this model could not be statistically distinguished from a more parsimonious model including only pathogen abundances (i.e., ΔAIC c <2 units), which suggests that the effect of tree size should be interpreted cautiously (Table 2).
Mean mortality across populations was 6.6% ± 1.7% and ranged from 0% to 27.9%. Across the distribution range of Q. suber, tree mortality was associated with population density and pathogen abundance but not with its interaction (Table 2). Specifically, those populations with higher tree density or higher pathogen abundance had also a higher proportion of dead trees (Figure 4). Our sampling did not detect a clear latitudinal pattern in tree defoliation and mortality, despite these two variables being positively correlated (R 2 = 0.40; p = 0.004).

| D ISCUSS I ON
We present here the results of a field study assessing the variations in demographic patterns across the complete latitudinal distribution of Q. suber, a tree species with an outstanding ecological and economic importance. Across the distribution of this species, we de-

| Reproductive investment across the distribution range
Q. suber is a masting species, with high intrapopulation synchrony in the production of large seed crops every two to three years as a mechanism

| Variations in the regeneration ability
A demographic structure dominated by old individuals and low recruitment is a reliable indicator of declining populations by recruitment failure, whereas populations characterized by a high proportion of seedlings and saplings compared to old individuals are considered to have a healthy status with high potential for population maintenance or expansion (Halpin & Lorimer, 2017). Thus, early detection of demographic changes in population structure can help us to identify potential declining areas before decline is irreversible.
In our study, we detected strong variation in demographic structure across the latitudinal distribution of Q. suber. Although we found populations with high regeneration capacity (over 1,000 seedlings ha −1 ) occurring at both extremes of the distribution, we identified populations where seedling density was low (ALC and TAZ sites with less than 100 seedlings ha −1 ) or very low (15 seedlings ha −1 at VAL site), clearly indicating a recruitment failure and therefore an early signal of population decline. We identified aridity as an important factor influencing seedling density across the species' distribution.
This result means that if the projected decrease in precipitation for the Mediterranean region is maintained (Giorgi & Lionello, 2008), it is very likely that recruitment failure would be the main constraint  (Avila, et al., 2017;Gómez-Aparicio et al., 2012). This pattern has been already described for some of the Spanish locations studied (Ibáñez et al., 2017), and we can now extend it to the southern distribution of the species. Our results also suggest far-reaching consequences for the distribution of the species on the mid-to long term, since the low recruitment, aged demographic structure and lower understorey dominance of Q. suber at its southern range are indicative of ongoing range contraction (extinction debt, Jackson & Sax, 2010), while the northern edge of the species evidence the colonization credit to expand its distribution on the coming decades F I G U R E 4 Univariate relationships between tree density (panel a) and Phytophthora cinnamomi abundance (expressed in colony forming units per gram of soil, panel b) on tree mortality. Panel c represents the variations in P. cinnamomi abundance across the latitudinal distribution of Q. suber (Matías & Jump, 2015), jointly resulting in a potential range displacement of the species.

| Defoliation and mortality intensities and their causes
Defoliation is a primary indicator of tree decline and has typically been attributed to the single or combined effects of climatic extremes such as winter frost or summer drought, defoliating insects and pathogens (Galiano, Martínez-Vilalta, Sabaté, & Lloret, 2012;Gómez-Aparicio et al., 2012;Sánchez, Caetano, Ferraz, & Trapero, 2002). Our results show that defoliation intensity was positively related to tree size across the complete latitudinal distribution of Q. suber. Bigger trees have a more complex hydraulic system and, in consequence, are more susceptible to embolism and hydraulic failure (Zhang et al., 2009), producing defoliation as a first step. Besides this structural cause, we also detected evidence for a positive effect of pathogen abundance on defoliation. The association between P. cinnamomi abundance and crown defoliation has been widely described for this species, causing loss of fine roots and therefore a reduction in the ability to uptake water from the soil (Sánchez et al., 2002;Serrano et al., 2015). This pathogen destroys fine roots causing a reduction in the ability to uptake water from the soil, therefore leading to similar water shortage as abiotic droughts, and eventually causing leaf shedding as a way to prevent xylem embolism through segmentation of the vascular system (Tyree, Cochard, Cruizat, Sinclair, & Ameglio, 1993). Yet, it is worth noting that defoliation typically shows high interannual variability (EU-ECE Forest Health Inventory, 2009), and a longer data set would be necessary to accurately determine the combination of factors causing tree defoliation.
The mortality levels detected across the Q. suber distribution in this study raise concern about the current risks for the species, reaching values as high as one fourth of the population at some sites (e.g., MER site). Our results indicate that mortality is determined by two main factors across the species' distribution: stand density and P. cinnamomi abundance. It is commonly assumed that denser populations increase competition for resources, inducing decreased survival or aggravating the effects of extreme dry events on population dynamics (Ruiz-Benito et al., 2013). Under a scenario of increasing aridity as the forecasted for the Mediterranean basin (Giorgi & Lionello, 2008), we could expect the effects of competition and drought to gain in importance over the following decades, particularly in dense populations. We already know that competition for resources disproportionately affects smaller size classes in dense stands (Wang et al., 2016), so this factor could also exacerbate the problem of recruitment failure detected at the southern distribution of Q. suber.
The pathogen abundance levels detected across the study populations were below the threshold for expression of disease in Q. suber seedlings under control conditions (Serrano et al., 2015), probably because populations were randomly selected without regard for decline signs. However, we found a clear positive relationship between P. cinnamomi abundance and Q. suber mortality. The strong association between P. cinnamomi and the decline and mortality of Q. suber has been known and studied for more than two decades (Brasier et al. 1993;Gómez-Aparicio et al., 2012;Sánchez et al., 2002), but this is the first report on the abundance of the patho-  Piou, & Vannini, 2006). Therefore, further research explicitly looking at the interactive effect of mortality causes is urgently needed for a deeper understanding of their consequences.

| CON CLUS IONS
Given the important ecological and economical roles of Q. suber in the Mediterranean basin, the detection of the factors affecting the population decline becomes especially relevant. Our results indicate potential for range displacement of the species driven by the recruitment failure at the southern edge of the distribution and a simultaneous high potential for range expansion at northern populations. The ongoing change in climate is altering the regeneration of Q. suber at the southern part of its distribution, primarily through hampered transition from sapling to adult. This pattern, together with the loss of dominance at the sapling stage, the high mortality rates in dense populations and the negative effect of P. cinnamomi abundance, highlights that the already threatened status of the species is explained by several drivers acting in combination. The decline of this important tree species at the southern part of the distribution might induce far-reaching changes in ecosystem functioning and provision of ecosystem services, from the loss of cultural and economic value to the impact on biological communities through cascading trophic effects (Acha & Newing, 2015;Correia, Haskell, Gill, Palmeirim, & Franco, 2015).
In order to turn back the current declining trends of Q. suber populations, we recommend the use of active restoration techniques enhancing the establishment of saplings at the southern extreme of the distribution through fencing and plantations certified to be free of P. cinnamomi. Moreover, it will be desirable to conduct management practices aimed at reducing stand densities as a way to diminish intraspecific competition for resources, to minimize adult mortality and to increase drought resistance (Ameztegui, Cabon, Cáceres, & Coll, 2017).
Equally important would be to promote control techniques aimed at reducing the abundance and dispersion of P. cinnamomi such as a strict control of the pathogen in forest nurseries (Jung et al., 2016), or the use of mineral amendments (Serrano, Leal, Vita, Fernández-Rebollo, & Sánchez, 2014) or biofumigation (Ríos et al., 2016). Our results show that conducting these actions is paramount at central-latitude populations to reduce the spread of the disease. Only by ameliorating the causes of adult mortality and improving regeneration would we be able to assure the permanence of these valuable forest ecosystems.

ACK N OWLED G EM ENTS
We

DATA ACCE SS I B I LIT Y
Data will be available from the Dryad Digital Repository upon acceptance.