INTRODUCTION
The role of climatic factors in shaping distribution patterns of living organisms is now unquestionable (Woodward & Williams Reference Woodward and Williams1987; Huntley et al. Reference Huntley, Bartlein and Prentice1989; Huntley Reference Huntley, Maltby, Holdgate, Acreman and Weir1999; Hughes Reference Hughes2000; Davis & Shaw Reference Davis and Shaw2001; Walther et al. Reference Walther, Post, Convey, Menzel, Parmesan, Beebee, Fromentin, Hoeg-Guldberg and Bairlein2002; Hamann & Wang Reference Hamann and Wang2006). While past and current climates have acted as a major influence on species distribution patterns, early signs of shifts in climate patterns challenge existing approaches to biodiversity conservation that inform laws and policies (Cliquet et al. Reference Cliquet, Backes, Harris and Howsam2009; Heller & Zavaleta Reference Heller and Zavaleta2009).
The current focus on protecting the ‘existing values’ demonstrated by species associations occurring within protected areas has been critiqued as a rather static approach (Dodd et al. Reference Dodd, Hardiman, Jennings and Williams2010). This approach ignores evidence of the individualistic nature of species’ response to environmental changes, which is likely to lead to communities that are compositionally unlike those found currently (Brubaker Reference Brubaker, West, Shugart and Botkin1981; Williams & Jackson Reference Williams and Jackson2007). Furthermore, local natural and anthropogenic disturbances (such as wildfires, grazing, management or land-use change) operate in synergy with global and regional shifts in climate, and play a major role in shaping patterns of vegetation distribution at regional and local scales (Suffling Reference Suffling1995; Forbes Reference Forbes1999; Weber et al. Reference Weber, Moloney and Jeltsch2000; Franklin et al. Reference Franklin, Syphard, Mladenoff, He, Simons, Martin, Deutschman and O'Leary2001).
Spatially explicit landscape scenarios of potential habitat modification that consider both regional climate change and local disturbances are needed to quantify changes in proxies (for example habitat loss, fragmentation, degradation or biomass) of ecosystem services (such as biodiversity or carbon sequestration) and can help to advance understanding of their effects on current biodiversity conservation approaches, at multiple scales. Scenario maps can be obtained by means of multiscale and multiprocess spatially explicit forest landscape simulation models (Scheller & Mladenoff Reference Scheller and Mladenoff2007; Xi et al. Reference Xi, Coulson, Birt, Shang, Waldron, Charles, Lafon, Cairns, Tchakerian and Klepzig2009). These models, considered ‘the essence of landscape ecological methods’ (Baker Reference Baker1989, p. 121), can simulate the long-term effects of ecosystem disturbances (for example forest harvesting and management, winds, wildfires, pests and diseases), and are suitable to multiple scale applications (Xi et al. Reference Xi, Coulson, Birt, Shang, Waldron, Charles, Lafon, Cairns, Tchakerian and Klepzig2009). Among models adapting stochastic methods and using integrated physical simulation methods to model key ecological processes, material cycles and energy flow (He Reference He2008), the LANDIS model has been used successfully in multiple contexts (He et al. Reference He, Mladenoff and Crow1999; Mladenoff Reference Mladenoff2004). This grid-based spatially explicit model can be used to simulate landscape-level forest changes, and natural (such as wildfires and windthrow events) and anthropogenic disturbance (for example logging) impacts. LANDIS-II includes variable time steps for different ecological processes, for example woody biomass and carbon dynamics (Sheller & Mladenoff Reference Scheller and Mladenoff2007), providing a platform to explicitly incorporate the role of global and regional climate change as a driving force, and to examine their coupled effects with local anthropogenic disturbances.
This study focuses on semi-natural grasslands (EEA [European Environment Agency] 2010), a globally threatened habitat type, and on the Mediterranean region. Grasslands harbour some of the most species-rich plant communities in Europe, hosting a remarkable set of endemic and orchid species (Wilson et al. Reference Wilson, Peet, Dengler and Pärtel2012). The Mediterranean region is known for its exceptionally high number (13000) of endemic plants (Myers et al. Reference Myers, Mittelmeier, Mittelmeier, Da Fonseca and Kent2000), yet constitutes a threatened conservation region where greater temperature increases, precipitation decreases and extreme events (such as floods, drought and heat waves) are anticipated (IPCC [Intergovernmental Panel on Climate Change] 2007).
The study area was in the Puglia region (southern Italy; Appendix 1, Fig S1, see supplementary material at Journals.cambridge.org/ENC), where agriculture imprints the landscape, forests cover only 9.25% of the land, and scattered non-forest semi-natural habitats, often contiguous or close to forests, occur within a dense network of protected areas. All these vegetation types are either of conservation concern (European Habitats Directive 2007) or protected by law. Local pressures have historically led to loss and fragmentation of grasslands, and climate change is expected to enhance this trend by inducing modifications of other habitat types.
In this study, we aimed to use the outputs of a LANDIS-II application to examine coupled impacts of local pressures and climate change on grasslands in connection with the distribution of tree species. The final goal was to derive suggestions for feasible grassland and forest conservation strategies considering vegetation dynamics at a landscape level, incorporating adaptation (IPCC 2007) into integrated grassland and forest management plans (Splittlehouse Reference Spittlehouse2005; Millar et al. Reference Millar, Stephenson and Stephens2007; Malmsheimer et al. Reference Malmsheimer, Heffernan, Brink, Crandall, Deneke, Galik, Gee, Helms, McClure, Mortimer, Ruddell, Smith and Stewart2008; Seppälä et al. Reference Seppälä, Buck and Katila2009).
METHODS
Study area
The study area (23800 ha, 16º 30′ 5.83″–16º 45′ 40.69″ W and 40º 48′ 13.81″–41º 0′ 5.68″ N, 380–520 m above sea level) is representative of the main vegetation types and associated threats of the Natura 2000 (N2K) site of community importance (SCI)/special protection area (SPA) Murgia Alta IT9120007 (European Habitats Directive 2007). The N2K site, also partly designated and managed as a national park (NP), spans over a wide karstic plateau (125 880 ha) surrounded by sedimetary plains under agriculture, and is subject to a Mediterranean macro-bioclimate with a pluviseasonal-oceanic bioclimate, ombrotypes ranging from dry to subhumid, and a meso-Mediterranean thermotype (Forte et al. Reference Forte, Perrino and Terzi2005). The N2K site constitutes a training site within the European Union (EU) Seventh Framework Programme (FP7/2007–2013) project BIO_SOS (www.biosos.eu), aimed at developing a pre-operational system for cost- and time-effective monitoring of land cover and habitat changes within and at the fringe of protected areas.
Semi-natural grassland (29800 ha, 24% of the total area), the ecosystem type of main conservation concern at this N2K site, was developed and maintained through a mix of anthropogenic and natural processes, including long periods of grazing by domestic stock, cutting and regimes of deliberate burning. From the 80000 ha that existed at the beginning of the 20th century, by the early 1990s, grassland extent had dropped to 40000 ha, due to transitions to other land uses (largely arable crops). Substantial losses further occurred between 1990 and 2000, mainly due to EU incentives promoting durum wheat production (Appendix 1, Fig S2, see supplementary material at Journals.cambridge.org/ENC) (Boccaccio et al. Reference Boccaccio, Labadessa, Leronni, Mairota, Calabrese, Pacucci, Occhialini and Russo2013). Pastures have also become prone to woody vegetation encroachment from neighbouring forest vegetation types, after the decline of traditional sheep grazing.
Other sub- and semi- natural vegetation types, some of which are also of conservation concern, include garrigues and scrubs, downy oak (Quercus pubescens Willd.) and kermes oak (Q. coccifera L.) coppice woodlands and conifer plantations. Forests are subject to an altered (anthropogenic rather than natural) fire regime, directly dependent on human intervention. A forest fire suppression policy is applied throughout the region, and prescribed burning is not in use.
More than one hundred animal species listed in European Directives or other red lists are associated with these vegetation types, which are cast within an agricultural matrix dominated by cereal crops and, to a lesser extent, comprised of olive groves, almond orchards and vineyards (Table 1).
Land cover/land use types and habitats of conservation concern in the study area. Sources: 1Regione Puglia (Reference Regione2008), 2 European Habitats Directive (2007) and standard data forms for the ‘Murgia Alta’ IT9120007 SCI/SPA Natura 2000 site.
Modelling
The reference LANDIS-II application (Leronni Reference Leronni2013) used for this study included six contrasting scenarios (Table 2), built under specified conditions of climate, tree-growth and dispersal, and disturbance regime, to assess the dynamics of the above ground biomass (AGB) of a set of Mediterranean forest tree species. In the first two scenarios, AGB trends were modelled under both current and projected climatic conditions (NCC_B versus CC_B). Antropogenic disturbances, namely fires (NCC_B_F versus CC_B_F) and management (NCC_B_F_M versus CC_B_F_M), were progressively introduced.
Scenarios built using LANDIS-II (Leronni Reference Leronni2013).
The initial forest structure and composition were derived from a combination of geographic information system (GIS) analysis and data collected in a 2011–2012 field survey of the study area based on stratified random sampling (n = 53). Four forest types (coniferous plantations, deciduous/semi-deciduous oak woodlands, evergreen oak woodlands, and shrubs) were identified consisting of 12 stand types (initial communities) with a unique combination of species-age cohorts (Leronni Reference Leronni2013).
Rates of establishment probabilities (SEP) and above ground net primary productivity (ANPP) for five tree species and shrub species, were computed through the ecosystem process model PnET-II (Aber & Federer Reference Aber and Federer1992; Aber et al. Reference Aber, Ollinger and Driscoll1997; Xu et al. Reference Xu, Gertner and Scheller2009). The values of the several individual species life-trait and ecophysiological variables were derived from the literature (Baldocchi et al. Reference Baldocchi, Hicks and Meyers1988; Ellenberg Reference Ellenberg1988; Bugmann Reference Bugmann1994; Aber et al. Reference Aber, Ollinger and Driscoll1997; Panaiotis et al. Reference Panaiotis, Carcaillet and M'Hamedi1997; Goodale et al.1998; Alifragis et al. Reference Alifragis, Smiris, Maris, Kavvadias, Konstantinidou and Stamou2001; Mediavilla & Escudero Reference Mediavilla and Escudero2003 a, b; Laurent et al. Reference Laurent, Bar-Hen, Francois, Ghislain and Cheddadi2004; Tapias et al. Reference Tapias, Climent, Pardos and Gil2004; Wright et al. Reference Wright, Reich, Westoby, Ackerly, Baruch, Bongers, Cavender-Bares, Chapin, Cornelissen, Diemer, Flexas, Garnier, Groom, Gulias, Hikosaka, Lamont, Lee, Lee, Lusk, Midgley, Navas, Niinemets, Oleksyn, Osada, Poorter, Poot, Prior, Pyankov, Roumet, Thomas, Tjoelker, Veneklaas and Villar2004; Niinemets & Valladares Reference Niinemets and Valladares2006; Fyllas et al. Reference Fyllas, Phillips, Kunin, Matsinos and Troumbis2007; Daas et al. Reference Daas, Montpied, Hanchi and Dreyer2008; Querejeta et al. Reference Querejeta, Barbera, Granados and Castillo2008; Paula et al. Reference Paula, Arianoutsou, Kazanis, Tavsanoglu, Lloret, Buhk, Ojeda, Luna, Moreno, Rodrigo, Espelta, Palacio, Fernández-Santos, Fernandes and Pausas2009; Xu et al. Reference Xu, Gertner and Scheller2009; Xi & Xu Reference Xi and Xu2010), while data relevant to the soil nutrient regime were derived from field surveys (n = 30). The study area was treated as a homogeneous ecoregion, based on the main chemical-physical characteristics of the soil and geomorphological features.
Historical monthly averages of climatic variables for the 1960–2010 period were derived from the records of four thermo-pluviometric stations close to the study area.
Climate change projections (2010–2060) were derived from the general circulation model GCM ECHAM5 SRES A1B (Roeckner et al. Reference Roeckner, Baeuml, Bonventura, Brokopf, Esch, Giorgetta, Hagemann, Kirchner, Kornblueh, Manzini, Rhodin, Schlese, Schulzweida and Tompkins2003), downscaled at RCM (regional circulation model, resolution 0.25°) for the Mediterranean basin (Christensen et al. Reference Christensen, Rummukainen, Lenderink, van der Linden and Mitchell2009), with atmospheric CO2 levels stabilizing at c.700 ppm by the year 2100, leading to a 1.6°C increase in mean annual temperature and a 1.3 mm decrease in total annual precipitation for the study area. Photosynthetically active radiation (PAR) intensity was indirectly estimated according to Hargreaves and Samani (Reference Hargreaves and Samani1982) and Weiss and Norman (Reference Weiss and Norman1985). These were extended to 2160, assuming the last time step (2055–2060) of PnET-II outputs were constant for the remaining part of the simulation.
The biomass succession extension was calibrated by comparing model outputs to the real biomass of each initial community, assessed by field survey (Scheller Reference Scheller2012).
The fire disturbance was simulated using the base fire extension (Scheller & Domingo Reference Scheller and Domingo2012), with the current fire regime specified by a mean observed fire rotation period (FRP) of 75 years, characterizing a moderate fire cycle in other Mediterranean vegetation types (Franklin et al. Reference Franklin, Syphard, Mladenoff, He, Simons, Martin, Deutschman and O'Leary2001). The calibration was performed through a systematic manual optimization of specified parameters, and comparison between observed average and modelled FRP values (Syphard et al. Reference Syphard, Yang, Franklin, He and Keeley2007).
The forest management was simulated through the biomass harvest extension (Scheller & Domingo Reference Scheller and Domingo2013), parameterized to mimic typical harvesting in the region based on forest management prescriptions. A moderate periodical thinning (10–20% biomass, 15 years rotation period) and typical coppicing with release (18–20 years cycle) were applied to coniferous plantations and deciduous oak woods, respectively. The prescription parameter ‘minimum time since damage’ was set to 10 years, according to current regulation, to prevent cutting in areas damaged by fires. Forest management included the elimination of encroaching vegetation on agricultural and other non-forest areas.
To assess within-scenario stochasticity, each scenario was replicated three times with a different random seed number. As replications yielded very low variances (2–3%) of response variables, scenario significance testing was not performed, and one replication was selected at random.
From this modelling (Appendix 1, Fig S3, see supplementary material at Journals.cambridge.org/ENC) the two most extreme scenarios (CC_B_F and CC_B_F_M, hereafter and in figures F and FM) were selected for a pairwise quantitative comparison of the spatially explicit dynamics of the biomass of a sub-set of the modelled Mediterranean forest species, namely Aleppo pine (Pinus halepensis, Miller), (Ph), downy oak, (Qp) and kermes oak (Qc) within grasslands.
The LANDIS-II output maps for the F and FM scenarios were overlaid in a GIS with a current land-cover/land-use map. Forty quasi-binary biomass (total and individual species) maps (grasslands versus encroachment) at the same grain and extent were obtained for the two scenarios, at five different points in time (year 2030 = t30, 2060 = t60, 2090 = t90, 2120 = t120 and 2150 = t150). These were used to quantify the area and rate of woody encroachment over the next 150 years. The ten total biomass maps were used to compute a set of landscape pattern indices (LPIs) considered suitable to capture the four spatial effects of fragmentation on habitat pattern: reduction in habitat amount, increase in number of habitat patches, decrease in sizes of habitat patches and increase in isolation of patches (Fahrig Reference Fahrig2003). These effects were assessed respectively by means of the percentage of landscape (PLAND), number of patches (NP), effective patch size (MESH) and patch cohesion index (COHESION), computed using FRAGSTATS 3.3 (McGarigal et al. Reference McGarigal, Cushman, Neel and Ene2002). To these, edge density (ED) and fractal index distribution (FRAC_AM) were added to evaluate effects of encroachment on grasslands habitat configuration in terms of the complexity of patch shape and edge length, also related to fragmentation. To explore the relationships between the habitat loss and the configuration components of fragmentation (McGarigal & McComb Reference McGarigal and McComb1995; Mairota et al. Reference Mairota, Cafarelli, Boccaccio, Leronni, Labadessa, Kosmidou and Nagendra2013) as well as to summarize the set of LPIs, we applied a principal components analysis (PCA).
LANDIS-II output tables were used to quantify encroaching biomass, density of colonizing vegetation (Mg ha−1) and annual average biomass encroaching rates (Mg ha−1 yr−1).
RESULTS
A progressive increase in forest area (from 16% to 35% of the landscape for F scenarios and to 23% of the landscape area for FM scenarios) and a loss to encroachment of almost one-third of grassland habitat (from 22% to 16%) (Fig. 1 a; Appendix 1, Fig S4, see supplementary material at Journals.cambridge.org/ENC) were observed over the simulation period. No significant differences (t-test, α = 0.05) were observed between the two scenarios in terms of changes in grassland habitat area nor in the behaviour of the LPIs describing grasslands fragmentation (Fig. 1).
Trends in (a) percentage of landscape (PLAND); (b) number of patches (NP); (c) effective mesh size (MESH) and patch cohesion (COHESION); (d) edge density (ED) and fractal index distribution (FRAC_AM) for grasslands. Climate change and wildfires (F), climate change scenario, wildfires and management (FM) scenario.
The first two PCA components (PC1, PC2) explaining 98% of the variance were considered for both scenarios (Table 3). PC1 is associated with the LPIs representing the four spatial effects of the process of fragmentation (PLAND, NP, MESH and COHESION), and PC2 with the two ancillary LPIs (ED and FRAC_AM). Increasing fragmentation of the grassland habitat in the landscape at different time steps was indicated by an increase in the degree of the association (PCA scores for cases; Table 3) of such landscapes to PC1 (r2 = 0.759 and r2 = 0.669, respectively, for the two scenarios), with PLAND, MESH and COHESION decreasing and NP increasing. A similar increase in the association between PC2 and both ED and FRAC_AM was only observed from the second time step (t30–t60) onwards, with the values of these two LPIs progressively decreasing. This seems to indicate a change in phase (sensu Forman Reference Forman1995 and Jaeger Reference Jaeger2000) of the fragmentation process from time step 2 (t30–t60). Thus, the three landscapes where habitat is most fragmented (corresponding to t90, t120 and t150) show higher values of association with both PC1 and PC2.
Summary of the principal components analysis by variable for both the climate change and fire scenario and the climate change, fire and management scenario.
The densities of encroaching woody biomass were rather low (3.09 Mg ha−1 and 1.55 Mg ha−1, respectively, for the two scenarios) but differed (t-test, α = 0.01) (Fig. 2 a, b).
(a) Trends and (b) variability in the total amount of woody biomass encroachment per unit area; (c) trends and (d) variability in average annual rates of total amount of woody biomass encroachment and area of woody encroachment in grasslands. F = climate change and wildfires scenario, FM = climate change, wildfires and management scenario.
Encroachment proceeded at different rates in terms of area in the two scenarios. In both cases a decreasing trend was observed (Fig. 2 c, d), however the biomass trend was different (t-test, α = 0.01). The average rate of encroachment of woody biomass increased for the F scenario, while it decreased under the FM scenario, associated with a lower variability in the average rate of biomass encroachment.
Different species showed different colonization trends, amounts and rates of both encroached area and biomass. Moreover, a spatially segregated encroachment by species was observed, and different species varied in spatial encroachment patterns, being either clumped or relatively dispersed (Appendix 1, Fig S5, see supplementary material at Journals.cambridge.org/ENC). This was particularly observed for Ph and Qp.
Overall, the area encroached by these two species was much lower than the area encroached by Qc (Fig. 3 a, b). This holds for both scenarios, and is also reflected in the magnitude of the average annual encroachment rates (Fig. 3 c, d). The area encroached by Ph decreased with time, and this species was clumped close to the margins of plantations. The area encroached by Qp increased with time, and the species was progressively less clumped, occurring at a relatively greater distance from existing oak coppices. The area encroached by Qc increased with time, and progressively occupied the central part of the study area. Management did not seem to affect the area's average annual encroachment rates for any species.
Trends in (a, b) overall area of woody encroachment and (c, d) average annual rates of change in area of woody encroachment for selected tree species. Ph = Aleppo pine, Qp = downy oak, Qc = kermes oak. F = climate change and wildfires scenario, FM = climate change, wildfires and management scenario.
The contribution of Ph and Qp to encroaching biomass was much lower than that contributed by Qc. This holds for both scenarios, and is also reflected in the magnitude of the average annual encroachment rates (Fig. 4 a, b). Management significantly affected encroaching biomass for all species but Ph (Fig. 4 c, d), as this tended to be lower in the FM scenario. However, contrasting trends were observed for different species. Ph exhibited a decreasing trend in both scenarios, while Qc depicted a positive trend. Qp showed a more complex pattern of response, decreasing in the F scenario and increasing in the FM scenario. Management seemed also to contrastingly affect biomass average annual encroachment rates (Fig. 4 e, f). While no effect was observed in the case of Ph, for Qc and Qp an opposite trend was shown, with the former increasing and the latter decreasing.
Trends in (a, b) biomass of woody encroachment per unit area, (c, d) variability of encroached woody biomass per unit area, and (e, f) variability of average annual rates of encroached woody biomass, for selected tree species. Ph = Aleppo pine, Qp = downy oak, Qc = kermes oak. F = climate change and wildfires scenario, FM = climate change, wildfires and management scenario.
DISCUSSION
Woody encroachment into grasslands has constituted an important global threat for this habitat type in many parts of the world over the past 150 years (Chopping et al. Reference Chopping, Su, Rango, Martonchik, Peters and Laliberte2008) altering the capacity of grassland ecosystems to provide essential ecological, hydrological and socioeconomic functions (Hudak & Wessman Reference Hudak and Wessman1998; Moleele et al. Reference Moleele, Ringrose, Matheson and Vanderpost2002; Huxman et al. Reference Huxman, Wilcox, Breshears, Scott, Snyder, Small, Hultine, William, Pockman and Jackson2005). The simulated changes confirm that, at the local scale, too major impacts might occur at both habitat and landscape levels.
At the habitat level, with regard to biodiversity, close to one-third of the current grassland habitat appears threatened by encroachment. Moreover, as indicated by LPIs, fragmentation of habitat pattern might change in character from an earlier to a later phase (for example, from the dissection to the dissipation phase, sensu Jaeger Reference Jaeger2000). Simulated changes also indicate that part of the grassland habitat is vulnerable to qualitative changes capable of affecting ecosystem functioning, as changes in plant community structure might lead to scrubby or semi-wooded vegetated habitat types. Modifications in vegetation composition and structure due to woody encroachment are likely to determine changes in habitat quality and render it progressively unsuitable to Italian feather (Stipa austroitalica Martinovsky), a priority species for conservation, and to a number of other grassland plant species, such as endemic orchids (such as the Apulian spider-orchid, Ophrys holosericea subsp. apulica [Danesh & Danesh] Buttler). Moreover, grassland habitat deterioration and loss could have cascade effects across food webs at the landscape level, which would impact many animal species that are locally, regionally or globally threatened, including butterflies (such as the Italian marbled white, Melanargia arge), grasshoppers (like the spiked magician, Saga pedo) and birds (such as the lesser kestrel, Falco naumanni and calandra lark, Melanocorypha calandra).
Ecosystem level changes (for example carbon sequestration efficiency) might also occur at the habitat level. Grasslands in Europe appear to sequester 3–4 times more carbon in the soil than forests (Ciais et al. Reference Ciais, Soussana, Vuichard, Luyssaert, Don, Janssens, Piao, Dechow, Lathière, Maignan, Wattenbach, Smith, Ammann, Freibauer and Schulze2010). Even though losses in soil organic carbon during secondary succession seem to be greater with increasing precipitation (Alberti et al. Reference Alberti, Leronni, Piazzi, Petrella, Mairota, Peressotti, Piussi, Valentini, Gristina, La Mantia, Novara and Rühl2011 and literature cited therein), a negative relationship between carbon soil stock and temperature cannot be excluded in areas with a meso-Mediterranean thermotype (sensu Rivas-Martínez Reference Rivas-Martínez2009), although this is contingent on other factors related to successional phases (such as land-use/cover legacies and litter quality) (La Mantia et al. Reference La Mantia, Gristina, Rivaldo, Pasta, Novara and Rühl2013).
At the landscape level, the overall mosaic pattern might be altered progressively as grassland decrease is paralleled by a forest cover increase. This might be possibly associated with compositional changes, as the gradual spread of Qc observed in encroached grasslands at the expense of the other two species could occur in other forests. The decline of Qp was expected, whereas given different life traits, seed dispersal and fire resistence strategies between Qc and Ph, a greater spatial spread and a more intense colonization of the latter could have been envisaged. LANDIS-II might not have been sensitive to sexual maturity, and effective seed dispersal distance as the same values for the two species were derived from the literature (Leronni Reference Leronni2013). However, a competitive advantage might have occurred for Qc because traits and parameters were reported having as higher values than those for Ph, such as fire and shade tolerance. Shade tolerance in particular plays a very important role in both real and LANDIS-II simulated forest succession, as it affects the ability of a tree to reproduce via seeds or resprouting (He & Mladenoff Reference He and Mladenoff1999 a, b). The probability of Qc establishment also appears to be 1.6 times higher than that of Ph (Leronni Reference Leronni2013). The pattern observed for Qc is compatible with other findings in Mediterranean-type ecosystems (Franklin et al. Reference Franklin, Syphard, Mladenoff, He, Simons, Martin, Deutschman and O'Leary2001; Acácio et al. Reference Acácio, Holmgren, Rego, Moreira and Mohren2009). This is consistent with the life strategies of this species, which appear likely to withstand climate change-induced increases in the intensity and duration of the water stress period. Qc typically concludes twig growth and bud formation between April and May, and it might therefore more successfully escape the drought (and related heat waves) which strikes particularly around the end of June and beginning of July, when Ph is still growing and pollinating. Qc can also take advantage of a second growth season at the onset of the wet season, when seeds also ripen (Tsiouvaras et al. Reference Tsiouvaras1987). The decline of Ph, despite its capabilitities as an early seeder, could be explained in relation to the fire regime, as there is evidence that many fire-adapted species may be extirpated by particular fire regimes (Keeley et al. Reference Keeley, Bond, Bradstock, Pausas and Rundel2012). Ph is one such species, and plantations of this species are regarded as fire-sensitive because, due to high fire frequency and incidence of crown fire, they may not be able to recover a sufficient canopy seed bank for population regeneration after fire (Pausas et al. Reference Pausas, Llovet, Rodrigo and Vallejo2008).
Thus the question of the effectiveness of the protected area emerges for the conservation of at least two categories of ecosystem services provided by grasslands, biodiversity and carbon sequestration efficiency. The results also bring about other longer-term strategical issues related to the maintenance of the overall landscape arrangement and to the conservation of the currently protected Qp habitat in changing environmental conditions. A growing body of research casts doubts on the capacity of existing protected areas and networks to fulfil conservation objectives grounded on climatic and biogeographic stability in an era of global change, and a range of possible options are being proposed to adapt conservation to such changes (Lemieux & Scott Reference Lemieux and Scott2005; Hannah Reference Hannah2008; Heller & Zavaleta Reference Heller and Zavaleta2009; Wiens et al. Reference Wiens, Seavy and Jongsomjit2011). In this study, it seems that a site-scale conservation strategy would be more appropriate than other approaches (such as regional planning and landscape connectivity) considering that biogeographic controls, existing regulatory constraints and socioeconomic conditions hinder the departure from the fixed-boundary protected area strategy. The same constraints indicate that a combined short/middle term (resistence/resilience, sensu Millar et al. Reference Millar, Stephenson and Stephens2007) strategy to vegetation changes might be appropriate to the twofold aim of maintaining grasslands and buffering other landscape level modifications. For the study area, particularly in the context of climate change (Leronni Reference Leronni2013), fire regime appears to be a major driver of change. Current forest management, however, whilst not effective in buffering changes in habitat spatial pattern, significantly influences the degree and the rate of qualitative change in terms of both vegetation composition and structure. Therefore, such local-scale strategies should integrate fire control, improved silviculture and direct encroachment control practices.
Fire prevention in particular might be effective in increasing the fire rotation period. Besides improved patrolling, fire breaks could limit the spread of fire into forests from croplands and grasslands, where burning, however illegal, is traditionally applied during the fire season as a means of soil fertility and pasture improvement. Fire breaks management could be combined with current direct encroachment control practices and the resumption of sheep grazing. Fire prevention silviculture (for example pruning associated with moderate thinning; Notarnicola Reference Notarnicola2012) could be implemented in older pine plantations, where the primary carrier of crown fires are dead branches accumulated in the lower parts of the canopy. Even though an attenuated fire regime might favour Ph, a non-native species, at the expense of the native oaks as observed, its encroaching spatial potential is limited to the surroundings of plantations. This renders direct control easier. Qp's encroachment potential seems to be effectively limited by current coppicing. However, this kind of silviculture regime, in the context of climate change, also tends to eliminate Qp from existing woodlands if Qc were to become dominant (Leronni Reference Leronni2013). This does not happen in the F scenario, consistent with the empirical findings of Fabbio et al. (Reference Fabbio, Merlo and Tosi2003) who showed higher biomass per unit area in non-managed coppice woodlands, which tended to become fully stocked and multistoreyed. Therefore, the conversion of coppice woodlands to tall forests, at least in higher fertility stands, could balance the need to control encroachment and and to enhance the resilience of Qp woodlands in the landscape.
From the perspective of adaptive management, further modelling simulating the alternative management options outlined and the scenario analysis proposed for the assessment of their implications, might greatly aid the design of a feasible strategy to manage environmental change at the local scale (Millar et al. Reference Millar, Stephenson and Stephens2007; Shang et al. Reference Shang, He, Xi, Shifley and Palik2012). However, such potential is presently hampered in the study area by the lack of biometric data for model parameterization and calibration for the different silvicultural options. Moreover, with reference to the specific case of grassland encroachment, LANDIS-II does not currently allow for the modelling of grass vegetation. Therefore, competition between trees and grasses was not modelled in this study, nor were concurrent changes in grassland species composition caused by the encroachment of woody species. Similarly, other changes that are likely to be induced by climate change could not be evaluated, for example those expected to derive from range shifts of the herbaceous species in the present biotic communities and of other species that might enter and possibly become dominant or invasive species (Bakkenes et al. Reference Bakkenes, Alkemade, Ihle, Leemans and Latour2002; Martinez-Meyer Reference Martinez-Meyer2005).
Such complex dynamics represent a challenge for future research. On the one hand, developments of the LANDIS-II model could be aimed at the integration of grassland ecosystem dynamic simulation models (see Thornley Reference Thornley1998) or at exploring the potential for the use of LANDIS-II within the framework of hybrid approaches (Midgley et al. Reference Midgley, Davies, Albert, Alfwegg, Hannah, Hughes, O'Halloran, Seo, Thorne and Thuiller2010), which integrate landscape models and both spatially explicit species distribution models (Elith & Leathwick Reference Elith and Leathwick2009) and population viability models (Akçakaya et al. Reference Akçakaya, Burgman, Kindvall, Wood, Sjögren-Gulve, Hatfield and McCarthy2004). On the other hand, LANDIS-II could effectively be used in combination with remote sensing (RS) based landscape change detection technologies, suitable for investigations relevant to grasslands habitat quality (Muldavin et al. Reference Muldavin, Neville and Harper2001; Tong et al. Reference Tong, Wu, Yong, Yang and Yong2004), especially when combined with spatial modelling (Bradley & Mustard Reference Bradley and Mustard2006), and with continuous approaches for landscape pattern analysis (Pearson Reference Pearson2002; Hofer et al. Reference Hofer, Wagner, Herzog and Edwards2007). Targeted field management actions would also be more effective in minimizing the impact of environmental changes if coupled with RS monitoring aimed at providing early signals of environmental changes (Jump et al. Reference Jump, Cavin and Hunter2010). RS data and techniques have proved very useful in helping managers identify early warning signs of climate change at regional (Alcaraz-Segura et al. Reference Alcaraz-Segura, Cabello, Paruelo and Delibes2009; Silang et al. Reference Silang, Ronggao and Yang2010; Altamirano et al. Reference Altamirano, Field, Cayuela, Aplin, Lara and Rey-Benayas2010) and at local (Lucas et al. Reference Lucas, Bunting, Paterson and Chisholm2008) scales, as they are based on early identification of changes in plant physiology and phenology (Nagendra et al. Reference Nagendra, Lucas, Honrado, Jongman, Tarantino, Adamo and Mairota2013). In particular, mapping through RS, including LiDAR (light radar technology), very high spatial resolution datasets and advanced image processing techniques (Nagendra et al. Reference Nagendra, Lucas, Honrado, Jongman, Tarantino, Adamo and Mairota2013), can be extremely important in providing early indication of woody encroachment in grasslands (Chopping et al. Reference Chopping, Su, Rango, Martonchik, Peters and Laliberte2008).
CONCLUSION
Comparative scenario studies are required to facilitate adaptive management, to tackle the combined challenges of climate change and anthropogenic pressure that are being increasingly felt by protected areas across the world (Thompson et al. Reference Thompson, Wiek, Swanson, Carpenter, Fresco, Hollingsworth and Foster2012). This study demonstrates that LANDIS-II outputs can be effectively turned into scenario maps of land-cover/land-use useful to assess how human interventions interact with ecological dynamics in complex systems at local scales. These can also support larger-scale spatial modelling and planning (Hannah Reference Hannah2008, Rose & Burton Reference Rose and Burton2009; Wiens et al. Reference Wiens, Seavy and Jongsomjit2011).
The coupled effects of anthropogenic forcings on an ecological system operate at different spatial scales and seem likely to produce legacy effects at different temporal scales (such as different species biomass encroachment rates), thus adding to the complexity of the dynamics of human-ecological systems. To further understanding of the potential trends of these coupled systems, and to map areas susceptible to future change, current spatially explicit modelling tools capable of handling coupled forcings must be improved and approaches integrating these with RS developed (Liu et al. Reference Liu, Dietz, Carpenter, Folke, Alberti, Redman, Schneider, Ostrom, Pell, Lubchenco, Taylor, Ouyang, Deadman, Kratz and Provencher2007). RS, in particular, would provide a way to map spatial legacies, as well as early warnings of current change (Hudak & Wessman Reference Hudak and Wessman1998; Chopping et al. Reference Chopping, Su, Rango, Martonchik, Peters and Laliberte2008; Nagendra et al. Reference Nagendra, Lucas, Honrado, Jongman, Tarantino, Adamo and Mairota2013).
ACKNOWLEDGEMENTS
This study was funded partly by the PhD programme of the University of Bari Aldo Moro, Italy, and by the European Union Seventh Framework Programme (FP7/2007–2013) under project BIO_SOS, grant no. 263435, and was carried out in collaboration with the Department of Forest and Wildlife Ecology, University of Wisconsin, USA. We thank the Editor-in-Chief, the Associate Editor and three anonymous reviewers for their valuable comments on the manuscript. P. Mairota wishes to dedicate this work to the memory of Boštjan Anko, Professor of Landscape Ecology, University of Ljubliana, Slovenja, first President of the IUFRO-LE working group.
