Assessment of the effects of biophysical and anthropogenic factors on woody plant encroachment in dense and sparse mountain grasslands based on remote sensing data

Abstract

Land abandonment exacerbated by climate change has led to increased woody plant encroachment of mountain grasslands in many regions of the world. The present study assessed woody plant encroachment below potential tree line in the Central Pyrenees of Spain and the association of this encroachment with changes in land use. Remote sensing data from Landsat-5 Thematic Mapper (TM) from the mid-1980s and mid-2000s were analyzed by supervised classification for identification of land cover types. The transition matrix indicated that shrublands were the most dynamic plant communities. Consequently, 21% of cultivated areas, 19% of dense grasslands, and 24% of sparse grasslands became shrublands during the period analyzed, and 35% of shrublands became forest. Generalized additive mixed models (GAMMs) were used to identify biophysical and anthropogenic factors that were significantly correlated with woody plant encroachment of dense and sparse grasslands. Distance to the nearest woody plant habitat (shrub or forest) was the most strongly correlated factor with woody plant encroachment of both types of grassland. This factor explained 69% and 71% of the variance in models of dense and sparse grasslands, respectively. Besides this factor, anthropogenic factors had larger effects on woody plant encroachment of dense grasslands, regions that were more productive and accessible. However, biophysical and especially topographic factors had slightly greater effects on woody plant encroachment of sparse grasslands, regions that were less productive and accessible. The changes in land cover that we observed indicated that land cover has become more homogeneous. There have been reductions in the variety, functions, and services of grasslands, particularly in areas below the potential tree line that are vulnerable to the development of woody plant habitats.

Keywords

land abandonment land cover change mountain grassland remote sensing Spanish Pyrenees woody plant encroachment

I Introduction

Woody plant encroachment of natural and semi-natural mountain grasslands is a widespread and acute problem throughout the world (Brandt et al., 2013; Fernández-Giménez and Fillat, 2012b; Lasanta-Martínez et al., 2005; Molinillo et al., 1997; Naito and Cairns, 2011; Ratajczak et al., 2012; Sitko and Troll, 2008; Stueve et al., 2011; Vicente-Serrano et al., 2005). Most of these habitats, especially mountain grasslands below the potential tree line, were created by the removal of shrubs and forests, and their persistence depends heavily on traditional land use practices such as livestock grazing and woody plant removal by fire or mechanical means (Didier, 2001; Holtmeier and Broll, 2007; Sankey et al., 2006; Sitko and Troll, 2008; Wehn et al., 2011). In the absence of these activities, these man-made habitats undergo ecological succession toward shrub- or tree-dominated vegetation (Lasanta-Martínez et al., 2005; Schirmel et al., 2011), leading to a reduction in the amount of mountain grasslands. Recent changes in human demographics and agricultural policies have led to these changes (Pôças et al., 2011), and global warming has exacerbated this problem (Batllori and Gutiérrez, 2008; Bolli et al., 2007; Holtmeier and Broll, 2007; Vicente-Serrano et al., 2005). Thus, a combination of social, economic, and natural processes, which operate from local to global scales, has altered the environments of mountainous regions (Gerard et al., 2010; Liverman and Cuesta, 2008; Mottet et al., 2006; Roura-Pascual et al., 2005). More studies on the extent of these changes and the factors that drive them are vital because these regions contain some of the best preserved natural communities and are protected by the European Union Habitats Directive.

Woody plant encroachment of grasslands reduces biodiversity (Ratajczak et al., 2012), increases the risk of fire (Lloret et al., 2002; Moreira et al., 2001), affects ecosystem functions by reducing the amount of productive grazing habitats for domestic and wild herbivores (Fernández-Giménez and Fillat, 2012b), and alters the culturally significant features and aesthetics of the landscape (Antrop, 2005). Nevertheless, woody plant encroachment can minimize landslides and soil erosion (García-Ruiz et al., 1996) and, in semi-arid grasslands, can reverse desertification (Maestre et al., 2009). In addition, woody plant encroachment of areas with high annual precipitation can increase carbon storage due to the increased net primary production (Knapp et al., 2008).

Since the Middle Ages, humans in the Pyrenees have increased the amount of cultivable land and summer grasslands by removal of forest and shrublands (Daumas, 1976; Kouba and Alados, 2012; Montserrat-Martí, 1992; Roura-Pascual et al., 2005). These man-made and natural mountain ecosystems provided a variety of resources, and, for centuries, high-elevation grasslands were used for extensive livestock grazing in summer. However, after the 1930s, livestock (mainly sheep) densities declined coinciding with an important decrease in the human population (IAE, 2012; Lasanta-Martínez et al., 2005). The Pyrenees experienced a substantial reduction in the number of agropastoral workers (mainly shepherds), many of whom migrated to large industrialized cities. Human depopulation, a reduction in livestock density, and changes in land use led to the abandonment of many areas, and this has led to changes in the composition of the vegetation. In particular, native shrubs have encroached the grasslands and cultivated areas and this process was followed by afforestation of the shrublands (Bartolome et al., 2005; Molinillo et al., 1997; Roura-Pascual et al., 2005). Previous studies have assessed the biophysical and anthropogenic factors that contribute to woody plant encroachment, especially the encroachment of cultivated areas and mountain grasslands (Brandt et al., 2013; Gellrich et al., 2007; Grau et al., 2008; Kouba and Alados, 2012; Lasanta-Martínez and Vicente-Serrano, 2007; Millington et al., 2007; Schulz et al., 2011; Sitko and Troll, 2008). However, few studies have compared woody plant encroachment in different types of mountain grasslands (Montané et al., 2010), such as the sparse and dense grasslands of the Pyrenees.

Remote sensing data are a valuable and cost-effective source for quantification of spatial-temporal changes in vegetation cover (Boyd and Danson, 2005; Stow et al., 2004). In particular, Landsat imagery has been used to measure changes in land cover over the last decades (Brandt et al., 2013; Gautam et al., 2003; Grau et al., 2008; Pôças et al., 2011; Schulz et al., 2011; Shalaby and Tateishi, 2007; Stow et al., 2004; Tømmervik et al., 2004).

The general objective of our research is to identify factors that promoted woody plant encroachment of dense and sparse mountain grasslands in Aragón, Central Pyrenees, Spain. The two grassland types under study differ in plant composition, biomass production, and geographic location, so are expected to respond differently on woody plant encroachment to biophysical and anthropogenic factors. Dense grasslands occur on low slopes and have deep soils, leading to high grass coverage. Conversely, sparse grasslands occur on steep slopes and have shallow soils, leading to low grass coverage (Fillat et al., 2008). Sheep and goats are the primary grazers of sparse grasslands and cattle and horses of dense grasslands, which are more accessible and have higher biomass production (Fernández-Giménez and Fillat, 2012a).

Our specific purposes are: (1) to use remote sensing data to identify the types of land cover in the mountainous areas in Aragón, Central Pyrenees, and quantify the spatial and temporal changes in these habitats, especially the loss of dense and sparse grasslands due to woody plant encroachment; and (2) to identify biophysical and anthropogenic factors that have led to woody plant encroachment of dense and sparse grasslands in Aragón, Central Pyrenees.

II Methods

1 Study area

The study area (138,363 ha) was in Aragón, Central Pyrenees, Spain, at 42º36’N, 0º00’E (Figure 1), which is within the Alpine Mountain region (600–3340 m asl). The area has a mountain climate that is characterized by equinoctial precipitation and a continental influence strongly affected by elevation. At 2200 m, near the potential tree line, the average annual precipitation is 1657 mm and the average daily maximum and minimum temperatures are 8.6ºC and 1.4ºC, respectively (AEMET, 2012; Cuadrat et al., 2007).

Figure 1.

Study area in the Central Pyrenees of Aragón, Spain, where the potential tree line is around 2100 m asl.

In the study area, climate has a strong influence on vegetation physiognomy along an elevation gradient (Gartzia et al., 2013). At high elevations (>2100 m), the vegetation is mostly sparse grasslands that are dominated by Festuca gautieri. Plants in these grasslands are adapted to harsh climatic conditions, the slopes are steep, the soil is very shallow, rocky outcrops are common, and vegetation cover is less than 50% (Fillat et al., 2008). In flatter areas, where the soils are deeper, there are dense grasslands that have grass coverage of more than 50%. These grasslands include various communities of a variety of pasture grasses, sedges, and leguminous species, and Festuca nigrescens and Agrostis capillaris, or Nardus stricta, Festuca eskia, or F. paniculata, are the dominant species (Fillat et al., 2008). Dense grasslands are adjacent to sparse grasslands at high elevations and adjacent to shrublands and forests at low elevations.

Typically, sheep and goats graze the sparse grasslands and cattle and horses graze the dense grasslands (Fernández-Giménez and Fillat, 2012a). Ecological succession from dense grassland to sparse grassland or vice versa is uncommon because of topographical limitations, which do not change. However, succession from grassland to woody community is common when anthropogenic pressures are reduced, especially below the potential tree line (2100 m asl) (Camarero and Gutiérrez, 2004), where forest is almost always the potential vegetation. In the shrublands, almost all of which are in transition to forest, Buxus sempervirens is predominant, along with Echinospartum horridum on basic soils and Juniperus communis or Rhododendron ferrugineum on acidic soils. Coniferous trees – mostly Pinus sylvestris or P. uncinata –, beech (Fagus sylvatica), or a variety of oaks (Quercus spp.) are predominant in the forests. At the bottoms of valleys, where the slopes are shallow and the soils are deep, the land was cultivated for crop production, particularly potatoes and cereals. However, these regions are currently meadows (Fernández-Giménez and Fillat, 2012a). Many cultivated areas have been abandoned in recent decades, and this has led to a transition to Genista scorpius or Rosa spp. shrub communities and forests.

2 Supervised classification of land cover using remote sensing data

a Remote sensing data acquisition and processing

Remote sensing data can be used to quantify spatial and temporal changes in the landscape (Boyd and Danson, 2005; Stow et al., 2004). Multi-temporal data sets, such as those from Landsat-5 TM, from 1984 until 2011, have allowed the study of changes in land cover in mountain areas (Brandt et al., 2013; Gautam et al., 2003; Grau et al., 2008; Pôças et al., 2011; Schulz et al., 2011; Shalaby and Tateishi, 2007; Stow et al., 2004; Tømmervik et al., 2004).

We added phenological information to increase the accuracy of the land cover classification by using multi-date imagery in the supervised classification (Brandt et al., 2013; Kuemmerle et al., 2006). We used cloud-free Landsat-5 TM imagery acquired on 3 July 1984 (early summer) and 14 September 1987 (late summer) for the first period and 4 August 2007 and 5 September 2007 (early-late summer) for the second period to identify and quantify the types of land cover (forest, shrubland, dense grassland, sparse grassland, and cultivated areas) in the study area. The imagery were provided by the USGS at level 1T (ortho-rectified image) and correspond to Landsat-5 TM scene 199-30. Prior to classification we used the ATCOR-2/3 (v 9.3), an atmospheric correction algorithm created by Dr Richter of the German Aerospace Center, to convert the digital value into reflectance (expressed as a percentage) and to remove atmospheric and topographic effects (based on the Lambertian assumption).

b Land cover supervised classification approach

We used 21 layers in each time period for the supervised classification of land cover (Figure 2) (Gartzia et al., 2013). There were six bands from the Landsat-5 TM imagery, except the thermal band (Kuemmerle et al., 2006), acquired in early summer, and another six bands in late summer, and there were two vegetation indices per image derived from these imagery: Normalized Difference Vegetation Index (NDVI) and the Soil-Adjusted Vegetation Index (SAVI) (Baret and Guyot, 1991). There were five topographical features (elevation, slope, north–south aspect, east–west aspect, insolation) that were derived from the 30 m resolution digital elevation model (DEM) (SITAR, 2012) using ArcGis 9.3 software. North–south and east–west aspects were cosine- and sine-transformed, respectively (Felicísimo, 1994).

Figure 2.

Protocol used for supervised classification of land cover in the study area (Aragón, Central Pyrenees), based on satellite images taken during the mid-1980s and mid-2000s (two images per period), two neocanals, and topographic data derived from the digital elevation model (DEM). Random Forest (RF) was used as a classifier.

The training data used in the supervised classification were based on 1000 randomly selected pixels (Gartzia et al., 2013) from imagery acquired in the mid-2000s. Each of the pixels was assigned to a land cover category based on a vegetation map in the area, field inspections, and ortho-photographs acquired in 1997 and 2006 (SITAR, 2012). The supervised classification was based on the non-parametric random forest (RF) classifier (Gartzia et al., 2013; Rodriguez-Galiano et al., 2012) using the randomForest package in the R software environment (Liaw and Wiener, 2002; R Development Core Team, 2011). The RF is not constrained by parametric restrictions and is not sensitive to collinearity and overfitting of the data (Breiman, 2001), allowing the use of layers with collinearity issues. This is useful for mountainous habitats, where topographic factors tend to be highly collinear (De la Riva, 1997; Gartzia et al., 2013; Shi et al., 2009).

Dense grasslands and cultivated areas were difficult to distinguish with the method used for supervised classification, so we initially classified both land cover types as dense grasslands. After the supervised classification, the crop and land use maps of 1978 and 2008 (MAGRAMA, 2010) were used to distinguish cultivated areas from dense grasslands.

The accuracy of the classified map from the mid-2000s was tested using another 1000 randomly selected pixels. These pixels are identified using the same methodology described in the training data identification. We calculated Cohen’s Kappa Index of Agreement (Cohen, 1960), a standard method in remote sensing research (Congalton and Green, 2009). After confirming the adequate accuracy of the supervised classification of the mid-2000s land cover map (Cohen’s Kappa = 0.91), mid-1980s imagery was classified based on the random forest parameters derived from the mid-2000s supervised classification (Figure 2) using 1980s remote sensing data and ancillary data.

3 Quantification of changes in land covers below the potential tree line (2100 m asl)

The primary objective of this study was to assess changes in land cover due to the ecological succession of vegetation, specifically the extent of woody plant encroachment of dense and sparse grasslands and cultivated areas, and the afforestation of shrublands. To minimize the potentially confounding effects of climate change, our analysis was restricted to below the potential tree line (Gellrich et al., 2007), where changes in land use are believed to have had a greater impact than global warming (Bryn, 2008; Gehrig-Fasel et al., 2007).

We used Chi-square test to measure the significant difference in the land covers between the two periods. Then, we used a transition probability matrix based on the two supervised classification maps to quantify the extent and direction of changes in vegetation (Brandt et al., 2013; Kouba and Alados, 2012; Mottet et al., 2006; Roura-Pascual et al., 2005). The probability (Pij) of transition from one cover class (i) to another cover class (j) was calculated as: Pij = (Aj 2000s/Ai 1980s), where Aj 2000s is the area of class j in the mid-2000s that was class i in the mid-1980s, and Ai 1980s is the total area of cover class i in the mid-1980s.

4 Identification of the main factors affecting woody plant encroachment in dense and sparse grasslands below the potential tree line (2100 m asl)

Biophysical factors and anthropogenic factors were used to identify the causes of changes of dense and sparse grasslands cover (Table 1):

Table 1.

Biophysical and anthropogenic factors included in the GAMMs to assess the causes of woody plant encroachment of mountain grasslands of the Central Pyrenees, Aragón, Spain. IGN = Spanish National Geographic Institute (IGN, 2012); IAE = Aragón Institute of Statistics (IAE, 2012); SITAR = Land Information System of Aragón (SITAR, 2012); CAA = Climatic Atlas of Aragón (Cuadrat et al., 2007).

	Explanatory factors used in GAMM	Factor type	Description	Spatial resolution (source)
Biophysical factors
Topographic
	Elevation	Not linear	Elevation	30 m (DEM)
	Slope	Linear	Slope in degrees	30 m (DEM)
*	North–south aspect	Linear	Cosine of the aspect	30 m (DEM)
	East–west aspect	Linear	Sine of the aspect	30 m (DEM)
*	Insolation	Linear	Solar Radiation Tool	30 m (DEM)
Abiotic
*	Lithology	Categorical	Basic, acid, quaternary mat.	1:50,000 (SITAR)
*	Resistance to the erosion	Categorical	Low/high resist. to the erosion	1:50,000 (SITAR)
*	Distance to rivers	Linear	Euclidean distance	1:25,000 (IGN)
Biotic
	Distance to woody habitats present in mid-1980s	Linear	Euclidean distance	30 m (Land Cover mid-1980s)
Climatic
*	Precipitation (mm/year)	Linear	Climatic Atlas of Aragon	100 m (CAA)
*	Mean of the minimum temperatures (ºC/day)	Linear	Climatic Atlas of Aragon	100 m (CAA)
Anthropogenic factors
*	Distance to main road	Linear	Euclidean distance	1:25,000 (IGN)
	Distance to passable road	Linear	Euclidean distance	1:25,000 (IGN)
	Distance to towns	Linear	Euclidean distance	1:25,000 (IGN)
	Distance to buildings	Linear	Euclidean distance	1:25,000 (IGN)
	Depopulation in historical municipality (HM)	Categorical	Depopulation <50% / >50%	HM level (IAE)
	Owner type of the land	Categorical	Private or public land	1:25,000 (SITAR)

*Factors not used in the model because of high collinearity (>0.7) with other factors or not significant in the model.

Biophysical factors. (a) Topographical factors derived from the DEM (SITAR, 2012): elevation (m), slope (º), north–south aspect (cosine of the aspect), east–west aspect (sine of the aspect), and insolation. (b) Abiotic factors: lithology (basic, acidic, or quaternary materials) (SITAR, 2012), resistance to erosion of the geological materials (SITAR, 2012), and distance to the nearest river (SITAR, 2012). (c) Biotic factors: distance to the nearest woody plant habitat that was present in the mid-1980s (based on the map prepared in this study). (d) Climatic factors: precipitation and temperature (Cuadrat et al., 2007).

Anthropogenic factors. (a) Distance to the nearest main road and any passable road (IGN, 2012; SITAR, 2012). (b) Distance to the nearest building or town (IGN, 2012; SITAR, 2012). (c) Ownership of the land (private or public) (SITAR, 2012). (d) Depopulation in each historical municipality between the 1930s and the 1980s (municipalities with more than 50% depopulation or municipalities with less than 50% depopulation) (IAE, 2012).

Generalized additive mixed models (GAMMs) were used to identify factors that significantly correlated with woody plant encroachment of dense and sparse grasslands. The GAMM models can accommodate spatial autocorrelation and non-linearity in the data (Zuur et al., 2009). In most cases, there was spatial autocorrelation in the landscape data and some of the explanatory factors included in the models had non-linear responses. Pixel identity was included as a random factor to avoid the effects of spatial autocorrelation. Examination of the residuals from the GAMMs (Zuur et al., 2009) confirmed that the assumptions of the model were met (normality, independence, and heteroscedasticity of the variance of the residuals). The GAMMs were assessed using the mgcv and nlme packages in R software environment (Pinheiro et al., 2011; R Development Core Team, 2011; Wood, 2011). The adequacy of the model was based on calculation of the Akaike Information Criteria (AIC), which considers goodness-of-fit and model complexity (Zuur et al., 2009). To selected the optimal model, we choose the one with lowest AIC, always deleting the not significant factors from the models.

The models included 500 randomly selected pixels from areas that were grasslands in mid-1980s and in the mid-2000s (absence of change in the grasslands), and 500 randomly selected pixels from areas that had changed from grasslands to woody plant communities (presence of the change in the grasslands). The data for dense and sparse grasslands were analyzed separately. Only factors with low collinearity (Pearson r < 0.7; Schulz et al., 2011) were included in the models.

III Results

1 Supervised classification maps and land cover changes in the Central Pyrenees

In the mid-1980s, the land cover above 2100 m asl (24% of the total study area) was 25% dense grassland and 75% sparse grassland. The land cover below 2100 m asl (76% of the total study area) was 28% forest, 38% shrubland, 23% dense grassland, 8% sparse grassland, and 3% cultivated areas (Figure 3a). Between the mid-1980s and the mid-2000s, the amount of forest increased by 13,904 ha (Figure 3). Over the same time, shrubland decreased by 6708 ha, dense grassland decreased by 4100 ha, sparse grassland decreased by 2438 ha, and cultivated areas decreased by 658 ha.

Figure 3.

Supervised classification of the main land covers in the study area (Central Pyrenees, Aragón, Spain) in the mid-1980s (a) and the mid-2000s (b). The percent of cover type refers to land cover below 2100 m asl. The ‘mask’ category includes all areas that did not fall within any of the prescribed land cover types used in the study (e.g. lakes, main rivers, main roads, towns, and shaded areas). (c) Grasslands and cultivated areas that had encroachment by woody plants.

The results of the chi-square test for the comparison of land cover types between the mid-1980s and the mid-2000s periods showed significant differences between both periods (chi-squared = 9.82, p-value = 0.04). The 96% of the changes in land cover occurred below 2100 m asl. Shrubs colonized grasslands and cultivated areas after they were abandoned (Figure 3c), and the amount of forest increased, particularly in shrublands. Between these two periods, 35% of the shrublands became forest, and only 5% of dense grasslands and 8% of cultivated areas were afforested (Figure 4). In sparse grasslands, afforestation was very limited, but shrub encroachment was 24%. In dense grasslands shrub encroachment was 19% and 21% in cultivated areas (Figure 4). Overall, woody plant encroachment, from grasslands or cultivated areas to shrublands or forests, was greater in cultivated areas (29%) than in dense and sparse grasslands (24% each).

Figure 4.

The probability (expressed as a percentage) that a land cover type underwent ecological succession to another type of vegetation. These values were measured for the mid-1980s to the mid-2000s in the Central Pyrenees, Aragón, Spain, below 2100 m asl.

2 Factors influencing woody plant encroachment in dense and sparse grasslands

In the GAMMs, distance to the nearest woody plant habitat (shrub or forest) was the most strongly correlated factor with the extent of woody plant encroachment in dense grasslands (explained 69% of the variance in the model) and sparse grasslands (explained 71% of the variance) (Figure 5; Table 2). In dense and sparse grasslands, 90% and 94% of the woody plant encroachment, respectively, occurred in the first 90 m from the shrub or forest habitat that existed in the mid-1980s. Beyond 330 m from the nearest woody plant habitat, less than 1% of the grasslands had been encroached from the mid-1980s to the mid-2000s. Slope explained 7% and 10% of the variance in the models of woody plant encroachment in dense and sparse grasslands, respectively (Table 2). About 75% of the woody plant encroachment occurred in areas where the slope was 15–35º in dense grasslands and 15–40º in sparse grasslands (Figure 6a). In areas where the slope was less than 5º, 60% of the sparse grasslands had been encroached (Figure 6b), corresponding to regions at elevations of 800–900 m asl (Figure 6d). These low-elevation areas which were woody encroached were close to rivers.

Figure 5.

Spatial representation of woody plant encroachment in mountain grasslands in the Central Pyrenees, Aragón, Spain, from the mid-1980s to the mid-2000s. This image shows a small portion of the total study area and indicates that woody plant encroachment mostly occurred in grasslands near woody plant habitats that were present in the mid-1980s, at mid-elevation regions, and near roads and towns. The ‘mask’ represents cultivated areas, lakes, main rivers and roads, towns, and shaded areas.

Figure 6.

Woody plant encroachment as a function of slope (a and b) and elevation (c and d) in regions below 2100 m asl in the Central Pyrenees, Aragón, Spain. Percentage encroachment is given as a proportion of total area encroached (a and c) and as a proportion at each slope or elevation range (b and d).

Table 2.

Factors that explained a significant amount of the variance in the GAMMs for woody plant encroachment of dense grasslands (top) and sparse grasslands (bottom) from the mid-1980s to the mid-2000s in the Central Pyrenees, Aragón, Spain at elevations below 2100 m asl.

			Woody encroachment in dense grasslands
Factor group	Factor	Factor type	Effect	Estimate	Std error	t value	F value (ANOVA)	% variance (of F)	Pr (>/t/)
	intercept			2.1	0.3	7.5			***
Biophysical	Elevation	Non-linear					10.6	5.7	***
	Slope	Linear	–	–0.057	0.016	–3.6	12.7	6.9	***
	East–west aspect	Linear	–	–0.23	0.12	–1.9	3.7	2.0	.
	Distance to woody habitat	Linear	–	–0.024	0.002	–11.3	127	68.6	***
Anthropogenic	Distance to passable road	Linear	–	–0.00045	0.00018	–2.5	6	3.2	*
	Depopulation >50%	Categorical	>50% +	0.66	0.18	3.7	13.8	7.5	***
	Land ownership	Categorical	Private +	0.75	0.26	2.9	8.3	4.5	**
	Distance to buildings	Linear	–	–0.00036	0.0002	–1.76	3	1.6	.

Significant codes: *** = 0.001; ** = 0.01; * = 0.05; . = 0.1. Significant positive (+) or negative (–) effect.

			Woody encroachment in sparse grasslands
Factor group	Factor	Factor type	Effect	Estimate	Std error	t value	F value (ANOVA)	% variance (of F)	Pr (>/t/)
	intercept			2.1	0.27	7.9			***
Biophysical	Elevation	Non-linear					8	7.0	**
	Slope	Linear	–	−0.034	0.01	−3.4	11.8	10.3	***
	Distance to woody habitat	Linear	–	−0.027	0.003	−9.1	82	71.4	***
Anthropogenic	Depopulation >50%	Categorical	>50% +	0.46	0.18	2.6	6.7	5.8	**
	Distance to towns	Linear	–	−0.00014	0.00005	−2.5	6.3	5.5	*

Significant codes: *** = 0.001; ** = 0.01; * = 0.05. Significant positive (+) or negative (–) effect.

The other factors that correlated with woody plant encroachment each explained less than 10% of the variance in the models (Table 2). In both types of grasslands, elevation correlated with woody plant encroachment in a non-linear manner. The extent of encroachment varied with elevation, in that 75% of the woody plant encroachment occurred at 1400–1900 m asl in dense grasslands, and at 1500–2100 m asl in sparse grasslands (Figure 6c). At lower elevations, however, the percentage of grasslands that had been encroached was higher. In dense grasslands, encroachment was more than 50% at 1000–1500 m asl, 25% at 1500–1900 m asl, and less than 3% near 2100 m asl (Figure 6d). In sparse grasslands, encroachment was more than 60% at riverbanks and eroded hillsides at 800–900 m, 40% at 1100–1400 m, and 4% near 2100 m. In addition, a westerly aspect was positively correlated with woody plant encroachment in dense grasslands.

Although anthropic factors do not explain a high proportion of the variance in the models, there are anthropic variables with significant correlation with the woody plant encroachment, particularly in dense grasslands (Table 2). There was extensive depopulation of the study area from the 1930s to the 1980s, and municipalities with high rates of depopulation (>50%) had the greatest amount of woody plant encroachment. In dense and sparse grasslands, 73% and 63% of the woody plant encroachment, respectively, occurred in municipalities that had high depopulation. In those areas, 22% and 18% of the dense and sparse grasslands, respectively, had woody plant encroachment. However, in areas that experienced less severe depopulation, 11% and 14% of the dense and sparse grasslands had transformed into woody plant habitats.

Finally, encroachment of dense grasslands was more likely to happen in private than public lands (Table 2). In addition, the extent of woody plant encroachment and distance to the nearest passable road or building were negatively correlated. In sparse grasslands, the distance to the nearest town explained a small but significant amount of the variance in the model. Distance to the nearest river, lithology, and resistance to erosion of the geological material did not explain significant amounts of variance in the GAMMs for woody plant encroachment of dense or sparse grasslands.

IV Discussion

The analysis showed that there has been a change in land cover types in the Central Pyrenees since the mid-1980s. A high woody plant encroachment of dense and sparse mountain grasslands below the potential tree line is occurring in grasslands created by human activities several centuries ago in the area. This ecological succession in the grasslands is a natural process and, without human intervention, will continue until these grasslands become forests or shrublands (Lasanta-Martínez and Vicente-Serrano, 2007; Lasanta-Martínez et al., 2005). Although these grasslands are semi-natural habitats, they are highly valued and the need for their conservation is widely recognized. The loss of grasslands began in the 1930s, when rural abandonment of the Pyrenees began (IAE, 2012).

Our results indicate that proximity to a woody cover is the main factor associated with succession of grasslands to shrublands or forests. In grasslands where the nearest woody habitat was more than 330 m away, the probability of woody plant encroachment was less than 1%, even if human activity had ceased decades earlier. Similar phenomena have been observed in other mountainous areas, such as the Alaska Range, the Ukrainian Carpathian Mountains, and the Swiss Alps (Bolli et al., 2007; Gellrich et al., 2007; Sitko and Troll, 2008; Stueve et al., 2011). In the Swiss Alps (Bolli et al., 2007) and Tierra de Fuego, Chile (Cuevas, 2000), the seed-rain and micro-environmental factors were important in the expansion of woody plant communities that were near pre-existing shrublands or forest habitats. This is consistent with ‘shrub autocatalysis’, a feedback mechanism in which shrub encroachment is greatest in areas where woody plant vegetation is already present (Brandt et al., 2013).

In our study area, woody plant encroachment progressed from low elevations to high elevations, especially in patches of grassland that were surrounded by woody plant habitats. Similar phenomena have occurred in the Swiss Alps and the Himalayas (Brandt et al., 2013; Gehrig-Fasel et al., 2007; Gellrich et al., 2007). In our study and in other studies, slope is the topographic factor that had the greatest effect on woody plant encroachment, particularly on mean slopes (Molinillo et al., 1997; Schulz et al., 2011), which are not suitable for human or livestock activities (Schulz et al., 2011). However, these areas have been used when human and livestock densities were higher. In addition, mean steep slopes are more easily colonized by shrubs than shallow slopes, where there is competition with grassland species (Guerrero et al., 1999; Komac et al., 2011). In contrast, steep slopes are less favorable for woody plant species (Gellrich et al., 2007), and this slows the rate at which woody plants can invade these habitats. In the dense grasslands of the Central Pyrenees, woody plant encroachment was most strongly correlated with a westerly aspect, a topographic factor that was strongly associated with human activities in the past (Poyatos et al., 2003). Areas facing the west were most productive because they did not have the harsh climate that is characteristic of north-facing slopes. But many of the west-facing slopes have been abandoned (Molinillo et al., 1997), and this might explain the more substantial woody plant encroachment of these areas.

Dense grasslands and sparse grasslands have been managed differently for many centuries, and these different management practices have influenced the nature and extent of woody plant encroachment. In particular, anthropogenic factors have had a greater influence in dense grasslands than sparse grasslands, as our results demonstrated. Humans have used dense grasslands more intensively because these areas offered the best topography, soil, and plant productivity for herding animals. Recent declines in human activities have contributed to an increase in woody plant encroachment of these dense grasslands. In sparse grasslands, however, slope and elevation have had a greater influence on woody plant encroachment than anthropogenic factors, as our results reveal. Thus, there are different patterns of woody plant encroachment in different grassland habitats (Montané et al., 2010).

In the Central Pyrenees, dense grasslands on private lands were more likely to have suffered from woody plant encroachment than those on public lands. This phenomenon was also observed on savannas in North America (Archer et al., 1989). Livestock pressure on private lands has usually been low to favor high productivity, which has led to an increased risk of woody plant encroachment. On the contrary, the communal pastures are intensively used, which are often overgrazed because individual self-interest can lead to misuse of a commonly held resource (Hardin, 1968). In addition, public lands have experienced less abandonment. Even when there was no local livestock within the municipality to grazing the public lands, the pastures were rented to others who brought livestock from outside the municipality.

The negative correlations that we observed between woody plant encroachment and the distances to the nearest building, town, or passable road suggest that human pressure on grasslands near these infrastructures has decreased and this has led to woody plant encroachment. This phenomenon has been observed elsewhere (Brandt et al., 2013; Gellrich et al., 2007; Grau et al., 2008; Schulz et al., 2011). For example, use of wood as a source of energy for cooking and heating was much more common in the past (Roura-Pascual et al., 2005), when the demand for wood led to harvesting of shrubs and trees in the most accessible areas. Nevertheless, the negative correlations that we observed are contrary to our initial expectation that the most remote sites were more likely to experience shrub encroachment and afforestation (Kouba and Alados, 2012; Schulz et al., 2011).

We found a relevant relationship between rural abandonment and woody plant encroachment of sparse and dense grasslands. Municipalities that had a greater woody plant encroachment were those that experienced the greatest depopulation between the 1930s and the 1980s. Since the 1980s, the population of this region has actually grown slowly (IAE, 2012), but agro-pastoral activity has continued to decline because the agricultural economy has become secondary to the tourism economy. Similar demographic changes occurred in the Italian Alps (Motta et al., 2006), and an increase in immigration to the Swiss Alps did not retard the trend of increased woody plant encroachment because farms continued to be abandoned (Gellrich et al., 2007).

V Conclusions

Significant woody plant encroachment has occurred in the mountain grasslands of the Central Pyrenees since the mid-1980s. This has been particularly severe below the potential tree line, where one-quarter of the dense and sparse grasslands has been lost due to woody plant encroachment. Proximity to shrublands or forest habitats was the most significant factor associated with woody plant encroachment of both types of grasslands. In the dense and sparse grasslands which have been managed differently for many centuries, the anthropogenic factors had a greater influence on woody plant encroachment of dense grasslands, while topographic factors had more influence on woody plant encroachment of sparse grasslands.

These grasslands are under the protection of the European Union Habitats Directive and woody plant encroachment of these grasslands cannot be reversed naturally, and human interventions, such as mechanical clearing and fires, have not always produced satisfactory results. Furthermore, the reversion of a woody plant system to grassland takes much more time than succession from grassland to shrubland. Thus, we expect that in the current socio-economic and natural circumstances, the woody plant encroachment of grasslands will continue. Therefore, if the objectives are to preserve the landscape, plant communities, biodiversity, and ecosystem functions and services, prompt action should be taken to preserve mountain grasslands. Our results provide a scientific basis for making political decisions on the management of these important regions.

Footnotes

Acknowledgements

We thank S. Benitez, H. Saiz, F. Fillat, and A. Aldezabal for assistance with the bibliographic search, for providing information about the study area, and for provided suggestions for improving earlier versions of the manuscript. Bruce MacWhirter and Scott Butler provided helpful suggestions on the penultimate version of manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

This research was funded within the framework of the research project CGL2011-27259 (Spanish Ministry of Economy and Competitiveness and Innovation), co-financed by the FEDER, project DIPA 125/2010 MMAMRM (Spanish National Park organization), and project FW7 ENV.2009.2.1.3.2 - LEDDRA (European Community).

References

AEMET (2012) Spanish Meteorological Agency website. Available at : http://www.aemet.es.

Antrop

(2005) Why landscapes of the past are important for the future. Landscape and Urban Planning 70: 21–34.

Archer

Hoffman

Danckwerts

(1989) How economic are the farming units of Leliefontein, Namaqualand? Journal of the Grassland Society of Southern Africa 6: 211–215.

Baret

Guyot

(1991) Potentials and limits of vegetation indices for LAI and APAR assessment. Remote Sensing of Environment 35: 161–173.

Bartolome

Plaixats

Fanlo

Boada

(2005) Conservation of isolated Atlantic heathlands in the Mediterranean region: Effects of land-use changes in the Montseny biosphere reserve (Spain). Biological Conservation 122: 81–88.

Batllori

Gutiérrez

(2008) Regional tree line dynamics in response to global change in the Pyrenees. Journal of Ecology 96: 1275–1288.

Bolli

Rigling

Bugmann

(2007) The influence of changes in climate and land-use on regeneration dynamics of Norway spruce at the treeline in the Swiss Alps. Silva Fennica 41: 55–70.

Boyd

Danson

(2005) Satellite remote sensing of forest resources: Three decades of research development. Progress in Physical Geography 29(1): 1–26.

Brandt

Haynes

Kuemmerle

. (2013) Regime shift on the roof of the world: Alpine meadows converting to shrublands in the southern Himalayas. Biological Conservation 158: 116–127.

10.

Breiman

(2001) Random forests. Machine Learning 45: 5–32.

11.

Bryn

(2008) Recent forest limit changes in south-east Norway: Effects of climate change or regrowth after abandoned utilisation? Norsk Geografisk Tidsskrift – Norwegian Journal of Geography 62: 251–270.

12.

Camarero

Gutiérrez

(2004) Pace and pattern of recent treeline dynamics: Response of ecotones to climatic variability in the Spanish Pyrenees. Climatic Change 63: 181–200.

13.

Cohen

(1960) A coefficient of agreement for nominal scales. Educational and Psychological Measurement 20: 37–46.

14.

Congalton

Green

(2009) Assessing the Accuracy of Remotely Sensed Data: Principles and Practices. Boca Raton, FL: CRC Press.

15.

Cuadrat

Saz

Vicente-Serrano

(2007) Atlas Climático de Aragón. Zaragoza: Gobierno de Aragón.

16.

Cuevas

(2000) Tree recruitment at the Nothofagus pumilio alpine timberline in Tierra del Fuego, Chile. Journal of Ecology 88: 840–855.

17.

Daumas

(1976) La vie rurale dans le haut Aragón oriental. Instituto de Estudios Oscenses y de Geografía Aplicada. Madrid: CSIC.

18.

De la Riva

(1997) Los Montes de la Jacetania: Caracterización Física y Explotación Forestal. Zaragoza: Consejo de Protección de la Naturaleza de Aragón.

19.

Didier

(2001) Invasion patterns of European larch and Swiss stone pine in subalpine pastures in the French Alps. Forest Ecology and Management 145: 67–77.

20.

Felicísimo

(1994) Modelos Digitales del Terreno: Introducción y Aplicaciones en las Ciencias Ambientales. Oviedo: Pentalfa Ediciones.

21.

Fernández-Giménez

Fillat

(2012a) Pyrenean pastoralists’ ecological knowledge: Documentation and application to natural resource management and adaptation. Human Ecology 40: 287–300.

22.

Fernández-Giménez

Fillat

(2012b) Pyrenean pastoralists’ observations of environmental change: An exploratory study in Los Valles Occidentales of Aragón. Pirineos, Revista de Ecología de Montaña 167: 143–163.

23.

Fillat

García-González

Gómez

Reiné

(2008) Pastos del Pirineo. Madrid: Departamento de Publicaciones, CSIC.

24.

García-Ruiz

Lasanta

Ruiz-Flano

. (1996) Land-use changes and sustainable development in mountain areas: A case study in the Spanish Pyrenees. Landscape Ecology 11: 267–277.

25.

Gartzia

Alados

Pérez-Cabello

Bueno

(2013) Improving the accuracy of vegetation classifications in mountainous areas. Mountain Research and Development 33(1): 63–74.

26.

Gautam

Webb

Shivakoti

Zoebisch

(2003) Land use dynamics and landscape change pattern in a mountain watershed in Nepal. Agriculture, Ecosystems and Environment 99: 83–96.

27.

Gehrig-Fasel

Guisan

Zimmermann

(2007) Tree line shifts in the Swiss Alps: Climate change or land abandonment? Journal of Vegetation Science 18: 571–582.

28.

Gellrich

Baur

Koch

Zimmermann

(2007) Agricultural land abandonment and natural forest re-growth in the Swiss mountains: A spatially explicit economic analysis. Agriculture, Ecosystems and Environment 118: 93–108.

29.

Gerard

Petit

Smith

. (2010) Land cover change in Europe between 1950 and 2000 determined employing aerial photography. Progress in Physical Geography 34: 183–205.

30.

Grau

Hernández

Gutierrez

. (2008) A peri-urban neotropical forest transition and its consequences for environmental services. Ecology and Society 13: 35.

31.

Guerrero

Alberto

Hodgson

. (1999) Plant community patterns in a gypsum area of NE Spain: Interactions with topographic factors and soil erosion. Journal of Arid Environments 41: 401–410.

32.

Hardin

(1968) The tragedy of the commons. Science 162: 1243–1248.

33.

Holtmeier

Broll

(2007) Treeline advance – driving processes and adverse factors. Landscape Online 1: 1–33.

34.

IAE (2012) Aragón Institute of Statistics website. Available at : http://www.aragon.es.

35.

IGN (2012) Spanish National Geographic Institute website. Available at : http://www.ign.es.

36.

Knapp

Briggs

Collins

. (2008) Shrub encroachment in North American grasslands: Shifts in growth form dominance rapidly alters control of ecosystem carbon inputs. Global Change Biology 14: 615–623.

37.

Komac

Alados

Camarero

(2011) Influence of topography on the colonization of subalpine grasslands by Echinospartum horridum . Arctic, Antarctic, and Alpine Research 43: 601–611.

38.

Kouba

Alados

(2012) Spatio-temporal dynamics of Quercus faginea forests in the Spanish Central Pre-Pyrenees. European Journal of Forest Research 131: 369–379.

39.

Kuemmerle

Radeloff

Perzanowski

Hostert

(2006) Cross-border comparison of land cover and landscape pattern in Eastern Europe using a hybrid classification technique. Remote Sensing of Environment 103: 449–464.

40.

Lasanta-Martínez

Vicente-Serrano

(2007) Cambios en la cubierta vegetal en el Pirineo aragonés en los últimos 50 años. Pirineos, Revista de Ecología de Montaña 162: 125–154.

41.

Lasanta-Martínez

Vicente-Serrano

Cuadrat-Prats

(2005) Mountain Mediterranean landscape evolution caused by the abandonment of traditional primary activities: A study of the Spanish Central Pyrenees. Applied Geography 25: 47–65.

42.

Liaw

Wiener

(2002) Classification and regression by randomForest. R News 2(3): 18–22.

43.

Liverman

Cuesta

RMR

(2008) Human interactions with the Earth system: People and pixels revisited. Earth Surface Processes and Landforms 33: 1458–1471.

44.

Lloret

Calvo

Pons

Díaz-Delgado

(2002) Wildfires and landscape pattern in the Eastern Iberian Peninsula. Landscape Ecology 17: 745–759.

45.

Maestre

Puche

Bowker

. (2009) Shrub encroachment can reverse desertification in semi-arid Mediterranean grasslands. Ecology Letters 12: 930–941.

46.

MAGRAMA (2010) Crop and land use maps. Available at: http://www.magrama.gob.es.

47.

Millington

Perry

Romero-Calcerrada

(2007) Regression techniques for examining land use/cover change: A case study of a Mediterranean landscape. Ecosystems 10: 562–578.

48.

Molinillo

Lasanta

García-Ruiz

(1997) Managing mountainous degraded landscapes after farmland abandonment in the Central Spanish Pyrenees. Environmental Management 21: 587–598.

49.

Montané

Casals

Taull

. (2010) Spatial patterns of shrub encroachment in neighbouring grassland communities in the Pyrenees: Floristic composition heterogeneity drives shrub proliferation rates. Plant Ecology 211: 267–278.

50.

Montserrat-Martí

(1992) Evolución glaciar y postglaciar del clima y la vegetación en la vertiente sur del Pirineo: Estudio palinológico (Evolution of glacial and postglacial climate in the meridional Pyrenees). Zaragoza: Monografías del Instituto Pirenaico de Ecología, CSIC.

51.

Moreira

Rego

Ferreira

(2001) Temporal (1958–1995) pattern of change in a cultural landscape of northwestern Portugal: Implications for fire occurrence. Landscape Ecology 16: 557–567.

52.

Motta

Morales

Nola

(2006) Human land-use, forest dynamics and tree growth at the treeline in the western Italian Alps. Annals of Forest Science 63: 739–747.

53.

Mottet

Ladet

Coqué

Gibon

(2006) Agricultural land-use change and its drivers in mountain landscapes: A case study in the Pyrenees. Agriculture, Ecosystems and Environment 114: 296–310.

54.

Naito

Cairns

(2011) Patterns and processes of global shrub expansion. Progress in Physical Geography 35(4): 423–442.

55.

Pinheiro

Bates

DebRoy

. (2011) nlme: Linear and nonlinear mixed effects models. R package version 3.1–102.

56.

Pôças

Cunha

Marcal

ARS

Pereira

(2011) An evaluation of changes in a mountainous rural landscape of northeast Portugal using remotely sensed data. Landscape and Urban Planning 101: 253–261.

57.

Poyatos

Latron

Llorens

(2003) Land use and land cover change after farmland abandonment: The case of a mediterranean mountain area (Catalan Pre-Pyrenees). Mountain Research and Development 23: 362–368.

58.

Ratajczak

Nippert

Collins

(2012) Woody encroachment decreases diversity across North American grasslands and savannas. Ecology 93(4): 697–703.

59.

R Development Core Team (2011) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available at: http://www.R-project.org.

60.

Rodriguez-Galiano

Ghimire

Rogan

. (2012) An assessment of the effectiveness of a random forest classifier for land-cover classification. ISPRS Journal of Photogrammetry and Remote Sensing 67: 93–104.

61.

Roura-Pascual

Pons

Etienne

Lambert

(2005) Transformation of a rural landscape in the eastern Pyrenees between 1953 and 2000. Mountain Research and Development 25: 252–261.

62.

Sankey

Montagne

Graumlich

. (2006) Lower forest–grassland ecotones and 20th century livestock herbivory effects in northern Mongolia. Forest Ecology and Management 233: 36–44.

63.

Schirmel

Mantilla-Contreras

Blindow

Fartmann

(2011) Impacts of succession and grass encroachment on heathland Orthoptera. Journal of Insect Conservation 15: 633–642.

64.

Schulz

Cayuela

Rey-Benayas

Schröder

(2011) Factors influencing vegetation cover change in Mediterranean Central Chile (1975–2008). Applied Vegetation Science 14: 571–582.

65.

Shalaby

Tateishi

(2007) Remote sensing and GIS for mapping and monitoring land cover and land-use changes in the northwestern coastal zone of Egypt. Applied Geography 27: 28–41.

66.

Shi

Zhang

Shi

Gao

(2009) Investigations on methods of land cover classification of TM image in mountain area. In: Zhang

Liu

. (eds) SPIE Proceedings Volume 7471: Second International Conference on Earth Observation for Global Changes, 25–29 May, China. Bellingham, WA: SPIE.

67.

SITAR (2012) Land Information System of Aragón. Available at : http://sitar.aragon.es.

68.

Sitko

Troll

(2008) Timberline changes in relation to summer farming in the western Chornohora (Ukrainian Carpathians). Mountain Research and Development 28: 263–271.

69.

Stow

Hope

McGuire

. (2004) Remote sensing of vegetation and land-cover change in Arctic tundra ecosystems. Remote Sensing of Environment 89: 281–308.

70.

Stueve

Isaacs

Tyrrell

Densmore

(2011) Spatial variability of biotic and abiotic tree establishment constraints across a treeline ecotone in the Alaska Range. Ecology 92: 496–506.

71.

Tømmervik

Johansen

Tombre

. (2004) Vegetation changes in the nordic mountain birch forest: The influence of grazing and climate change. Arctic, Antarctic, and Alpine Research 36: 323–332.

72.

Vicente-Serrano

Lasanta

Romo

(2005) Analysis of spatial and temporal evolution of vegetation cover in the Spanish Central Pyrenees: Role of human management. Environmental Management 34: 802–818.

73.

Wehn

Pedersen

Hanssen

(2011) A comparison of influences of cattle, goat, sheep and reindeer on vegetation changes in mountain cultural landscapes in Norway. Landscape and Urban Planning 102: 177–187.

74.

Wood

(2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1): 3–36.

75.

Zuur

Leno

Walker

. (2009) Mixed Effects Models and Extensions in Ecology with R. Berlin: Springer.