A new method for quantifying treeline-ecotone change based on multiple spatial pattern dimensions

Treeline-ecotone spatial patterns and their dynamics reflect underlying processes. Changes in ecotone pattern may reflect changes in natural drivers or land-use practices. However, characterizing these dynamics presents a major challenge, limiting our ability to map, understand and predict changes in the upper limits of mountain forests. This paper proposes a new method using multiple pattern dimensions to describe treeline-ecotone spatial pattern shifts. This standardized protocol should be able to (i) distinguish different types of treeline-ecotone patterns within a large study area, (ii) characterize temporal pattern shifts in spatial pattern between two or more dates. We mapped alpine treeline ecotones (ATE) at 648 sites in the eastern French Pyrenees using aerial images from ~ 1955 and ~ 2015, identifying forest and non-forest areas at the hillslope scale. Extracted patch metrics were summarized using a Principle Component Analysis (PCA) and spatial pattern change was quantified from the shift in the PCA space and compared to elevational shifts. Three clusters of patterns were distinguished: diffuse, discrete and island-forming ATEs. Between 1955 and 2015, about half of the sites changed from one pattern cluster to another. Shifts into discrete ATEs were associated with smaller and negative elevational shifts, while shifts into diffuse ATEs coincided with the highest positive elevational shifts. The proposed method allows a standardized and repeatable quantification of vegetation pattern change in alpine treeline ecotones based on historical aerial imagery. Seeing the importance of treeline-ecotone shifts for alpine biodiversity, we encourage the use of this protocol to better understand treeline dynamics at treelines globally.


Introduction
Vegetation transitions such as treeline ecotones (in mountains and in subarctic areas) are widely recognized as being sensitive to environmental changes in combination with anthropogenic changes (e.g.Camarero and Gutiérrez 2004;Holtmeier and Broll 2005).This sensitivity results in ecotonal dynamics that have many implications (e.g. for biodiversity conservation, ecosystem restoration, and gravitational hazard management) and that have been evidenced for a long time (see Holtmeier and Broll 2019 for an historical overview).Spatio-temporal ecotonal dynamics can manifest in different ways.In most studies, treelineecotone dynamics are characterized by elevational (or latitudinal) shifts (Ameztegui et al. 2016;Bonanomi et al. 2018;Elliott and Cowell 2015).However, not only the position but also the spatial pattern of treeline ecotones is subject to change.Bader et al. (2021) recently proposed a typology of treeline-ecotone spatial patterns, distinguishing, from a perpendicular top-down perspective, diffuse, discrete, and island patterns at the hillslope scale.They pointed out that these different patterns should reflect underlying processes (e.g.tree mortality, environmental heterogeneity, wind and snow effects or even seed dispersal) and thus can be used to understand treeline dynamics.
Yet treeline-ecotone spatial patterns or forms can change over time (e.g.densification or upward scattering of treeline ecotone; see Ameztegui et al. 2021) and recent studies have shown that such changes in patterns may not be systematically associated with an elevational shift (Feuillet et al. 2020;Morley et al. 2020;Treml and Chuman 2015).Whereas temperature is the major determinant of the position of the potential treeline, which will thus shift predictably with global warming (Körner and Paulsen 2004), the spatio-temporal heterogeneity observed in treeline ecotones implies that climate change alone is not sufficient to predict their dynamics.Clearly, other processes interact locally to modulate treeline-ecotone position and dynamics.These scale-dependent underlying processes can be revealed by treeline spatial patterns (Bader et al. 2021).As Harsch and Bader (2011) point out, these spatial patterns could be "a potential key to understanding treeline dynamics".Therefore, characterizing treeline-ecotone spatial patterns and their time-dependent changes is a prerequisite to disentangling the role of different environmental and anthropogenic factors in treeline-ecotone dynamics worldwide (see Holtmeier and Broll 2007 for an overview of the biotic and abiotic factors potentially involved).However, the quantitative operationalization of this characterization presents a major challenge, due to the multi-dimensionality and regional variability of the spatial patterns found in alpine treeline ecotones.
To date, no standard protocol for quantifying the temporal change in treeline-ecotone spatial pattern has been proposed.The consequent methodological heterogeneity makes results among studies hard to compare and may explain some of the inconsistencies reported worldwide.The different methods used include, for example, neighborhood analysis (Beloiu and Beierkuhnlein 2019;Humphries et al. 2008;Camarero et al. 2015), computation of forest-cover density shifts (Sigdel et al. 2020) or landscape metrics (Ameztegui et al. 2021;Malandra et al. 2019), andvisual interpretation (Treml andVeblen 2017).Landscape-and more specifically treeline ecotone spatial patterns-can be described by the configuration of cover patches, which can be characterized by aspects such as patch density, shape, complexity, and the degrees of clumping and connectivity.To highlight spatial pattern shifts in terms of patch configurations, these must therefore be characterized in a reproducible way, integrating all these aspects.
In this study we propose an easily replicable methodology for characterizing treeline-ecotone pattern shifts based on patch-based techniques for quantifying landscape characteristics (Hesselbarth et al. 2019;McGarigal et al. 2009;Kupfer 2012).We set out to quantify treeline ecotone Vol.: (0123456789) spatial pattern changes at the hillslope scale, replicating sites at a regional scale, integrating the multiple dimensions of spatial patterns.To achieve this objective, we characterized treeline spatial patterns and changes in these patterns within a multidimensional space.We aimed to (i) characterize treeline ecotone patterns within a large study area, (ii) characterize temporal pattern shifts, and (iii) relate pattern shifts to elevational shifts.

Data and methods
Treeline-shaping ecological processes and their hierarchy could evolve over time, resulting in changes in spatial patterns or the pathways of ongoing pattern development.The protocol developed here aims to evaluate, in a standardized way, changes in alpine treeline ecotone (ATE) spatial pattern between two dates based on repeat aerial photography.The protocol is easily replicable, the workflow being implemented entirely using R software and open source packages (R core Team, 2022;Hesselbarth et al. 2019).Figure 1 graphically sums up the proposed methodological protocol described below.

Study area
Analyses were applied to the ATE of the French eastern Pyrenees (Fig. 2), covering the Pyrénées-Orientales and Ariège departments.The study area includes some major massifs: the Canigou (2,785 m a.s.l.), Puigmal (2,910 m a.s.l.), Carlit (2,921 m a.s.l.) and Madrès (2,469 m a.s.l.) in the Pyrénées-Orientales department, and the Trois-Seigneurs (2,199 m a.s.l.), Montcalm (3,077 m a.s.l.), and Mont Valier (2,838 m a.s.l.) in the Ariège department.It includes various climatic, topo-geomorphologic and anthropogenic conditions, with for example: (i) a climate gradient with Mediterranean influence increasing towards the Pyrénées-Orientales, with a concurrent increase in summer droughts and decrease in oceanity; (ii) lithological variation with crystalline (e.g.Canigou and Carlit) and sedimentary (e.g.Puigmal) those from ~ 1955.Supplementary individuals (supp.ind.purple dots) refers to patterns for which coordinates are predicted based on the PCA constructed from active individuals, but who did not contribute to the PCA building Vol:.( 1234567890) massifs; and (iii) a more marked abandonment of both the agro-silvo-patoral system and the extraction of metals in the Pyrénées-Orientales from the midnineteenth century (Milian et al. 2012;Saulnier et al. 2020;Métailié and Paegelow 2004;Vacquié 2015).Mountain pine (Pinus uncinata Ramond ex.DC.) is the most prevalent species at the forest upper limit, along with various deciduous tree species like Fagus sylvatica, mainly in the western part of the study area (Ariège).Pinus uncinata, which presents an ecological optimum in the upper-montane forest, is the main tree species found at Pyrenean ATE (Feuillet et al. 2020;Améztegui et al. 2010;Cantegrel 2019;Cantegrel 1986).Since the alpine vegetation is generally dominated by herbaceous species, the spectral contrast with the dark-colored P. uncinata trees is strong and classification of panchromatic and visible-color aerial imagery was feasible in our study area.

Data acquisition
Because tree growth and colonization processes tend to be slow in treeline ecotones, we aimed at covering the longest possible time span available in aerial photographs for the study area and were lucky to have high-quality images available already from the 1950s.Spatial patterns at the upper forest limit in the study area were therefore delimited from panchromatic historical (1956 and 1953 respectively for Ariège and Pyrénées-Orientales) and recent colored RGB (2016 and 2015 respectively for Ariège and Pyrénées-Orientales) 0.5-m resolution orthophotographs.These data were provided by the National Institute of Geographic and Forest Information -IGN (https:// www.ign.fr/ insti tut/ ident itycard).Radiometric pre-processing and ortho-rectification with digital terrain models were carried out by the data provider (IGN; https:// geose rvices.ign.fr/ docum entat ion/ donne es/ ortho/ bdort ho).We assessed the spatial fit of the two image dates by measuring the distance between 144 distinctive landscape features (e.g.rocks) on the two overlain images.The average shift was 2.3 m.

Image classification
For each site, orthophotographs (with a 0.5-m digital resolution) were first segmented using the segment mean shift algorithm (ArcMap v.10.5.1) to group contiguous close mean value pixels in the same segment.A supervised classification based on the random forest algorithm was then applied to these segmented images using the EnMAP-Box toolbox (Breiman 2001;Van der Linden et al. 2015;EnMAP-Box Developers, 2019) based on the scikit-learn Python library (Pedregosa et al. 2011).Accuracy values were computed based on confusion matrices from a sample of 500 random points for each year and each department (Table 1).Historic and recent forest/nonforest binary classifications were obtained for each site set (Fig. 1).
These classifications were visually checked and manually corrected.Finally, the classifications were sifted to keep only the forest patches formed of more than 100 pixels (25 m²).This threshold was determined visually to exclude the smallest forest patches and small isolated trees as we work at the resolution of small clusters of trees and to focus only on the lower part of the ATE (near the timberline).Excluding these small tree-covered patches has important implications for the interpretation, as discussed below ("Discussion" section).This threshold was also determined to exclude artifacts, such as boulders classified as trees, which were sometimes unavoidable in the classification of the panchromatic images.

Sampling design
Sets of two (one for each year) rectangular plots (300 m x 200 m, following the study of Dearborn and Danby 2020) were arranged along the upper forest limit parallel to the main slope centered on the upper forest limits visually identify as the transition zone between forest and alpine grassland of the respective year and spaced about 500 m apart on the classified orthophotographs.Thus each plot incorporates approximately 50% forest and 50% grassland so as to be centered on the ecotone.Each plot was manually arranged by the same operator to control for sampling bias.A total of 648 plots per year was delineated, 457 in Ariège and 191 in Pyrénées-Orientales (Fig. 2).Each site was therefore described by two (recent and historical) plots, more or less overlapping depending on the elevational change of the treeline ecotone between the two dates (Fig. 2).

Spatial pattern characterization by landscape metrics
The ATE spatial pattern was characterized within each rectangular plot from the computation of landscape metrics at the class level (forest / non-forest) using the "landscapemetrics" package (Hesselbarth et al. 2019) in the R software (version 4.0.5;R Core Team, 2022).We used 19 landscape metrics (Table 2) in four categories (aggregation, area and edge, core area, and spatial pattern) to describe the spatial patterns.Apart from the mean of the metric for each plot, we also calculated dispersion parameters (variation coefficient and standard deviation) for those metrics that were calculated at the patch  (Hesselbarth et al. 2019;McGarigal 2015).For those metrics that are computed as a mean of the values per patch, we also used the variation between patches as a separate metric in the PCA analysis, considering both the coefficient of variation (cv) and the standard deviation (sd).Statistical computations were performed only for the forest class, except from the contiguity index, which was calculated taking into account forest and non-forest pixels   Vol:.( 1234567890) level, to quantify heterogeneity within each plot.This resulted in a total of 33 (19 + 14) metrics that were included in the principal component analysis (PCA; see below).

Information reduction by PCA
A reduction of the information by principal component analysis (PCA) was carried out based on the metrics calculated for the historical ecotone plots, while adding recent ecotones as supplementary individuals (note that in factor analysis, supplementary individuals do not contribute to the component building but do have factorial coordinates and can then be positioned in the factorial spaces).This first step thus allowed to synthesize with a reduced number of variables (the so-called principal components) the information contained in the initial 33 variables.This analysis served as a basis for clustering plots in groups with similar characteristics, for quantifying profiles of pattern change, as described below, and for identifying the dominant dimensions distinguishing spatial patterns in our study area.The number of PCA dimensions for subsequent analysis was chosen based on three main constraints: PCA eigenvalues; explained variance; and multifactorial analysis interpretability.The PCA was performed using the FactoMineR R package (Le et al. 2008) and visualized using the Factoextra R package (Kassambara and Mundt 2020).

Clustering and qualitative assessment of spatial pattern change
A hierarchical classification on principal components (HCPC) of historical ecotones was performed, also using the FactoMineR R package (Le et al. 2008).Since a HCPC based on the first four principal components showed that the first two components provided an almost perfect distinction of the clusters, we then used only the two first principal components to build clusters, facilitating visualization and the comparison between the years.The number of clusters was derived from the tree partition that maximized the relative loss of within-cluster inertia.The factorial plane was then divided into subspaces corresponding to the cluster map based on a Voronoi tessellation with an aggregation of polygons per cluster (Fig. 4).Recent ecotones were then projected in these cluster-based subspaces in order to count the shifts of ecotones from one cluster to another.The "v.test" statistics (Escofier and Pagès 2008;Feuillet et al. 2012) was used to describe what patch metrics differ significantly between the clusters and thus define them.The resulting clusters were associated with the pattern typology described by Bader et al. (2021) and presented graphically by the most representative ecotone plots, i.e. those located closest to the centroid of the cluster.Shifts between clusters, i.e. from one treelineecotone type to another, between the two dates were recorded as a qualitative measure of pattern change.

Quantitative assessment of spatial pattern change
Changes in pattern between the two dates were quantified through a spatio-temporal PCA (Cossart et al. 2020).The patch metrics of the recent plots were projected onto the factorial plane defined by the historical metrics, and coordinates on the first and second axes were then extracted for each site both for the historical and the recent plot.The distance between these pairs of plots was calculated from these coordinates.This distance then defined the degree of change in the ecotone spatial pattern within each site: the greater the distance, the more the ecotone pattern had shifted.Additionally, we calculated a directional pattern shift by extracting the absolute shift along the first PCA axis, along which the three pattern types were most strongly distinguished.Here, positive distances indicated a shift towards more diffuse and negative distances a shift towards more discrete treelineecotone types.The differences in elevational shift between groups of sites with qualitatively different spatial pattern shifts was determined by performing Kruskall-Wallis tests.Finally, a pairwise comparison Tukey's test was performed to identify significant differences between the treeline-ecotone plot clusters.

Descriptive statistics
The two sets of 648 plots were located at an average of 1850 and 1885 m a.s.l. for plots centered on historic and recent treelines, respectively.The elevational change in ATE position, as defined by the five uppermost forest pixels, was 39.4 m on average across the entire study area with a range from − 212.0 m to + 307.8 m.

Information reduction by PCA
The first four PCA axes had eigenvalues of 14.2, 6.0, 3.2 and 2.0.The first and second axes accounted for 65.2% of the total variance with respectively 45.8% and 19.4% (Fig. 3), while the third and fourth accounted for respectively 10.4% and 6.6%.Due to the strong reduction in information between the first two and the third and fourth dimensions, we selected only the first two dimensions for subsequent analyses.
The first PCA axis reflects 13 area and edge, shape, core area and aggregation metrics, which contribute more than 70% to its construction.Among these metrics, 11 are strongly positively correlated with it (r > 0.75; p < 0.05) and two are negatively correlated with it (r < − 0.75; p < 0.05; Table S1).Thus, an increase in this first dimension is associated with an increase in total core area ("tca"), core area ("ca"), percentage of landscape ("pland"), largest patch area ("lpi"), variability of radius of gyration ("gyrate_sd" and "gyrate_cv"), core (core_sd), total area (area_sd), aggregation, cohesion and clumpiness ("clumpy"), and a decrease of division index ("division") and normalized landscape shape index ("nlsi").This dimension thus describes the aggregation level and the size of the patches forming the forest landscape with positive values indicating ecotones with few large and aggregated patches, and negative values indicating many dispersed small patches.
For the second component, the six most correlated metrics (r > |0.75|; p < 0.05; Table S1) contribute more than 65% to its formation.An increase in this second dimension is associated with an increase of contiguity and core area indices ("contig_cv", "contig_sd", "cai_cv" and "cai_sd") and a decrease in the average of these indexes ("contig_mn" and "cai "mn").This dimension thus describes the level of connectedness within the ecotone.Positive values in this dimension are associated with ecotones formed by heterogeneous connections between small patches and negative values with well-connected and large patches.

Clustering of historical and recent treeline ecotone plots
Three clusters were distinguished, and these separated out most strongly along the first PCA axis.The first cluster is located to the left, the second cluster in the middle towards the top (high values on the second  2 for abbreviations of the metrics Vol.: (0123456789) axis), and the third one to the bottom right (Fig. 4, Table S2).

Characterization of the clusters
The majority of the patch metrics differ between the three clusters (Fig. 5).The three historical (~ 1955) clusters made it possible to distinguish three types of forest treeline-ecotone plots: diffuse, island and discrete (Fig. 6).These were characterized by an aggregation and increasing patch size, and a decrease in the number of patches from diffuse to discrete, through island.The transition between the forest below and non-forest vegetation above became increasingly distinct from cluster 1 to cluster 3, as illustrated in Fig. 6.These terms are based on the terminology suggested in Bader et al. ( 2021), although we apply the term "diffuse" more broadly, including patterns with small patches rather than just those with spreading single trees.

Profiles of change
The distance calculated between recent and historical plot coordinates projected in the historical PCA (Fig. 4c) constitutes the change in the spatial pattern of the treeline-ecotone site at the two dates.Thus, the greater the distance between two points of a plot pair, the greater the change in spatial pattern has been.
Of the 648 ecotone sites, 52% (n = 337) experienced changes in spatial patterns between ~ 1955 and ~ 2015: changing clusters from diffuse to island or discrete (n = 163), from island to diffuse or discrete (n = 125), or from discrete to diffuse or island (n = 49).The average distance between the treelineecotone plots of 1955 and 2015 in the factorial plane is 5.0, with a range from 0.3 to 16.8.
Overall, there was a shift towards more discrete patterns, with the number of diffuse treelines decreasing and the number of island and discrete treelines increasing over time (with 333, 220 and 95 plots in 1955, and 240, 262 and 146 plots in 2015, respectively for cluster 1, 2 and 3; Table S2).More than half of the 333 plots with a diffuse spatial pattern in ~ 1955 remained in this cluster in 2015 (51%, n = 170), while 40% (n = 134) evolved into an island form and 9% (n = 29) evolved towards a discrete form (Fig. 7).Of the 220 plots with an island spatial pattern, 43% (n = 95) remained in this cluster in 2015, 25% (n = 54) shifted to a diffuse form and 32% (n = 71) to a discrete form.Of the 95 plots with a discrete spatial pattern, 48% remained in this cluster in 2015 (n = 46), 35% shifted to an island form (n = 33) and 17% to a diffuse form (n = 16).As expected, the shifts in the factorial plane were smallest for the sites not changing cluster and largest for the sites changing between cluster 1 and 3 or vice versa (Fig. 7).

Comparison between pattern shifts and elevational shift
The Kendall correlation between elevational shifts and extent of the spatial pattern changes was not significant (p = 0.9) and the LOESS regression indicated the absence of any relationship between elevational shifts and the distance of pattern shift in the PCA plane (Fig. 8a).This result indicates that elevational shifts can take place without large shifts in the pattern and vice versa.However, these two types of shifts are not independent, which becomes clear when considering the quality of the pattern shift rather than just the distance.Shifts into discrete ATE types tended to coincide with small or downward elevational shifts, while changes into diffuse ATE types coincided with larger upward elevational shifts (Fig. 8b).Directional change along the first PCA axis showed a negative correlation with elevational shift (Pearson r = − 0.23, p < 0.05).This confirms that ATE becoming more diffuse had a larger upward shift while those becoming more discrete moved the least or downward in elevation.The observed elevational and pattern change differed between the clusters, i.e. between ecotone spatial types.Historically as well as, more clearly, recently diffuse patterns changed the least in terms of the pattern metrics (Fig. 9a and b), although overall their elevational change was similar to that of the other forms (Fig. 9c and d).They were situated at the highest elevation compared to the other spatial patterns (Fig. 9e and f).Those ecotones that were already diffuse in 1955 and remained diffuse until 2015 were located at the highest elevations of all, while most ecotones that transitioned from being diffuse to discrete were situated at lower elevations, similar to ecotones that were discrete already in 1955 (Fig. 9f).

Discussion
The new method used demonstrated that the studied alpine-treeline ecotones showed variable and partly strong shifts in their spatial pattern, and that these shifts were not necessarily accompanied by a strong elevational shift.Likewise, strong elevational shifts (up to 200 m upwards or 100 m downwards) not necessarily implied strong changes in the spatial pattern, as also suggested by Feuillet et al. (2020) based on forest cover density as a proxy of spatial pattern.So the intensity of the change was not correlated, but the directions of the shifts were clearly related: shifts towards discrete ecotones were associated with small or downwards elevational shifts while shifts towards diffuse patterns were associated with the largest upwards shifts.These associations align with theoretical predictions, since discrete treelines, especially when also abrupt (i.e. an abrupt change in vegetation height from tallish trees to low-stature alpine vegetation; in contrast to discreteness, which refers only to the change in forest cover, as seen from above) are thought to be caused by disturbances (i.e.mortality above but not inside the ecotone1 ), while diffuse treelines may indicate an active colonization by trees of the alpine zone (Bader et al. 2021).
Since treeline elevational shift and pattern shift are not necessarily correlated (as shown here and previously in Feuillet et al. 2020), but pattern shifts can have important consequences for ecosystem functioning in the treeline ecotone, quantifying directional pattern change is important.As shown here, a reduction of the information held in different pattern metrics via PCA can be a useful way of identifying the main dimensions of spatial pattern and defining the direction of pattern change.In contrast to most previous studies on treeline-ecotone change, which have focused mainly on infilling processes, quantified as increases in tree cover or density within the ecotone (e.g.Treml and Chuman 2015; Batllori et al. 2010), our method explicitly uses landscape metrics that describe ecologically meaningful pattern properties like clumping, connectivity and core-edge ratios.On the one hand, these patterns describe the landscape available to e.g.animals and establishing plants.On the other hand, just like point patterns based on elaborate field sampling of tree locations (e.g.Batllori et al. 2010;Wang et al. 2016;Camarero et al. 2000), they indicate which processes are acting to form and perhaps shift the treeline ecotone (Bader et al. 2021).
Since the PCA is sensitive to the metrics included, with both the variability within metrics and the number of correlated metrics describing certain pattern dimensions (e.g.connectivity) determining the principal components, the array of metrics included needs to be determined with care.For specific questions about ecotone pattern change, metrics could be considered individually or in tailored combinations (e.g. if connectivity of either forest or alpine habitats for animal movement within the ecotone is of interest).Likewise, such specific studies could determine cluster boundaries and define clusters (ecotone types) based on specific criteria taking into account functional properties of these metrics.However, for more general and explorative studies like the one presented here, the use of the full array of pattern dimensions, reduced to manageable dimensions allowing the interpretation of shift directions by PCA, holds the advantage of capturing a maximum of information about the patterns, reducing the risk of missing subtle changes.Likewise, the statistical approach to defining clusters used avoids subjective decisions on cluster boundaries and is suitable to describe change trends.Until the scientific community agrees on the quantitative criteria for defining treeline-ecotone spatial types, such a statistical approach is the best we can do.
Since the contrast between the P. uncinata canopies and the soil and vegetation background in our study area was strong, it was possible to identify trees even on panchromatic aerial imagery from the 1950s.However, even under these favorable conditions, and even on the more modern multispectral images, an automated identification of single trees, especially smaller trees as found towards the treeline, is challenging at best on imagery with a 0.5-m resolution.To avoid spurious results, we decided to exclude the smallest segments, which may represent single trees but also rocks and their shadows, especially on the panchromatic images.As a result, the diffuseness and possibly the elevational advance of the studied treeline ecotones were probably underestimated.This does not invalidate the method, however, which can also be applied in cases where single trees can be identified unequivocally thanks to very-high resolution spectral imagery and/or the additional use of Lidar data to allow a clearer demarcation of tree cover.If sizes of individual trees can be derived, e.g. from Lidar data (e.g.Coops et al. 2013) or by correlating horizontal canopy sizes or shade lengths (Dial et al. 2022) to tree height, we recommend including additional metrics describing the treeline ecotone in three dimensions, such as the height of the uppermost trees, the rate of height change with elevation, and the distribution of height-to-width ratios of the trees (to detect Krummholz forms).In most cases, such metrics will not be available for historical ecotones, and the spatial resolution may not allow the same level of detail as for recent images.Since for detecting change the same resolution should be used for historical and recent data, the historical data set the standard.However, additionally having access to more detailed metrics for recent ecotones is useful both for future change detection and for better understanding the processes that have led to the recent pattern.For example, whether advancing trees are stunted to Krummholz tells us whether they suffer from dieback, e.g.due to strong winds or frost desiccation, which would limit further advance.And whether a stable discrete treeline is composed of tall trees or not tells us whether these trees are growth limited or not -if they Vol.: (0123456789) are not, disturbance is likely maintaining the discrete form.
Our treeline ecotones were relatively simple to work with, not only because the trees were easily distinguished, but also because they are composed on a single dominant tree species among herbs and small shrubs.This contrasts strongly with e.g.tropical treeline ecotones, where species richness of trees and tall shrubs can be very high and although many ecotones are very abrupt and discrete, those that are not tend to show no clear binary distinction between forest and alpine vegetation due to the presence of tall shrub vegetation (e.g.Bader et al. 2007;Hofstede et al. 2014).This does not preclude the use of the proposed methodology though.More complex vegetation patterns (i.e. with more than two classes) may also be quantified by landscape metrics.Although this remains to be tested, we expect that these can then be used to distinguish pattern types and to detect pattern shifts in the same way as described here.
We developed the method with this case study based on images from two dates 60 years apart.Shorter and multiple time intervals would have allowed a finer matching of changes to potential drivers, especially for changes due to disturbances.However, we consider that generally for alpine treelines, where processes tend to be slow (disturbance is fast but the resulting patterns will last until dispersal and growth have permitted a recolonization by trees), this long time period is suitable to describe change, if only two dates area going to be compared.On the other hand, the method could also be applied to other spatial objects for which spatial pattern shifts may be much faster, such as urban sprawl patterns, urban ecosystem services (Schwarz 2010; Grafius et al. 2018), the effects of natural disaster such as typhoons on landscapes (Lin et al. 2006), or ecological connectivity (Herrera et al. 2018).In such cases data points should clearly be chosen at shorter time intervals.Analyzing more than two points in time can be done, as long as data are available and can be processed, using the same methodology and is useful to refine the description and interpretation of spatial pattern shifts.
The fact that diffuse treelines occurred at higher elevations than discrete ones and that ecotones changing from diffuse to abrupt ended up located at much lower elevations than those remaining diffuse, suggests that the discrete ecotones were, in the majority of cases, formed by the removal of forest patches, likely by disturbances.In addition, (re)colonization of the area above the closed forest is either not yet completed (trees being smaller than our lower size threshold) or hindered completely by repeated disturbances, either natural or anthropogenic, or limitations to seedling establishment imposed by the substrate (e.g. a lack of soil).For some tree species, shade-dependence or sensitivity to frost at the seedlings stage may also limit tree establishment outside the protection of forest (this mechanism probably explains the maintenance of discrete abrupt treelines in southernhemisphere treelines composed by Nothofagaceae, Wardle 2008).Although such a mechanism may play a role for the beech (Fagus sylvatica) treelines in the western part of our study area, it cannot explain the variety of patterns observed in the eastern part of the study area, where one single tree species, Pinus uncinata, composed all the ecotones.Although our results cannot exclude the possibility that the more diffuse forms are also shaped by disturbances, especially land use practices like cutting and grazing, it appears that in the Eastern Pyrenees, the more diffuse treelines represent the least disturbed.
Interestingly, the spatial pattern appeared to have only a weak predictive value for future elevational or pattern shifts, while it reflected past pattern shifts but not elevational shifts within the last 60 years.This may appear to contradict theoretical expectations, where diffuse treelines are predicted to be more dynamic and more responsive to climate change than discrete treelines (Bader et al. 2021).However, this prediction holds only if the processes leading to certain treeline forms are not changed.If patterns are caused by disturbance (e.g.due to human land-use), as is often (though not always) the case with abrupt discrete treelines, and this disturbance stops, then the dynamics will change (in this case, trees will likely colonize the space above the closed forest) and the pattern loses predictive power.Also, in treeline ecotones that are all composed of the same tree species (like in the eastern part of our study area), the spatial patterns are less likely to reflect fundamentally different ecological processes between sites, but rather reflect disturbance history.Disturbances are predictable only to a small extent.Lastly, excluding the smaller single trees from the classification and pattern descriptions prevented the detection of recent colonization patterns, important indicators of ongoing Vol:.( 1234567890) shifts in pattern and elevation.Thereby "recent" can sometimes comprise even the 60 years covered by this study, since trees tend to grow very slowly when approaching treeline.As mentioned above, a higher resolution of aerial imagery and auxiliary spatial data could at least partially lift this limitation.
As demonstrated, our protocol can be applied to characterize spatial pattern changes at the regional scale.However, it presents a limitation inherent to geomatics, for our analysis allows us to characterize treeline-ecotone pattern change only in a two-dimensional (2D) way, as seen from above.Thus, some treeline-ecotone characteristics remain invisible from this viewpoint and at this scale (e.g.tree dimensions, species composition, age and growth form), whereas a field survey could reveal them (Treml and Veblen 2017;Camarero et al. 2000Camarero et al. , 2004)).Conversely, a field survey does not allow as many sites to be compared as a remote-sensing-based analysis does.In order to articulate the advantages of both approaches, we advise first carrying out a geomatic analysis and then a field survey, for a detailed study of a subset of sites.This dual work (geomatic and field) thus allows one to target sites that are often difficult to access, and to carry out a representative sampling of the diversity of spatial pattern change.

Conclusion
The study of treeline-ecotone spatial patterns is key to understanding the underlying mechanisms at work, in order to disentangle the multiple interactions associated with its dynamics.A first step consists in characterizing the spatial patterns and their temporal change.However, as spatial pattern change is understudied compared to elevational shift, no standardized methodology has been proposed so far.The methodology presented here takes into account the multidimensional nature of spatial patterns.Our change analysis emphasizes the variety in the evolution of treeline-ecotone spatial patterns.This suggests a variety of associated underlying processes, and thus the importance of capturing spatial pattern temporal shifts when analyzing treeline-ecotone dynamics.Finally, to interpret the observed changes in spatial pattern, we recommend exploring spatial associations between these changes and potential drivers (e.g.climate and anthropogenic land use change, topogeomorphology).Combined with such an analysis, the methodology presented here allows another step forward in understanding what drives the dynamics of alpine treeline ecotones and how mountain ecosystems will respond to global change in all its multiplicity.

Fig. 1
Fig. 1 Graphic summary of the proposed protocol to quantify alpine-treeline ecotone spatial pattern change, applied in the eastern part of the French Pyrenees with aerial images from ~ 1955 and ~ 2015.Active individuals (active ind., green dots) refers to patterns (in plots) used to construct the PCA, i.e.

Fig. 2
Fig. 2 Study area location map and sampling design.a Distribution of the sites (n = 648) in the Eastern Pyrenees.Shades of blue and green indicate 1-km elevational bands from lowland (blue) to > 3000 m (yellowish green), the dark blue to the east is the Mediterranean Sea.b Location and spatial configuration of the example sites shown in c and d. c Historical aerial image showing two plots with the binary forest/non-forest clas- percentage of each patch in core area, i.e. pixels that have only neighbors in the same class.values in 3 × 3-pixel windows as a fraction of the maximum value, i.e. for a landscape of one class.Central pixel counts as 1 and sets the class, horizontal and vertical neighbours of the same class count as 2 and diagonal same-class pixels count as 1. Background pixels are count as coefficients) were estimated to evaluate the association between quantifications of the spatial patterns change (distance in factorial plane) and the elevational treeline shift.Finally, LOESS regression was estimated for testing the ability to predict quantitative spatial pattern change based on the elevational shift indicator.

Fig. 3
Fig. 3 Correlation of the landscape metrics describing 648 alpine treelineecotone sites in the French Pyrenees with the first two PCA axes based on the historical (~1955) landscape metrics.The length and color of the arrows indicate their contribution in percent to the construction of the two main dimensions (see color scale).See table 2 for abbreviations

Fig. 4 a
Fig. 4 a Clusters of treeline-ecotone plots based on a hierarchical classification in a PCA-space based on patch metrics from historical (~1955) aerial imagery.b Subdivision of the factorial plane into subspaces based on a Voronoi tessellation and assignment of each polygon to the cluster map.c Recent

Fig. 6
Fig. 6 Treeline-ecotone plots closest to the center of each cluster, as defined on the first two PCA axes based on landscape metrics of forest cover, for historical (~ 1955) treeline ecotones in the French Pyrenees

Fig. 7
Fig. 7 Number of Pyrenean treeline-ecotone sites (in color) changing from one ecotone type to another between ~1955 and ~2015 (total n = 648 sites), and the average distance moved in the PCA plane by this groups, shown as numbers (unitless)

Fig. 9
Fig. 9 Violin plots comparing different treeline-ecotone spatial types according to a clustering of landscape metrics, either historical ~1955 (a and c) or recent ~2015 (b and d) in terms of spatial pattern shifts (a and b), elevational shifts (c and d), or recent elevation (e).f Violin plots comparing the recent ele-

Table 1
Accuracy assessment (n = 500) of each classification performed on historical(~ 1955)and recent (~ 2015) orthophotographs for each French departement (Ariège and Pyrénées-Orientales) before visual control and manual correc-tion.P_accuracy = Producer accuracy (refers to false negative segments); U_accuracy = User accuracy (refers to false positive segments); Kappa = Overall agreement between classified pixels and ground truth. 1 = forest area, 0 = non-forest area

Table 2
Description of the landscape metrics computed from algorithms provided by the "landscapemetrics" package