Enhancing the Resilience of a Mediterranean Forest to Extreme Drought Events and Climate Change: Pinus — Tetraclinis Forests in Europe

: The southeast Iberian Peninsula is the only place in the European Community where Tetraclinis articulata (Vahl) Masters populations are native. In this area, the optimal ecological niche for this species is occupied by Pinus halepensis (Miller). The increasing intensity of extreme drought events induced by climate change causes severe declines in pine forests, while providing expansion opportunities for established Tetraclinis populations. Within the framework of the LIFE-TETRACLINIS project, a study has been designed to simulate the pine forest decline effects on the population dynamics of this protected species. This work investigates the effects of decreasing competition on T. articulata specimens with limited reproductive activity. To induce the reproductive activity of these specimens through increasing the availability of light, the surrounding pines were removed within a 15 m radius. Increased light availability was modelled using “Light Detection And Ranging” (LiDAR) data, and changes in the main reproductive parameters were registered throughout the study period. A signiﬁcant increase in the reproductive population was achieved, as well as the cones produced per specimen and the recruitment. Findings obtained are promising for the habitat management in continental Europe and enhancing this forest system’s resilience to extreme drought events and climate change.


Introduction
Drought-induced forest decline is a global phenomenon potentially driven by climate change [1]. Over the past decade, some authors have estimated that the frequency, duration, and severity of drought periods will increase in Mediterranean countries during the 21st century [2,3]. In the southeast of the Iberian Peninsula, these events have induced severe decline episodes in forests mostly dominated by Pinus halepensis Miller (Aleppo pine) [4,5]. This is especially true for species located at their limits of distribution and whose climatic suitability may be compromised by these events [6]. Forest management in Iberian Peninsula needs to take into account the predicted impacts in the context of climate change [7].
The southeast of the Iberian Peninsula is a key area for assessing climate change effects on ecosystem biodiversity. This region is an ecotone between the Mediterranean biome and the shrublands of semiarid areas of subtropical character [8][9][10]. Its coastal Pinus halepensis Miller-Tetraclinis articulata (Vahl) Masters mixed forest stands are particularly responsive to effects of climate change. In this regard, the hydraulic architecture of Tetraclinis provides a better resilience against drought events [11]. Indeed, regional research based on species distribution models shows opposite trends for both species. These models suggest an expansion of T. articulata habitat [12] in combination with a decrease in in the Aleppo pine area due to the increase of severe drought events recurrence [5]. However, declining pine forests do not have perfect spatial coincidence with the current populations of Tetraclinis. Therefore, a deeper understanding of these substitution mechanisms is required in order to facilitate a transition into a resilient forest to climate change.
Traditional practices used to expand T. articulata population areas have been exclusively focused on reforestation projects [13], which have not yet been evaluated in terms of cost-effectiveness in comparison with other potential alternatives. Until now, no initiatives have been implemented to increase recruitment in Tetraclinis populations with limited demographic dynamics.

Solar Irradiation and Reproduction
The role of solar irradiation in conifers reproduction has been studied extensively. Within the same population, trees located in areas of high light exposure tend to produce more cones than others [14,15]. Whereas, those located in shaded areas or closed stands may have low or even no production [16,17]. The relationship between the initiation of reproductive activity and tree size is modulated by environmental conditions, including the availability of light [18]. In fact, environmental conditions that promote vegetative growth could delay the beginning of reproductive activity [18][19][20]. Besides reproductive maturity, light availability also plays an essential function in the germination of seeds and recruitment rates. There is abundant documentation about environmental factors related to the seed germination process in the Pinus genus over the Mediterranean [21][22][23][24][25][26][27]. However, research on this topic for T. articulata is more recent [28][29][30]. Recruitment dynamics of the European populations for this species are strongly influenced by the sunlight exposure [31,32].

Interspecific Competition
In a Pinus halepensis-Tetraclinis articulata mixed forest, both species are in strong competition for light and other resources as they are at the lower limit of their distributional ranges [28,31]. Besides light availability, water availability caused by aspect and drainage flow is another limiting factor for both species in the southeast Iberian Peninsula [31]. In these locations the T. articulata ecological niche is partially occupied by P. halepensis [32]. The northern slopes are dominated by Aleppo pine and Tetraclinis is displaced to the south and southeast. The eastern slopes are where both species establish a most balanced competition. Shaded areas with a dense pine canopy contain a number of isolated T. articulata trees. Most of these specimens remain reproductively inactive due to the limited solar irradiation received.
We developed an experiment within the framework of the LIFE13 NAT/ES/00436 project (Conservation of the habitat 9570* Tetraclinis articulata forest in the European continent) to trigger reproductive behavior of T. articulata specimens through the creation of gaps in the forest canopy. Targeted specimens were in an inactive reproductive phase despite their size. The objective was to investigate the sequential response caused by the increase of light irradiation on the reproductive dynamics of this species (namely, reproductively active status, cone production, and recruitment). For this purpose, the results are evaluated across a range of light intensity and pine canopy.

Experimental Design
The study area is located at the Regional Park of "Calblanque, Monte de las Cenizas y Peña del Águila", in the southeast of the Region of Murcia in Spain ( Figure 1). We monitored a total of 29 plots (10 control and 19 experimental). In 2016, circular plots of 15 m radius were established centered on individual T. articulata specimens located in a matrix of Aleppo pine. The specimens were selected due to their limited or absent reproductive activity (i.e., absence of cones and nearby recruitment). In the fall of 2017, most of the P. halepensis specimens were removed from within the experimental plots, retaining 0%-50% pine canopy cover ( Table 1). The Aleppo pine specimens were logged and debarked, subsequently spreading the wood pieces over the plot. halepensis specimens were removed from within the experimental plots, retaining 0-50% pine canopy cover ( Table 1). The Aleppo pine specimens were logged and debarked, subsequently spreading the wood pieces over the plot.  Table 1. Pine canopy cover ranges in the intervened plots before (2016) and after (2018) the experiment.

Data Collection
Over the period 2016-2019, we collected diameter, reproductive status, cone production, and recruitment data for all the T. articulata specimens that were present inside the plots. Specifically, 2016 and 2019 data were collected from June to September and 2018 data from February to March. The diameter of the primary stem was measured at 20 cm above the ground or immediately at ground level when the tree size was less than 20 cm tall [13]. Cone production was estimated by visual counting. Total number of new T. articulata recruits were annotated and the percentage of pine canopy for each plot was fieldcalculated by linear transects before and after the experiment.

Data Collection
Over the period 2016-2019, we collected diameter, reproductive status, cone production, and recruitment data for all the T. articulata specimens that were present inside the plots. Specifically, 2016 and 2019 data were collected from June to September and 2018 data from February to March. The diameter of the primary stem was measured at 20 cm above the ground or immediately at ground level when the tree size was less than 20 cm tall [13]. Cone production was estimated by visual counting. Total number of new T. articulata recruits were annotated and the percentage of pine canopy for each plot was field-calculated by linear transects before and after the experiment.

LiDAR-Based Light Model
Due to the many difficulties encountered in reproducing light measurements (presence of clouds, day of the month, and exact time of day), a LiDAR-based light model was applied instead. For this purpose, we used the 2015-2016 LiDAR data available from the Spanish National Geographic Information Center (https://www.cnig.es/home, accessed on 2 September 2020) and GRASS GIS r.sun.daily model [33] in combination with the Laser Penetration Index (LPI) based on LiDAR data [34,35]. The 2016-point cloud was processed into a bare-ground and canopy Digital Elevation Models (DEMs) using FUSION/LDV v.4.20 [36]. The final DEMs grid size was 4 m × 4 m in order to decrease the data processing effort. Due to the lack of LiDAR data for 2019, we decided to change the canopy raster values of the plots using the available orthophoto images. The DEMs were used to calculate annual light models before and after the experiment. The light models were corrected by using LPI in order to estimate the annual direct and diffuse irradiation at ground level received through the forest canopy.

BACI Analysis Design
The BACI analysis [37] was designed to compare differences between 10 control and 10 impact (experimental) plots with similar initial biometric characteristics (Control diameter = 15.3 ± 4.3 cm, Impact diameter = 15.7 ± 6.2 cm) and the LiDAR-based light model's data. The other 9 plots were excluded due to significant differences in the diameter of the specimens (31.2 ± 6.4 cm). The analysis was applied to direct solar irradiation models and the number of reproductive specimens by using an ANOVA approach using R 'car' package [38] in order to explore significant differences between the groups before (2016) and after (2018) the intervention.

Reproductive Status Variation
Variation in reproductive status of the 19 experimental plots has been analyzed using binomial GLM [39]. The reproductive status was used as a binomial variable (active or inactive) versus the diameter (measured in cm) of all the specimens in each plot in 2016 and 2018. Both obtained models were used to estimate the difference in the reproductive status of the involved specimens after the experiment.

Demographic Response Observed in the Experimental Plots
In an attempt to explain the effects caused on the demographic dynamics of T. articulata during this study, we fitted negative binomial GLMs from 'MASS' R package [40]. Due to the low recruitment observed before the experiment (2016), the total number of specimens was used as alternative. We used the recruitment of 2019 as dependent variable. Explanatory variables with a direct effect on the species population dynamics were selected: two models of annual solar irradiation (direct and diffuse), total number of terrain's drainage flow lines in each plot, total number of reproductive specimens, previous year's cone production, and P. halepensis canopy for each survey plot. The number of drainage flow lines were calculated with SAGA module [41] using a 4 m × 4 m spatial resolution DTM. Furthermore, we used these variables in multiple regression models as a more complex approach to their interaction with demographic dynamics. These models were selected by using the minimum Akaike information criterion method (AIC) [42].

BACI Analysis
No significant differences were observed in the initial conditions between control and experimental plots. Statistics obtained by ANOVA tests are summarized in Table 2 and Figure 2. There are significant differences for both annual direct solar radiation and the number of reproductive specimens between the control and impact plots after the intervention. Differences are also observed for these parameters among the impact plots. The experiment resulted in a significant increase of direct solar irradiation and also in the number of reproductive specimens. The experiment resulted in a significant increase of direct solar irradiation and also in the number of reproductive specimens. Impact before-Impact after Direct solar irradiation 13.038 0.0019 ** Reproductive specimens 6.6977 0.0186 * Figure 2. Boxplot of (a) annual direct solar irradiation and (b) total number of reproductive specimens (red point represents the mean value). Cb: control before, Ca: control after, Ib: impact before, Ia: impact after.

Observed Variations in the Reproductive Status
Total number of T. articulata specimens at the experimental plots was 93 in 2016, 134 in 2018, and 545 in 2019. There were some important changes in the reproductive parameters of the specimens located on managed plots (Table 3). An increase of ten new specimens in the reproductive population was combined with a higher average cone production rate per reproductive specimens. The recruitment numbers in subsequent years suggest better availability of germination microsites and a great increase in the seed availability.

Observed Variations in the Reproductive Status
Total number of T. articulata specimens at the experimental plots was 93 in 2016, 134 in 2018, and 545 in 2019. There were some important changes in the reproductive parameters of the specimens located on managed plots (Table 3). An increase of ten new specimens in the reproductive population was combined with a higher average cone production rate per reproductive specimens. The recruitment numbers in subsequent years suggest better availability of germination microsites and a great increase in the seed availability. In 2016 there were 18 reproductive specimens of T. articulata in the experimental plots and this number increased to a total of 28 in 2018. Obtained GLM models relating maturity and diameter of the specimens are presented in the Table 4. To achieve a 50% chance of being reproductively active in 2016, the specimens needed a diameter of at least 15 cm (Figure 3a). After the experiment, this minimum diameter was reduced to about 10 cm. The observed difference (Figure 3b) between both curves denotes the diameter ranges whose probability of reaching reproductive maturity increased substantially. The change was greater in specimens with a diameter of 10-15 cm.
In 2016 there were 18 reproductive specimens of T. articulata in the experimental plots and this number increased to a total of 28 in 2018. Obtained GLM models relating maturity and diameter of the specimens are presented in the Table 4. To achieve a 50% chance of being reproductively active in 2016, the specimens needed a diameter of at least 15 cm (Figure 3a). After the experiment, this minimum diameter was reduced to about 10 cm. The observed difference (Figure 3b) between both curves denotes the diameter ranges whose probability of reaching reproductive maturity increased substantially. The change was greater in specimens with a diameter of 10-15 cm. The variation in the 2018 cone production rate was related with two main factors: the addition of new reproductive specimens and the greater response observed in those that were already active. From 2018's total production, only 13.55% was due to new reproductive specimens. The 2018 average production ratio was 2721 cones per previously reproductive specimen and 597 cones per new reproductive specimen. This represented an increase by a factor of 1.88 over 2016 observed levels.

Demographic Response Analysis
Pearson's correlation indices between the variables used in negative binomial GLM are shown in Table 5. The same predictor variables for 2016 and 2019 (solar irradiations, pine canopy cover, reproductive status, and cone production) were used for these models. As response variables, the total number of specimens (2016) and the observed recruitment (2019) were used. The variation in the 2018 cone production rate was related with two main factors: the addition of new reproductive specimens and the greater response observed in those that were already active. From 2018's total production, only 13.55% was due to new reproductive specimens. The 2018 average production ratio was 2721 cones per previously reproductive specimen and 597 cones per new reproductive specimen. This represented an increase by a factor of 1.88 over 2016 observed levels.

Demographic Response Analysis
Pearson's correlation indices between the variables used in negative binomial GLM are shown in Table 5. The same predictor variables for 2016 and 2019 (solar irradiations, pine canopy cover, reproductive status, and cone production) were used for these models. As response variables, the total number of specimens (2016) and the observed recruitment (2019) were used.
The most obvious are the linear correlations between both light intensities (direct and diffuse) and their negative association with pine forest canopy. The initial number of T. articulata specimens in 2016 is correlated with these three variables and with the number of previous reproductive specimens. The association of 2019 reproductive variables to the light irradiations and pine canopy seems to be more complex, with a progressive reduction in their linear correlation for the reproductive specimens (higher), cone production, and recruitment (lower). In the 2016 GLM analysis, the Tetraclinis specimen number showed a significant relation with both solar irradiation models (Table 6). In contrast, Aleppo pine forest canopy is negatively correlated. Regarding the recruitment observed in 2019, the cone production of the previous year is the most significant variable and direct solar irradiation effect disappears. The only variables maintained during the period evaluated were diffuse solar irradiation and P. halepensis canopy. After the experiment, response was greatest under 17%-18% Aleppo pine canopy cover (Figure 4). A rapid exponential response of T. articulata recruitment relative to the cones produced in the previous year is observed in 2019.    Multiple regression models (Table 7) indicate an explained deviance value of 48.15% for the total number of T. articulata specimens (2016) and 52.97% and for the recruitment after the experiment (2019). In the first case, the variables were the direct solar irradiation and the number of reproductive specimens. In 2019, the variables were the cone production of the previous year and pine canopy. Diffuse irradiation could have been incorporated to this model since it was significant (p-value <0.05), but it would have a negative effect by increasing the variance inflation factors (VIF). Table 7. Multiple regression negative binomial GLM models obtained for the total number of specimens (2016) and the observed recruitment (2019). Int.: intercept, Pr.a, b: predictor a or b, Pr. 2 : quadratic predictor. Significance codes: * < 0.05, ** < 0.01, *** < 0.001. Recruitment data collected in 2018 were used to study the early effects of the experiment. The only significant explanatory variable for the observed recruitment in that year was the number of drainage flow lines ( Table 8). The importance of this variable is very significant for this first year's recruitment. In addition, it shows preferences for two environmental scenarios: (i) convex areas of limited ecological competition with pine (low number of flow lines) and (ii) moderately concave areas (medium to high number of flow lines).

Discussion
This experiment increased the light availability for Tetraclinis articulata through the creation of gaps in a forest matrix dominated by Pinus halepensis. Increased direct solar irradiation received by Tetraclinis specimens triggered their reproductive and recruitment dynamics. These are characteristic responses of forest species in ecological competition under favorable changing conditions [43,44].
Three events are triggered after the reduction of pine canopy: (i) reproductively active T. articulata population increases, (ii) cone production and thus the number of available seeds rises, and (iii) population recruitment increases. These events begin when competition for direct sunlight decreases and they follow a sequential pattern.
First response to the decreased light competition is the increase of reproductive specimens during the first year after the experiment. Results of the BACI analysis suggest that increasing of solar irradiation seems to have direct repercussions in reproductive dynamics of T. articulata. Indeed, the number of reproductive specimens significantly increased only one year after the experiment. Increased production is a pattern previously observed in conifers exposed to direct solar irradiation [14,15]. Observed variations in the reproductive behavior curve after the experiment suggest that specimens with a diameter of 10-15 cm were activated from reproductive dormancy. The behavior of these specimens was comparable to those populations without Pinus halepensis competition [32].
It is known that competitive relationships can lead to morphological and biochemical responses in plants [45,46]. Under low direct solar irradiation conditions (shaded areas), T. articulata specimens tend to respond unfavorably to Aleppo pine presence [28,31], even affecting their reproductive maturation. Delay or absence of reproductive activity in trees located in shaded areas has been previously highlighted in the past [16][17][18][19][20]. Though there was no previous clear evidence in this regard for T. articulata, increased direct solar irradiation received by the specimens has immediately stimulated their reproductive and recruitment dynamics.
The second response is the increase in cone production, which is already noticeable in 2018 but is reflected demographically in 2019's observed recruitment. The effect of increased light on 2018 cone production results in those specimens previously active almost doubling their production rate per specimens (from 1445 to 2721). The new specimens that joined the reproductive cycle in that same year also showed a similar increase in 2019 (from 597 to 950 production rate per specimen). It was in this last year when the effects of masting became evident [13], with a reduction in the number of cone-producing specimens.
The third response is increased recruitment. This occurs already during the first year after the intervention, but it is greater after two years (in 2019) and is related to the increase in cone production observed during the previous year. Increased recruitment has already been observed occurring in 2018, without a significant change in seed production. In this case, the increased recruitment would be associated with an increased light availability. Tetraclinis articulata seedlings are considered strongly heliophilous [32], so the observed number of recruits in 2018 would correspond to higher availability of solar light. In contrast to the findings of previous authors [28], no adverse effects on recruitment caused by pine litter have been detected [47]. This is probably caused by the limited accumulation of pine litter in the experimental plots.
In relation to the recruitment behavior observed the first year after the experiment in relation to drainage flow lines, T. articulata presents two different trends. In shaded areas in competition with Aleppo pine, the specimens are relegated to areas of low number of flow lines with a high number of fissured rocks. In areas influenced by stronger solar irradiation, the species needs to compensate for evapotranspiration and requires greater water availability (a higher number of flow lines in the plots). This differential species behavior has been previously described by other authors [31].
Several interesting responses have been found in the single predictor models obtained in order to explain the observed recruitment and the initial number of T. articulata in the experimental plots. As expected, the set of correlated variables associated with incident radiance and tree canopy retained a very important weight in the number of total specimens (2016) and recruitment (2019). The role played by this set of variables on the reproductive activity becomes complex after the experiment. In 2019, direct solar irradiation had a strong positive correlation with the reproductive activity initiation and cone production, but its linear correlation with the recruitment is negative. Diffuse irradiance is apparently more relevant for cone production and recruitment, although in the latter case the correlation is not linear. The pine canopy is highly correlated to the solar irradiations (particularly with diffuse), although it may have an added role in the interception and redistribution of precipitation at the germination microsites [48]. Pine canopy is highly relevant during the whole study period. This variable shows a maximized value for recruitment in 2019, with a maximum level of 17%-18%. Without pine cover the recruitment of T. articulata is frequent, while with 40%-50% of pine canopy the recruitment is null. This information could be very useful in the management of these mixed-forests and also is consistent with previous studies about this species [28,30]. Recruitment shows an exponential response to the cone production, predominantly explained the observed recruitment in 2019. These results are consistent with previous studies done in the same area [32].
The evolution of the variables in multiple predictor models is particularly interesting. In the 2016 multiple regression model, the number of specimens depends on the reproductive specimens and direct solar irradiation. This behavior is consistent with previous T. articulata research [28] and demographic models of the southeastern Iberian populations [32], where the number of adult specimens is the primary factor driving the recruitment process. The linearity of the radiation and pine canopy cause its correlation to be sufficient for only one of them to be included in the model. The 2019 recruitment model explanation is absorbed by the previous year cone production and the pine cover. Diffuse solar irradiation also appears to play a minor role following the entry of pine canopy in the model. The quadratic behavior of these last two variables suggests them as modulators of the recruitment response against the number of available seed.
Obtained results highlight the capacity of T. articulata populations to thrive when gaps are created in the pine canopy. This behavior is consistent with other studies that report the species as a weak competitor compared to Aleppo pines [28,31]. The expected decline of the Aleppo pines due to extreme drought events related to climate change may provide a potential expansion opportunity for Tetraclinis populations in the Iberian southeast [5,12]. This experiment evidenced the solar light competition relation between both species, but the recruitment of T. articulata improves in plots with reduced canopy cover of P. halepensis (around 18%). Therefore, the initial competitive interaction between these species turns into facilitation when the pine canopy is reduced. A trade-off between available light and the water balance of the recruited seedlings may be the reason for this. The greater role played by diffuse compared to direct irradiation during the recruitment process seems to support this idea.

Conclusions
This experiment proves that competition for solar irradiation is a major factor conditioning the reproductive activity of T. articulata. Interventions aimed at decreasing competition with Aleppo pines triggered a reproductive response in only two years, increasing the effective number of reproductive specimens and also the cone production. Abundant recruitment of Tetraclinis has been observed with pine canopy cover around 18%, while at 50% it is practically absent. Therefore, Pinus halepensis behavior changes from a competitor to a facilitator species, increasing Tetraclinis articulata recruitment. After occupying the new microsites available during the first year, the recruitment of new specimens appears to be dependent on the number of available seeds. This experiment illustrates the high reproductive potential of T. articulata once gaps are created in the pine canopy, an expected scenario according to available climate change model predictions. However, gaps in the pine canopy may not overlap with Tetraclinis specimens, so this substitution process should be facilitated through appropriate forest management strategies. These findings may be very useful for the development of management and conservation measures for the European populations of this species in the southeast of the Iberian Peninsula, as well as improving the resilience of these forests to drought events and climate change. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/jmmp83/Forest_Tetraclinis_data, accessed on 1 March 2021.