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

Déborah Birre, Thierry Feuillet, Roman Lagalis, Johan Milian, Frédéric Alexandre, David Sheeren, Roberto Serrano-Notivoli, Matthieu Vignal, Maaike Y Bader

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Déborah Birre, Thierry Feuillet, Roman Lagalis, Johan Milian, Frédéric Alexandre, et al.. A new method for quantifying treeline-ecotone change based on multiple spatial pattern dimensions. Landscape Ecology, 2023, 38, pp.18. 10.1007/s10980-022-01589-4. hal-03942369

HAL Id: hal-03942369

https://hal.inrae.fr/hal-03942369

Submitted on 17 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A new method for quantifying treeline-ecotone change based on multiple spatial pattern dimensions

Déborah Birre · Thierry Feuillet · Roman Lagalis · Johan Milian · Frédéric Alexandre · David Sheeren · Roberto Serrano-Notivoli · Matthieu Vignal · Maaike Y. Bader

Received: 5 October 2022 / Accepted: 23 December 2022
© The Author(s), under exclusive licence to Springer Nature B.V. 2023

## Abstract

*Context* 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.

*Objective* 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.

*Method* 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.

*Results* 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

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10980-022-01589-4.

D. Birre (✉) · F. Alexandre
Université Sorbonne Paris Nord, Pléiade, UR 7338, F-93430 Villetaneuse, France
e-mail: deborah.birre@univ-paris13.fr

T. Feuillet
University of Caen, UMR 6266 IDEES CNRS, Caen, France

R. Lagalis
University of Paris Cité, Paris, France

J. Milian
University of Paris 8, UMR LADYSS CNRS, Saint-Denis, France

D. Sheeren
University of Toulouse, INRAE, UMR DYNAFOR, Castanet-Tolosan, France

R. Serrano-Notivoli
University of Zaragoza, Zaragoza, Spain

M. Vignal
University of Avignon, UMR 7300 ESPACE CNRS, 74 Rue Louis Pasteur, 84000 Avignon, France

M. Y. Bader
Faculty of Geography, Ecological Plant Geography, University of Marburg, Marburg, Germany
ATEs coincided with the highest positive elevational shifts.

*Conclusion* 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.

**Keywords** Alpine treeline ecotone · Landscape metrics · Structural shift · Spatio-temporal PCA · French Pyrenees

## 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, treeline-ecotone 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 tree-line 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; Mory 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 tree-line 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), and visual interpretation (Treml and Veblen 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
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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

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

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 inogeneity; (ii) lithological variation with crystalline (e.g. Canigou and Carlit) and sedimentary (e.g. Puigmal)
<figure><img src="image_4.png" /><figcaption>**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 classification (~1955). **d** The same sites as in c on recent imagery showing the respective plots and binary classification (~2015). Image source: BD-ORTHO/BD-ORTHO Histo IGN-f (0.5-m resolution). Source for a: From "Contour des departements français issus d'OpenStreetMap", © the contributors of Open-StreetMap under license from ODbL (02/01/2018) and ©Euro-Geographics for the administrative boundaries</figcaption></figure>

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 mid-nineteenth 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/institut/identity-card). Radiometric pre-processing and ortho-recognition with digital terrain models were carried out by the data provider (IGN; https://geoservices.ign.fr/documentation/donnees/ortho/bdortho). We assessed the spatial fit of the two image dates by measuring the distance between 144 distinctive
<caption>**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 correction.</caption>

Pendorf accuracy (refers to false negative segments); Uendorf accuracy (refers to false positive segments); Kappa = Overall agreement between classified pixels and ground truth. 1 = forest area, 0 = non-forest area

<table><thead><tr><th rowspan="2">Year</th><th rowspan="2">Department</th><th colspan="2">Pendorf accuracy</th><th colspan="2">Uendorf accuracy</th><th rowspan="2">Kappa</th></tr><tr><th>1</th><th>0</th><th>1</th><th>0</th></tr></thead><tbody><tr><td rowspan="2">~1955</td><td>Ariège</td><td>0.88</td><td>0.97</td><td>0.92</td><td>0.96</td><td>0.87</td></tr><tr><td>PO</td><td>0.72</td><td>0.99</td><td>0.86</td><td>0.97</td><td>0.76</td></tr><tr><td rowspan="2">~2015</td><td>Ariège</td><td>0.85</td><td>0.93</td><td>0.91</td><td>0.88</td><td>0.79</td></tr><tr><td>PO</td><td>0.83</td><td>0.98</td><td>0.94</td><td>0.93</td><td>0.83</td></tr></tbody></table>

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/non-forest 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 identified 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
<figure><img src="image_5.png" /><figcaption>**Table 2** Description of the landscape metrics computed from algorithms provided by the “landscapemetrics” package (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</figcaption></figure>

<table><thead><tr><th>Name in package</th><th>Abbreviation</th><th>Description</th><th>Values</th><th>Units</th><th>%</th><th>Category</th></tr></thead><tbody><tr><td>Aggregation index</td><td>Ai</td><td>Percent of same-class neighbor pairs. 100 : same-class neighbors pixels maximally aggregated.</td><td>[0-100]</td><td></td><td></td><td>Aggregation</td></tr><tr><td>Clumpiness index</td><td>Clumpy</td><td>Deviation from a random distribution of same-class pixels.<br/>- 1: same-class pixels disaggregated.<br/>0: same-class pixels randomly arranged.<br/>1: same-class pixels aggregated.</td><td>[-1-1]</td><td>/</td><td></td><td>Aggregation</td></tr><tr><td>Patch cohesion index</td><td>cohesion</td><td>Patch perimeter along same-class neighbors as percentage of maximum patch perimeter.<br/>0: same-class pixels disaggregated and subdivided, 100: all same-class pixels connected.</td><td>[0-100]</td><td>/</td><td></td><td>Aggregation</td></tr><tr><td>Division index</td><td>Division</td><td>Probability that two pixels randomly selected are not in the same patch.<br/>0: One big patch<br/>1: Only 1-pixel patches</td><td>[0-1]</td><td>Probability</td><td></td><td>Aggregation</td></tr><tr><td>Normalized landscape shape index</td><td>Nlsi</td><td>Total patch-edge length of same-class pixels rescaled to the minimum and maximum same-class patch length possible.<br/>0: landscape is constituted by a single square patch<br/>1: landscape is maximally disaggregated</td><td>[0-1]</td><td>/</td><td></td><td>Aggregation</td></tr><tr><td>Patch density</td><td>pd</td><td>Number of patches per area</td><td>>0</td><td># patches / 100 ha</td><td>Aggregation</td><td>Aggregation</td></tr><tr><td>Splitting index</td><td>split</td><td>Number of patches obtained by dividing the patches of the same class into patches of the same size<br/>1: One patch<br/>Increases with same-class patch area reduction and landscape division.</td><td>[1; total landscape pixels]</td><td>/</td><td></td><td>Aggregation</td></tr><tr><td>Patch area</td><td>area</td><td>Mean same-class patch area.</td><td>>0</td><td>ha</td><td></td><td>area and edge</td></tr><tr><td>Total class area</td><td>ca.</td><td>Sum of the area of all patches of the same class.</td><td>>0</td><td>Ha</td><td></td><td>area and edge</td></tr></tbody></table>
<table><caption>**Table 2** (continued)</caption><thead><tr><th>Name in package</th><th>Abbreviation</th><th>Description</th><th>Values</th><th>Units</th><th>Category</th></tr></thead><tbody><tr><td>Edge density</td><td>ed</td><td>Sum of the patch edge length (in meter) of the same class over the total class area (in hectares).</td><td>≥0</td><td>m/ha</td><td>area and edge</td></tr><tr><td>Radius of gyration</td><td>gyrate</td><td>Mean distance between each pixel and the centroid of its patch.</td><td>≥0</td><td>m</td><td>area and edge</td></tr><tr><td>Largest patch index</td><td>lpi</td><td>Percent of landscape included in the largest patch.</td><td>[0-100]</td><td>%</td><td>area and edge</td></tr><tr><td>Percentage of landscape</td><td>pland</td><td>Addition of all patches area over total landscape area converted in percent.</td><td>[0-100]</td><td>%</td><td>area and edge</td></tr><tr><td>Total edge</td><td>te</td><td>Edge length of all patches of the same class.</td><td>≥0</td><td>m</td><td>area and edge</td></tr><tr><td>Core area index</td><td>cai</td><td>Mean of the percentage of each patch in core area, i.e. pixels that have only neighbors in the same class.</td><td>[0-100]</td><td>%</td><td>core area</td></tr><tr><td>Core area</td><td>core</td><td>Mean central area within each patch</td><td>≥0</td><td>ha</td><td>core area</td></tr><tr><td>Total core area</td><td>tea</td><td>Total central area within the plot</td><td>≥0</td><td>ha</td><td>core area</td></tr><tr><td>Contiguity index</td><td>contig</td><td>Sum of contiguity values in 3×3-pixel windows as a fraction of the maximum value, i.e. for a landscape of one class.<br>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 0.</td><td>[0;1]</td><td>/</td><td>Shape</td></tr><tr><td>Shape index</td><td>shape</td><td>Shape complexity (compared to a square)</td><td>≥1</td><td>/</td><td>Shape</td></tr></tbody></table>

<figure><img src="image_6.png" /></figure>
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 treeline-ecotone 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 treeline-ecotone types.

## Characterization of elevational shift and comparison with pattern shift

Elevational shift was determined using a digital elevation model (1 m in resolution, provided by the National Institute of Geographic and Forest Information - IGN) by computing the difference in elevation of the five highest forest pixels at the two dates following the method developed in Feuillet et al. (2020), based on the two previously described supervised classifications. The relationship between spatial pattern change and elevational shifts was then determined by performing bivariate descriptive analyses. Correlations (Kendall rank correlation
**Fig. 3** Correlation of the landscape metrics describing 648 alpine treeline-ecotone 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

<figure><img src="image_7.png" /></figure>

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.

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.

# Results

## 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
<figure><img src="image_8.png" /><figcaption>**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</figcaption></figure>

(~2015; red dots) plots projected on the historical cluster-based subspaces and comparison with historical treeline-ecotone (black dots). Cluster 1 (blue) = diffuse, cluster 2 (yellow) = island and cluster 3 (grey) = discrete treeline ecotone

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Identification of the variables significantly (p<0.05) associated with each historical (~ 1955) treeline-ecotone clusters. No dot is shown for non-significant associations. Positive and negative v.test values indicate high and low metrics values respectively, i.e. those that are over- or under-represented within each cluster. Colors indicate the general category of the metric. See Table 2 for abbreviations of the metrics</figcaption></figure>

("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
<figure><img src="image_10.png" /><figcaption>**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</figcaption></figure>

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 treeline-ecotone 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).
<figure><img src="image_11.png" /><figcaption>**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)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8** a LOESS regression of pattern change against elevational shift. Pattern change is defined as the distance between historical (~1955) and recent (~2015) landscape metrics in the historical factorial PCA plane, irrespective of the direction of change. Kendall's tau is indicated at the top left. n.s.: non-significant. The colors of the points correspond to the pattern</figcaption></figure>

<figure><img src="image_13.png" /></figure>

type in year. b Comparison of elevational shifts between different pattern transitions between historical (~1955) and recent (~2015) ecotones in the historical factorial plane. Black points indicate the means. The dotted red line refers to the overall mean.

## 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.
<figure><img src="image_14.png" /></figure>

Historical landscapes (~1955)

pwc: Tukey; p-adjust: Bonferroni

<figure><img src="image_15.png" /></figure>

Historical landscapes (~1955)

pwc: Tukey; p-adjust: Bonferroni

<figure><img src="image_16.png" /></figure>

Recent landscapes (~2015)

<figure><img src="image_17.png" /></figure>

Recent landscapes (~2015)

pwc: Tukey; p-adjust: Bonferroni

<figure><img src="image_18.png" /></figure>

Recent landscapes (~2015)

pwc: Tukey; p-adjust: Bonferroni

<figure><img src="image_19.png" /></figure>

Historical landscapes (~1955)

Diffuse

Island

Discrete

Discrete

**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-

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

vation of ecotones at which different transitions between tree-line-ecotone spatial types have taken place. Black points refer to average values. Dashed red line corresponds to the mean of all sites. Pairwise p-values (in a–d) are based Tukey's pairwise comparison tests

(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 ecotone<fnref n="1" />), while diffuse

<footnote n="1">¹ Mortality outside the ecotone or outside of existing forest patches could also be due to shade dependence or frost sensitivity of treeline tree species, i.e. stress rather than disturbance (Bader et al. 2021). However, in our study area all treeline ecotones were composed of *P. uncinata*, and climatic gradients were mild, so that different patterns were unlikely caused by such species-specific sensitivities and more likely by disturbance history.</footnote>
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 tree-line, 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
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 dis-
tinguished, but also because they are composed on a
single dominant tree species among herbs and small
shrubs. This contrasts strongly with e.g. tropical tree-
line 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 pat-
terns (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 driv-
ers, especially for changes due to disturbances. How-
ever, 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 eco-
system services (Schwarz 2010; Grafius et al. 2018),
the effects of natural disaster such as typhoons on
landscapes (Lin et al. 2006), or ecological connec-
tivity (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 chang-
ing from diffuse to abrupt ended up located at much
lower elevations than those remaining diffuse, sug-
gests 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 com-
pleted (trees being smaller than our lower size thresh-
old) 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-depend-
ence or sensitivity to frost at the seedlings stage
may also limit tree establishment outside the protec-
tion of forest (this mechanism probably explains the
maintenance of discrete abrupt treelines in southern-
hemisphere 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 unci-
nata*, 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 expecta-
tions, 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 cer-
tain 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 eco-
tones 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 dif-
ferent ecological processes between sites, but rather
reflect disturbance history. Disturbances are predict-
able only to a small extent. Lastly, excluding the
smaller single trees from the classification and pat-
tern descriptions prevented the detection of recent
colonization patterns, important indicators of ongoing
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-dimen-
sional (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. 2000, 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 interac-
tions associated with its dynamics. A first step con-
sists 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 evolu-
tion of treeline-ecotone spatial patterns. This suggests
a variety of associated underlying processes, and thus
the importance of capturing spatial pattern tempo-
ral 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, topo-
geomorphology). Combined with such an analysis,
the methodology presented here allows another step
forward in understanding what drives the dynam-
ics of alpine treeline ecotones and how mountain
ecosystems will respond to global change in all its
multiplicity.

**Acknowledgements** This study was funded by the National Research Agency (ANR SpatialTreeP - https://anr.fr/Projet-ANR-21-CE03-0002) coordinated by Thierry Feuillet. We also benefited from the support of University of Paris 13 and the Pléiade laboratory.

**Funding** The authors have not disclosed any funding.

Declarations

**Conflict of interest** The authors have not disclosed any conflict of interests.

References

Ameztegui A, Brotons L, Coll L (2010) Land use changes as major drivers of mountain pine (*Pinus uncinata* Ram.) Expansion in the Pyrenees. Glob Ecol Biogeogr 19(5):632–641

Ameztegui A, Coll L, Brotons L, Ninot JM (2016) Land use legacies rather than climate change are driving the recent upward shift of the mountain tree line in the Pyrenees. Glob Ecol Biogeogr 25(3):263–273

Ameztegui A, Morán-Ordóñez A, Márquez A, Blázquez-Cas-
ado Á, Pla M, Villero D et al (2021) Forest expansion in
mountain protected areas: trends and consequences for the
landscape. Landsc Urban Plann 216:104240

Bader MY, Rietkerk M, Bregt AK (2007) Vegetation structure and temperature regimes of tropical alpine treelines. Arct Antarct Alp Res 39(3):353–364

Bader MY, Llambí LD, Case BS, Buckley HL, Toivonen JM, Camarero JJ et al (2021) A global framework for linking alpine-treeline ecotone patterns to underlying processes. Ecography 44(2):265–292

Batllori E, Camarero JJ, Gutiérrez E (2010) Current regenera-
tion patterns at the tree line in the Pyrenees indicate simi-
lar recruitment processes irrespective of the past distur-
bance regime. J Biogeogr 37(10):1938–1950

Beloiu M, Beierkuhnlein C (2019) Differences in the spatial structure of two *Pinus cembra* L. populations in the Carpathian Mountains. Forests 10(4):326

Bonanomi G, Rita A, Allevato E, Cesarano G, Saulino L, Di Pasquale G et al (2018) Anthropogenic and environmental factors affect the tree line position of *Fagus sylvatica* along the Apennines (Italy). J Biogeogr 45(11):2595–2608

Breiman L (2001) Random forests. Mach Learn 45(1):5–32

Camarero JJ, Gutiérrez E (2004) Pace and pattern of recent treeline dynamics: response of ecotones to climatic variability in the Spanish Pyrenees. Clim Change 63:181–200
Camarero JJ, Gutiérrez E, Fortin MJ (2000) Spatial pattern of subalpine forest-alpine grassland ecotones in the Spanish Central Pyrenees. For Ecol Manag 134(1-3):1-16

Camarero JJ, García-Ruiz JM, Sangüesa-Barreda G, Galván JD, Alla AQ, Sanjuán Y et al (2015) Recent and intense dynamics in a formerly static Pyrenean treeline. Arct Ant-arct Alp Res 47(4):773-783

Cantegrel R (1986) Seuils écologiques et organisation biologique des pineraies alticoles de Pinus uncinata Ram: nouvelles perspectives en Néouvielle (Pyrénées occidentales)&nbsp;Bulletin de la Société Botanique de France. Actualités Botaniques. 133(1):193-208

Cantegrel R (2019) De l'origine des Pins de montagne euro-péens. Revue forestière française 71(3):249-268

Coops NC, Morsdorf F, Schaepman ME, Zimmermann NE (2013) Characterization of an alpine tree line using airborne LiDAR data and physiological modeling. Glob Change Biol 19(12):3808-3821

Cossart E, Pic J, Le Guen Y, Fressard M (2020) Spatial patterns of vineyard abandonment and related land use transitions in Beaujolais (France): a multiscale approach. Sustainability 12(11):4695

Dearborn KD, Danby RK (2020) Spatial analysis of forest-tundra ecotones reveals the influence of topography and vegetation on alpine treeline patterns in the subarctic. Ann Am Assoc Geogr 110(1):18-35

Dial RJ, Maher CT, Hewitt RE, Sullivan PF (2022) Sufficient conditions for rapid range expansion of a boreal conifer. Nature. https://doi.org/10.1038/s41586-022-05093-2

Elliott GP, Cowell CM (2015) Slope aspect mediates fine-scale tree establishment patterns at upper treeline during wet and dry periods of the 20th century. Arct Antarct Alp Res 47(4):681-692

EnMAP-Box Developers (2019) : EnMAP-Box 3 - A QGIS Plugin to process and visualize hyperspectral remote sensing data. https://enmap-box.readthedocs.io

Escofier B, Pagès J(2008) Analyses factorielles simples et multiples: Objectifs, méthodes, et interprétation, Malakoff. France: Dunod

Feuillet T, Mercier D, Decaulne A, Cossart E (2012) Classification of sorted patterned ground areas based on their environmental characteristics (Skagafjörður, Northern Iceland. Geomorphology 139:577-587

Feuillet T, Birre D, Milian J, Godard V, Clauzel C, Serano Notivoli R (2020) Spatial dynamics of alpine tree lines under global warming: what explains the mismatch between tree densification and elevational upward shifts at the tree line ecotone? J Biogeogr 47(5):1056-1068

Grafius DR, Corstanje R, Harris JA (2018) Linking ecosystem services, urban form and green space configuration using multivariate landscape metric analysis. Landscape Ecol 33(4):557-573

Harsch MA, Bader MY (2011) Treeline form-a potential key to understanding treeline dynamics. Glob Ecol Biogeogr 20(4):582-596

Herrera JM, Alagador D, Salgueiro P, Mira A (2018) A distribution-oriented approach to support landscape connectivity for ecologically distinct bird species. PLoS ONE 13(4):e0194848

Hesselbarth MHK, Sciaini M, With KA, Wiegand K, Nowosad J (2019) Landscapemetrics: an open-source R tool to calculate landscape metrics. Ecography 42:1648-1657

Hofstede RGM, Dickinson KJM, Mark AF, Narváez E (2014) A broad transition from cloud forest to páramo characterizes an undisturbed treeline in Parque Nacional Llanganates, Ecuador. Arct Antactic Alp Res 46(4):1-12

Holtmeier FK, Broll G (2005) Sensitivity and response of northern hemisphere altitudinal and polar treelines to environmental change at landscape and local scales. Glob Ecol Biogeogr 14(5):395-410

Holtmeier FK, Broll GE (2007) Treeline advance-driving processes and adverse factors. Landsc Online 1:1-33

Holtmeier FK, Broll G (2019) Treeline research—from the roots of the past to present time. A review. Forests 11(1):38

Humphries HC, Bourgeron PS, Mujica-Crapanzano LR (2008) Tree spatial patterns and environmental relationships in the forest-alpine tundra ecotone at Niwot Ridge. Colo USA Ecol Res 23(3):589-605

Kassambara A, Mundt F(2020) factoextra: Extract and visualize the results of multivariate data analyses. R package version 1.0.7. https://CRAN.R-project.org/package=factoextra

Körner C (2012) Alpine treelines: functional ecology of the global high elevation tree limits. Springer Science & Business Media

Körner C, Paulsen J (2004) A world wide study of high altitude treeline temperatures. J Biogeogr 31(5):713-732

Kupfer JA (2012) Landscape ecology and biogeography rethinking landscape metrics in a post-FRAGSTATS landscape. Prog Phys Geogr 36(3):400-420

Le S, Josse J, Husson F (2008) FactoMineR: an R package for multivariate analysis. J Stat Softw 25(1):1-18

Lin YP, Chang TK, Wu CF, Chiang TC, Lin SH (2006) Assessing impacts of typhoons and the Chi-Chi earthquake on Chenyulan watershed landscape pattern in central Taiwan using landscape metrics. Environ Manage 38(1):108-125

Malandra F, Vitali A, Urbinati C, Weisberg PJ, Garbarino M (2019) Patterns and drivers of forest landscape change in the Apennines range, Italy. Reg Environ Chang 19(7):1973-1985

McGarigal K, Tagil S, Cushman SA (2009) Surface metrics: an alternative to patch metrics for the quantification of landscape structure. Landscape Ecol 24(3):433-450

McGarigal K(2015) FRAGSTATS help. University of Massachusetts: Amherst, MA, USA, 182

Métailié JP, Paegelow M(2004) Land abandonment and the spreading of the forest in the eastern French Pyrenees in the nineteenth to twentieth&nbsp;

Milian J, Eychenne C, Barthe L (2012) Du symbole à l'action: agriculture et pastoralisme dans la politique du Parc naturel régional des Pyrénées Ariégeoises.&nbsp;Sud-Ouest Européen. Revue Géographique des Pyrénées et du Sud-Ouest 34:51-64

Morley PJ, Donoghue DN, Chen JC, Jump AS (2020) Montane forest expansion at high elevations drives rapid reduction in non forest area, despite no change in mean forest elevation. J Biogeogr 47(11):2405-2416
Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O et al (2011) Scikit-learn: machine learning in Python. J Mach Learn Res 12:2825-2830

R Core Team (2022) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. &nbsp; &nbsp; https://www.R-project.org/

Saulnier M, Cunill Artigas R, Foumou LF, Buscaino S, Métailié JP, Galop D, Py-Saragaglia V (2020) A study of late Holocene local vegetation dynamics and responses to land use changes in an ancient charcoal making woodland in the central Pyrenees (Ariège, France), using pedoanthracology. Veg History Archaeobotany 29(2):241-258

Schwarz N (2010) Urban form revisited—selecting indicators for characterising european cities. Landsc urban Plann 96(1):29-47

Sigdel SR, Liang E, Wang Y, Dawadi B, Camarero JJ (2020) Tree to tree interactions slow down himalayan treeline shifts as inferred from tree spatial patterns. J Biogeogr 47(8):1816-1826

Treml V, Chuman T (2015) Ecotonal dynamics of the altitudinal forest limit are affected by terrain and vegetation structure variables: an example from the Sudetes Mountains in Central Europe. Arct Antarct Alp Res 47(1):133-146

Treml V, Veblen TT (2017) Does tree growth sensitivity to warming trends vary according to treeline form? J Biogeogr 44(7):1469-1480

Vacquié L (2015) Modélisation prospective et échelles spatiales en montagne: application aux Pyrénées françaises (Doctoral dissertation, Université Toulouse le Mirail-Toulouse II)

Van der Linden S, Rabe A, Held M, Jakimow B, Leitão PJ, Okujeni A et al (2015) The EnMAP-Box—A toolbox and application programming interface for EnMAP data processing. Remote Sens 7(9):11249-11266

Wang Y, Pederson N, Ellison AM, Buckley HL, Case BS, Liang E, Camarero J (2016) Increased stem density and competition may diminish the positive effects of warming at alpine treeline. Ecology 97(7):1668-1679

Wardle P (2008) New Zealand forest to alpine transitions in global context. Arct Antarct Alp Res 40(1):240-249

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

<figure><img src="image_20.png" /></figure>