# Beyond mean fitness: Demographic stochasticity and resilience matter at tree species climatic edges

Arnaud Guyennon, Björn Reineking, Roberto Salguero-Gomez, Jonas

Dahlgren, Aleksi Lehtonen, Sophia Ratcliffe, Paloma Ruiz-Benito, Miguel

Zavala, Georges Kunstler

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

► To cite this version:

Arnaud Guyennon, Björn Reineking, Roberto Salguero-Gomez, Jonas Dahlgren, Aleksi Lehtonen, et al.. Beyond mean fitness: Demographic stochasticity and resilience matter at tree species climatic edges. Global Ecology and Biogeography, 2023, 32 (4), pp.573-585. 10.1111/geb.13640. hal-04032620

HAL Id: hal-04032620

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

Submitted on 29 Mar 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>
# Beyond mean fitness: Demographic stochasticity and resilience matter at tree species climatic edges

Arnaud Guyennon<fnref n="1" />

Björn Reineking<fnref n="1" /> ID

Roberto Salguero-Gomez<fnref n="2" /> ID

Jonas Dahlgren<fnref n="3" />

Aleksi Lehtonen<fnref n="4" />

Sophia Ratcliffe<fnref n="5" /> ID

Paloma Ruiz-Benito<fnref n="6" />,<fnref n="7" /> ID

Miguel A. Zavala<fnref n="7" />

Georges Kunstler<fnref n="1" /> ID

**Correspondence**

Georges Kunstler, Univ. Grenoble Alpes, INRAE, LESSEM, 2 rue de la Papeterie - BP 76 F-38402, St-Martin-d'Hères, France.

Email: georges.kunstler@inrae.fr

**Funding information**

Agence Nationale de la Recherche, Grant/Award Number: ANR-16-SUMF-0002, ANR-20-CE32-0005-01 and ANR-20-EBI5-0005-03; EU FP7 ERA-NET, Grant/Award Number: REFORCE; H2020 Environment, Grant/Award Number: 101000574

Handling Editor: Margaret Mayfield

**Abstract**

**Aim:** Linking local population dynamics and species distributions is crucial to predicting the impacts of climate change. Although many studies focus on the mean fitness of populations, theory shows that species distributions can be shaped by demographic stochasticity or population resilience. Here, we examine how mean fitness (measured by invasion rate), demographic stochasticity and resilience (measured by the ability to recover from disturbance) constrain populations at the edges compared with the climatic centre.

**Location:** Europe: Spain, France, Germany, Finland and Sweden.

**Period:** Forest inventory data used for fitting the models cover the period from 1985 to 2013.

**Major taxa:** Dominant European tree species; angiosperms and gymnosperms.

**Methods:** We developed dynamic population models covering the entire life cycle of 25 European tree species with climatically dependent recruitment models fitted to forest inventory data. We then ran simulations using integral projection and individual-based models to test how invasion rates, risk of stochastic extinction and ability to recover from stochastic disturbances differ between the centre and edges of the climatic niches of species.

**Results:** Results varied among species, but in general, demographic constraints were stronger at warm edges and for species in harsher climates. Conversely, recovery was more limiting at cold edges. In addition, we found that for several species, constraints at the edges were attributable to demographic stochasticity and capacity for recovery rather than mean fitness.

**Main conclusions:** Our results highlight that mean fitness is not the only mechanism at play at the edges; demographic stochasticity and population capacity to recover also matter for European tree species. To understand how climate change will drive species range shifts, future studies will need to analyse the interplay between population mean growth rate and stochastic demographic processes in addition to disturbances.

**KEYWORDS**

European tree species, integral projection model, population dynamics, recovery, species range, stochasticity

<footnote n="1">¹Univ. Grenoble Alpes, INRAE, LESSEM, 2 rue de la Papeterie, St-Martin-d'Hères, France</footnote>

<footnote n="2">²Department of Zoology, University of Oxford, Oxford, UK</footnote>

<footnote n="3">³Department of Forest Resource Management, Swedish University of Agricultural Sciences, Umeå, Sweden</footnote>

<footnote n="4">⁴Natural Resources Institute Finland (Luke), Helsinki, Finland</footnote>

<footnote n="5">⁵NBN Trust, Unit F14-18 St Mary's Gate, Nottingham, UK</footnote>

<footnote n="6">⁶Departamento de Biología y Geología, Física y Química Inorgánica, Escuela Superior de Ciencias Experimentales y Tecnología, Universidad Rey Juan Carlos, Móstoles (Madrid), Spain</footnote>

<footnote n="7">⁷Grupo de Ecología y Restauracion Forestal, Departamento de Ciencias de la Vida, Universidad de Alcalá, Edificio de Ciencias, Campus Universitario, Madrid, Spain</footnote>
# 1 | INTRODUCTION

Given the magnitude of the projected climate changes, the distribution of tree species across Europe is likely to change significantly (Cheaib et al., 2012). Understanding how local population dynamics control tree species range limits is crucial to predict range shifts (Schurr et al., 2012). However, we still have a very crude understanding of this relationship.

The Hutchinsonian niche concept states that species range limits occur in the environmental conditions where population performance does not allow them to persist (Godsoe et al., 2017; Hutchinson, 1978). This concept is at the core of a long-standing pos-
tulate in biogeography known as the “abundant centre hypothesis” (Brown, 1984; Pironon et al., 2017), which proposes that local pop-
ulation performance decreases towards the range limits. Although species range limits could also be influenced by other processes, such as dispersal and non-equilibrium dynamics (Holt et al., 2005), the focus on the local population dynamics level proposed by this hypothesis has guided empirical studies on range limits. Most of these studies have focused on the importance of a decline of tree population growth rate towards the species range limits (Csergo et al., 2017; Le Squin et al., 2021; Purves, 2009). However, theoreti-
cal studies have demonstrated that the links between species range limits and local population dynamics could be more complex than merely an effect on mean population growth rate (Holt et al., 2005; Sexton et al., 2009).

Holt et al. (2005) proposed three mechanisms that could lead to stable range limits. The first mechanism is based on the classical idea that species are present where their mean population growth rate allows their presence to be maintained. Previous studies gener-
ally used density-independent models and were thus estimating mean finite population growth rate (Csergo et al., 2017). However, for populations with strong density dependence, such as trees, inva-
sion rate is more appropriate than population growth rate (Le Squin et al., 2021; Pagel et al., 2020; Purves, 2009). The second mecha-
nism, demographic stochasticity, describes the random fluctuations in population size attributable to probabilistic discrete events of in-
dividual tree recruitment and death (quantified by the demographic variance; for an in-depth definition, see Melbourne, 2012), which might ultimately result in local extinction (Grimm & Wissel, 2004). Extinction risk increases when demographic variance increases or when the number of individuals decreases (Engen et al., 2001). The third mechanism, environmental stochasticity, assumes that temporal variations in extrinsic environmental conditions, such as climate extremes or disturbances, might affect population per-
sistence and thus species distribution (Holt et al., 2005; Ovaskainen & Meerson, 2010). In forest ecosystems, the ability of the population to recover from external disturbance is crucial (Seidl et al., 2017).

When exploring changes in local population dynamics at spe-
cies range limits, it is, therefore, crucial to analyse these mecha-
nisms jointly. Demographic and environmental stochasticity remains
largely underexplored but could explain why populations experience
local extinctions at their climatic range limits, even if their population

growth rate is positive (Holt et al., 2005). In addition, these three
mechanisms are not mutually exclusive, but their respective roles
are likely to vary between limits with different physiological con-
straints. Kunstler et al. (2021) recently showed that constraints on
tree demographic rates were different between drought- and cold-
limited climatic limits, because tolerance to different abiotic stresses
requires different adaptive strategies. Finally, even if species were
adapted to their environment, the response of local population dy-
namics might be expected to depend on the climatic optimum of the
species, with stronger responses at range limits in more extreme cli-
matic conditions owing to the disproportionate effects of climatic
extremes (Kunstler et al., 2021).

Recently, several studies have assessed how population dynamics
drive tree species distributions using National Forest Inventories (NFIs)
(Kunstler et al., 2021; Le Squin et al., 2021; Purves, 2009; Schultz
et al., 2022; Thuiller et al., 2014; Yang et al., 2022). However, to our
knowledge, there have been no systematic tests of the respective
roles of demographic and environmental stochasticity for range limits
of tree species (but for shrub response to fire disturbance in South
Africa, see Pagel et al., 2020), probably because most studies either
ignored recruitment or assumed that it was independent of climate
(Kunstler et al., 2021; Le Squin et al., 2021; but see Purves, 2009).
Recruitment, however, is a key stage of the life cycle to explore the
role of stochastic processes properly (Grubb, 1977; Holt et al., 2005).

Here, we assessed the relative importance of the three mech-
anisms presented above on the continental distributions of 25
European tree species (9 gymnosperms and 16 angiosperms). We
extended an integral projection model (IPM) recently developed
for European tree species (Kunstler et al., 2021) by adding species-
specific climate- and density-dependent recruitment models. The
IPMs developed here describe the full life cycle of each species.
As such, they allow us to estimate metrics of population perfor-
mance representative of the three mechanisms proposed by Holt
et al. (2005), then to test how they differ between the centre and
the edges of the species climatic niches (for an overview of the met-
RICS and the tests, see Figure 1). For each species, we differentiated
a "hot and dry" and a "cold and wet" climatic edge. More specifically,
we tested the following hypotheses:

(H₁) Mean population performance, measured by the invasion rate, decreases at the edge relative to the centre (Brown, 1984).

(H₂) The risk of stochastic extinction increases at the edge rela-
tive to the centre because of a higher demographic stochasticity
and/or a smaller tree density at equilibrium (Holt et al., 2005).

(H₃) The ability to recover from stochastic disturbances de-
creases at the edge compared with the centre.

(H₄) We also tested whether the role of the three mechanisms (H₁–H₃) differs between the hot and dry edge versus the cold and wet edge.

(H₅) Finally, we tested whether the strength of limitation at the edge is stronger for species edges in the extremes of European climate (hot edges of hot-distributed species and cold edges of cold-distributed species).
<figure><img src="image_3.png" /></figure>

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

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

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

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

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

FIGURE 1 Conceptual figure illustrating (a-c) the three groups of mechanisms that could limit species distribution at their edges proposed by Holt et al. (2005) and (d-f) the approach to test their responses at the hot and dry or the cold and wet edges, and their variation depending on species climatic centre. (a) Mean fitness is estimated by the invasion rate as the ability of the population to grow when rare (black lines represent two different population trajectories of invaders and the red arrows their estimated invasion rates). (b) Demographic stochasticity is measured as the variability of tree density solely attributable to stochasticity of vital rates and its effect on the risk of stochastic extinction (lines represents stochastic tree density variations in small populations, which results in extinction for the continuous line at the red cross). (c) Resilience to disturbance is measured as the recovery time of a tree population after disturbance (represented by the red arrows). (d) Values of population performances, *m*, for three species at their climatic niche centre (black circle), hot and dry edge (red circle) and cold and wet edge (blue circle) (values along the x-axis represent their positions on the climatic gradient). (e) Index of response at the edge in comparison to the centre ($\Omega^{edge}$) for the three species. Filled points represent significant species responses, and the horizontal line represents the overall effect allowing us to test whether there is a general response of $\Omega$ across all species. (f) Variations of $\Omega$ with the climatic niche centre of species (i.e., median of their positions along the climatic axis). Panels (d-f) presents the expected results according to our hypotheses: population performance decline at the edges, which is equivalent to $\Omega < 0$ at each edge, and $\Omega$ decrease is stronger at the hot and dry edges for species occurring in a hot climate and at the cold and wet edge for species occurring in a cold climate. See Section 2 Materials and Methods for a full description of the metrics, $\Omega$, and the statistical tests.

2 | MATERIALS AND METHODS

2.1 | Forest inventory and climatic data

2.1.1 | National Forest Inventory dataset

To fit vital rate functions (growth, survival and recruitment), we used the European forest inventory data compiled in the FunDivEUROPE project (Baeten et al., 2013). The dataset contains information on individual trees in 91,528 plots across Spain, France, Germany, Sweden and Finland, with records of species identity, diameter at breast height (d.b.h.) and status (alive, dead or harvested) at two surveys. These data allow both the tracking of individual growth and survival and the description of local competition. We focused our analysis on 25 species with >2000 individuals and >500 plots in the data, in order to ensure a good coverage of their range and good estimation of their demographic rates. The mean number of individuals and plots per species were 38,531 and 5568, respectively, and the census time ranged from 5 to 20 years, with a median of 11 years

(for more detail, see Supporting Information Table S2). The minimum d.b.h. of trees was 10 cm, and no data were available on either seed production by conspecific adult trees or seedling and sapling growth/survival. We thus did not disentangle the different stages leading to the ingrowth of a tree with 10 cm d.b.h. (i.e., trees that grew larger than the 10 cm d.b.h. threshold between two surveys).

Survey design varied between countries, but plots were generally circular, with variable radii depending on tree size (largest radius ranging from 10 to 25 m; see protocols in Supporting Information S1). We excluded from the analyses all plots with records of harvesting operations or disturbances between the two surveys, which would otherwise influence our estimation of local competition.

2.1.2 | Climate variables

Following Kunstler et al. (2021), we used two climatic variables known to control the physiological performance of trees to fit our vital rates functions (see Supporting Information 3.1): the sum
of degree days >5.5°C (hereafter sgdd) and the water availability index (wai). The value of wai is calculated as (P - PET)/PET (Ratcliffe et al., 2017), where P is the annual precipitation and PET the potential evapotranspiration. Daily temperature and P were extracted from the study by Moreno and Hasenauer (2016) and PET from the Climatic Research Unit data (Harris et al., 2014). Climate variables were averaged over the period covered between 2 years before the first survey and the year of the last survey, in order to account for the fact that climate might have a delayed effect on demographic rates.

## 2.2 | Integral projection models

An IPM predicts the size distribution, $n(z', t+1)$, of a population at time $t+1$ from its size distribution at time $t$, $n(z, t)$, based on a kernel $K(z', z)$ (with $z$ and $z'$ the size at time $t$ and $t+1$) (Easterling et al., 2000). Here, we consider size as the diameter at breast height (d.b.h.). The kernel, $K(z', z)$, can be split into the survival and radial growth kernel, $P(z', z)$, and the fecundity kernel, $F(z', z)$, as follows: $K(z', z) = P(z', z) + F(z', z)$. The survival and radial growth (hereafter growth) kernel, $P(z', z)$, is defined as $P(z', z) = s(z) \times G(z', z)$, with $s$ being the survival function and $G$ the growth kernel. The fecundity kernel, $F(z', z)$, gives the size distribution of newly recruited trees at time $t+1$ as a function of the size distribution at time $t$. More details on the IPM are provided in Supporting Information 2.1.

Below, we describe the fitting of the recruitment, growth and survival functions (more details are provided in Supporting Information 2.2 and 2.3). Each of these vital rate functions were fitted separately for each species. The impact of climate on vital rates was modelled through two potential alternative shapes: asymptotic or quadratic polynomial. This allowed us to capture alternative climate responses, such as increasing, decreasing or bell shaped. To account for uncertainty in the climatic response shape, for each species, we fitted 100 models to 70% of resampled data and selected each time the best climatic response model based on the Akaike information criterion (i.e., lowest AIC; see Burnham & Anderson, 2002). Then, we evaluated the goodness-of-fit on the remaining 30% of the data (see Supporting Information 2.2). In the remaining analysis, we used the 100 models to translate the uncertainty in the vital rate functions into the metrics of population dynamics.

### 2.2.1 | Recruitment function

We developed a recruitment model that accounted for two main processes: fecundity of the conspecific trees (represented by a power function of the basal area of conspecifics); and the competitive effect of heterospecific and conspecific trees (represented by an exponential function of their basal area; see Supporting Information 2.3). Basal area is the sum of cross-sectional area of trees at 130 cm above ground per plot area (in metres squared per hectare) and is a good surrogate of tree stand biomass. After thorough exploration of different

distributions for the number of recruited trees, we fitted, for each species, a model with a negative binomial distribution using the approach presented above for the climate response. Given that the angle count sampling method used in the German NFI makes recruitment analysis difficult, we excluded this country from the recruitment analysis. We used country-specific intercepts to account for variance owing to national specificities (e.g., differences in protocols between NFIs), and an offset for the different number of years between surveys.

Finally, in the IPM, we included a delay in tree recruitment to account for the time it takes for a sapling to reach the minimum d.b.h., meaning that a newly recruited tree is integrated into the population only after 10 years (see Supporting Information 2.4).

### 2.2.2 | Growth and survival functions

The radial growth and survival were modelled as functions of d.b.h., basal area of competitors and climatic conditions (sgdd and wai), in addition to country-specific intercepts (as for recruitment). A normal random plot effect accounting for unexplained variation at the plot level was included in the growth model. No random plot intercept was included in the survival model, because in most plots no individuals died between the surveys, making the estimation of a random plot effect difficult. Growth models were fitted with a log-normal distribution. Survival models were fitted with a generalized linear model with a binomial error and a complementary log-log link with an offset representing the number of years between the two surveys to account for variable survey times between plots (Morris et al., 2013). Models with interactions between the climate variables and both size and competition were also tested, in order to allow trees to have different climatic responses depending on their size or their competitive environment. Equations are presented in Supporting Information 2.2, and more details are given by Kunstler et al. (2021).

Harvesting is present in all populations and probably leads to a lower natural mortality rate compared with unmanaged forests. Thus, a fixed harvest rate was added to natural mortality in the kernel P. We chose to use the mean annual probability of harvesting over the entire dataset and not to include variability in the harvest rate, because we are focusing on the climatic drivers of species distribution and not on the effect of management.

## 2.3 | Simulations of population dynamics

We simulated dynamics of discretized size distribution $X_t$ [number of individuals per size class, corresponding to integration of $n(z, t)$ over each size class] with a matrix formulation of the IPM as follows:

$$X_{t+1} = (P(BA_t^{tot}) + F(BA_t^{het}, BA_t^{con})) \times X_t, \quad (1)$$

with P and F being the matrices representing the kernel P and F with the d.b.h. range divided into 700 bins (for the numerical integration,
see Supporting Information 2.4 and 2.5). Owing to the density dependence of growth, survival and recruitment rates, the matrix *P* depends on the basal area of competitors at time *t*, $BA_t^{tot}$, and the matrix *F* on heterospecific and conspecific basal area, $BA_t^{het}$ and $BA_t^{con}$, respectively.

To explore the effect of demographic stochasticity on the dynamics of small populations, we also developed an individual-based model (IBM) based on the same vital rate functions as the IPM (see Supporting Information 2.7). In contrast to the IPM, the IBM allows us to simulate discrete finite populations, which is key to exploration of demographic stochasticity in small populations. For each species, we ran 100 IBM and IPM simulations using the 100 resampled vital rate functions to represent their uncertainty.

### 2.3.1 | Equilibrium

All population metrics, with the exception of invasion rate, were computed starting from equilibrium, because observed tree distributions were highly variable in the forest inventory data. We identified the size distribution at equilibrium $X_e$ for each species and climatic conditions by running simulations with various random initial states until the variations in $X_t$ were negligible. There is no direct analytical solution of the equilibrium for density-dependent IPMs. Nonetheless, we checked that our simulations matched the analytical solution for IPMs with a constant transition matrix, *P*, calculated at the equilibrium basal area (as proposed by Rebarber et al., 2012; Townley et al., 2012; see Supporting Information 2.6).

For a small number of species and models, simulations did not reach equilibrium because they predicted a continuous increase in basal area. We discarded models that continued to increase above 200 m² ha⁻¹ of basal area at the end of the simulation (the observed maximum basal area in our dataset is 126 m² ha⁻¹). Given that simulations work on continuous population abundances, there is no strict extinction. However, there might be very low tree density, which will make the computation of some metrics numerically unstable (recovery from perturbations, for example). In the simulations, we defined a lower limit for basal area of 1 m² ha⁻¹ (corresponding to one tree of 10 cm in a circle of 10 m diameter) under which populations were not analysed.

If a simulation did not lead to demographic equilibrium (i.e., basal area <1 m² ha⁻¹ or increasing to >200 m² ha⁻¹), the simulation was discarded from further analysis. If, in a set of simulations, <50% of the models showed demographic equilibrium, we excluded the full set of simulations (see Supporting Information Table S8). In total, only c. 9% of the resampled models did not lead to an equilibrium, and 5% of species' edges were excluded (see Supporting Information 2.6).

## 2.4 | Population metrics

### 2.4.1 | Invasion rate

Invasion rate was used to evaluate mean fitness. In size-structured populations, the invasion rate is measured by the net reproductive rate, $R_0$, of a rare invader (Falster et al., 2017). In our density-dependent

IPM, we estimated $R_0$ by assuming that the basal area of the invader was small and had no density-dependent effects on the matrices *F* and *P*. Doing so allowed us to use the same equation as for density-independent IPMs (the dominant eigenvalue of the matrix *F* ($I - P^{-1}$); see Supporting Information 2.8 and Ellner et al., 2016). As we considered that the invader was rare, we set the conspecific basal area to a low value of 0.1 m² ha⁻¹ in *F*. We computed $R_0$ for two conditions of heterospecific competition: no heterospecific competition (where $BA^{het} = 0$); and a high level of heterospecific competition (where $BA^{het} = 60$ m² ha⁻¹, corresponding to a dense closed forest in our data).

### 2.4.2 | Demographic stochasticity

To evaluate the effect of demographic stochasticity, we derived the time to extinction for finite populations with 250 IBM simulations for each species, climatic condition and resampled model. We initiated simulations by randomly sampling a finite number of trees from the distribution at equilibrium $X_e$ and for a surface of 100 m² and ran the simulation for 1000 years. Following Grimm and Wissel (2004), we extracted the parameter $T_m$ from these simulations, which corresponds to the intrinsic mean time to extinction. Although this provides estimates of time to extinction for a very small population that are likely to be much shorter than for large populations in the field, it has the advantage of providing a tool for comparing stochastic extinction between the edge and the centre.

Then, we derived two metrics that drive time to extinction: the density at equilibrium and the demographic variance. Density at equilibrium was computed from long-term simulations, as presented above. We computed the demographic variance from time series of the total reproductive values (Engen et al., 2009; Jaffré & Le Galliard, 2016) estimated with long-term IBM simulations (3000 years), on a plot area of 1 ha (see Supporting Information 2.8).

### 2.4.3 | Ability to recover from population disturbance

We used damping time to test the ability of a population to recover from disturbances, and two metrics related to short-term responses. Damping time (i.e., the time to converge to a stable size structure after a disturbance) is independent of the size structure of perturbations (see computation in Supporting Information 2.8). This metric, however, does not account for short-term transient evolution of the size distribution after a disturbance, because it is computed around equilibrium. Analytical metrics that characterize population transient dynamics cannot be used with density-dependent models (Capdevila et al., 2020). We, thus, used simulations to derive two other metrics: (1) $T_0$, the time for the first return to equilibrium density (regardless of the tree size distribution); and (2) $T_{half}$, the time after which at least half of the initial tree reduction has been recovered permanently. For each species, we disturbed its population at demographic equilibrium by reducing the density of the largest
trees (above the diameter 66th percentile) by half, then simulated its dynamics for 1000 years. We extracted the two metrics from these simulations. For systems that do not present oscillations (i.e., low damping time), these two metrics will be highly correlated.

## 2.5 | Response at the edge

### 2.5.1 | Niche centre and edge definitions

Owing to the high correlation between the two climate variables, we defined the climatic position of each species along a single climatic gradient, the first axis of the principal component analysis of the two climatic variables (as Kunstler et al., 2021). For each species, the niche centre was the median of the first axis, the hot edge was the 5th percentile and the cold edge the 95th percentile. To ensure that the edges for each species corresponded to valid borders of species distribution, we excluded species and edges where occurrence probability did not decline. Occurrence probability was computed using BIOMOD2 (Thuiller et al., 2009), using presence/absence data from the study by Mauri et al. (2017) (see Supporting Information 3.2 and 3.3).

### 2.5.2 | Tests of response at the edge

Using the 100 resampled species-specific IPMs, we predicted the seven metrics at the climatic centre ($M_{centre}$), at the hot and dry climatic edge ($M_{hot}$) and at the cold and wet climatic edge ($M_{cold}$). We then measured the relative response at the edges as $\Omega_{hot} = \log(M_{hot}/M_{centre})$ and $\Omega_{cold} = \log(M_{cold}/M_{centre})$. For each metric and edge type, we tested whether $\Omega = \log(M_{edge}/M_{centre})$ was significantly different from zero ($H_1-H_4$) using a mixed model with edge type effect (hot or cold) as a fixed effect and a random species effect:

$$ \Omega_{sp,i}^{\text{edge}} = K_{sp}^{\text{edge}} + \sigma_{sp}^{\text{edge}} + \sigma, \quad (2) $$

where $K_{sp}^{\text{edge}}$ is the edge effect.

To test whether the mean climatic position of the species influenced its response at the edge ($H_5$), we analysed, for each metric and edge type, the relationship between $\Omega$ and the species climatic centre conditions. We performed a regression between the value of the median climatic condition and $\Omega$ for each edge (taking into account the variance of the 100 resampled models).

As in any comparative analysis, phylogenetic proximity might introduce bias in our analysis, because we might expect similar population metrics between closely related species. That is why we also tested the robustness of the relationships to the phylogenetic proximity of the tree species with a phylogenetic generalized least squares regression (see Supporting Information 5; Symonds & Blomberg, 2014).

All analysis were conducted in R CRAN (R Core Team, 2021); vital rates were estimated using lme4 (Bates et al., 2015) and glmmTMB (Brooks et al., 2017).

Of the 27 tree species analysed, two were fully discarded: *Acer campestre* owing to the absence of equilibrium; and *Prunus padus* owing to the absence of decline in its prevalence at niche borders.

# 3 | RESULTS

## 3.1 | Metrics of performance at edges relative to the centre

### 3.1.1 | Invasion rate

The invasion rate was generally lower at the hot edge than at the centre, both in the absence of heterospecific competition and at high levels of heterospecific competition (Figure 2). However, the overall effect across all species was stronger at high levels of heterospecific competition. In the absence of competition, 55% of the studied species had a significantly negative relative response $\Omega$, whereas at high levels of competition 64% had significantly negative $\Omega$ (see Supporting Information 4, Table S7). Relative differences at the cold edge were not significant, with or without competition.

### 3.1.2 | Demographic stochasticity

The time to extinction was lower at the hot edge compared with the centre, with $\Omega$ significantly negative for 10 of 22 species (45%; Figure 2). Of the two potential drivers of time to extinction, only $\Omega$ of tree density was also significantly less than zero at the hot edge (11 of 22 species, 50%). The $\Omega$ of the demographic variance was not significantly positive at the hot edge. No significant effects were detected at the cold edge.

### 3.1.3 | Ability of population to recover from disturbance

Of the three metrics used to study recovery from disturbance, we found a significant effect across all species only for damping time; the damping time tended to be longer at the cold edge compared with the centre, indicating slower recovery (8 of 15 species, 53%; Figure 2). There was no difference in damping time between the hot edge and the climatic centre. Lastly, we found no differences at either edge type for the time to reach equilibrium density or the time until the perturbation intensity was permanently halved across all species; we found as many species with positive as negative responses.

## 3.2 | High species variability in response at the edge

There was a high variability in species response at the edge, with several species showing no effect or even a higher mean performance
<figure><img src="image_9.png" /><figcaption>**FIGURE 2** Relative metrics, Ω, by edge. Each symbol represents a species. The error bar is the range of 5th and 95th percentiles. Relative metrics significantly different from zero (see text) are represented by filled circles, otherwise by open circles. Coloured thick horizontal lines represent the edge effect on relative metrics over all species (variable K in Equation 2). The significance of relative metrics over all species (see main text) is shown with a symbol (ns or *).</figcaption></figure>

at the edge rather than a decrease. This was particularly true for invasion rate without competition; see, for example, *Abies alba* or *Picea abies* in Figure 3. Interestingly, for these species, stochastic processes might compensate for the lack of effect on mean fitness at the species level.

At the hot edge, among the nine species that did not show a decline of mean fitness or had a contrasted mean fitness response (one metric decreased and the other increased), three were constrained by the extinction time (see, for example, *Juniperus thurifera* or *Quercus petraea* in Figure 3). At the cold edge, it was the case

for 3 of 10 species (see, for example, *Pinus uncinata* or *Pinus nigra* in Figure 3).

### 3.2.1 | Species responses vary with their climatic centre

Part of the variability in species response was related to the position of the climatic centre of the species. Several metrics of response at the edge were more severely constrained for species with niche
centres in more extreme climates, **Figure 4**. At the hot edge, Ω for the invasion rate without competition, tree density, T₀ and Thalf were significantly more strongly reduced for species with mean climatic positions in hotter and dryer conditions. At the cold edge, only the invasion rate without competition showed a significant trend, with a stronger reduction in species with a mean climate in colder conditions. These results were robust to the inclusion of phylogenetic structure in the residuals; it only affected the relationship of invasion rates at the cold edge, and damping time at the hot edge (see **Supporting Information 5**).

## 4 | DISCUSSION

Despite considerable variation across species, our results show both a consistent decrease in invasion rate and an increase in extinction risk at the hot edge across all species. These patterns were

not observed at the cold edge, where only species occurring in extremely cold climates showed a reduction in these two metrics. In contrast, we found a decrease in resilience to perturbation at the cold edge in most species.

### 4.1 | Several demographic processes drive species distribution

#### 4.1.1 | Invasion rate (H₁)

Our results demonstrate a limitation in the invasion rate at the hot edge. This limitation is exacerbated in species that occur in extremely hot and dry climates. These results are consistent with those of a previous study, which found that life span decreased at the hot edge (Kunstler et al., 2021). The reduced invasion rate at the hot edge is probably the result of this shorter life span, but also of

<figure><img src="image_10.png" /><figcaption>**FIGURE 3** Direction and significance of relative differences of the population performance metrics between edges and the climatic centre (Ω) for each analysed species at the hot and dry edge and the cold and wet edge. Species are ordered from the one showing a significant reduction of invasion rate on the left to the one with the opposite response on the right. Green indicates significant constraints on the metric in agreement with the expected direction (expected direction is indicated by minus and plus signs on the left, decrease or increase in the metric). Red indicates a significant effect in the opposite direction. Blue indicates a non-significant response.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 4 Relative differences of the population performance metrics between edges and the climatic centre (Ω) along the first principal component axis of species mean climatic conditions. Regression lines are plotted when significant (p-value <5%). Species relative metrics significantly differ from zero (see main text) are represented by filled circles, otherwise by open circles.</figcaption></figure>

the lower recruitment at this edge (observed for species occurring in hot climate; see *Supporting Information 4.1*). Previous studies have proposed that competition could strengthen the limitation on mean fitness at the edge (see Louthan et al., 2015). Here, we found only weak evidence for this, because competition increased the number of species with reduced invasion rate at the edge only from 12 to 14 (see *Quercus petraea* and *Abies alba*; Figure 3).

We found no clear evidence of a general limitation in invasion rate at the cold edge. Only species distributed at the cold extreme of the gradient showed signs of reduced invasion rate (such as *Betula* and *Pinus uncinata*). This insight emerged in our new model covering the full life cycle, but not in study by Kunstler et al. (2021), who did not report a decrease in survival or life span for these species. This might be driven by a low recruitment rate at the cold edge for species in extremely cold climates (*Supporting Information Figures S15 and S16*).

We found that the invasion rate in the absence of competition was more strongly constrained in species from harsh extremes of European climate (hot and dry or cold and wet; H₅). In contrast, invasion rate was not limited in climatic conditions typical of temperate regions, where productivity is high (Jung et al., 2007).

Direct comparison with previous studies is difficult because they differ in their way of representing species distribution and computing mean fitness (Le Squin et al., 2021; Purves, 2009; Thuiller et al., 2014). However, even if the structure of the model is different from ours, it is interesting that Purves (2009), in a study on East North American tree species, found a significant decrease in invasion rate at the northern edge but not at the southern edge. The fact that decreases in fitness occur at the opposite edges for tree species in East North America and in Europe might be related to differences in climatic space between these continents. The European southern edge corresponds to a hot and dry climate, whereas the southern edge of tree species in East North America is not limited by drought (Zhu et al., 2015).

### 4.1.2 | Demographic stochasticity (H₂)

The mean time to extinction represents an integrative metric of the demographic stochasticity, which increases when tree density decreases and demographic variance increases (Ovaskainen & Meerson, 2010). At the hot edge among the 10 species showing a
shorter time to extinction, this decline could be related to a change in either demographic variance or tree density, or both. This suggests that these two processes could reinforce each other to result in a stronger reduction in the time to extinction. Our results are also interesting in light of the abundant-centre hypothesis (Brown, 1984), which postulates a decrease in tree density at the edge of a species range. Indeed, our analysis of tree abundance at long-term equilibrium showed that this hypothesis is far from being supported for all edges and species, even with the abundance predicted at equilibrium. This is in agreement with previous large-scale analyses of observed tree abundance by Dallas et al. (2017), Pironon et al. (2017) and Sagarin and Gaines (2002).

### 4.1.3 | Population recovery after disturbances ($H_3$)

Holt et al. (2005) stated that increases in environmental variability can explain range limits despite the absence of a decrease in mean fitness. Here, we explored the role of the time to recovery from disturbances. Disturbance is a key component of environmental variability for tree species. We found an overall significant increase of damping time at the cold edge. This change in a metric of long-term recovery might be connected to the slower tree growth at this edge reported by Kunstler et al. (2021). In contrast, the effect for the short- and medium-term metrics of population recovery ($T_0$ and $T_{half}$) were extremely variable between species (yet, 7 of 14 species showed a longer time to return to the equilibrium density). It is noteworthy that these metrics are extracted from simulations that might lead to a higher variability than the analytical approach used for damping time. At the hot and dry edge, species variability was extremely large. We found evidence for an increase of short- and medium-term metrics only for species with a climatic centre in extreme hot and dry conditions.

A key limitation of our approach to disturbance is that we explored only a single type of abstract disturbance, whereas the real disturbance regime might vary across the species range and play a role in setting distribution limits (Schultz et al., 2022; Senf & Seidl, 2021; Sheil, 2016). In addition, it would be crucial also to explore how interannual variability in climatic conditions, another key component of environmental variability, affects population dynamics. This interannual variability can have a negative, null or positive effect on the population performance depending on the response of demographic rates (Le Cœur et al., 2022). Estimating how natural disturbances and interannual climatic variability might affect tree vital rates and population dynamics at the continental scale remains, however, challenging.

### 4.2 | Stronger constraints at the hot edge ($H_4$)

At the hot and dry edges, we found that the invasion rate was constrained, and we observed increased stochastic risk of extinction for numerous species. Conversely, constraints at the cold edge were

less clear, with an indication of a lower resilience in general and a reduced invasion rate only for species in extreme cold conditions.

These differences might emerge if drought results directly in an increased mortality and higher extinction risk, whereas cold stress could reduce vital rates and population dynamics and thus mechanically increase its response time to disturbances. These differences might also reflect a degree of disequilibrium between the current and potential distribution (Guyennon et al., 2022). Global warming might lead to an increase in drought pressure at the hot edge (Carnicer et al., 2011) and, in contrast, a decrease in mortality at the cold edge in Europe (Neumann et al., 2017). Such disequilibrium in the southern and northern edges has already been identified in East North America (see Talluto et al. (2017)).

### 4.3 | On the challenges of connecting population dynamics and species distribution

It is striking that most studies (including this one) found limited concurrence between mean fitness and species distribution (Kunstler et al., 2021; Le Squin et al., 2021; Purves, 2009; Thuiller et al., 2014). The novelty of our study is that we show that when mean fitness is not constrained at the edge, stochastic processes can play a key role. Nonetheless, there is still a large variability in species responses, with several species having no clear indication of performance constraints or even better performances at edges (Salix caprea and Larix decidua at the hot edge and Juniperus thurifera at the cold edge). Several factors might explain the results for these species. First, we explored species distribution in climate space using only two key climatic variables. Even if these variables discriminate the distribution of the 25 tree species in Europe well (see Supporting Information Figure S14), species distribution might be influenced by other climatic variables or by other abiotic factors, such as soil variables. Second, besides environmental space, species distribution can also be analysed in geographical space (see Pironon et al., 2017). In geographical space, dispersal limitation and a decrease in suitable habitat availability can also explain species range limits in a metapopulation framework (Holt & Keitt, 2000). Third, species distributions are not necessarily in equilibrium with current climatic conditions. Highly managed species (such as Pinus pinaster or Picea abies) can be planted outside their native range. Svenning and Skov (2004) also argued that tree species might still be in the process of slow recolonization since the last glacial age. Here, by initiating our simulations at equilibrium, we effectively removed all legacy effects. Then, as discussed above, climate change could modify constraints at edges (Clark et al., 2021) and explain the difference observed at the hot versus the cold edge. Finally, our ability to capture the complex population dynamics of long-lived organisms, such as trees, is still limited and might explain the poor match with the distribution. For instance, our models do not consider potential variability in seed production, juvenile growth or survival, which could, however, also constrain species ranges (Clark et al., 2021). Indeed, we modelled recruitment with a constant lag to account for the time to grow from a seed to a tree of 10 cm
d.b.h. across species and climate. This shows that our knowledge on the recruitment phase is still extremely limited. In addition, we explored the role of competition in a relatively crude way, considering the competitor effect only through basal area and ignoring the complexity of multispecies interactions. A full exploration of its role would require analysis of how the stochastic dynamics of a multispecies community constrains species range (Godsoe et al., 2017).

## 5 | CONCLUSION

Our study is one of the first to tease apart several mechanisms that could lead to species range limits using field data across a large set of species at the continental scale. Our results show that the mean fitness might not be the only mechanism at play at the edge; demographic stochasticity and the ability of a population to recover also matter for European tree species. Thus, to understand how climate change will drive species range shifts, we encourage ecologists to analyse the full life cycle of trees and explore how the average population growth rate interacts with stochastic processes and recovery from disturbances in driving species ranges.

### AUTHOR CONTRIBUTIONS

A.G. and G.K. conceived the ideas and designed methodology with the help of S.R. and R.S.-G.; S.R. formatted the forest inventory and climatic data with the help of P.R.-B., M.A.Z., G.K. and A.L.; G.K. and A.G. analysed the data and developed the IPM with the help of B.R. and R.S.-G.; A.G. and G.K. led the writing of the manuscript with important inputs from B.R., S.R. and R.S.-G.

### ACKNOWLEDGEMENTS

This work was funded by the REFORCE - EU FP7 ERA-NET Sumforest 2016 through the call "Sustainable forests for the society of the future", with the ANR-French National Research Agency—as national funding agency (grant ANR-16-SUMF-0002). G.K. was funded by the ANR DECLIC (grant ANR-20-CE32-0005-01). G.K. and B.R. were funded by BiodivERsA ERA-NET Cofund project FUNPOTENTIAL (ANR funding grant number ANR-20-EBI5-0005-03) and RESONATE H2020 project (grant number 101000574). The NFI data synthesis was conducted within the FunDivEUROPE project funded by the European Union's Seventh Programme (FP7/2007-2013) under grant agreement number 265171. We thank Gerald Kandler (Forest Research Institute Baden-Württemberg) for his help in building the German data. We thank the MITECO ("Ministerio para la Transición Ecológica y Reto Demográfico"), the Johann Heinrich von Thunen-Institut, the Natural Resources Institute Finland (LUKE), the Swedish University of Agricultural Sciences and the French Forest Inventory (IGN) for making NFI data available. G.K., S.R. and R.S.-G. initiated this work in the working group sAPROPOS (Analysis of PROjections of POpulationS), which was supported by sDiv (Synthesis Centre of the German Centre for Integrative Biodiversity Research, iDiv), funded by the German Research Foundation (FZT 118). M.A.Z. was supported by grant DARE (RTI2018-096884-B-C32). We are

grateful to Fabian Roger for providing code to build the species phylogeny. We are grateful to Laura Touzot for her comments on an earlier version of the manuscript.

### CONFLICT OF INTEREST

The authors declare that there is no conflict of interest regarding the publication of this article.

### DATA AVAILABILITY STATEMENT

No new data are associated with this work. Data are available from Spanish National Forest Inventory (https://www.miteco.gob.es/en/biodiversidad/servicios/banco-datos-naturaleza/informacion-disponible/ifn3_bbdd_descargas.htm.aspx and https://www.miteco.gob.es/en/biodiversidad/servicios/banco-datos-naturaleza/information-disponible/ifn2_descargas.aspx), French National Forest Inventory (https://inventaire-forestier.ign.fr/spip.php?article532) and German National Forest Inventory (https://bwi.info/Download/de/BWI-Basisdaten/ACCESS2003/). The Swedish and Finnish data with blurred geographical coordinates and the climatic data used in the analysis are available (Ratcliffe et al., 2020) from the Dryad Digital Repository (https://doi.org/10.5061/dryad.wm37pvmkw). R code, fitted growth, survival and recruitment models and estimated metrics for each species and climatic conditions are provided in a Zenodo repository (https://zenodo.org/record/6424490).

### ORCID

Björn Reineking https://orcid.org/0000-0001-5277-9181

Roberto Salguero-Gomez https://orcid.org/0000-0002-6085-4433

Sophia Ratcliffe https://orcid.org/0000-0001-9284-7900

Paloma Ruiz-Benito https://orcid.org/0000-0002-2781-5870

Georges Kunstler https://orcid.org/0000-0002-2544-1940

### REFERENCES

Baeten, L., Verheyen, K., Wirth, C., Bruelheide, H., Bussotti, F., Finér, L., Jaroszewicz, B., Selvi, F., Valladaresh, F., Allan, E., Ampoorter, E., Auge, H., Avacariei, D., Barbaro, L., Barnoaiea, I., Bastias, C. C., Bauhus, J., Beinhoff, C., Benavides, R., ... Scherer-Lorenzen, M. (2013). A novel comparative research platform designed to determine the functional significance of tree species diversity in European forests. *Perspectives in Plant Ecology, Evolution and Systematics*, 15(5), 281–291. https://doi.org/10.1016/j.ppees.2013.07.002

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. *Journal of Statistical Software*, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., & Bolker, B. M. (2017). *glmmTMB* balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. *The R Journal*, 9(2), 378–400. Retrieved from https://journal.r-project.org/archive/2017/RJ-2017-066/index.html

Brown, J. H. (1984). On the relationship between abundance and distribution of species. *The American Naturalist*, 124(2), 255–279.

Burnham, K. P., & Anderson, D. R. (2002). *Model selection and multimodel inference, a practical information-theoretic approach*. Springer.

Capdevila, P., Stott, I., Beger, M., & Salguero-Gómez, R. (2020). Towards a comparative framework of demographic resilience. *Trends in Ecology & Evolution*, 35, 776–786.
Carnicer, J., Coll, M., Ninyerola, M., Pons, X., Sanchez, G., & Penuelas, J. (2011). Widespread crown condition decline, food web disruption, and amplified tree mortality with increased climate change-type drought. *Proceedings of the National Academy of Sciences*, 108(4), 1474–1478.

Cheaib, A., Badeau, V., Boe, J., Chuine, I., Delire, C., Dufrêne, E., François, C., Gritti, E. S., Legay, M., Pagé, C., Thuiller, W., Viovy, N., & Leadley, P. (2012). Climate change impacts on tree ranges: Model intercomparison facilitates understanding and quantification of uncertainty. *Ecology Letters*, 15(6), 533–544. https://doi.org/10.1111/j.1461-0248.2012.01764.x

Clark, J. S., Andrus, R., Aubry-Kientz, M., Bergeron, Y., Bogdziewicz, M., Bragg, D. C., Brockway, D., Cleavitt, N. L., Cohen, S., Courbaud, B., Daley, R., Das, A. J., Dietze, M., Fahey, T. J., Fer, I., Franklin, J. F., Gehring, C. A., Gilbert, G. S., Greenberg, C. H., ... Zlotin, R. (2021). Continent-wide tree fecundity driven by indirect climate effects. *Nature Communications*, 12(1), 1242. https://doi.org/10.1038/s41467-020-20836-3

Csergo, A. M., Salguero-Gomez, R., Broennimann, O., Coutts, S. R., Guisan, A., Angert, A. L., Welk, E., Stott, I., Enquist, B. J., McGill, B., Svenning, J.-C., Violle, C., & Buckley, Y. M. (2017). Less favourable climates constrain demographic strategies in plants. *Ecology Letters*, 20, 969–980. https://doi.org/10.1111/ele.12794

Dallas, T., Decker, R. R., & Hastings, A. (2017). Species are not most abundant in the Centre of their geographic range or climatic niche. *Ecology Letters*, 20(12), 1526–1533.

Easterling, M. R., Ellner, S. P., & Dixon, P. M. (2000). Size-specific sensitivity: Applying a new structured population model. *Ecology*, 81(3), 694–708.

Ellner, S. P., Childs, D. Z., & Rees, M. (2016). *Data-driven modelling of structured populations*. Springer.

Engen, S., Lande, R., Sæther, B.-E., & Dobson, F. S. (2009). Reproductive value and the stochastic demography of age-structured populations. *The American Naturalist*, 174(6), 795–804.

Engen, S., Sæther, B.-E., & Møller, A. P. (2001). Stochastic population dynamics and time to extinction of a declining population of barn swallows. *The Journal of Animal Ecology*, 70(5), 789–797.

Falster, D. S., Brännström, Å., Westoby, M., & Dieckmann, U. (2017). Multitrait successional forest dynamics enable diverse competitive coexistence. *Proceedings of the National Academy of Sciences*, 114(13), E2719–E2728.

Godsoe, W., Jankowski, J., Holt, R. D., & Gravel, D. (2017). Integrating biogeography with contemporary niche theory. *Trends in Ecology & Evolution*, 32, 488–499. https://doi.org/10.1016/j.tree.2017.03.008

Grimm, V., & Wissel, C. (2004). The intrinsic mean time to extinction: A unifying approach to analysing persistence and viability of populations. *Oikos*, 105(3), 501–511.

Grubb, P. J. (1977). The maintenance of species-richness in plant communities: The importance of the regeneration niche. *Biological Reviews*, 52, 107–145.

Guyennon, A., Reineking, B., Dahlgren, J., Lehtonen, A., Ratcliffe, S., Ruiz-Benito, P., Zavala, M. A., & Kunstler, G. (2022). Colonization and extinction dynamics and their link to the distribution of european trees at the continental scale. *Journal of Biogeography*, 49(1), 117–129.

Harris, I., Jones, P., Osborn, T., & Lister, D. (2014). Updated high-resolution grids of monthly climatic observations—the cru ts3. 10 dataset. *International Journal of Climatology*, 34(3), 623–642.

Holt, R. D., & Keitt, T. H. (2000). Alternative causes for range limits: A metapopulation perspective. *Ecology Letters*, 3(1), 41–47.

Holt, R. D., Keitt, T. H., Lewis, M. A., Maurer, B. A., & Taper, M. L. (2005). Theoretical models of species' borders: Singles species approaches. *Oikos*, 108(1), 18–27. https://doi.org/10.1111/j.0030-1299.2005.13147.x

Hutchinson, G. E. (1978). *An introduction to population ecology*. Yale University Press.

Jaffré, M., & Le Galliard, J.-F. (2016). Population viability analysis of plant and animal populations with stochastic integral projection models. *Oecologia*, 182(4), 1031–1043.

Jung, M., Vetter, M., Herold, M., Churkina, G., Reichstein, M., Zaehle, S., Ciais, P., Viovy, N., Bondeau, A., Chen, Y., Trusilova, K., & Feser, F. (2007). Uncertainties of modeling gross primary productivity over europe: A systematic study on the effects of using different drivers and terrestrial biosphere models. *Global Biogeochemical Cycles*, 21(4), GB4021.

Kunstler, G., Guyennon, A., Ratcliffe, S., Rüger, N., Ruiz-Benito, P., Childs, D. Z., Dahlgren, J., Lehtonen, A., Thuiller, W., Wirth, C., Zavala, M. A., & Salguero-Gomez, R. (2021). Demographic performance of European tree species at their hot and cold climatic edges. *Journal of Ecology*, 109(2), 1041–1054. https://doi.org/10.1111/1365-2745.13533

Le Coeur, C., Yoccoz, N. G., Salguero-Gómez, R., & Vindenes, Y. (2022). Life history adaptations to fluctuating environments: Combined effects of demographic buffering and lability. *Ecology Letters*, 25(10), 2107–2119.

Le Squin, A., Boulangeat, I., & Gravel, D. (2021). Climate-induced variation in the demography of 14 tree species is not sufficient to explain their distribution in eastern north america. *Global Ecology and Biogeography*, 30, 352–369.

Louthan, A. M., Doak, D. F., & Angert, A. L. (2015). Where and when do species interactions set range limits? *Trends in Ecology & Evolution*, 30(12), 780–792. https://doi.org/10.1016/j.tree.2015.09.011

Mauri, A., Strona, G., & San-Miguel-Ayanz, J. (2017). Eu-forest, a high-resolution tree occurrence dataset for europe. *Scientific Data*, 4(1), 1–8.

Melbourne, B. A. (2012). Demographic stochasticity. In J. Gross (Ed.), *Sourcebook in theoretical ecology*. University of California Press.

Moreno, A., & Hasenauer, H. (2016). Spatial downscaling of european climate data. *International Journal of Climatology*, 36(3), 1444–1458.

Morris, W. K., Vesk, P. A., & McCarthy, M. A. (2013). Profiting from pilot studies: Analysing mortality using bayesian models with informative priors. *Basic and Applied Ecology*, 14(1), 81–89.

Neumann, M., Mues, V., Moreno, A., Hasenauer, H., & Seidl, R. (2017). Climate variability drives recent tree mortality in europe. *Global Change Biology*, 23(11), 4788–4797.

Ovaskainen, O., & Meerson, B. (2010). Stochastic models of population extinction. *Trends in Ecology & Evolution*, 25(11), 643–652.

Pagel, J., Treurnicht, M., Bond, W. J., Kraaij, T., Nottebrock, H., Schutte-Vlok, A., Tonnabel, J., Esler, K. J., & Schurr, F. M. (2020). Mismatches between demographic niches and geographic distributions are strongest in poorly dispersed and highly persistent plant species. *Proceedings of the National Academy of Sciences*, 117(7), 3663–3669. https://doi.org/10.1073/pnas.1908684117

Pironon, S., Papuga, G., Villellas, J., Angert, A. L., García, M. B., & Thompson, J. D. (2017). Geographic variation in genetic and demographic performance: New insights from an old biogeographical paradigm: The Centre-periphery hypothesis. *Biological Reviews*, 92(4), 1877–1909. https://doi.org/10.1111/brv.12313

Purves, D. W. (2009). The demography of range boundaries versus range cores in eastern US tree species. *Proceedings of the Royal Society B: Biological Sciences*, 276(1661), 1477–1484. https://doi.org/10.1098/rspb.2008.1241

R Core Team. (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Retrieved from https://www.R-project.org/

Ratcliffe, S., Dahlgren, J., Lehtonen, A., Adonteng, C., Ruiz-Benito, P., Zavala, M. A., Kaendler, G., Mäkipää, R., & Kunstler, G. (2020). Forest inventory data from Finland and Sweden for: Demographic performance of european tree species at their hot and cold climatic edges, plus ancillary climate data. *Dryad Dataset*. https://doi.org/10.5061/dryad.wm37pvmkw

Ratcliffe, S., Wirth, C., Jucker, T., van der Plas, F., Scherer-Lorenzen, M., Verheyen, K., Allan, E., Benavides, R., Bruelheide, H., Ohse, B., Paquette, A., Ampoorter, E., Bastias, C. C., Bauhus, J., Bonal, D., Bouriaud, O., Bussotti, F., Carnol, M., Castagneyrol, B., ... Vesterdal, L. (2017). Biodiversity and ecosystem functioning relations in
european forests depend on environmental context. *Ecology Letters*, 20(11), 1414-1426.

Rebarber, R., Tenhumberg, B., & Townley, S. (2012). Global asymptotic stability of density dependent integral population projection models. *Theoretical Population Biology*, 81(1), 81-87.

Sagarin, R. D., & Gaines, S. D. (2002). The 'abundant Centre' distribution: To what extent is it a biogeographical rule? *Ecology Letters*, 5(1), 137-147.

Schultz, E. L., Hülsmann, L., Pillet, M. D., Hartig, F., Breshears, D. D., Record, S., Shaw, J. D., DeRose, R. J., Zuidema, P. A., & Evans, M. E. K. (2022). Climate-driven, but dynamic and complex? A reconciliation of competing hypotheses for species' distributions. *Ecology Letters*, 25(1), 38-51.

Schurr, F. M., Pagel, J., Cabral, J. S., Groeneveld, J., Bykova, O., O'Hara, R. B., Hartig, F., Kissling, W. D., Linder, H. P., Midgley, G. F., Schröder, B., & Singer, A. (2012). How to understand species' niches and range dynamics: A demographic research agenda for biogeography. *Journal of Biogeography*, 39(12), 2146-2162.

Seidl, R., Thom, D., Kautz, M., Martin-Benito, D., Peltoniemi, M., Vacchiano, G., Wild, J., Ascoli, D., Petr, M., Honkaniemi, J., Lerer, M. J., Trotsiuk, V., Mairota, P., Svoboda, M., Fabrika, M., Nagel, T. A., & Reyer, C. P. O. (2017). Forest disturbances under climate change. *Nature Climate Change*, 7(6), 395-402. https://doi.org/10.1038/nclimate3303

Senf, C., & Seidl, R. (2021). Storm and fire disturbances in Europe: Distribution and trends. *Global Change Biology*, 27(15), 3605-3619. https://doi.org/10.1111/gcb.15679

Sexton, J. P., McIntyre, P. J., Angert, A. L., & Rice, K. J. (2009). Evolution and ecology of species range limits. *Annual Review of Ecology, Evolution, and Systematics*, 40(1), 415-436. https://doi.org/10.1146/annurev.ecolsys.110308.120317

Sheil, D. (2016). Disturbance and distributions: Avoiding exclusion in a warming world. *Ecology and Society*, 21(1), 10.

Svenning, J.-C., & Skov, F. (2004). Limited filling of the potential range in European tree species. *Ecology Letters*, 7, 565-573.

Symonds, M. R., & Blomberg, S. P. (2014). A primer on phylogenetic generalised least squares. In L. Z. Garamszegi (Ed.), *Modern phylogenetic comparative methods and their application in evolutionary biology* (pp. 105-130). Springer.

Talluto, M. V., Boulangeat, I., Vissault, S., Thuiller, W., & Gravel, D. (2017). Extinction debt and colonization credit delay range shifts of eastern north american trees. *Nature Ecology and Evolution*, 1(7), 1-6.

Thuiller, W., Lafourcade, B., Engler, R., & Araújo, M. B. (2009). Biomod-a platform for ensemble forecasting of species distributions. *Ecography*, 32(3), 369-373.

Thuiller, W., Münkemüller, T., Schiffers, K. H., Georges, D., Dullinger, S., Eckhart, V. M., Edwards, T. C., Jr., Gravel, D., Kunstler, G., Merow, C., Moore, K., Piedallu, C., Vissault, S., Zimmermann, N. E., & Zurell, D. (2014). Does probability of occurrence relate to population dynamics? *Ecography*, 37(12), 1155-1166. https://doi.org/10.1111/ecog.00836

Townley, S., Rebarber, R., & Tenhumberg, B. (2012). Feedback control systems analysis of density dependent population dynamics. *Systems and Control Letters*, 61(2), 309-315.

Yang, X., Angert, A. L., Zuidema, P. A., He, F., Huang, S., Li, S., Li, S.-L., Chardon, N. I., & Zhang, J. (2022). The role of demographic compensation in stabilising marginal tree populations in North America. *Ecology Letters*, 25(7), 1676-1689. https://doi.org/10.1111/ele.14028

Zhu, K., Woodall, C. W., Monteiro, J. V., & Clark, J. S. (2015). Prevalence and strength of density-dependent tree recruitment. *Ecology*, 96(9), 2319-2327.

## BIOSKETCH

Arnaud Guyennon is interested in large-scale species distributions, in particular the link between local population/metapopulation dynamics and global range. His work focuses on European forest species.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Guyennon, A., Reineking, B., Salguero-Gomez, R., Dahlgren, J., Lehtonen, A., Ratcliffe, S., Ruiz-Benito, P., Zavala, M. A., & Kunstler, G. (2023). Beyond mean fitness: Demographic stochasticity and resilience matter at tree species climatic edges. *Global Ecology and Biogeography*, 32, 573-585. https://doi.org/10.1111/geb.13640