Effects of climate change on the distribution of Pinus sylvestris L . stands in Spain . A phytoclimatic approach to defining management alternatives

This paper presents some contributions on the possible effects of climatic change on the distribution of Pinus sylvestris L. stands in Spain. We studied the phytoclimatic status of Scots pine in current climate conditions (period 1951-1999) and in projected future climate conditions (2050). The phytoclimatic diagnosis followed a modified version of the Allué-Andrade phytoclimatic system. This calculation determines potential areas for Scots pine totalling 8,444,700 ha in current climate conditions and only 1,269,100 ha in 2050 conditions. The shrinkage of the area is especially pronounced in the southern half of Spain, where the model predicts that the species will practically disappear from the Baetic mountain ranges and from a major part of the Central System. In phytoclimatic terms, the maximum values of suitability correspond chiefly to the oroborealoid transitional to nemoral subtype VIII(VI)1 in current conditions, but in projected future conditions the maximum values correspond to areas currently assigned to oroarticoid transitional to oroborealoid subtypes X(VIII). The lowest viability scores are found in the southern half of the Iberian Peninsula and in Mediterranean transitional and nemoromediterranean subtypes. The results also suggest that stands of Pinus sylvestris will migrate upwards and will encounter a serious limitation in the scant availability of high mountain areas other than in the large northern massifs like the Pyrenees and the Cantabrian Cordillera, which accounts for their scant capacity to colonize new areas (2,134 km as opposed to extinction over 73,890 km).


Introduction
The use of various predictive modelling techniques to forecast the future effects that climate change may produce in the habitats of forest species is a field of growing interest to researchers and forestry managers, especially in the Mediterranean area (Hansen et al., 2001;Bakkeness et al., 2002;Pearson and Dawson, 2003;Gracia et al., 2005;Petit et al., 2005;Del Barrio et al., 2006;Benito et al., 2008;Franklin, 2009).As a result of the efforts made in modelling habitat suitability, predictive techniques have become more numerous and have been improved in recent years (Guisan and Zimmerman, 2000;Pearson, 2007;Elith and Leathwick, 2009;Elith and Graham, 2009).
Forests are particularly sensitive to climate change, because the long life-span of trees does not allow for rapid adaptation to environmental changes (Davis et al., 2005;Kremer, 2007;Lindner et al., 2010).Species that lie at the limits of their natural range of distribution appear to be particularly sensitive to the effects of climate change (Crawford, 2008;Holtmeier, 2003).One such case is Scots pine (Pinus sylvestris L.), the pine with the largest natural area and the most widespread in Europe and Asia, stretching from east to west, from eastern Siberia to Galicia, and north to south, from Scandinavia to the south of the Iberian Peninsula (Sierra Nevada) where the most southerly specimens are found (Ruiz de la Torre, 2006).The area covered by Scots pine in Spain is 1,280,000 ha (Montero et al., 2008).About 605,200 ha are natural stands and 674,800 ha planted (Cañellas et al., 2000).
A naturally-occurring species in taiga forests in northern Europe and Asia, its northern area of distribution is more or less continuous, on mesetas and plains, whereas in the south, in the Mediterranean basin, it is becoming increasingly fragmented and confined to mountain areas (Barbéro et al., 1998).Because of the considerable economic and ecological wealth of this species in Spain (Montero, 1994;Cañellas et al., 2000) the vulnerability of natural and artificial stands to climate change is a matter of particular interest and concern, especially considering that the mountain ecosystems of southern Europe may be among the worst affected by such change (Cubash et al., 1996;Watson et al., 1997).
Species distribution models (SDMs) are empirical models relating species occurrence to environmental variables based on statistical or other response surfaces.SDMs have become an important component of conservation planning in recent years and predicte global change impacts on plant species distributions (Thuiller et al., 2008).A wide variety of modeling techniques have been developed for this purpose (Guisan and Thuiller, 2005;Elith and Leathwick, 2009).Some of the earliest SDMs defined the hyper-rectangle that bounds species records in multi-dimensional environmental space, weighting each predictor equally (Box, 1981).A number of alternative modeling algorithms have been applied to classify the probability of species' presence (and absence) as a function of a set of environmental variables (Pearson, 2007): Gower Metric (Carpenter et al., 1993), Ecological Niche Factor Analysis (Hirzel et al., 2002), Maximum Entropy (Phillips et al., 2006), Genetic algorithm (Stockwell andPeters, 1999), Artificial Neural Network (ANN) (Pearson et al., 2002;Benito et al., 2006), Regression (e.g., Lehman et al., 2002;Elith and Graham, 2006;Leathwick et al., 2006;Elith and Leathwick, 2007) such as generalized linear model (GLM), generalized additive model (GAM), boosted regression trees (BRT), multivariate adaptive regression splines (MARS), and Multiple methods BIOMOD (Thuiller, 2003).
This paper reports the assay of a methodology to evaluate the possible effects of climate change on the distribution of Pinus sylvestris stands in the Iberian Peninsula.In particular the authors have attempted to answer a number of questions, such as the likely future area of distribution, what its altitudinal range of occurrence will be, what part of the existing potential area will continue to be potential for this species in the climate conditions of the future, what part will cease to be so and what currently unsuitable areas will become suitable in the future.

Material and methods
The phytoclimatic system used is based on the models of Allué-Andrade (1990 and1997) modified by García-López and Allué (2003).This phytoclimatic system was chosen for the present study because it makes it possible not only to f it a station into a previouslydefined phytoclimatic category in qualitative terms but also to quantify the adjustment of the station to that category or phytoclimatic type, and likewise all the other types in the system, using relative position coordinates and phytoclimatic distances, between these and relative to factorial phytoclimatic ambits.
Possibly one of the most original features of the cited SDMs is that they do not make direct use of an environmental space formed by the predictors (factorial space) but a Euclidean phytoclimatic scalar space, whose chief advantages over the former are the possibility of defining phytoclimatic distances in it as dual measures of proximity and potentiality, by first defining a phytoclimatic adjustment function with physiological rather than probabilistic implications.The fact that the model is more a biological than a statistical one means that it is possible to work with less transformations of the initial variables and to minimize the wellknown black box effect of some SDMs whose complex algorithms make it difficult to understand their functioning and interpret their results in biological terms.
The fact that absence data are not required (Franklin, 2009) is an advantage in such spoiled environments as Mediterranean ones, where the absence of a species may very probably be due human intervention.Not only the initial predictors but also the specific values of the predictors present different weights in the calculation through their characterizing powers, and their contribution to the f inal result can be readily assessed in the model's calculation matrix.By numerically quantifying the phytoclimatic distance of a point from the set of ambits of real existence of the species considered in the calibration of the model it is possible to continually assess all their levels of phytoclimatic suitability as fuzzy sets (Zadeh, 1965).
The model used in this study is therefore much like a SDM with a physiological rather than a statistical response surface (Liu et al., 2009).The fact of not depending on the frequency with which the target species occurs is another advantage, since according to Walter (1970), the physiological and ecological optima (of frequency) do not usually coincide, and hence a species' distribution in terms of frequency or density is not a good indicator of its suitability -in this case phytoclimatic.
Pinus sylvestris was assigned an autoecological factorial ambit determined from the 4,942 sampling points in the II National Forestry Inventory (NFI) and the corresponding provinces now available from the III NFI (see Villanueva, 1990, andBravo et al., 2002, for details about NFI).Sampling points were selected setting apart all registers in which pine occurred naturally as the first dominant tree in the formation.The appropriate values of the factors in Table 1 were obtained by means of cross-referencing with the regionalized climatic data base of Gonzalo ( 2008) for the compendium period 1951-1999.
According to García-López and Allué (2003), the borderline of each ambit can be defined in very close correspondence with the cluster of points in 12-dimensional factorial hyperspace by calculating a convex hull that will convert it to a hyperpolyhedron and can be  projected on to planes formed by pairs of factors in order to perform the specif ic calculations for the phytoclimatic model.The 12 factors used (Table 1) are the original ones from the model (Allué-Andrade, 1990 and1997) except for PV, which replaces the probable duration of frosts and offers greater diagnostic power than the latter (calculated as the number of months in which TMMFi ≤ 0°C) in terms of mean characterizing power (García-López and Allué, 2008).
By running diagnostics on a point with the autoecological phytoclimatic system constructed using the programme CLIMATFOREST 1.0 (García-López and Allué, 2009c), in each case we can generate a scalar e psy for that station's adjustment to the factorial ambit of Pinus sylvestris, which is defined by the convex hull and includes the point analysed, with e Psy ≥ 0 and ≤ 1.Each scalar functions as an index of relative phytoclimatic suitability of a forest species with respect to the optimum.
For the purposes of this article, «phytoclimatic suitability» (Allué, 1996) means the degree to which a site is suited to host certain taxa or syntaxa, principally in terms of staying power (self-regenerating capacity), ability to compete with other species and resistance to diseases.
Any approach to ecological modelling has little merit if the predictions cannot be, or are not, assessed for their accuracy using independent data (Verbyla and Litaitis 1989).It is generally accepted that robust measures of prediction success make use of independent data, i.e. data not used to develop the prediction model (Pearson, 2007).The simplest, and most widely used, measure of prediction accuracy is the number of correctly classified cases (Fielding and Bell, 1997;Anderson et al., 2003).
If data are partitioned the size of training set must decrease and this can reduce model accuracy.There is a trade-off between having a large test set that gives a good assessment of the classifier's performance and a small training set which is likely to result in a poor classifier.In our case a random sample constituting 30% of the initial sampling parcels was reserved and not used in the construction of the model, following Huberty (1994).These parcels were used exclusively to determine the model's diagnostic capacity in terms of percentage of correct predictions of the eventual dominant species in the forest formation.The results were highly satisfactory, with correct determination of 97.5% of the components of the validation sample.
The phytoclimatic system thus constructed was applied to two factorial data bases for the Iberian Pe-ninsula at a resolution of 1kmx1km, using a specific large-scale calculation model derived from CLIMAT-FOREST 1.0.The first of them, representing the current phytoclimatic conditions, was calculated on the basis of regionalized climatic variables (Gonzalo, 2008) for the period 1951-1999.The second factorial data base, representing future climatic conditions, was calculated by taking the first base and applying the relative coefficients of variation assigned for the year 2050, again at a resolution of 1 km × 1 km to the Iberian Peninsula, derived from the Global Climate Model CCCMA (Canadian Centre for Climate Modelling and Analysis) for the A2 emissions scenario (Hijmans et al., 2005).
The IPCC Fourth Assessment Report (Solomon et al., 2007) applied different scenarios based on the Special Report on Emission Scenarios (Nakicenovic et al., 2000).The SRES scenarios are grouped into four scenario families (A1, A2, B1, and B2) that explore alternative development pathways, covering a wide range of demographic, economic, and technological driving forces and resulting GHG emissions.The emission projections are widely used in the assessments of future climate change, and their underlying assumptions with respect to socioeconomic, demographic and technological change serve as inputs to many recent climate change vulnerability and impact assessments.
The A2 scenario family describes a very heterogeneous world.The underlying theme is self-reliance and preservation of local identities.Fertility patterns across regions converge very slowly, which results in continuously increasing global population.Economic development is primarily regionally oriented and per capita economic growth and technological changes are more fragmented and slower than in other storylines.This scenario was chosen for our study because it predicts medium-to-high warming over all the available scenarios.A2 describes an increasing global population, but with slower economic growth than in the other scenarios.
In the case of the future factorial data base, the points inside the present ambit of existence of Pinus sylvestris, previously defined as a convex hull, were selected in order to find the species' potential phytoclimatic area of existence.

Results
From the phytoclimatic diagnosis of the factorial base of current climatic conditions and the factorial base of future climatic conditions, each containing around 500,000 geographical points, it was possible to generate in the first place two phytoclimatic diagnostic data bases in which each point was assigned an indicator of its internal (genuine) position, proximal external (analogous) position or remote external (disparate) position with respect to the factorial ambit of existence of Pinus sylvestris (convex hull in the 12-dimensional factorial space) and a corresponding Suitability Index (e Psy ).
According to these results, the current potential phytoclimatic area of Pinus sylvestris, considered as the set of points inside the convex hull formed in the 12-dimensional factorial hyperspace is 84,447 km 2 (Fig. 1a).The potential area of high viability, that is the area presenting values of e Psy > 0.7, is 43,972 km 2 (Fig. 1b).
In the future climate scenario considered, this potential area is reduced to 12,691 km 2 (Fig. 2).
In phytogeographic terms, Table 2 shows the results for the Mean Suitability Index (e Psy ) of Pinus sylvestris as a function of the altitudinal range considered for the current phytoclimatic scenario and for the future phytoclimatic scenario.According to these results, highlysuitable pine stands (e Psy > 0.70) are currently found from an altitude of 1,200 m upwards, while in the future scenario it is from 1,500 m upwards.In addition, the altitudinal range of maximum phytoclimatic suitability is currently 1,500-1,800 m, while in the future scenario it is 1,900-2,300 m.Above 1,800 m altitude the future area of distribution will be more suitable for Scots pine than the present one.
In phytoclimatic terms, Table 3 shows the results for the Mean Suitability Index (e Psy ) of Pinus sylvestris as a function of the phytoclimatic subtype to which they belong, for the current phytoclimatic scenario and for the future phytoclimatic scenario.In the latter case we have given the results for the present subtype at each point and the future subtype at the same point.
According to these results, stands presenting high suitability (e Psy > 0.7) are currently situated in subtypes VI(VII), VI, VIII(VII), VIII(VI) 1 , VIII(VI) 2 and X(VIII), while in the future scenario they are situated at stations which are currently classified in subtypes VIII(VI) 1 , VIII(VI) 2 and X(VIII) but which in the future scenario would also be classif ied in subtypes VI(VII), VI, VIII(VII), VIII(VI) 1 and VIII(VI) 2 .In addition, the subtype of maximum phytoclimatic suitability at present is VIII(VI) 1 , whereas in the future scenario the stations of maximum suitability would be ones that are currently classified in subtype X(VIII) but in the forecasted future conditions would also enter subtype VIII(VI) 1 .Note also that pine stands situated in sub-types classified as less suitable, such as IV(VI) 1 and IV(VII), would no longer appear there in the future scenario.Table 4 attempts to answer a number of questions, such as how much of the current potential area of Pinus sylvestris will continue to be potentially suitable for this species in future climatic conditions, how much will cease to be so and what currently unsuitable areas will become suitable in the future.According to these results, 87% of the potential area (73,890 km 2 ) will cease to be suitable in the future and only the remaining 10,557 km 2 will continue to be suitable for these pine stands.Only 2,134 km 2 not currently suitable for Scots pine will be suitable for it in 2050.
The altitude and suitability data in Table 4 show that the areas lost in 2050 are generally the ones at lower mean altitudes within the potential phytoclimatic distribution area of the species (1,101 m), those at intermediate altitudes remain stable (1,448 m) and areas at high mean altitudes increase (16,41 m).Also, the areas that remain stable present greater average suitability than those that lose or gain.

Discussion and conclusions
The results indicate that climate change may have a strong impact on the distribution and phytoclimatic potentialities of Scots pine stands in Spain if the numeric predictions of the future climate scenario used materialize, given that the potential area for the species Our results appear to fall within a similar range to those of the few other studies of current potential areas for Scots pine in the Iberian Peninsula, particularly the most recent one (Benito et al., 2006), which estimates that potential area at between 32,300 km 2 and 103,800 km 2 using three predictive techniques based on automatic learning of physiographic and phytoclimatic data (learning machine).
In the case of future potential areas, Benito et al. (2008) estimate the area at 3,907 km 2 in 2050, well below our own figure of 12,691 km 2 .This difference in predictions is surely a result of the different models and techniques used to generate future scenarios; while in the case cited they considered the A2 emission scenario, they did so on the basis of the CSIRO and HadCM3 models, the second of which at least is more pessimistic than our model in its predictions of warming.The same authors predict a 99% reduction from the present area by 2080.
The shrinkage of the area is especially pronounced in the southern half of Spain, where the model predicts that the species will practically disappear from the Baetic mountain ranges and from a major part of the Central System.In fact it appears that some of the pine stands situated in the south-eastern Iberian Peninsula have already begun to show clear symptoms of decay, such as severe defoliation (Navarro et al., 2007), and our own results suggest that this trend will continue.The possible disappearance of Iberian Scots pine populations is particularly worrying considering that these populations are genetically distinct from the European varieties (Prus-Glowacki & Stephan, 1994;Prus-Glowacki et al., 2003;Robledo-Arnuncio et al., 2005), which means that the potential loss of genetic diversity is incalculable.Symptoms of decay in Pinus sylvestris stands also have been observed in the middleeast, in the Iberian mountains (Alquézar et al., 2008).
In practice the model corroborates the results reported by other authors concerning the general rise in the altitudinal limit of mountain treelines (e.g., Grace et al., 2002;Hotmeier and Broll, 2007;Bake and Moseley, 2007;Klanderud and Birks, 2003;Dullinder et al., 2004).But in any case these predictions should be treated with caution since the climatic causes of upland migration and enrichment of tree species may overlap with other causes such as abandonment of alpine livestock grazing (Gehring-Fasel et al., 2007).
The results also suggest that stands of Pinus sylvestris will migrate upwards and will encounter a serious limitation in the scant availability of high mountain areas other than in the large northern massifs like the Pyrenees and the Cantabrian Cordillera, which accounts for their scant capacity to colonize new areas (2,134 km 2 as opposed to extinction over 73,890 km 2 ).In García-López and Allué (2009b) we evaluated the enhancement of the competitive ability of fagaceae in medium-high altitudinal areas currently occupied by Scots pine in Castilla y León.Altitudinal migration will occur to a considerable extent at the cost of the area now identified as oroborealoid subtype X(VIII), which underlines how important such humid mountain areas will become as refuges for Scots pine populations, whose main direction of colonisation will be vertical.In any case, it must be taken into acount the complex and heterogeneous response of the high mountain ecosystems (Malanson, 2001).
In addition to these limitations we should cite the reservations expressed by Davis et al. (1998) or Pearson and Dawson (2003) regarding the construction of factorial ambits of existence in studies of this kind since all the possible situations for each species are not taken into account when constructing factorial spaces, particularly in such degraded forest environments as those found in the Mediterranean basin.It is worth noting in this respect that both the present and the future scenarios predict extensive potential phytoclimatic areas in the Cantabrian Cordillera, where the species today is only found in small residual enclaves, although there is a considerable amount of palaeobiogeographic evidence that its disappearance is relatively recent and due to human activity, chiefly the use of fire for purposes of animal husbandry (Costa et al., 1990;García-Antón et al., 2001;Rubiales et al., 2010) When dealing with limitations we cannot ignore the existence of strong biological factors which will surely modulate future predictions, for instance the capacity of species to spread, the connectivity among ecosystems fragmented by climate change, variation in patterns of inter-specific competition or the ability of species to adapt to the rate of change (Woodward, 1990;Thomas et al., 2001).
Another important issue is intrinsic genetic adaptative capacity of individuals (individual heterozygosity, acclimation, epigenetic response), communities or species, which can make changes necessary in the hulls as originally calculated (Bradshaw and McNeilly, 1991;Kremer, 2007;Lindner et al., 2010).
Also, predictions relating to the area of potential rather than real phytoclimatic distribution of Pinus sylvestris refer to an idealization of the species' capacity to survive long-term at a given location.Even if climate data for the year 2050 are used to make the predictions, these must be treated with caution as they do not reflect the probable state of those forests for that year but the state that may be expected if these new conditions persist over the long term.The mitigating and homeostatic properties of forest ecosystems will hold back change unless those condition should persist for long periods.
Also, in future we must take into account factors that are as yet little-known such as inter-species competition.Changes in competitive relations in Pinus sylvestris forests mean that this species will very probably behave differently from what was originally predicted.One aspect of particular interest is the question of invasion of some low-and medium-altitude pine forests by deciduous broadleaf species, which will assuredly alter the relationship of the pines with the surrounding environment (Amarasekare, 2003).This will affect the shapes of the convex hulls; in this connection Hutchinson (1957) defined the fundamental niche as the n-dimensional hypervolume where a species, in the absence of competition, is able to persist indefinitely (Araújo and Guisan, 2006).
As well as these limitations we must consider the limitations proper to future climate scenarios.We would stress the importance of further elaborating this methodology in the future so that the model takes into consideration more climate scenarios and thus provides a more comprehensive picture than the one offered in this paper.An important issue when considering adaptation responses to climate change is the uncertainty in the predictions of future climate.With more scenarios it would be possible to make progress in separating the variability in the predictions associated with the use of the model from the variability associated with uncertainty as to the future course of climate change.
Another practical limitation of the model derives from the territorial scope within which it is constructed.The area of distribution of Pinus sylvestris extends far beyond the Iberian Peninsula, and therefore the climatic hulls that have been constructed do not reflect the entire potential phytoclimatic ambit of the species (Petit et al., 2005).Improvement of the definition of the convex hulls will be a crucial issue in the future for enhancing the predictive efficiency of the model.
And another source of uncertainty is the future possibility of combinations of as-yet unknown factorial values, so that the model may determine that an ambit is not suitable for the target species when at a future date they may be shown to be compatible (Fitzpatrick and Hargrove, 2009).
In short, the thrust of the ideas set out in the foregoing paragraphs is that the shrinkages in area predicted in this study will most probably prove to be smaller in the event.
Among the manifold implications that the results of this study may have for management of Pinus sylvestris forests in Spain, we should particularly highlight those relating to future reafforestation projects.Our findings point to a need to conserve the considerable genetic wealth of Spanish Scots pine forests (Agúndez et al., 1992), especially the most southerly ones (Pinus sylvestris var nevadensis) in the context of a steep latitudinal gradient in the species (Oleksyn et al., 1998).The reasons for this are twofold: firstly, the fact that conditions in southern Spanish forests are more xeric than in northern forests and regeneration is more difficult (Castro et al., 1999 and2002) renders them more vulnerable to climate change.And secondly, these southern forms are better adapted to drought conditions (Alía et al., 2001) and hence are especially valuable when planning reafforestation initiatives in more northerly areas where future climatic conditions may come to resemble those prevailing today in southern massifs (Galera and Albertos, 1990).Thought should be given to the use of material from regions of origin other than those where the forestation is planned (Catalán, 1991).
. Calculated on the basis of the quotient As/Ah, where Ah is the humid area of the climodiagram (Pi curve above the Ti curve, i.e., 2Ti < Pi) and As is the dry area of the climodiagram (Pi curve below the Ti curve, i.e., 2Ti > Pi).Ti and Pi are the mean temperature and precipitation of the month i Duration of aridity in the sense of GAUSSEN, that is the number of months in which the Ti curve is above the Pi curve, i.e., minima of the month with the lowest mean temperature Average of the maxima of the month with the highest mean temperature Certainty of frost.Calculated as the number of months in which Ti ≤ 4°C Period of free plant activity, calculated as the number of months in which Ti ≥ 7.5°C, not counting periods where A > 0 Thermal oscillation.Calculated as TMC-

Table 1 .
Climate change and distribution of Pinus sylvestris in Spain331 Phytoclimatic factors used

Table 2 .
Altitudinal distribution of the mean Phytoclimatic Suitability Index (e Psy ) of Pinus sylvestris for the present and future climatic scenarios.Std indicates the Standard Deviation.Cursive letter indicates high-suitability altitudinal ranges (ePsy > 0.7) and bold letter indicates maximumsuitability altitudinal ranges

Table 3 .
Distribution by phytoclimatic subtype of the mean Phytoclimatic Suitability Index (e Psy ) of Pinus sylvestris for the present and future climatic scenarios.Std indicates the Standard Deviation.Cursive letter indicates high-suitability subtypes (e Psy > 0.7) and bold letter indicates maximum-suitability subtypes

Table 4 .
Climate change and distribution of Pinus sylvestris in Spain335 Shift of areas capable of sustaining Scots pines, between present and future climatic conditions