Regeneration Dynamics in Fragmented Landscapes at the Leading Edge of Distribution: Quercus Suber L. as A Study Case

Aims We studied the regeneration dynamics of woodlands and abandoned old elds in a landscape dominated by Quercus suber in its lower limits for rainfall and temperature. Two hypotheses were established: (1) recruitment of Q. suber is restricted more by abiotic variations than other species adapted to more extreme Mediterranean conditions; and (2) decreases in precipitation reduce growth, but temperature positively affects growth in the leading cold edge of this species distribution area. Methods We selected nine sites containing forest stands and old elds with and without tree remnants, and analyzed stand structure, soil parameters and tree growth. Results and September.


Introduction
Human activity has in uenced the structure and dynamics of forests in the Mediterranean region over centuries (Barbero et al. 1990; Scarascia-Mugnozza et al. 2000; Chauchard et al. 2007; Camisón et al. 2015). This is the case of forests and woodlands dominated by the cork oak, Q. suber L., managed systems that are protected by the European Union (Habitat directive 92/43/EEC). Q. suber is an evergreen tree species present in the western Mediterranean region that extends through the Iberian Peninsula to the western rim of the Italian peninsula, as well as some Mediterranean islands (Corsica, Sardinia, Sicilia and Balearic Islands) and coastal plains and hilly areas in North Africa, ranging from Morocco to Tunisia (Blanco et al. 1997; Magri et al. 2007; Pausas et al. 2009a). Covering about 2.2 million hectares, the current distribution area of this species is rather patchy, suggesting that in the past there was a more continuous distribution and that much of what we see today is relictual (Blanco et al. 1997, Carrión el al. 2000, Sánchez-Palomares et al. 2007, Pausas et al. 2009a, Jovellar et al. 2010). Additionally, intensi cation of ongoing climate change is expected to increase temperatures and the length of dry spells in the Mediterranean basin (Kovats et al. 2014). Consequently, areas suitable for cork oak are predicted to become generally reduced over the twenty-rst century, owing also to intermediate and high CO 2 emission scenarios Vesella et al. 2017).
Cork oak can form closed woodlands, but in most parts of its geographic distribution in the Iberian Peninsula and North Africa it usually appears in open woodlands with other dominating or co-dominating Quercus (Q.ilex , Q.faginea) or Pinus (P. pinea, P.pinaster) species (Blanco et al. 1997; Bugalho et al. 2009). Moreover, these open woodlands have been traditionally managed as agro-silvopastoral systems providing goods such as pasture for livestock, acorns for the high-quality pork industry, cereal crops, rewood, cork and game species. However, these woodlands also provide services in the form of recreational tourism, as cork oak systems represent an important cultural heritage in the Mediterranean region (Bugalho et al. 2009;Ovando et al. 2009;Vallejo et al. 2009). In addition, in many parts of Europe, since the second half of the twentieth century, the Iberian Peninsula has been subjected to substantial land abandonment due to socioeconomic factors. These factors have led to the decline and/or disappearance of traditional management systems, crop abandonment, shrub encroachment and an increased risk of re (Bugalho et al. 2009;Vallejo et al. 2009). The recruitment of oaks after land abandonment is usually a slow process, especially in shrublands, where succession is strongly delayed or almost halted in a condition known as arrested succession ( In this work, we took advantage of the heterogeneity of the landscape generated by anthropogenic activities in the lower ecological limits of rainfall and temperature supported by Q. suber as a representative case to study the regeneration dynamics and growth at the leading distribution edge. This area is expected to be highly vulnerable to changes in traditional management and can be useful for detecting early warnings signs. As such, we have established two main hypotheses. This rst hypothesis states that recruitment of Q. suber is more restricted by the abiotic variations induced by the dominant elements of vegetation (nursing effects) than other species more adapted to extreme Mediterranean conditions (e.g.: Q. ilex). Therefore, we expected that Q. suber recruitment would respond more drastically to changes in vegetation structure through the landscape, performing worse than Q. ilex in harsh environments like open old elds. To test this hypothesis, stand structure, soil parameters and spatial patterns of trees and saplings were analyzed in three types of landscape use. The second hypothesis states that tree growth is driven by local climatic conditions (mainly precipitation and temperature) affected by climate change, and that a decrease in precipitation reduces tree growth in water limited regions. We therefore predicted that temperature would positively affect tree growth in the leading cold edge of its distribution.

Study site
The study area is located in the north subdivision of the Central plateau of the Iberian Peninsula ( Fig. 1a; 41° 07' N, 5° 47'W; 800-850 m a.s.l.). The mean annual precipitation is around 380 mm, with a typical Mediterranean period of low precipitation during July and August. The mean annual temperature is around 12°C, with mean temperatures between 3-4°C and 20-21°C during the coldest and the warmest months, respectively. This area contains the biggest Q. suber woodland in the northwest part of the Iberian Peninsula, with an extension of almost 2,000 ha (Guerra Velasco 2015). The study area is located in the distributional edge of this species, as Q. suber typically occurs between 0 and 800 m a.s.l., it requires an annual precipitation between 600 and 1000 mm and an average temperature around 15°C (Blanco et al. 1997 Plot establishment and data collection Nine sites were selected encompassing three forest stands dominated by Q. suber, three old elds with tree remnants and three old elds without tree cover. Site selection was based on changes in landscape use determined by comparing aerial photographs of the Spanish Inter-ministry ight from 1973-1986 (https://fototeca.cnig.es) and aerial photographs taken from Google Earth from 2015 (Fig. 1bc). The sites located in old elds, with and without tree remnants, had been cultivated in the 70s and 80s and abandoned in the late 1990s. The forest stands were non-cultivated areas managed for traditional cork extraction, which is still harvested today. In each site, we established a 40 x 40 m plot (1600 m 2 ) where the origin point of each one was randomized and their sides were oriented in the directions of the cardinal points. All data were collected between October 2015 and June 2017.
All trees, both live and dead, and the saplings and seedlings in each plot were recorded. Trees were de ned as individuals with a diameter at breast height (dbh) ≥5 cm and saplings as individuals with a dbh < 5 cm and height > 200 cm. Small seedlings were de ned as individuals < 50 cm in height and large seedlings were de ned as individuals ≥ 50 and ≤ 200 cm in height. Resprouting individuals were not considered, and in some doubtful cases the soil surrounding each plant was excavated to check for independence from the nearest plants. The positions of the trees and saplings were located to the nearest centimeter using measuring tapes that were aligned with the sides of the plots, providing X and Y coordinates in a Cartesian plane. Shrub cover (%) was measured sampling ve line transects (5 m in length) per plot.

Dendrochronological analysis
Increment cores from all trees were extracted using Pressler increment borers (Häglof, Sweden) at 0.3-0.4 m above ground level to obtain the most accurate age for each tree (Veblen 1992), and at 0.6 m when the tree centers were rotten. Increment cores were mounted and sanded following the procedure established by Stokes and Smiley (1968), and the annual rings were counted using a stereomicroscope (SMZ800, Nikon, Japan). When the cores did not reach the pith, the number of rings to the center was estimated using the geometric procedure described by Duncan (1989). If the center was rotten, the rings counted in the non-rotten section of the core were considered as the minimum age for that tree. Cores were scanned at 2000 dpi resolution (Perfection V550, Epson, Japan), and tree-ring widths were measured with a 0.01 mm resolution on the scanned JPG images using the software CooRecorder 7.6 (Cybis, Sweden). The visual and statistical cross-dating of the tree-ring width series was done and checked using the software CDendro 7.6 (Cybis, Sweden) and Cofecha (Holmes 1983), respectively. The tree-ring width series were detrended to obtain mean residual chronologies of the tree-ring width indices. First, a spline function was tted to each tree-ring width series to obtain standardized indices which preserve the high-frequency variability potentially related to climate. Second, an autoregressive model was applied to remove the rstorder temporal autocorrelation in the detrended series and generate residual indices. Third, a biweight robust mean was computed to produce residual mean chronologies for each species. Detrending residual chronologies were obtained with the package dplR (Bunn 2008

Soil analysis
Five soil subsamples per plot were taken using a soil core sampler at depths between 0 and 30 cm. One subsample was taken in the center of the plot and four subsamples were taken 20 m apart from the center towards each corner of the plot. The subsamples were pooled in a single sample for structural and chemical soil analysis. The samples were air-dried and passed through a 2-mm sieve before laboratory analysis.
The pH was determined in distilled water (in a ratio 1:2.5) using a CRISON micropH 2001 pH-meter. Soil texture was quanti ed by the Robinson's pipette method after previous dispersion with sodium hexametaphosphate (Loveland and Whalley 1991). Organic carbon was determined by wet oxidation with a dichromate-sulphuric acid mixture (Walkley 1947). Residual dichromate was back titrated using ferrous sulphate. The difference in added FeSO 4 compared with a blank titration determined the amount of easy oxidizable organic carbon (Walkley 1947). We used a conversion factor of 1.72 to convert organic carbon to organic matter (Nelson and Sommers 1996). Available Ca, Mg and K were extracted using ammonium acetate 1M and determined by plasma ICP-MS. P was determined by the Bray I method, modi ed from Bray and Kurtz (1945).

Statistical analysis
To analyze how seedling regeneration (small seedlings, large seedlings and saplings) could be explained by soil and structural stand variables, a Redundant Analysis (RDA) was performed (Borcard et al. 2018). Previously, seedling data was transformed using a hellinger transformation, as suggested by Legendre and Gallagher (2001) for abundance data. To avoid collinearity within sets of soil and structural explanatory variables, we examined the correlation between variables and groups of variables using a hierarchical cluster analysis approach (complete linkage method). When the correlation between variables or groups of variables was greater than 0.6, only one variable was selected, which resulted in variance in ation factors (VIF) always less than ve units within the set of variables (Quinn and Keough 2002). For the soil variables we selected organic matter, P (%), pH and sand content (%), and for the structural stand variables we selected basal area and tree density (Table S2, Figs S1 and S2). Then, we performed a variation partitioning analysis (9999 random permutations) to analyze how much of the variance in seedling regeneration was explained by these variables (Borcard et al. 2018). Variation partitioning and RDA were analyzed with the vegan library under the R environment (R Development Core Team 2018), whereas HH library was used to evaluate variance in ation. Bootstrapped Pearson's correlation functions were calculated between the mean residual chronologies of each tree species and monthly climatic variables (mean, minimum and maximum temperature, precipitation). Con dence envelopes were obtained by calculating 1000 bootstrap samples and tested for signi cance using the 95% percentile range method (Dixon 2001). The climatic window for these analyses spanned from the previous September (year t-1), i.e. prior to the year of tree-ring formation, up to the following September (year t). This analysis was performed with the function dcc implemented in the package bootRes (Zang and Biondi 2013) under R environment (R Development Core Team 2018). The trends in monthly or seasonal climatic data were analyzed using the Kendall τ statistic with one-tail tests, as we were interested in decreasing rainfall and increasing temperature patterns.

Results
Woodlands showed tree densities of 95.8± 29.2 trees ha -1 and basal areas of 11.84±1.34 m 2 ha -1 (mean± se), where Q. suber was the dominant tree species. Mean tree dbh was 29.9 cm (con dence interval of 24.7-36.2 cm at 95% level). Quercus species showed a broad and multiage structure, but with most of the trees recruiting between 1930 and 1980 (Fig. 2a). Old elds with tree remnants showed lower tree density and basal area than woodlands (62.5± 12.5 trees ha -1 , 5.20±1.95 m 2 basal area ha -1 , mean± se), and with most of the trees (66.7% of total trees) recruiting after 1970 (Fig. 2b). Mean tree dbh was 19.6 cm (con dence interval of 16.1-23.2 cm at 95% level). Q. suber was also the main species in the old elds with tree remnants, although Q. ilex dominated tree recruitment after 1990 (Fig. 2b).
In relation to seedling recruitment, oak recruitment was quite abundant in the woodlands (2410.4± 719.8 seedlings ha -1 , mean± se) and dominated by Q. suber, especially in the shortest height classes (Fig. 3a). Signi cant differences in the mean recruitment of seedlings were observed between species and height classes, with clear differences in the shortest categories but insigni cant ones in the highest categories ( Table 2). In the old elds with tree remnants, oak recruitment was moderately abundant (1266.7± 445.6 seedlings ha -1 , mean± se, Fig. 3b) and seedling recruitment varied between species and height class ( Table 2). By contrast, Q. ilex was the only species that recruited in the old elds without tree cover, although with low densities (33.3± 10.4 seedlings ha -1 , mean± se, Fig. 3c). Sapling density of Quercus species was low in woodlands (66.7± 29.6 saplings ha -1 ) and old elds with tree remnants (43.7± 64.9 saplings ha -1 ) and zero in old elds without tree remnants. Differences in sapling density were found between forest types (df= 2, F= 7.06, p= 0.005) but no differences were found between species (df= 2, F= 0.62, p= 0.56).
In woodlands and old elds with tree remnants, small seedlings of Q. suber showed a clear clumped spatial pattern at distances up to 8-12 m and were associated to Quercus trees, with more seedlings than expected under a random process up to distances of 5-7 m from each tree (Table 1). Small seedlings of Q. ilex presented a clumped pattern, but at greater distances in old elds with tree remnants than in woodlands. Q. ilex small seedlings were also associated with trees but at lower distances than Q.suber, especially in woodlands (Table 1). By contrast, large seedlings of both Q. suber and Q. ilex mostly showed a random pattern distribution and were independently distributed from trees in woodlands whereas in old elds with tree remnants only large seedlings of Q. suber showed association with trees at distance of 4-7 m (Table 1).
In the RDA, the rst axis represented a contrast between old elds without tree cover at the right, where only small seedlings of Q. ilex were present, and woodlands at the left, where small and large seedlings of Q. suber dominated tree regeneration (Fig. 4). Basal area, pH and soil organic matter showed a strong correlation with the rst axis, with higher values of these three variables associated with woodlands. The second axis was related to the abundance of Q. faginea with regard to tree regeneration, with the sites with the higher scores in the upper panel showing a higher abundance for this species (Fig. 4).
Variance partition analyses showed that soil and structural variables explained most (81.0%) of the variance in the abundance of regeneration across sites (Fig. 5). Shared explained variation for both sets of variables was 31.6%, with soil variables contributing slightly more to the unshared explained variation in seedling abundances than structural variables (soil variables: 27.1%, F= 3.85, p= 0.014, permutations= 999; structural stand variables: 22.3%, F= 5.67, p= 0.004, permutations= 999).
Ring growth was positively correlated with rainfall in January and June (Fig. 6a), and September temperatures (mean, maximum and minimum) and minimum winter-early spring temperatures also positively in uenced ring growth (Fig. 6bcd). Monthly rainfall in January and June did not show a decreasing trend over time (Table S1), although annual rainfall signi cantly decreased (Table S1, Fig.   S3). All mean, maximum and minimum temperatures in January, February, March and September did not show temperature increases over time, except for the minimum temperature in September (Table S1).

Discussion a) Regeneration dynamics of Quercus
Succession was considerably arrested in the plots without tree remnants two to three decades after cultivation abandonment. Our results indicate that remnant oak trees are great accelerators of forest recovery after cultivation abandonment. Also, old elds with tree remnants boost oak regeneration almost 40-fold in comparison to old elds without trees, and woodlands show a 2-fold increase in oak regeneration in comparison to old elds with tree remnants. Tree cover is considered to play a fundamental role in Quercus recruitment through several processes. First, acorns are heavy fruits with limited dispersal, where most of the acorns produced only reach areas situated close to the parent trees ( Our results are also consistent with these patterns, as tree basal area and density were strongly correlated with the rst RDA ordination axis that ordered sites with a strong gradient of seedling density from right to left (Fig. 4). Additionally, the spatial analysis of small seedlings of Q. suber and Q.ilex, which accounted for most of the seedling abundance (71.6%) showed a clustered pattern associated with Quercus trees in both woodlands and old elds with tree remnants, although at shorter distances in Q. ilex than in Q. suber. By contrast, Quercus large seedlings showed a random spatial pattern independent of the trees, highlighting that tree cover is only a limiting factor during seedling establishment. These changes in the spatial pattern with ontogenetic development support previous ndings that suggest that the positive effects of shaded microhabitats are reversed for seedling development (Pérez-Ramos et al 2011, Pausas et al. 2009b). In Q. suber and Q. ilex, seedling establishment and survival are improved under shade (Espelta et al. 1995, Pausas et al. 2009b, Smit et al. 2009, Pérez-Ramos et al. 2012). However, it has been also shown that low light suppresses growth, that these species establish 'seedling banks' under dense tree cover and that more open canopy conditions are required for saplings and trees to develop (Espelta et al. 1995, Pausas et al. 2009b, Pérez-Ramos et al. 2010, 2012. In addition, we also found signi cant differences in the regeneration densities of Q. suber,Q. ilex, Q. faginea, with the abundance of Q. suber strongly correlated with the basal area and tree density of stands. In woodlands, Q. suber dominated seedling regeneration in the categories with the shortest height, but the differences in seedling abundance between Q. ilex and Q. suber disappeared in the tallest categories and showed similar densities. In the old elds with tree remnants, Q. ilex and Q. suber showed similar seedling densities, and although the interaction was not signi cant, the abundance of Q. ilex large seedlings was 5-6 fold higher than that of Q. suber. These results indicate that Q. ilex has a higher seedling survival rate and more likely to reach the sapling stage, probably due to their higher tolerance to shade (Sevilla 2008), and in particular their greater tolerance to hydric stress during the summer (Plieninger et al. 2010, González-Rodríguez et al. 2011, San-Eufrasio et al. 2020). Thus, although Q. suber trees dominate the landscape of our study area, the results indicate that Q. ilex can produce a greater number of young trees as compared to Q. suber (Fig. 2b), owing to the better performance of Q. ilex seedlings. Only in woodlands can Q. suber partially compensate for their lower performance and survival, producing a higher number of seedlings. Conversely, the deciduous Q. faginea present the lowest abundance, especially in the old elds (with or without tree remnants), which is in line with its lower tolerance to drought compared to evergreen species ( The presence of tree remnants showed a strong effect over soil parameters with higher concentrations of organic matter, N, exchangeable cations (K + , Ca 2+ , Mg 2+ ) and slightly less acidic soils than in old elds without tree remnants (Table S2). The deep root system of Q. suber trees uptakes and pumps basic cations, especially Ca 2+ , from the deep to upper soil layers throughout litterfall production and decomposition, which signi cantly improves soil nutrient availability (Serrasolses et al. 2009, Rossetti et al. 2015. The results of the variance partitioning showed that soil and stand structural parameters explain a signi cant amount of the shared variation in the amount of regeneration of Quercus species (Fig. 5), which makes sense since the in uence of trees on soil characteristics are spatially correlated.
However, soil variables also explain a signi cant amount of the variance (27.1%) not explained by structural stand variables and are important for understanding the differences found in Quercus regeneration between plots. Among soil variables, organic matter content summarizes most information on nutrient availability with which is strongly correlated (Table S2, Fig. S1), but it is also directly involved in improving soil water retention in abandoned cultivars, especially in soils with high sand content (Rawls , in our study site, the mean and minimum temperatures during winter and/or spring and the mean and maximum temperatures during September exerted a positive in uence on tree growth. This is consistent with the cold temperature conditions occurring in winter and early spring in our study area, located at the leading temperature edge of its distribution limits in the Iberian Peninsula. The large vessels in oak trees are very sensitive to winter embolisms caused by freezing temperatures, and in spring the reactivation of growth is greatly dependent on hydraulic conductivity recovery (Hacke and Sauter 1996, Lebourgeois et al. 2004). This highlights the importance of mild winter-early spring temperatures on Q. suber tree growth in the "cold leading edge" of its distribution. However, although winter temperatures are a limiting factor and warming has room for net positive effects on tree growth (Sánchez-Salguero et al., 2015), climatic data did not reveal a signi cant temperature increase during winter and early spring in this study site. On the other hand, minimum temperatures during September have signi cantly increased during the last decades and have had a positive impact on tree growth by most likely extending the favorable weather (Marqués et al. 2018). However, the positive effect of an extended growing season in early autumn can be counteracted by a decrease in rainfall, especially in the study area where annual rainfall is also in the lower edge of Q. suber distribution limits, as has also been recently observed in Mediterranean forest at the limits of species distribution (Madrigal-González et al. 2018, Marqués et al. 2018). Although no signi cant decreasing trends were found in the more critical months (Fig. 6a), there is a decreasing trend in mean annual rainfall (Fig. S3), so more detailed studies will be needed to understand the combined effects of temperature and rainfall on tree growth within future climate change scenarios.

Conclusions
Regeneration dynamics were strongly modeled by the presence and density of tree cover in the fragmented landscape in the cold and dry edge of Q. suber distribution. We have seen that tree cover affects seedling abundance differently through both direct (acorn availability and shading) and indirect (improving nutrient availability and water retention in soils) mechanisms. Consequently, the stress tolerant species Q. ilex was the only species found to recruit in open old elds decades after being abandoned. Also, the presence of isolated tree remnants in old elds allowed the recruitment of Q. suber, but Q. ilex had a higher abundance of large seedlings and young trees; only in woodlands did both species show similar recruitment success. These results indicate than Q. ilex is more abundant in this landscape compared to Q. suber, even though the latter species is the dominant tree in our study area.
In addition, the climate-tree growth relationship of Q. suber presents contrasting results. Higher temperatures were found to have a positive effect on tree ring-growth, especially during winter/early spring and early autumn. On the other hand, winter and spring rainfall recharges water soil reserves promoting tree growth during the favorable season. These opposing effects increase pose uncertainty in predicting Q. suber growth and productivity under climate change scenarios involving higher temperatures and less rainfall. However, when the demographic processes are considered, less water availability is likely to play a critical role in favoring Q. ilex recruitment in contrast to Q. suber, which could lead to changes in dominance hierarchies, especially in abandoned areas with scarce or absent tree cover. Table 2 Results of the two-way variance analysis testing for differences in seedling recruitment between species and height classes.