Predicting suitability of forest dynamics to future climatic conditions: the likely dominance of Holm oak [Quercus ilex subsp. ballota (Desf.) Samp.] and Aleppo pine (Pinus halepensis Mill.)

Abstract

Key message

Composite logistic regression models simulating the potential effect of global climate change on forests dynamics in the southern Iberian Peninsula identify Holm oak [ Quercus ilex subsp. ballota (Desf.) Samp.] and Aleppo pine ( Pinus halepensis Mill.) as the chief beneficiaries of the anticipated environmental shifts, whereas other oak species and conifers suffer a decline.

Context

The ten most important tree species (five oaks and five conifers) in Southern Spain were selected for the study. The study area, corresponding to the region of Andalusia, is located in an interesting position between Central European and North African climates. The territory also exhibits the most extreme patterns of rainfall in the Iberian Peninsula.

Aims

This study aims to model the potential distribution of the ten species in response to climate change, in several time periods, including the present and two future twenty-first century dates.

Methods

The potential distributions within the different scenarios were simulated using logistic regression techniques based on a set of 19 climate variables from the WorldClim 1.4 project. The scenarios were drawn from the RCP 2.6 and 6.0 in the CCSM4 Global Circulation Model. The resolution of the output maps was 30 arc-seconds.

Results

The simulation predicted increased distribution areas for Q. ilex and P. halepensis under the four future scenarios as compared to present. The eight remaining taxa suffered a severe retraction in potential distribution.

Conclusion

Global climate change is likely to have a significant impact on forest dynamics in southern Spain. Only two species would benefit to the detriment of the others. Logistic Regression is identified as a robust method for carrying out management and conservation programmes.

1 Introduction

The last few decades have seen an increase in multidisciplinary approaches to forest dynamics. Land use, paleoecology, modelling and genetic diversity are just some of the disciplines which have been brought to the study (Alix-Garcia et al. 2016; Conedera et al. 2016; Vacchiano et al. 2016; Feurdean et al. 2017). Forest dynamics depend on biotic (species interactions) and abiotic (climate, orography, pedology, soil chemistry, water content, etc.) processes, which modify stand structures on both temporal and spatial scales (Frelich 2002). In Europe, anthropogenic activity has been a significant driver of forest dynamics since before the Common Era (BCE) (Watson 1996), and this has led to a decline in species in some cases, although some pioneer and hardy species have been favoured (Nicolaci et al. 2014; Conedera et al. 2016).

Potential shifts in the distribution of tree species and their communities as a result of global climate change are of great concern to scientists and researchers. According to the IPCC (2014), mean temperatures are set to rise and rainfall patterns are likely to alter. One result of these changing conditions is that the Mediterranean Basin could undergo a decrease in rainfall and consequently suffer an increase in aridity (Schröter et al. 2005; Klausmeyer and Shaw 2009; Planton et al. 2012; Vitale et al. 2012). This scenario would make the area unsuitable for those species with a narrow ecological range, whilst those with a wider tolerance would benefit. Thermophilous/xerophilous vegetation assemblages are likely to increase and exotic species could establish.

Species Distribution Models (SDMs) are used to simulate future habitat suitability for the purposes of species management, conservation biology, biogeography and climate change research (Guisan and Zimmerman 2000), this latter focussing on species distribution in the context of global climate change (Morán-Ordóñez et al. 2016; Thorson et al. 2016).

The main goal of the work presented in this paper is to evaluate the forest dynamics of ten tree species in southern Spain using a Logistic Regression (LR) approach. In total, five scenarios across three separate time periods were investigated: the current period, as reference, and four forecasted scenarios from the Representative Concentration Pathways (RCP) in the CCSM4 Global Circulation Model (GCM), specifically RCP 2.6 and RCP 6.0, each over two periods; 2041–2060 and 2061–2080 (called 2050 and 2070 in the whole text, respectively).

2 Material and methods

2.1 Study area

The study area was the region of Andalusia (Fig. 1) in Southern Spain. It displays a wide range of environmental and topographical conditions, including both the driest and wettest sites in Spain. The collision between the African and European tectonic plates in the Tertiary caused the formation of the Baetic range, mainly made up of limestones from marine sediments (Rosenbaum et al. 2002). The average elevation is high, with Mulhacén, at 3479 m a.s.l., the highest peak in the Iberian Peninsula. Another significant range is the Sierra Morena, which extends along northern Andalusia, and was created by the Hercynian folding (Palero et al. 2003). It is characterised by acid soils, principally quartzite, schist and slate. The average elevation of this range is much lower than the Baetic range, rising no higher than 1300 m a.s.l. in the eastern extreme of the range. The two ranges are separated by the Guadalquivir depression, which is characterised by low elevation and a fertile smooth surface formed by marls (Wheeler 1996; Villa Díaz 2013). These traits result in a largely agricultural landscape with vestigial patches of vegetation. In view of these features, the study area is especially suitable for studying forest dynamics under the impact of expected climate change in the twenty-first century.

Fig. 1

Location of the study area

As can be seen in Fig. 1, there is a west-east gradient. Generally, the western area is dominated by oaks, while the more easterly area, with higher orography, is better suited to conifers like pines, and even one fir species (Abies pinsapo Boiss.). Extreme rainfall patterns are also peculiar to the study area. Both the driest (subdeserts) and the wettest sites in the Iberian Peninsula are found here (de Castro et al. 2005). The huge area (around 87,697 km²) and the abovementioned characteristics make this an interesting area to model and simulate potential species distribution.

2.2 Test species

The evergreen oaks Quercus ilex subsp. ballota (Desf.) Samp. and Q. suber L. are common species throughout the Iberian Peninsula. The former is widespread due to its high degree of ecological plasticity. It can grow in different substrates and tolerate xeric conditions (López González 2007). The latter is especially prevalent in the southwest of the Iberian Peninsula, as it prefers acid soils in thermic areas with a minimum annual rainfall (Costa Tenorio et al. 2005). They are most typically found scattered across savannah-like pasturelands known as Dehesas in Spanish and montados in Portuguese (Smit et al. 2009; Arosa et al. 2017). By contrast, deciduous oaks are ecologically stricter, although Quercus faginea Lam. can often be found with evergreen oaks in mixed forests in southern Spain (Muñoz Álvarez 2010). Finally, Quercus pyrenaica Willd. and Q. canariensis Willd. are the most water dependent and least common oaks (Valle 2004) in this area.

High mountain pines such as Pinus nigra Arnold subsp. salzmannii (Dunal) Franco and P. sylvestris L. var. nevadensis H. Christ can be found in the higher reaches of the Baetic range. Pinus nigra can withstand harsher xeric conditions than P. sylvestris (Pasho et al. 2012) and is mainly found inland, whereas the latter occurs nearer the coast. Abies pinsapo Boiss. is a relict fir which grows at a mean elevation slightly lower than the two pines above, although it can reach almost 1900 m a.s.l. (López-Tirado and Hidalgo 2014). The main area of distribution of this conifer is the Straits of Gibraltar, an area coming under the influence of the Atlantic rainfall fronts and humidity. Exploited by humans over the centuries, A. pinsapo has benefitted from conservation programmes but since 1950 which have increased its distribution (Ceballos and Ruiz de la Torre 1979; Ruiz de la Torre et al. 1994; Linares and Carreira 2009; Linares et al. 2009), and nowadays covers 36–50 km² (Fernández-Cancio et al. 2007; Linares et al. 2010). Finally, Pinus halepensis Mill. and P. pinaster Aiton grow at elevations below the high mountain conifers. The latter inhabits an intermediate belt, coming into contact with P. nigra and P. sylvestris above and P. halepensis—the most thermophilous pine—below (López González 2007).

2.3 Methodological remarks

Logistic Regression is a type of Generalised Linear Model (GLM) which uses the presence/absence of data to hypothesise a formula that can be applied to the whole study area for both present and forecasted periods (Hidalgo et al. 2008). In this study, five oak species (two evergreens and three deciduous) and five conifers (four pines and one fir) were modelled at high resolution (Table 1). The potential distribution maps obtained can be used to study the dominant areas of species and the overlapping areas between two or more of them (Siles et al. 2010; Bonthoux et al. 2013; Pottier et al. 2013; Mikolajczak et al. 2015; Chai et al. 2016).

Table 1 Studied taxa by group and abbreviation

Taxa were studied using SPSS Statistics 20 [IBM, 2011], following the forward conditional method. This is a very useful stepwise technique in which the most statistically significant variable is chosen in the first step, then the second most significant, and so on until the least significant—or correlated—ones are rejected. All 19 climate variables from WorldClim 1.4 project were used (Hijmans et al. 2005) at a resolution of 30 arc-seconds (Table 2). Tree species distribution was retrieved from the regional Environmental and Water Agency (Agencia de Medio Ambiente y Agua, AMAYA). Models were verified by the Area Under the Curve (AUC) technique. A total of five periods were studied. First, the current period was considered from a dataset which ranged from 1960 to 1990. Then, four forecasted scenarios were processed with the downscaled and bias-corrected CCSM4 GCM, downloaded from the WorldClim 1.4 project. The original data are available at the CMIP Phase 5 home page (http://cmip-pcmdi.llnl.gov/cmip5). Two emission scenarios based on RCP were each considered in two distinct time periods, as follows: (i) RCP 2.6 (2050), (ii) RCP 2.6 (2070), (iii) RCP 6.0 (2050) and (iv) RCP 6.0 (2070). These scenarios provide time- and space-dependent trajectories of concentrations of greenhouse gases and pollutants resulting from human activity. Changes in land use, but with different socioeconomic assumptions, are also included. The first two scenarios assume an impact of 2.6 W/m², whilst the latter two assume an impact of 6.0 W/m². For further details, the reader is directed to the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC 2014).

Table 2 Climate variables used in the Logistic Regression models

2.4 Management of the potential distribution maps

Raster files showing potentiality ranging from 0 to 1 (maximum potentiality) were standardised, using ArcGIS 10 to process the data and gather all the raster data into one vector point shapefile. This was managed with the selection by attributes tool. Specific records were selected by means of different formulae. Cells which had a value under 0.5 for all the studied taxa were dismissed. Dominance of a species was detected when its score was larger than the others by at least 0.2. Cells that showed codominance, i.e. two species with higher potentialities than the others but within 0.2 of each other, were also identified (Mikolajczak et al. 2015). The number of such combinations involving two codominant taxa (considering all 10 species included in this study) was 45. All of these were scrutinised and the 12 most representative combinations were displayed in the final map. A new field in the vector file database was used for registering each dominance or codominance detected. This procedure was carried out for all five periods.

3 Results

The single most significant finding from the collation of the potential distribution of the ten tree species is that Q. ilex and P. halepensis gradually spread over a wider territory during the five periods. As a result, the distribution maps become simpler as the time period progresses (Figs. 2, 3 and 4). This trend also affects the patterns of codominance among the taxa, which tend to decrease, with the exception of Q. ilex-P. halepensis in all forecasted scenarios, Q. pyrenaica-P. nigra in 2050 and P. nigra-P. halepensis in 2070. The verification of fitness of the models is displayed in Table 3. Suitability maps for each species are supplied in the appendixes.

Fig. 2

Dominant and codominant areas for the ten taxa in the current period

Fig. 3

Dominant and codominant areas for the ten taxa in the RCP 2.6 forecasted scenarios. Please see Fig. 2 for legend

Fig. 4

Dominant and codominant areas for the ten taxa in the RCP 6.0 forecasted scenarios. Please see Fig. 2 for legend

Table 3 AUC values for the tree species studied

Table 4 shows the predicted areas (km²) of dominance and codominance. The chief findings, expressed as a percentage of the study area, are as follows. With regard to oaks, Q. ilex almost doubles its area from 20.71% in the current period to 40.96% in the RCP 6.0 (2050) scenario. Other species remain dominant only in the east (the higher elevations of the Baetic mountains), on the southern coast (the Strait of Gibraltar and surrounding areas) and in the west (the Sierra Morena). It is in these latter two areas that the major stands of Q. suber are found, which, according to our results, suffers a gradual decrease from 7.71% in the current period to around 0% in the RCP 6.0 (2070) scenario. In between these two areas of dominance, there is a current Q. ilex-Q. suber codominance, which decreases from its present 19.25 to 0% in the RCP 6.0 (2070) scenario. Quercus faginea is dominant in relatively small zones within the study area. It is worth noting that this species exhibited a larger predicted area in codominance with Q. ilex in the current period. Nonetheless, in the future scenarios, Q. ilex-Q. faginea is drastically reduced. In the case of Q. pyrenaica, an irregular trend can be noted, which starts from a situation similar to that for Q. faginea. In the current period, Q. pyrenaica displays dominance in only 0.03% of the whole study area. In general, the predictions for Q. pyrenaica dominance show a smaller area, although codominance with Q. ilex, and especially P. nigra, actually increases in some of the forecasted scenarios. The final oak, Q. canariensis, shows 0.21% dominance and 4.29% Q. suber-Q. canariensis codominance for the total study area in the current period. However, its predicted area is gradually reduced to 0% in the RCP 6.0 scenario, making it was the most negatively impacted of all the species.

Table 4 Predicted area (km²) of dominance and codominance between target species. The data has been computed with values over 0.5 of probability, but unsuitable areas which show the surface under 0.5 of probability for all the species

With respect to the conifers, these were generally distributed in the eastern reaches of the study area. The endangered A. pinsapo shows good potential distribution in the current period, including codominance with oaks (0.96% of the whole area), chiefly in the southern coastal mountains. However, the forecast for RCP 2.6 (2070) predicts zero dominance and codominance for this fir. Pinus nigra also experiences a significant decrease in the RCP 6.0 (2070) scenario, especially notable in the inland north-eastern mountain ranges, whereby its 9.48% dominance in the current period is predicted to be reduced to 0.47%. On the other hand, P. nigra-P. halepensis codominance increases in both RCP (2070) scenarios. Pinus sylvestris was found to have limited dominance in the current period in comparison with the other species, with only a P. sylvestris-P. nigra codominance being worthy of note (2.55% of the study area). Again, forecasts for all scenarios predict the disappearance of all current stands. P. halepensis is the second most dominant conifer in the current period (1.10% of the total study area). The most significant finding with respect to this pine is its potential to increase its dominance to 24.08% in the RCP 6.0 (2070) scenario. Moreover, a high number of significant codominances is also predicted, specifically with the three target taxa of Q. ilex, P. nigra and P. pinaster. Finally, with regards to P. pinaster, almost no dominant areas were found over the five periods studied, with only codominances with P. nigra and P. halepensis showing any significance.

4 Discussion

4.1 Initial remarks

This paper describes an SDM-based distribution simulation of ten species of oak and conifers in southern Spain in the context of global climate change. Although oaks are largely distributed throughout Spain, conifers may become dominant especially in mountain areas such as the Baetic range (Costa Tenorio et al. 2005).

Human-induced disturbances in Spain have taken place over many centuries, especially in afforestation programmes involving pines (Montero 1997; Zavala and Zea 2004; Muñoz Álvarez 2010; González-Muñoz et al. 2014). Forest fires and the practice of grazing in pine forests hinder the establishment of late successional species like oaks (Naudiyal and Schmerbeck 2016). In terms of altitudinal range, oak communities are generally found beneath conifers (Costa-Tenorio et al. 2005). Nevertheless, the disturbances like the above can lead to unclear situations in the study area, such as Quercus faginea subsp. alpestris (Boiss.) Maire growing at higher altitudes than A. pinsapo (Blanca et al. 2009). Thus, it is important to know the degree of alteration induced by humans in a given territory.

Following Swets (1988), the AUC scores (see Table 3) provide a high degree of confidence. It is important to highlight that the results are based on the actual niches of the species for the study area. Although biotic processes, such as species interactions, constitute relevant factors in this kind of study, they are not easy to be obtained and applied. Bioclimatic variables have been used to predict suitability areas in both present and future scenarios.

4.2 Modelling oak forest dynamics

The taxa with the widest predicted distribution in the study area are Q. ilex and P. halepensis. The former is a species with great ecological plasticity, as testified by the many published papers detailing the adaptations and strategies of this versatile tree (Morales et al. 2002; Corcuera et al. 2004; Sardans et al. 2006a; Galván 2012; Leiva and Díaz-Maqueda 2016). The latter, P. halepensis, is the most thermophilous of the conifers studied, and is able to withstand xeric conditions and grow on gypsum areas. By contrast, the other species suffer a decrease in predicted area. They can be considered more vulnerable in terms of specific traits, and are thus unable to achieve the high degree of adaptation displayed by Q. ilex and P. halepensis. For instance, Q. suber is less adapted to the climatic conditions of southern Portugal than Q. ilex (Martins et al. 1999), and is less resistant to xeric conditions (Peguero-Pina et al. 2009). It favours thermophilous areas with high annual rainfall and acid soils (López González 2007), traits which account for Q. suber enduring only in two main patches in the western reaches of the study area under the environmental conditions displayed in our models.

In the current period, Q. faginea shows only a small area of predicted dominance, although in co-occurrence with Q. ilex it has a larger distribution. In both cases, dominance and codominance, the forecasted scenarios predict a reduction in area. Differences between seasonal growth in Q. ilex—late summer and autumn—and Q. faginea—spring—indicate a greater importance regarding hydraulic conductivity for the former (Albuixech et al. 2012). Q. ilex can grow in less favourable periods, as it is less water demanding than Q. faginea. This trait is also relevant to Q. pyrenaica and Q. canariensis (Muñoz Álvarez 2010). Quercus pyrenaica can be found scattered throughout the study area. It occupies terrain with high annual rainfall, generally on north-facing slopes. It is the oak that grows at the highest elevation in the study area—2100 m a.s.l. (López González 2007). According to our models, this species would undergo an upward migration of the treeline as a consequence of climate change. The Sierra Morena stands would be lost due to the lack of higher mountains, whereas in the Baetic Range upward migration might occur. At the higher sites, Q. pyrenaica would compete with P. nigra, as codominance is predicted by the model. By contrast, Q. canariensis is clustered in the area around the Straits of Gibraltar, an area considered one of the wettest in the Iberian Peninsula. It co-occurs with Q. suber, although Q. canariensis usually inhabits riverbanks and humid depressions (Urbieta 2008). Table 4 shows how the predicted area for Q. suber-Q. canariensis gradually decreases. The low response by Q. canariensis seedlings to shifts in water availability has been noted (Quero et al. 2008), and hence, the expected increase in aridity would lead to a significant decrease in this species. A sizeable area of dominance by Q. canariensis is noted only in the RCP 2.6 (2050) scenario, which can be accounted for by the corresponding decrease predicted for Q. suber.

4.3 Modelling conifer forest dynamics

Abies pinsapo is an endangered relict fir, which in the course of the twentieth century has suffered afforestation in some areas (Álvarez 2013), as the management of stands is leading to a process of mortality (Abellanas et al. 2016). Abies pinsapo shares with Q. canariensis both a high demand for water and a restricted distribution, and consequently follows a similar trend to the latter in our LR model, although with more drastic consequences. This is consistent with the sensitivity of marginal Iberian forests to variations in environmental conditions (Génova Fuster 2013).

The high mountain pines, P. nigra and P. sylvestris, suffer a gradual and a total decrease, respectively. Pinus nigra undergoes an upward migration in the same way as Q. pyrenaica. The potential distribution of P. nigra stands extends over mountain ranges up to 2400 m a.s.l. whilst the predicted areas for P. sylvestris reach higher ranges (up to 3479 m a.s.l.). By contrast, P. sylvestris would be unable to find suitable areas at the highest sites due to edaphic conditions not considered here. In both cases, P. sylvestris and P. nigra show patterns of adaptation to climate (Vizcaíno-Palomar et al. 2016), providing evidence of genetic variation that could smooth the potential decrease in distribution.

Of the four pines studied, P. halepensis is considered the most resistant to drought, undergoing the least defoliation under such stress, whilst P. sylvestris undergoes the greatest (Sánchez-Salguero et al. 2012). Our results highlight the extensive ability of P. halepensis to co-occur with other taxa, such as Q. ilex, P. nigra and P. pinaster. The Q. ilex-P. halepensis distribution is the one that achieves the largest distribution in comparison with the other codominant combinations over the periods studied (see Table 4), reinforcing the assertion that P. halepensis is a species which is tolerant of many climates and soils (Fernández-Ondoño et al. 2010). In a mixed Q. ilex-P. halepensis community, P. halepensis is more competitive, taking shallow water whilst Q. ilex uses deep-soil water (del Castillo et al. 2016), enabling the community to become, according to our results, the largest in future scenarios. However, Q. ilex can be more competitive than P. halepensis in certain traits, such as showing a greater ability to extract phosphorus from the soil, and having a higher root biomass and nutrient allocation (Sardans et al. 2006b). Finally, P. pinaster is autochthonous to the Iberian Peninsula and widely planted for timber production (Garcia-Gonzalo et al. 2011; Fernandes et al. 2016), such that its current distribution goes beyond its original extension (Aguiar et al. 2007). Most of the stands planted in the Sierra Morena have led to decrepit specimens, due to both unsuitable sites and poor management (Muñoz Álvarez 2010). The intermediate characteristics of this species with respect to the other conifers would seem to explain why P. pinaster exhibits almost zero dominance in the predicted scenarios, whilst showing codominance with P. nigra and P. halepensis at the upper and lower ranges of its elevation, respectively.

4.4 Conclusions

Each of the taxa studied, together with its specific floristic composition, constitutes the climax of the vegetation series in the study area (Rivas-Martínez 1987; Valle 2004). In this respect, SDM represents a powerful tool for studying single species and/or forest dynamics (Alba-Sánchez et al. 2010; Attorre et al. 2013; Bede-Fazekas et al. 2014; Vessella et al. 2015; Al-Qaddi et al. 2016; López-Tirado and Hidalgo 2016a, 2016b). The study presented in this paper used an LR model to predict the forest dynamics of the most significant trees in southern Spain at specific points in the twenty-first century. The main conclusion of this study is that Q. ilex and P. halepensis could be favoured by predicted environmental shifts in the coming decades, and in consequence, mixed stands of Q. ilex-P. halepensis codominance could spread across the study area. By contrast, the other eight taxa are likely to contract in terms of potential distribution. This is especially the case with Q. canariensis and A. pinsapo, whose current potential areas could be lost within a matter of years. It is also worth noting how the expected harsher environmental conditions cause an increase in unsuitable areas (cells with probability < 0.5) for all the species in the periods studied. In general, tree species cannot colonise beyond a range of 20 km over the course of a century. In this regard, the results reported here are highly relevant to management and conservation programmes, and in particular to the Spanish and Andalusian regional governments with respect to carrying out strategies in both the economic and ecological spheres.