Research area and species

The research area is situated in the humid part of the calcareous Majella massif (reaching 2793 m asl and 32 km from the Adriatic Sea; Figure 1), which has mean annual precipitation of about 1200 mm in Majella National Park (MNP) and 1470 mm at the maritime face of the massif (van Gils et al 2013). Winter precipitation falls as snow above 1400 m asl. The elevational bioclimatic gradient of the Majella massif ranges from arctic-alpine to temperate-montane (uninhabited beech forest belt from 1000 to 1700 m asl) through sub-Mediterranean (broadleaf deciduous oak, hop hornbeam, manna ash forest, and farmland mosaic) to the Mediterranean biome (olive orchards and evergreen oak forest patches) (van Gils et al 2013). P. mugo krummholz covers 14 km² in 3 major patches (Mt Rapina, Blockhaus, and Mt Ugni; Figure 1) at elevations of 1221–2640 m (Table 1), that is, below and above the beech forest treeline. Summer pastoralism in the central Apennines dates back millennia (Forni 1998). Pasture has expanded by cutting and burning P. mugo (Migliaccio 1966). However, empirical observations have shown that P. mugo krummholz is resilient against burning (Leys et al 2014); this is consistent with observations of recuperating krummholz in the Majella massif at the Colle Ingotto (2000–2200m) at the fire scar of the wildfire in 2000 (personal observation by Hein van Gils and Giampiero Ciaschetti; 4 October 2010). Pastoralism declined sharply after World War II and then gradually to the edge of extinction today (Palombo et al 2013; van Gils et al 2013). The herders observed during the field study have abandoned the research area since then (personal communication, Giampiero Ciaschetti; 16 October 2016). The human population has almost halved in the municipalities of the MNP since the 1950s (Comuni-Italiani 2010). Evidence of abandoned pastoral practices and structures is abundant in the MNP (Micati 2000; van Gils et al 2008; van Gils et al 2010; Manzi 2016).

FIGURE 1 

Location of Majella National Park, with P. mugo krummholz distribution (Blasi et al 1999) and the 3 research sites: Mt Rapina (14°4′E / 42°8′N ∼ 14°6′E / 42°6′N), Blockhaus (14°6′E / 42°9′N ∼ 14°9′E / 42°7′N), and Mt Ugni (14°09′E / 42°07′N). (Map by Li Dai and Hein van Gils)

TABLE 1 

Extent and elevation of P. mugo krummholz in 1954, expansion during 1954–2007 at the Mt Rapina and Blockhaus sites (Figure 1), and predicted expansion from 2007 to 2060 by the best-fitting model at probability 0.6 at Mt Rapina.

P. mugo is a prostrate (<5 m) multistemmed, invasive (Kasak et al 2015), woody plant. Its seeds are dispersed by gravity, wind, birds, and seed-caching rodents (Müller-Schneider 1986); layering is common. Seedlings are scarce at close proximity to a seed source (<20 m) and negligible at further distance (>100 m) (Dullinger et al 2004). P. mugo krummholz is widespread in the Alps and the Balkans but rare in the Apennines. When present, P. mugo dominates the treeline (Dullinger et al 2004; Boratynska et al 2005). It is a pioneer with a fragmented distribution (Ali et al 2006). One of its most geographically isolated locations is the Majella massif (Boratynska et al 2005). P. mugo invades avalanche chutes and abandoned subalpine pastureland. This invasion poses a threat to endemic herbaceous species and influences local natural resource management (Dullinger et al 2004; Hoch and Körner 2009).

Change detection

The color 2007 APs served as the reference for ortho-rectification of the scanned black-and-white 1954 APs (Table 2). Ground truth was obtained in September 2009 by visiting randomly generated points with at least 60-m intervals within a 30-m zone on both sides of hiking tracks passing through the P. mugo patches at Mt Rapina (n=32) and Blockhaus (n=79).

TABLE 2 

Input materials for change detection and expansion modeling.

We recorded 4 vegetation cover classes within a plot of 10-m radius at each point: beech forest (>30% canopy cover; tree height >5 m), P. mugo krummholz (>30% canopy cover), low vegetation (<30% woody cover; >30% herbaceous cover) and bare ground (<30% vegetation cover). A 5th class, shadow, was distinguished on the 2007 APs. Only 2 classes could be identified on the 1954 APs: forest (beech forest, P. mugo krummholz plus shadow) and nonforest (low vegetation plus bare ground).

An automated object-oriented classification computer program (eCognition) was applied to reduce the salt-and-pepper effect (Hay et al 2005) resulting from using high-resolution imagery (Pekkarinen 2002). The object-oriented procedure consisted of multiresolution segmentation and subsequent supervised land cover classification of the segments (Radoux and Defourny 2007).

The ground observations were used for the accuracy assessment of the classified 2007 APs. For the 1954 AP accuracy assessment, cover classes were assigned by visual inspection of the APs in an area with a 10-m radius for around 50 random points per site (Platt and Schoennagel 2009). Agreement between classified APs and ground truth was calculated in R (R-Development-Core-Team 2009). In addition, a weighted Kappa index of agreement was obtained by entering the proportion of each cover class as a weight in the confusion matrix calculation.

Cover changes from 1954 to 2007 were detected using a raster calculator in ArcMap and then intersected with the Digital Elevation Model (DEM) obtaining the upslope krummholz migration. The pixels changed from nonforest (1954) to P. mugo krummholz (2007) were assigned as krummholz expansion.

Expansion modeling

Mt Rapina was selected for modeling because of its regular bell-shaped distribution of P. mugo krummholz with elevation (Figure 2) and its superior cover-classification accuracies. A logistic regression model was developed using the generalized linear model method in the R environment (Fox 2002), using several relevant predictors (Table 2) according to the hypotheses. The generalized linear model allows a choice of link functions for the response variable; in this case a binomial response link was chosen. Neighborhoods of 5–100 m for krummholz, and 10–100 m for beech forest in the 1954 APs were tested as predictors, both as ground surface distance and as distance in 2-dimensional (2D) planar projection. In the krummholz expansion area, 2000 random points were generated and intersected with each potential predictor.

FIGURE 2 

Elevation distribution of P. mugo in 1954, in 2007, and expansion from 1954 to 2007 at Mt Rapina and Blockhaus. The dip at 1900–2000 m elevation in the 1954 distribution at Blockhaus is reflected in the 2007 curve and in the expansion curve. (Graphs by Li Dai)

The distribution of continuous predictors was tested for symmetry, and those with skewness > 0.5 were log-transformed. The variance inflation factor was calculated to assess colinearity of predictors. The significance of correlations between categorical predictors and krummholz expansion was calculated using the Pearson's chi-square test. After selecting predictors, models with single and multiple predictors were compared for goodness-of-fit using the area under the receiver operating characteristic curve (AUC) and parsimony.

The probability of expansion of krummholz at Mt Rapina for the period 2007–2060 was forecast with the best-fitting model, while substituting the residual grazing map used in modeling (1999) with its 2007 update and the 1954 APs with the 2007 APs. Grazing is the only component of pastoralism for which we have consistent spatiotemporal data at fine resolution; we therefore consider grazing to be a practical proxy of other dimensions of pastoralism, including trampling, de-bushing (eg van Gils et al 2014), fuelwood harvesting, and purposive burning.

Meteorology and dendrochronology

In the absence of meteorological stations at the research site, records from the only station at alpine elevation in the central Apennines (Campo Imperatore, 2137 m asl; 56 km northwest of Blockhaus) and from 2 nearer stations at lower elevations (Pescocostanzo, 1395 m asl; 26 km south of Blockhaus; Sant'Eufemia a Majella, 870 m asl; 26 km south of Blockhaus) were examined. Pearson's correlations between these 3 stations were calculated for the mean annual and monthly maximum, minimum, and mean temperatures and for monthly total precipitation. To investigate trends, a linear model in R was fit to seasonal temperature and summer precipitation. Seasonal temperature values used were winter (December–January–February), spring (March–April–May) and summer (June–July–August). Summer precipitation was computed as the sum of the July and August values. Because of snowmelt at the treeline in April and May (MNP 2016), moisture limitation for P. mugo recruitment seems unlikely in June and therefore we assessed July–August precipitation.

Dendrochronological analysis was conducted at Blockhaus (2000–2100 m) and Mt Ugni (2100–2200 m) in 3 plots (30 × 20 m) at 10-m intervals in each site (Figure 1). The 15–20 P. mugo specimens with the largest stem diameter were sampled by increment borer in each plot. Cores (n = 122) were dated by counting rings, and ring widths were measured using LINTAB equipment, a stereomicroscope coupled with the Time Series Analysis Program. Individual raw ring-width chronologies were cross-dated visually and checked statistically using the program COFECHA (Grissino-Mayer 2001). The successfully cross-dated individual raw chronologies (n = 90) were averaged to obtain a mean raw chronology per plot. Subsequently, Pearson's correlation was calculated for the mean raw chronologies of the 3 plots within each site.

The raw series were averaged per calendar year to generate a mean raw chronology per site, and Pearson's correlation was calculated for the mean raw chronologies of the 3 plots within each site. Successively the individual raw chronologies were averaged per calendar year to generate a mean raw chronology per site, and Pearson's correlation was calculated for the mean raw chronologies of the 3 plots within each site. Successively, the individual raw chronologies were standardized using ARSTAN (Cook and Holmes 1984) to remove variability due to tree age and nonclimatic site factors.

The variability (interannual trend) linked to climate was preserved by standardization of the raw data into indexed chronologies (mean value = 1.0) using smoothing splines, and then by dividing the ring width measurements by the values obtained from the fitted curve. By standardization, mean standard chronologies per plot were created and successively used for climate-growth analyses. In addition, mean standard chronologies per plot were created and then averaged per calendar year to generate master standard chronologies per site.

The mean standard chronologies of both sites were correlated with each other and with climate predictors using the linear model in R and examining the model success by the coefficient of determination (adjusted R²). Forward stepwise linear regression was then applied to the individual predictors with significant relations. The final model was matched with linear modeling assumptions. The relationships between standard dendrochronologies and the climatic data were also assessed by the Pearson's correlation; Student's t-test was used to determine its significance.

Change detection and krummholz expansion model

Overall accuracies of the AP classifications were satisfactory (>0.8) across classes, sites, and years. Krummholz expanded at Mt Rapina and Blockhaus (1.5% year^–1 and 1.4% year^–1 respectively), more than doubling its extent from 1954 to 2007 (Figure 3), largely within the 1700–2100 m belt (Figures 2; 3; Table 1). At Mt Rapina, the treeline migrated 35 m (2573 to 2608 m) and 65 m at Blockhaus (2442 to 2507 m) (Table 1). The elevation distribution and expansion curves of P. mugo krummholz at Mt Rapina and Blockhaus are bell-shaped with a dip from 1900–2000 m at the latter (Figure 2).

FIGURE 3 

P. mugo krummholz expansions from 1954 to 2007 at Mt Rapina and Blockhaus. (Maps by Li Dai)

A full model included solar radiation, which was not a significant predictor and thus dropped, while a reduced model was highly significant (P < 0.005). All neighborhood variables were associated with P. mugo expansion. Combining the best neighborhood predictors in the nested model yielded the proximity to krummholz (planar 2D) and the proximity to beech forest (surface distance) as the better distance predictors (Table 3). The remaining potential predictors showed no colinearity (variance inflation factor < 10). All categorical variables were significantly (P < 0.005) associated with P. mugo expansion.

TABLE 3 

Logistic regression coefficients for the significant predictors in the best-fitting model.

A single predictor, proximity to krummholz, well explained P. mugo expansion (AUC = 0.82). After stepwise elimination of predictors, the best-fitting and most parsimonious model (AUC = 0.87) showed that the probability of expansion significantly increased with the following conditions (Table 3): proximity to P. mugo krummholz, presence of residual grazing in 1999, decreasing elevation, a negative 10-m neighborhood effect of beech forest, and proximity to beech forest beyond this 10-m neighborhood. At the expansion probabilities of 0.7 and 0.6, the best-fitting model generated twice as many true positives as false positives and forecast that the P. mugo treeline will be 26–30 m higher in 2060 than in 2007.

Climate trends and dendrochronologies

The temperature and precipitation from Campo Imperatore, Pescocostanzo, and Sant'Eufemia were correlated significantly (P < 0.001). Pescoconstanzo was selected for further analysis due to its geographic proximity and relatively complete record. The pre–World War I records were discarded because of apparent instrumentation modification after the war. The spring (2.87°C century^–1; R² = 0.16) and winter (4.38°C century^–1; R² = 0.25) minima have increased since 1919. Summer minimum temperatures increased (3.17°C century^–1; R² = 0.16) and summer maximum temperatures decreased (–2.69°C century^–1; R² = 0.12). A long-term trend could not be detected in the remaining seasonal temperature values, mean annual temperature, or in summer rainfall. The extrapolated trend in spring minima is not predicted to reach the threshold temperature for growth (6 °C; Hoch and Körner 2009) until 2086 (Figure 4).

FIGURE 4 

Trend in spring mean minimum temperature from 1928 to 2007 and extrapolation to 2086 (solid red line) with confidence intervals (dashed red lines) and prediction intervals (dotted blue lines). The temperature values (Pescoconstanzo 1395 m) cannot be taken as representative for the krummholz sites above the treeline, only the trend. The actual temperature values may be lower by several degrees. The predicted minimum does not reach the plant-physiological threshold for growth of 6°C during the extrapolation period. The upper 95% interval reaches the threshold in 2011 and the upper 95% confidence interval in 2064. (Graph by David Rossiter)

The mean raw chronologies per plot at both sites were significantly (P < 0.01) correlated (Blockhaus: 0.45–0.81 and Mt Ugni: 0.22–0.62), allowing averaging per site. The P. mugo radial growths at Blockhaus and Mt Ugni were moderately (r = 0.58) and positively correlated for the same years, while showing no temporal pattern over the 20th century. A quarter of the P. mugo growth was predicted at both Blockhaus (R² = 0.25) and Mt. Ugni (R² = 0.22) from the spring maximum temperature and the summer rain (Figure 5). The warmer the spring or the wetter the summer, the larger the radial growth in the same calendar year.

FIGURE 5 

Master raw chronologies at Blockhaus (solid grey line) and Mt Ugni (solid black line) with their trend lines; the dotted lines show the number of core samples included in the record. (Graph by Caterina Palombo)

Krummholz expansion

The P. mugo krummholz expansion rate (1.4–1.5% year^–1) is similar to the beech forest expansion rate in the montane belt of MNP (1.2% year^–1) and forest expansion into abandoned farmland across Europe (van Gils et al 2008). The doubling of the krummholz extent in MNP over half a century corresponds to the statistics from Switzerland for all forests extending to the treeline over a 12-year period (Gehrig-Fasel et al 2007). The upslope migration of P. mugo krummholz in the central Apennines (6–11 m decade^–1) concurs with the positive trend detected by Palombo et al (2013) and is close to the krummholz expansion in the Alps (13 m decade^–1; Tattoni et al 2010). However, the range of the upslope treeline migration in the Alps is 25–50 m decade^–1 (Gehrig-Fasel et al 2007; Chauchard et al 2010).

The slower expansion of krummholz in the central Apennines may be explained by its higher elevation, given that the latter is a negative expansion predictor (Table 3). The single predictor model (proximity to P. mugo krummholz) may explain the expansion by the dispersal limitations of P. mugo. Further, the best-fitting multipredictor model (Table 3) suggests environmental limitations for P. mugo recruitment.

The positive influence of residual grazing in 1999 on P. mugo expansion conforms to findings that moderate livestock impact stimulates expansion of woody plants into pastures (van Gils et al 2008; Luo and Dai 2013), possibly in this case through livestock trampling creating a seedbed (Germino et al 2002) or by preclusion of a litter layer unsuitable as seedbed (Dullinger et al 2004). If we accept grazing as a proxy for pastoralism, associated de-bushing may also play a role (Piussi 1994; van Gils et al 2014).

The absence of P. mugo within close (10 m) proximity to beech forest may result from moisture depletion by beech roots (Ellenberg 1978) or the shade effect (Holtmeier 2009). At the same time, P. mugo expansion is more probable in the close proximity of beech forest beyond a 10 m beech forest neighborhood. This suggests that beech and P. mugo share site preferences. Proximity to P. mugo as an expansion predictor of the same species is found in the Alps (Dullinger et al 2004). The number of false (33) versus true (66) positives predicted by the best-fitting model suggests that P. mugo is slow to expand over suitable sites.

Our significant neighborhood and distance predictors (Table 3) indicate that the seed dispersal mode may have differential impacts on krummholz expansion. Gravity and wind dispersal may drive expansion at close proximity (<70 m). The weaker predictor “beyond 70 m krummholz neighborhood” may be associated with pine-seed dispersal by birds or seed-caching rodents. Further, both distance and neighborhood values appear beyond the range where layering may contribute to krummholz expansion.

Of the tested expansion predictors, 7 are redundant for the best-fitting model, namely stocking density in residual grazing, geology, terrain type, soil type, incoming solar radiation, slope, and aspect. The range in current stocking density (<1 versus >1 sheep ha^–1) seems too low for a differential impact on expansion. The substrate (geology, soil type) seems to be homogeneously shallow and calcareous in the study area. The topo-climate (incoming solar radiation, terrain type, slope, and aspect) and substrate variables were also found to be redundant for predicting beech forest expansion on the Majella massif (van Gils et al 2008). The redundancy of topo-climatic variables makes climate warming an unlikely driver of expansion. Instead, the successful single predictor and multipredictor models are primarily driven by proximity to a seed source, residual grazing in 1999, and avoidance of close (10 m) proximity to beech forest.

Our forecast of the krummholz expansion supports the claim that the regionally rare P. mugo krummholz is in no danger of extinction. However, the arctic-alpine vegetation invaded by the krummholz represents a plant biodiversity hotspot (Stanisci et al 2005; van Gils et al 2012). The forecasted further expansion of krummholz may threaten this botanical biodiversity and dependent fauna diversity (Kasak et al 2015).

Climate trends and dendrochronology

The trend and magnitude of rising minimum winter and spring temperatures found in this research is consistent with findings (+0.4°C decade^–1) at lower elevations in the Apennines (Izzo et al 2004; Kumar et al 2005). The declining summer maximum temperatures during the 20th century in the research area are similar to those measured elsewhere in the central Apennines and differ from trends in increased maximum temperatures from the northern Apennines. The lower elevation of the meteorological stations in the northern Apennines versus central Apennines (Pescoconstanzo 1395 m) may explain the difference, because at midlatitudes elevation positively associates with rising minimum temperatures in spring and winter (but see Hu et al 2014) as well as with less overall climate warming (Liu et al 2009; Lu et al 2010). The extrapolated trend of rising minimum spring temperatures does not reach the threshold for growth during the 21st century (Figure 4), not even in the beech forest belt at 1395 m. This suggests that the P. mugo growing season may not be extended during the period of the forecast.

The local long-term climate trend of the 20th century, higher winter and spring temperature minima, appears not to affect radial growth of P. mugo. The reason may be that the higher minima are still below freezing point in winter and well below the threshold for growth even in spring at the much lower elevation of Pescoconstanzo (1395 m). Consequently, the current trend of rising spring temperature minima may not result in increased radial growth in the 21st century (Figure 4). The stable climatic growth conditions for P. mugo over the 20th century and our forecast of local climate warming in terms of the plant physiological minimum temperature suggest that upslope migration of suitable krummholz temperature conditions is unlikely within the period of forecast.

Similar to our findings, the effects of current climate change on the elevational position of krummholz seem minor compared with the effect of land use change in midlatitudinal Eurasia (Palombo 2013; Geanta et al 2014; Schickhoff et al 2015). Moreover, the authors of the latter article conclude that krummholz is relatively insensitive to climate change.