Pinus mugo Krummholz Dynamics During Concomitant Change in Pastoralism and Climate in the Central Apennines

The dynamics of Pinus mugo krummholz during concomitant change in pastoral land use and climate in central Italy since the mid-20th century was investigated. Krummholz dynamics were detected using sequential aerial photography and fitted to a logistic regression model with elevation, grazing, proximity to beech forest, and proximity to krummholz as explanatory variables. Dendrochronological series were correlated with temperature and precipitation and fitted to a linear model. During this period krummholz doubled in extent and migrated 35–65 m upslope. Expansion was positively associated with krummholz proximity, residual pastoral grazing, and proximity to beech forest beyond 10 m and negatively associated with elevation and beech forest closer than 10 m. The logistic regression model forecasts krummholz migration by an additional 30 m upslope by 2060. During the 20th century, winter and spring minimum temperatures increased but did not result in increased radial stem growth of P. mugo. The combined evidence suggests that krummholz dynamics can be explained by the legacy of summer pastoralism and the dispersal limitations of P. mugo, rather than by climate change.


Introduction
The dynamics at the alpine treeline modify biodiversity (Stanisci et al 2005), alter landscapes, influence stakeholder benefits (Grace et al 2001), and have significant implications for the carbon cycle (van Gils et al 2008). From the Alps to the Himalayas the prevailing treeline in extra-tropical Eurasia is located below the climatic treeline due to centuries of downslope pasture extension (Ellenberg 1978;Blyakharchuk et al 2007;Miehe et al 2007;Fazlur-Rahman 2009;Tattoni et al 2010). Today's anthropogenic treelines include stunted Pinus mugo Turra, known as ''krummholz'' or ''crooked wood'' (Horvath et al 1974;Palombo et al 2013). Natural alpine treelines are rare in extra-tropical Eurasia (Gehrig-Fasel et al 2007;Dai et al 2013) but widespread on continents without a history of alpine pastoralism.
The contemporary dynamics of anthropogenic treelines is attributed to either abandonment of pastoralism (Dullinger et al 2004;Gehrig-Fasel et al 2007;Chauchard et al 2010;Tattoni et al 2010), climate change (Boden et al 2010;D ıaz-Varela et al 2010), or their combination (Motta et al 2006;Palombo et al 2014). An increase in ambient temperature raises the thermal boundary of vegetation growth and induces treeline dynamics (K€ orner 1998;Harsch et al 2009). Gottfried et al (2012) studied vegetation at the treelines of 60 summits in all major European mountain ranges and concluded that climate change gradually transformed mountain plant communities. Further, land use change, specifically abandonment of farmland and alpine pasture, has led to invasion by woody vegetation (Gehrig-Fasel et al 2007;van Gils et al, 2008;Palombo et al 2014).
Efforts to distinguish between land use and climatedriven dynamics at the treeline appear to be limited to locations of relatively low elevation (,2500 m asl) and high latitude (.458) as well as places where abandonment of pastoralism predates 20th-century climate change. In contrast, the P. mugo krummholz treeline in the central Apennines is at higher elevation (.2500 m asl) and lower latitude (428 N) -and is situated in the Mediterranean basin where a reduction in summer rainfall and aboveaverage warming is predicted (IPCC 2007) and where abandonment of pastoralism is concomitant with 20thcentury climate change (Palombo et al 2013). Treelines above 2500 m asl are common in midlatitudinal  Eurasia. The dynamics of anthropogenic treelines have been correlated to climate changes by spatial overlay (Gehrig-Fasel et al 2007;D ıaz-Varela et al 2010), regression (Dullinger et al 2004), and circumstantial evidence (Motta et al 2006;Chauchard et al 2010). These studies have extrapolated climate change from meteorological measurements outside the research sites at elevations below the treeline, generally without validation for the study sites. In cases where validation of the extrapolation has been attempted (Dullinger et al 2004), it only covers a fraction of the modeling time span. Moreover, recent research suggests that the pattern of climate change at midlatitudinal high elevations, such as our research site, differs from that at lower elevations (Liu et al 2009;Lu et al 2010). Finally, no previous research provides a spatially explicit forecast of krummholz migration direction and rate.
Remote sensing images, especially from spaceborne sensors, have not been widely used for the treeline studies due to their relatively recent availability (from late 1970s) and the coarse spatial resolution (approximately 80 m) of early imagery. Aerial photographs (APs) are appropriate for historical treeline change detection because of their longer availability (from the early 1950s) and finer spatial resolution (,1 m) (Okeke and Karnieli 2006). A successful method to relate tree growth to local climate at very high spatial and temporal resolutions is dendrochronological analysis (Hoch and K€ orner 2009); this is especially applicable where the historical climate record is unavailable at high elevations. Treeline change studies using both sequential APs and dendrochronological analysis have seldom been reported in the scientific literature.
Therefore, the aim of this study was to identify, model (second half of the 20th century), and forecast (first half of the 21st century) the dynamics of P. mugo krummholz in the central Apennines, with particular focus on the possible impacts of summer pastoralism, climate change, and neighborhood and other environmental predictors. Our research was based on spatiotemporally consistent krummholz and climate change measurements using the dendrochronology of P. mugo as an in situ proxy for climate changes. The following hypotheses were tested: 1. P. mugo krummholz expansion has been predominantly into abandoned pastures in the beech forest belt below today's treeline during the second half of the 20th century (H 0 ), or the krummholz extent below the treeline is stable (H 1 ). The null hypothesis is suggested by beech forest expansion into montane pastures in the central Apennines (van Gils et al 2008) and in the Alps (Dullinger et al 2004). 2. P. mugo krummholz has migrated upward in elevation over the second half of the 20th century (H 0 ) or is stable (H 1 ). Upslope migration is plausible given the upslope-migrating beech forest in the central Apennines (van Gils et al 2008). The alternative hypothesis is plausible in view of the stable treeline in the northern Apennines (Holtmeier 2009). 3. P. mugo krummholz expansion during the second half of the 20th century can (H 0 ), or cannot (H 1 ) be modeled using neighborhood and environmental predictors. The null hypothesis is plausible for neighborhood variables, but not for environmental predictors (Dullinger et al 2004;van Gils et al 2008). If a wellfitting model can be developed, a forecast will be attempted. 4. The high-elevation climate shows a warming trend (H 0 ) or is stable (H 1 ) over the 20th century. Global trends and meteorological records from the northern Apennines (Kumar et al 2005) support the null hypothesis. Findings from midlatitudinal high elevations (Liu et al 2009;Lu et al 2010) and from the central Apennines (Izzo et al 2004) suggest that spring and winters may be warmer, but mean annual and summer temperatures are stable (H 1 ). Any significant climate trends will be used for a forecast. 5. The radial stem growth of P. mugo shows that this species is affected by climate warming (H 0 ) or is stable (H 1 ) over the 20th century. The null hypothesis (H 0 ) is supported by findings on the same species from the French Alps and the Apennines (Rolland et al 1998;Palombo et al 2013), the alternative hypothesis (H 1 ) by dendrochronologies of Pinus cembra across the Alps (Carrer et al 2007).
By testing the hypotheses, the relationships among krummholz dynamics, climate change, environmental factors, and anthropic activities were analyzed. Results should be beneficial to policy-makers and local residents designing sustainable development strategies.

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 temperatemontane (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 2 in 3 major patches (Mt Rapina, Blockhaus, and Mt Ugni; Figure 1) at FIGURE 1 Location of Majella National Park, with P. mugo krummholz distribution (Blasi et al 1999) (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 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€ uller-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€ orner 2009).

Change detection
The color 2007 APs served as the reference for orthorectification 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).
An automated object-oriented classification computer program (eCognition) was applied to reduce the salt-and-  pepper effect (Hay et al 2005) resulting from using highresolution imagery (Pekkarinen 2002). The objectoriented 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) Dai) 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.
The distribution of continuous predictors was tested for symmetry, and those with skewness . 0.5 were logtransformed. 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)

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.
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 2 ). 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).
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.
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.878C century -1 ; R 2 ¼ 0.16) and winter (4.388C century -1 ; R 2 ¼ 0.25) minima have increased since 1919. Summer minimum temperatures increased (3.178C century -1 ; R 2 ¼ 0.16) and summer maximum temperatures decreased (- 2.698C century -1 ; R 2 ¼ 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 8C; Hoch and K€ orner 2009) until 2086 (Figure 4).
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 2 ¼ 0.25) and Mt. Ugni (R 2 ¼ 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. Shortened in text to ''proximity to krummholz.'' 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 debushing 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 topoclimate (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.48C 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.

Conclusions
The combination of change detection with sequential APs and dendrochronology provides a spatiotemporally consistent method to distinguish between drivers of krummholz changes above and below the beech forest treeline. The short distance migration in elevation of the P. mugo krummholz in the central Apennines is consistent with both climate warming and abandonment of pastoralism. The pattern of krummholz expansion is largely determined by proximity of P. mugo as seed source. In the krummholz belt, 20th-century climate change is expressed in higher winter and spring minimum temperatures. Although the dendrochronology of P. mugo is shown to be a suitable proxy for the occasional warmer springs or wetter summers, no climatic trend was embedded in radial stem growth of P. mugo. The declining maximum summer temperatures, stable summer rainfall, and mean annual temperatures at the research site conform to findings from high-elevation areas at midlatitudes in Eurasia but contrast with findings from the Alps and regional predictions.
The forecasts for both krummholz expansion and the local climate suggest that upslope migration of krummholz will be very modest for the foreseeable future. Only the inward expansion within the krummholz belt into abandoned pasture may be substantial. The combined spatial and temporal evidence suggests that the upslope migration of P. mugo krummholz since the mid-20th century can be explained by the legacy of summer pastoralism and the dispersal limitations of P. mugo. The regionally rare P. mugo krummholz is not threatened by extinction in the Majella massif; however, the endemicrich arctic-alpine vegetation may be compromised by the inward expansion of krummholz.