Assessing changes in species distribution from sequential large-scale forest inventories

Abstract

• Context

It is assumed that global change is already affecting the composition, structure and distribution of forest ecosystems; however, detailed evidences of altitudinal and latitudinal shifts are still scarce.

• Aims

To develop a method based on National Forest Inventory (NFI) to assess spatio-temporal changes in species distributions.

• Methods

We develop an approach based on universal kriging to compare species distribution models from the different NFI cycles and regardless of the differences in the sampling schemes used. Furthermore, a confidence interval approach is used to assess significant changes in species distribution. The approach is applied to some of the southernmost populations of Pinus sylvestris and Fagus sylvatica in the Western Pyrenees over the last 40 years.

• Results

An increase of the presence of the two species in the region was observed. Scots pine distribution has shifted about 1.5 km northwards over recent decades, whereas the European beech has extended its distribution southwards by about 2 km. Furthermore, the optimum altitude for both species has risen by about 200 m. As a result, the zone in which the two species coexist has been enlarged.

• Conclusions

This approach provides a useful tool to compare NFI data from different sampling schemes, quantifying and testing significant shifts in tree species distribution over recent decades across geographical gradients.

1 Introduction

The impact of global change on forests is emerging as a major concern for the twenty-first century society (FAO 2000). Forestry and conservation po1icy-makers need to understand how the distributions of species are affected by global change in order to tackle its effects on forests. Although it is assumed that global change is already affecting the composition, structure and distribution of forest ecosystems at different spatial and temporal scales, reported evidences of altitudinal and latitudinal shifts are still scarce (Peñuelas and Boada 2003).

Large-scale forest surveys such as the National Forest Inventory (NFI) can provide a valuable tool for monitoring changes in forest biodiversity and can be of great relevance in the conservation and management of natural resources. In the second half of the twentieth century, most developed countries started undertaking periodical NFIs covering the entire forest area. Although NFIs were primarily designed to estimate forest resources, they are increasingly being employed to assess the impact of global change on forest ecosystems (Thuiller et al. 2003). When permanent plots are inventoried sequentially, forest evolution can be analysed through direct comparison between inventories (Poljanec et al. 2010; Vilá-Cabrera et al. 2011); otherwise, species–environment relationships should be analysed from species distribution models (SDMs) (Guisan and Zimmermann 2000). One problem that arises when comparing species distributions is whether or not the changes are related to the different sampling schemes employed across successive NFIs. This shortcoming can be dealt with in a similar way to change-of-support problems by using block kriging techniques (Yoo and Trgovac 2011).

Environmental variables such as precipitation, temperature and elevation exhibit spatial dependence, which is partly responsible for the spatial pattern observed in vegetation distribution (Miller et al. 2007). In addition, spatial autocorrelation can also result from ecological processes involved in forest dynamics occurring at different scales, leading to spatial continuity in species distribution (Bellehumeur and Legendre 1998). As a consequence, spatial dependence in biogeographical data has recently been identified as an important area of SDM research, especially for species with a broad distribution, which are generally better modelled by including spatial autocorrelation in the model (Segurado and Araújo 2004). Universal kriging may be considered a type of spatial regression, incorporating spatial autocorrelation as well as relationships between environmental variables in the spatial prediction of species distribution (Montes and Ledo 2010). It may be important to consider this latter quality of universal kriging in SDMs as species distribution varies over environmental gradients. The altitudinal range where species with a broad latitudinal distribution appear varies as a result of a decline in temperature with increasing elevation and more northerly latitude. This variation in modelled relationships over space, known as non-stationarity, shows interdependence with spatial scale, so local and landscape approaches may reveal species–environment relationships that global models average out (Osborne et al. 2007; Miller and Hanham 2011).

The Iberian Peninsula, like other Mediterranean regions, constituted an important refuge for flora during the Quaternary glaciations. This fact, together with the current interaction of the Eurosiberian, Mediterranean and Alpine biogeographical zones, made the Iberian Peninsula a biodiversity hotspot where many species find the limits of their distribution; many of them highly fragmented in different mountain ranges. This is the case of Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.), two species widely distributed across Europe and which have some of their southernmost populations in the mountain ranges of the Iberian Peninsula. The most recent climatic records for the Mediterranean region reveal an increase in the number of heat waves and droughts during the last century (Alpert et al. 2008). This pattern, together with the important impact of human activity and the sensitivity of the ecosystems to contrasts in climatic conditions, makes the mountain ecosystems of the Iberian Peninsula particularly vulnerable to the effects of global change (Engler et al. 2011). It is expected that woody species will respond to changes in environmental conditions by migration towards suitable areas or through adaptation to the new conditions. Otherwise, they may decline or even become locally extinct (Aitken et al. 2008). Taking into account not only these premises but also the results of previous studies in European mountain forests (Lenoir et al. 2008; Poljanec et al. 2010), we hypothesize an upward shift in the distribution of Scots pine and European beech at regional scale over recent decades in the Western Pyrenees.

The aim of this study is to develop a new method based on universal kriging models to assess spatio-temporal changes in species distribution from large-scale sequential forest inventories regardless of the differences in the sampling schemes used. For testing the usefulness of this approach, we also aim to answer the following questions: (a) Is it possible to detect significant shifts in the spatial distribution of P. sylvestris and F. sylvatica in the Western Pyrenees over recent decades? (b) Can we quantify the extent of these shifts across geographical gradients?

2 Material and methods

2.1 Study area

This study focuses on the Western Pyrenees in the Navarra region of Spain between E.D. 50 Universal Transverse Mercator (UTM) 4.734 S–4.764 N and 642 E–685 W. The study area, which covers 741,60 km², is characterized by a steep altitudinal gradient, which ranges from approximately 600 masl at the southern limit to 2,400 masl at the northern limit. Due to its geographical location, where several biogeographical regions overlap, the region has a variable climate. At more northerly latitudes, where the highest altitudes are found and the oceanic influence of the Cantabrian Sea is less notable, the mean annual precipitation is 1,063 mm, and the mean annual temperature is about 8 °C, which reflects the typical continental characteristics of the Pyrenees. By contrast, in the mid latitudes, mean annual precipitation levels are higher (1,300 mm), as are the temperatures (10.5 °C). The southernmost area, however, is influenced by a sub-Mediterranean climate with moderate summer drought. Hence, the mean annual rainfall in this part of the study area is lower (1,044 mm), although the mean annual temperature is higher (11 °C) (Ninyerola et al. 2005).

Scots pine and beech form the naturally occurring uneven-aged forests that dominate the study area. In the Iberian Peninsula, the majority of Scots pine stands can be found in the more continental, south-central Pyrenees; becoming scarcer in the northwest (Costa et al. 1997), i.e. in the areas influenced by the Cantabrian Sea, where beech stands are more common. Beech and Scots pine occur in the montane altitudinal belt between 600 and 1,600 masl, with beech occupying higher and shadier slopes. In the lower montane altitudinal belt, P. sylvestris can be found alongside marcescent Quercus species like Quercus faginea Lam., Quercus pyrenaica Willd. and Quercus pubescens Willd., whereas above these forests, the sub-alpine belt (1,600–2,300 masl) is dominated by Pinus uncinata Ram. and the alpine belt (over 2,300 masl) is the domain of herbaceous, sometimes shrubby vegetation. Both species can co-occur in some zones, where European beech is found in more mesic conditions and Scots pine where more xeric conditions prevail.

2.2 NFI dataset

The Spanish NFI is a forest monitoring system covering the entire forested area of Spain. Inventories are undertaken approximately every 10 years with an intensity of one sampling point every 1 km². At each sampling point, trees are measured in concentric circular plots with increasing radii from 5 to 25 m. In the 5-m-radius plot, trees with a diameter at breast height (dbh) ≥ 75 mm and trees with 25 ≤ dbh < 75 mm (saplings) are measured, while the recruitment density is monitored (individuals of less than 1.30 m height or dbh < 25 mm). In the 10-m-radius plot, trees with dbh ≥ 125 mm are measured; in the 15-m-radius plot, trees with dbh ≥ 225 mm are measured; and in the 25-m-radius plot, trees with dbh ≥ 425 mm are measured. In the Navarra region, four cycles of the Spanish NFI have been completed. Although the NFI2 (1986–1995) and NFI3 (1997–2006) cycles have been used in numerous studies (Vilá-Cabrera et al. 2011), research involving the comparison of data from the four NFIs conducted since 1960 has never been undertaken. Since different sampling designs were used in the various NFI cycles, the number of plots and locations varied from one to another (Table 1).

Table 1 Spanish National Forest Inventory (NFI) cycles sampling periods, source of information used for forest area assessment, sampling period and number of plots measured in the study area in each NFI cycle

In order to analyse changes in the spatial distribution of species using the NFI databases, presence/absence indicator variables with a value of 1 if a species was present in a plot and 0 if it was absent were defined for both species. The presence of a species in a given plot was registered according to size, recruitment, saplings and trees with dbh ≥75 mm. Therefore, the presence/absence indicator variable is a regionalized variable that (for the xy location of the plot centre) will be 1 if there are individuals within the respective radius and 0 if not. The variable-radius sampling design of the NFI plots smoothes the sampling estimator surface and minimizes the variance of the estimator (Williams 2001).

2.3 Universal kriging

Universal kriging (UK) is a spatial regression procedure that incorporates the spatial autocorrelation in the estimation of a regionalized variable. In the universal kriging model, the value of the variable Z(s ₀) at location s ₀ is expressed as a polynomial of the auxiliary variables ${\displaystyle \sum_{k=0}^p{\beta}_k{f}_k\left({s}_0\right)}$ (i.e. mean function), which account for a spatial trend, and a spatially autocorrelated residual process δ(s ₀) (Matheron 1969):

$$ Z\left({s}_0\right)={\displaystyle \sum_{k=0}^p{\beta}_k{f}_k\left({s}_0\right)}+\delta \left({s}_0\right) $$

(1)

Universal kriging was used to model the relationships between species presence/absence indicator variables from the NFI dataset and the latitude (y, E.D. 50 UTM, in kilometres), elevation (h, in metres above sea level) and exposure (cos(φ), the cosine of the aspect azimuth) derived from a 25 m resolution digital elevation model (DEM) of Spain, for each inventory. In the case of the Scots pine presence/absence indicator variable, exploratory analyses indicated differences in the altitude at which presence prediction was maximum across the latitudinal gradient; therefore, a polynomial of order 2 including h ², h and the h × y interaction was used to model the latitudinal variation in altitudinal maximum in the mean function. The model performance was notably improved by incorporating the cosine of the aspect azimuth and the cos(φ) × y interaction to model the latitudinal variation of this variable. Therefore, the mean function for the Scots pine presence/absence indicator variable was the following:

$$ {\displaystyle \sum_{k=0}^p{\beta}_k{f}_k\left({s}_0\right)}={\beta}_0+{\beta}_1h+{\beta}_2{h}^2+{\beta}_3 hy+{\beta}_4 \cos \left(\varphi \right)+{\beta}_5 \cos \left(\varphi \right)y $$

(2)

In the case of the beech presence/absence indicator variable, preliminary analyses did not show any interaction between altitude and latitude. However, the altitudinal distribution clearly showed a maximum presence probability. In addition, a latitudinal trend was noted for cos(φ); thus, the mean function finally chosen was the following:

$$ {\displaystyle \sum_{k=0}^p{\beta}_k{f}_k\left({s}_0\right)}={\beta}_0+{\beta}_1h+{\beta}_2{h}^2+{\beta}_3y+{\beta}_4\kern0.5em \cos \left(\varphi \right)+{\beta}_5 \cos \left(\varphi \right)y $$

(3)

The universal kriging prediction p(Z,s ₀) of the regionalized variable Z(s ₀) can be interpreted as the probability of the presence of each species analysed in s ₀ (Goovaerts 1994) and is calculated as:

$$ \mathrm{p}\left(Z,{s}_0\right)={\displaystyle \sum_{i=1}^n{\lambda}_i\cdot Z\left({s}_i\right)} $$

(4)

under the unbiasedness conditions:

$$ \begin{array}{l}\begin{array}{ll}{\displaystyle \sum_{i=1}^n{\lambda}_i\cdot {f}_k}\left({s}_i,{t}_i\right)={f}_k\left({s}_0,{t}_0\right)\hfill & k=0\dots p\hfill \end{array}\hfill \\ {}{f}_0\left({s}_i,{t}_i\right)\equiv 1\kern3em \forall \kern2em i\hfill \end{array} $$

(5)

where Z(s _i ) is the value of the variable at n sampling points. The variance of the universal kriging prediction (σ _UK(s ₀)) was calculated as in Cressie (1993).

The variogram that describes the degree of spatial dependence of the regionalized variable was modelled using the spherical variogram defined by the following parameters: the nugget (the semivariance value at the origin), the autocorrelation range (the distance at which the semivariance stabilizes) and the sill (the semivariance for distances greater than the autocorrelation range). One of the crucial steps when using the universal kriging is the variogram parameters and β coefficient fitting. The variance least square (VLS) method used for the variogram and mean function estimation seems more suitable than maximum likelihood methods for the presence/absence indicator variable because it does not require multi-Gaussian distributional assumption (Montes and Ledo 2010).

The kriging prediction may give values outside the range [0,1], although these deviations are usually of small magnitude (Goovaerts 1994). There are different techniques used for constraining the indicator kriging prediction to [0,1] (Tolosana-Delgado et al. 2008). For the comparison between inventories pursued in this study, such constrains are not necessary; however, to map the occurrence of the species (Fig. 1), the kriging indicator predictions outside [0,1] were set to the nearest bound (upward–downward correction; Deutsch and Journel 1992).

Fig. 1

Universal kriging predictions for the distribution areas of P. sylvestris (black) and F. sylvatica (dark grey) and areas where both species coexist (light grey) from NFI1 (1971), NFI2 (1989–1990), NFI3 (1999–2000) and NFI4 (2009–2010). Kriging was performed over a 500 × 500-m grid covering the surveyed areas, and a threshold of 0.5 for the indicator variable value was used to define the presence/absence of each species at each prediction point

2.4 Cross-validation of the universal kriging model

A leaving-one-out cross-validation was used to determine the prediction bias and accuracy of the prediction variance estimates of the fitted model (Cressie 1993). Differences between the predicted and observed values (SEE, sum of estimation errors) were used to assess the bias.

To test the accuracy of the prediction variance estimation, the ratio of the variance of the residuals to the prediction variance was estimated through the fitted variogram (σ _UK(s ₀)). Variance of the standardized estimation errors (VSEE) was calculated. Values close to 1 indicated a similar distribution for the observation and prediction errors.

2.5 Testing the significance of changes in species distribution between inventories

To test the statistical significance of the changes in species distribution along the gradients of the analysed variables between inventories, a novel significance test based on block kriging variances (Isaaks and Srivastava 1989) was developed. The block kriging estimates the mean of the variable for a determined region (or block) through averaging the universal kriging weights λ _i for the points resulting from the discretizacion of the block. The block mean and the block prediction variance were estimated for the strata determined by the classification of the forest area (a) in ten altitude classes, (b) in seven latitude classes and (c) in nine exposure classes, which constitute the blocks in this case. The block mean of the presence/absence indicator variable can be interpreted as the mean probability of species presence in each stratum.

As kriging provides unbiased estimators, the block mean does not depend on the number of plots in each inventory, although the kriging variance depends on the spatial autocorrelation and the distance between samples, which varies with the plot density. To assess differences between inventories at 95 % confidence level, confidence intervals for the presence/absence block prediction at each altitude, latitude and exposure class were calculated from the kriging variance for each inventory date (Isaaks and Srivastava 1989):

$$ \left(\mathrm{p}\left(Z,{s}_0\right)-1.96{\sigma}_{\mathrm{UK}}\left({s}_0\right),p\left(Z,{s}_0\right)+1.96{\sigma}_{\mathrm{UK}}\left({s}_0\right)\right) $$

To determine the significance of the differences between the indicator variable kriging predictions of the different NFI cycles, these were compared with the size of the confidence intervals. The use of confidence intervals allows the differences to be presented as a range of magnitudes, rather than just assess the statistical significance (Katz 1992). The absence of overlapping between the confidence intervals at two different inventory cycles gives a conservative test at 95 % confidence for the difference between the predicted values. Note that the presence/absence prediction mean is estimated for each class, so bias due to differences in forest area between inventories is avoided.

2.6 Mapping changes in species distribution

To graphically assess the species distribution at each inventory date, a threshold probability indicating the species presence was set at 0.5 (Montes et al. 2005), checking that the proportion of predicted points above the indicator threshold was unbiased. The indicator variables of both species were inferred in a 500 × 500-m grid covering the forest area at each inventory cycle to build the prediction maps.

2.7 Software used

ArcGis 9.2 was used to interpolate 25-m resolution DEM from the topographical map of Spain and to derive latitude, altitude and cosine of the aspect azimuth at sample and prediction points. Geostatistical analyses were performed using a Microsoft® VisualBasic® application developed by the authors.

3 Results

3.1 Variogram models

The variogram models fitted for the Scots pine and beech differed considerably. The nugget effect was smaller for the Scots pine (Table 2), whereas the spatial autocorrelation range was larger for beech (Table 3). The autocorrelation range in the Scots pine model progressively diminished from approximately 3 km in 1971 to 0.86 km in 2010, which indicates a decrease in the spatial continuity of the Scots pine indicator variable in the study area. In contrast, the progressively larger autocorrelation ranges of beech (from 9 km in 1971 to 22 km in 2010) indicated increasing continuity in the spatial distribution. The analyses of cross-validation residuals showed that the universal kriging models fitted performed satisfactorily in terms of bias and kriging variance estimation in both cases (Tables 2 and 3).

Table 2 Universal kriging models, VLS estimations of variogram parameters and β coefficients for P. sylvestris and explanatory variables; elevation, square of the elevation, exposure (Expo) and the interaction between latitude, elevation and exposure. Cross-validation sum of residuals (SEE) and the residual variance/prediction variance ratio (VSEE) are also shown

Table 3 Universal kriging model, VLS estimations of variogram parameters and β coefficients for Fagus sylvatica and explanatory variables; elevation, the square of the elevation, latitude, exposure (Expo) and the latitudinal interaction with the exposure. The cross-validation sum of residuals (SEE) and the residual variance/prediction variance ratio (VSEE) are also shown

3.2 Mapping forest range changes

The most notable changes in the universal kriging prediction maps (Fig. 1) were the expansion of the Scots pine from 71 to 79 % of the plots from 1971 to 2010 and an increase from 53 to 72 % in the case of beech. In addition, there were changes in the distribution areas of Scots pine and beech. For instance, in the NFI1, the two species predominantly occupied distinct areas (coexisting in 31 % of the study area), whereas in the NFI3 and NFI4, the area in which the two species coexisted had increased to approximately 42 % of the study area.

3.3 Assessing species distribution changes between inventories

Differences for a range of magnitude between species distribution across the NFI are assessed through the confidence intervals of the block kriging mean predictions for each of the latitude, altitude and exposure classes considered. The NFI1 and NFI2 show wider confidence intervals due to the lower sampling density. The most notable differences for both species arose at mid latitudes (4,747–4,752 km), where the block mean kriging prediction for NFI3 and NFI4 falls outside the confidence intervals for the NFI1 and NFI2 (Fig. 2). These differences were significant at 95 % level (no overlapping between prediction confidence intervals) for beech between the NFI1 and the NFI4.

Fig. 2

Universal kriging block mean prediction for the presence/absence indicator variables by latitudinal (above), altitudinal (mid) and exposure (below) classes for P. sylvestris (a, c, e) and F. sylvatica (b, d, f) from the NFI1, NFI2, NFI3 and NFI4 (from light to dark grey) data. Error bar lines represent the 95 % confidence intervals ci-u upper confident interval, ci-l lower confident interval of the block mean for each inventory and species. The black arrow indicates the direction of the shifting trend of each tree species in the study area

3.3.1 Latitudinal shifts

The results revealed changes in the latitudinal distribution range of both species at the spatial scale and time frame studied. The block kriging prediction for the ten latitudinal classes suggests a major increase in the presence of Scots pine in the mid-northern latitudes of the study area from 1971 to 2000 and in the northern latitudes from 2000 to 2010 (Fig. 2a). This increase was concomitant with a decrease in Scots pine presence prediction in southern latitudes from 2000 to 2010, which indicates a possible shift in the distribution range of this species towards northern latitudes within the study area. In contrast, the presence of beech generally increased in the region both in northern latitudes (mainly from 1971 to 1990) and mid-southern latitudes (from 1990 to 2010) (Fig. 2b), where this species was scarce in the first inventories.

3.3.2 Shifts in elevation

Figure 2c, d reveals a general shift in the presence/distribution of both species towards higher elevations. In the case of Scots pine, the optimum altitude (corresponding to the maximum presence prediction values) was found at about 700 m in 1971 and 800–900 m in 2010. The maximum altitude at which the presence prediction was above the 0.50 threshold was 1,250 m in 1971 and 1,450 m in 2010 (Fig. 2c). In the case of beech, the optimum altitudinal prediction varied over the period 1971–2010, from 1,300–1,400 to 1,450–1,600 m. The maximum altitude at which the species was found shifted from 1,550 to 1,750 m (Fig. 2d).

3.3.3 Changes in exposure distribution

The block averages of the presence by exposure class show that this environmental factor is not a determinant for Scots pine distribution (Fig. 2e). In contrast, the presence prediction for beech has increased for northern exposures, where this species is more frequent (Fig. 2f).

3.4 Trade-off between altitude/exposure and latitude

The inclusion of altitude for latitude and exposure for latitude interactions in the universal kriging model provided a deeper insight into the species distribution at regional scale. In the study area, there was a south–north elevation gradient that influenced the altitudinal species distribution. However, Figs. 3a and 4 show that for both species, the largest elevation shifts took place at mid- and southern latitudes, where the sub-Mediterranean influence is greater and the elevations are lower. The presence increased for northern exposures and southern latitudes for both species (Fig. 3b), following the elevation gradient.

Fig. 3

Universal kriging block predictions of the mean distribution likelihood of each species and NFI, considering the interactions between altitude and latitude (above) as well as exposure (as the cosine of the aspect azimuth (cos(φ)) and latitude (below))

Fig. 4

Schematic reconstruction of P. sylvestris and F. sylvatica forest dynamics from 1971 to 2010 in the Western Pyrenees. The spatial expansion induced by latitudinal and altitudinal shifts of both species is shown based on the combination of universal kriging results for the two species in the NFI1 and NFI4. Light grey trees represent the sporadic occurrence of each species with a presence prediction value below 0.5

4 Discussion

4.1 A new method for assessing changes in species distribution from long-term forest inventories

In this study, we have presented a new method based on universal kriging for the early detection of changes in the distribution of forest species using long-term information derived from the NFI. Kriging techniques provide optimal and unbiased estimates of a regionalized variable in the presence of spatial autocorrelation, allowing the uncertainty of the estimates to be quantified (Mandallaz 2000). In our case study, spatial autocorrelation accounted for a large part of the semivariance of the presence/absence indicator variable for both P. sylvestris and F. sylvatica. This fact has also been reported for other species with widespread distribution area (Segurado and Araújo 2004). The spatial autocorrelation reflects the spatial continuity of the main factors involved in species dynamics (dispersal, mortality, environmental/physical barriers, historical biogeography and socioeconomic factors). Furthermore, the changes observed in the variogram of each species across the successive inventories provide relevant information for understanding the spatial dynamics and patterns of the species in the region.

The distribution of the analysed species is non-stationary at the scale of the study, whereas F. sylvatica dominates in the Atlantic influenced north-western forests, and P. sylvestris is distributed in the southernmost areas of the study region. Universal kriging, as opposed to other kriging techniques, is capable of explaining the non-stationarity of the species distribution in the study area through the mean function depending on the auxiliary variables included in the model (Cressie 1993). Our results reveal how local environmental variations can also constrain the species distribution, highlighting the benefit of using local and landscape approaches to identify species–environment relationships that are averaged out by global models (Osborne et al. 2007). Furthermore, the species–environment relationships can also display non-stationarity for broadly distributed species (Miller and Hanham 2011), such as the upward shift in the altitudinal distribution of the Scots pine across the latitudinal gradient in the study area or the preference of both species for less exposed locations as the Mediterranean influence increases at southern latitudes (Jump et al. 2009). These latitudinal changes in the species–environment relationships are incorporated through the altitude × latitude and exposure × latitude interactions in the universal kriging model. The species–environment interactions can be a key aspect of analysis when studying changes in the species distribution ranges in transitional bioclimatic regions.

Direct comparison of NFI plot data may lead to overestimation of changes in species distributions due to the fluctuating population dynamics, changes in forest inventory area and variations in sampling density across the successive NFIs (Loiselle et al. 2003). The block kriging prediction resolves this problem by providing an unbiased estimator of the mean for the blocks from the information contained in the entire data set. Block kriging incorporates the change of support from point to block, smoothing variations caused by the different sampling schemes used (Yoo and Trgovac 2011) and allowing direct comparison. Figure 2 shows the wider confidence intervals derived from the block kriging variance for the NFI1 and NFI2, which had lower sampling densities than the NFI3 and NFI4. The confidence interval approach was found to be a useful technique for identifying trends in the differences between inventories along environmental gradients. This approach provides information about the magnitude of the differences between inventories relative to the accuracy of the kriging prediction, which may be an alternative to hypothesis testing (Katz 1992) to assess the significance of the changes observed in species distribution.

This new method may provide the basis for more studies at a broader scale. Although in this study we have analysed presence/absence data to identify evidence of changes in the spatial distribution of species, other variables of biological interest derived from plot level information of NFIs could be analysed using these techniques. Based on the changes in species distribution observed using this methodological approach, future research could focus on disentangling the effect of land use changes and climate change on forest dynamics.

4.2 Changes in the distribution of forests

The methodological approach based on the universal kriging model allowed us to obtain evidence of some of the impacts of global change on species distribution over recent decades at regional scale. Firstly, the universal kriging model indicated a progressive increase in the occurrence of the two species in the study area, which agrees with previous studies reporting a global increase in forests in mountainous areas of other developed countries (Gellrich and Zimmermann 2007). Secondly, with regard to the altitudinal distribution of both Scots pine and beech, the approach based on the universal kriging model revealed a vertical shift of approximately 200 m towards higher elevations between 1971 and 2010; a finding supported by other research conducted at European scale (Lenoir et al. 2008).

As expected, the distribution of Scots pine in the study area was found to have shifted over recent decades towards northern latitudes and higher altitudes. This result is consistent with previous predictions for pine species in the Iberian Peninsula (Benito Garzón et al. 2008). In addition, the spatial autocorrelation range of Scots pine decreased between 1971 and 2010, implying a probable fragmentation and loss of continuity of its distribution area. This finding, together with the decreasing presence of the species over the last decade in the southernmost areas, where the sub-Mediterranean influence is higher, may reflect an increase in the presence of other tree species in the most recent NFI. In these zones, it may be that Scots pine is being replaced by broadleaved taxa like Q. pubescens, which is more adaptable to drier conditions, as reported in studies of other marginal populations of this species (Gimmi et al. 2010). The shift towards the northernmost latitudes and higher altitudes along with both the decreasing presence in the southernmost latitudes and the fragmentation of its distribution area may suggest a northward retraction of the Scots pine in the region.

Beech populations have shown the greatest expansion in the study area over the four decades covered by this study, increasing its presence at higher altitudes. This seems to confirm the ‘deciduous tree invasion’ in the sub-alpine belt predicted by Kräuchi and Kienast (1993). The progressive increase in the spatial autocorrelation range of the beech indicator variable reflects an increasingly continuous distribution of the species in the study area, as reported in other areas of Europe with similar scenarios (Poljanec et al. 2010). The presence of beech increased significantly between NFI1 and NFI2 in the northernmost mountains, where this species was already widely distributed. However, in the most recent NFIs, the presence of this species unexpectedly increased at mid latitudes, where it might be expected that the temperate, humid ecological requirements of this species would not be satisfied due to the weaker Atlantic influence. However, it must be considered that this increase mainly occurred on north facing slopes, reflecting the ability of beech to spread easily, although only to favourable biotopes.

The altitudinal and latitudinal shifts in the distribution range of both species have led to an extension of the area where the Scots pine and beech coexist (Figs. 1 and 4). These findings provide new evidence of changes in the forest composition. Scots pine and beech can be considered ‘engineer species’ because their presence substantially characterizes the environment of a site and can define habitats. Therefore, variations in the forest composition resulting from the expansion or retraction of these species populations might have profound implications on the diversity and distribution of other species associated with their forests. For example, endangered species of folivores such as the capercaillie (Tetrao urogallus), very sensitive to habitat fragmentation, find their westernmost Pyrenean habitat in these forests (Rodríguez and Obeso 2000). In the future, the biotic interactions established between the two species analysed in this study will play a crucial role in the persistence and stability of these forest areas. Moreover, these sensitive transitional forest areas will be central to monitoring the effects of global change and to the development of successful conservation strategies for the region.

4.3 Conclusion

This work provides a new tool for the early detection of changes in forest species and constitutes a first step towards tackling the effects of global change on forests. The methodological approach proposed allows the spatial and temporal shifts in species distribution over recent decades to be analysed using long-term NFI and regardless of differences in the sampling schemes used in each. Furthermore, the confidence interval approach proposed allows us to assess whether the differences detected in the distribution of species are significant.