Introduction

Climate change has already prompted species to shift their range toward higher latitudes and elevations (Parmesan 2006; Chen et al. 2011; Virkkala and Lehikoinen 2014). Species persistence in response to climate change may critically depend on their ability to expand their range and track suitable environments. While most current predictive approaches ignore important biological mechanisms such as demography, dispersal, and biotic interactions, these play key roles in species response to environmental change (Urban et al. 2016). In response, several calls have been made for models to incorporate the processes mediating species response (Fordham et al. 2013; Stralberg et al. 2015, 2019) and mechanistic approaches have been developed to improve the realism of projections. Some recently developed models (e.g., dynamic range models and forest landscape models) already improve projections, but more work is required to increase accuracy and usability as they remain rarely employed in conservation planning when compared to correlative species distribution models (SDMs, Guisan and Thuiller 2005; Franklin and Miller 2009; Guisan et al. 2013). The challenge now lies in the development of approaches that are accessible, customizable and integrate multiple processes and their interplay (Thuiller et al. 2013; Urban et al. 2016; McIntire et al. 2022). A strong theoretical background is necessary to take on this challenge and guide the development of approaches to balance complexity and tractability in species distribution modelling (Thuiller et al. 2013).

Explicit modelling of the processes that underlie distribution dynamics is challenging (Hefley et al. 2017; Briscoe et al. 2021). Dynamic range models provide a successful example of incorporating demographic processes and dispersal to improve the accuracy of species distribution projections (Briscoe et al. 2021). They are based on niche theory, assuming that species occur at locations where the environment allows positive growth rates (Hutchinson 1957; Godsoe et al. 2017). However, such models are often difficult to parameterize because measuring growth rate is challenging (McGill 2012) and requires very specific data on species response to abiotic conditions. Indeed, on top of being computationally intensive (Snell et al. 2014), the data required to parameterize these models are rarely available (Urban et al. 2016). Furthermore, local demography on its own may be insufficient to explain broad-scale species distribution, suggesting that processes at broader scales must also be considered (Le Squin, Boulangeat, and Gravel 2021), including dispersal limitations, disturbances, and biotic interactions (Urban et al. 2016; Stephan, Mora, and Alexander 2021).

Another approach recently proposed is derived from metapopulation theory. Metapopulations are expected to persist in heterogenous landscapes if colonization is sufficient to balance local extinctions (Hanski and Ovaskainen 2000). The environment may constrain these two processes and limit metapopulation persistence. Distribution limits eventually emerge over environmental gradients at this location where persistence becomes critical. Furthermore, distributions may be constrained by the amount of suitable conditions in a region. As a result, a species may be absent from a region, or a portion of a gradient, despite the occurrence of suitable conditions if these are not abundant enough or if extinction is too high relative to colonization. Metapopulation theory also makes an ideal framework to incorporate several elements of complexity such as landscape heterogeneity, dispersal, and biotic interactions. Realistic landscape structures can be represented with spatially explicit patch occupancy models (Hanski 1999b; Ramiadantsoa, Hanski, and Ovaskainen 2018). An incidence function is used to scale colonization to patch isolation and extinction risk to patch area (MacArthur and Wilson 1967; Hanski 1999b; Schnell et al. 2013; Huang, Pimm, and Giri 2019). Colonization and extinction rates can also be modulated to represent competitive, mutualistic or antagonistic interactions (Hanski 1999b; Gravel et al. 2011; Fordham et al. 2013; Vissault et al. 2020). The metapopulation framework may thus be understood as a flexible approach to integrate fundamental processes driving distribution dynamics.

Landscapes are highly heterogeneous and dynamic. They are continuously affected by changes that can be slow or fast. Disturbances, environmental changes, and biotic interactions are processes that may cause species distribution to be constantly out of equilibrium with their niche (Ovaskainen and Hanski 2002; Svenning et al. 2014; Boulangeat et al. 2018). Non-equilibrium dynamics are especially marked in plants that are limited by slow demography and restricted dispersal (Svenning and Sandel 2013; Savage and Vellend 2015; Talluto et al. 2017; Vissault et al. 2020). Representing this reality requires an adapted approach and metapopulation theory offers the opportunity to model non-equilibrium dynamics (Hanski and Simberloff 1997; Ovaskainen and Hanski 2002). Recent studies have indeed documented species distributions that do not match the distribution of their favourable climate and that present extinction debts and colonization credits (Savage and Vellend 2015; Talluto et al. 2017). Metapopulation models have shown the trailing edge of current tree distribution to be persisting despite unfavourable climatic conditions as slow demography delays the extinction of populations. At the leading edge, dispersal limitations and competition prevent trees from colonizing favourable habitats (Talluto et al. 2017). The ability of metapopulation models to study and describe dynamic landscapes therefore makes them particularly suitable to study persistence under changing climate.

In this paper, we show how metapopulation theory can be used to model and thus complement the interpretation of species distribution in a changing environment. In addition, we illustrate how metapopulation theory can be used to leverage species distribution models by accounting for the complexity arising from biotic interactions, demography, and landscape structure. To achieve this goal, we first review the theory to account for these key ecological processes in distribution modelling and present associated sources of complexity. We then illustrate the effect and importance of these processes on persistence and distribution dynamics using the Bicknell’s Thrush (Catharus bicknelli) as a case study. The Bicknell’s Thrush is a threatened species in Canada with restricted distribution. We conclude that metapopulation theory can improve the interpretation and the use of habitat projections, notably under increasing climate warming by accounting for the spatial arrangement of habitats.

Key Concepts Arising From Metapopulation Theory

We first review the theoretical framework to incorporate key ecological processes into a mechanistic approach of range dynamics. We frame these processes in the context of a bottom-up system where the distribution of a focal species (e.g., a predator, a habitat specialist or a mutualist) is contingent on the distribution of a trophically lower-level species (e.g., a prey, a vegetation type or a host). Our approach thus integrates dispersal, demography, and biotic interactions. We study distribution dynamics under climate warming and the associated sources of complexity brought by landscape structure. We then contextualize the various effects of climate change on persistence using a conceptual habitat specialist species as an example. The resulting changes to the species’ range support the emergence of distribution changes of greater complexity than predicted by correlative approaches and show that accounting for spatial arrangement of habitats is necessary to capture distribution changes.

Model description

The classic metapopulation model describes species distribution over a set of suitable patches of habitat connected by dispersal (Levins 1969, 1970). Regional dynamics are driven by colonization and extinction events, which corresponding rates depend on local environmental conditions. Together they define the species distribution limits. The dynamics may be complexified with the representation of several trophic levels, where high-trophic level species occur exclusively at locations occupied by lower-level species (Fordham et al. 2013). We adopt the specialist-habitat terminology throughout this study to lighten the text and fit the example, even if the results are more general and can apply to any bottom-up system (e.g., predator-prey or host-mutualist).

Consider a simple system composed of a specialist species tracking the spatial distribution of a dynamic favourable habitat, such as a particular type of vegetation patch providing shelter and food. The model represents the dynamics of the occupancy of three possible states: empty, occupied by the favourable habitat alone (\(H\)) or in co-occurrence with the specialist (\(S\)). The landscape is heterogeneous and each local patch is characterized by the abiotic environmental condition (\(E\)). Dynamics of occupancy are given by the following system of differential equations:

\[ \frac{dH(E)}{dt} = c_{H}(E)H(1-H) - e_{H}(E)H \]

\[ \frac{dS(E)}{dt} = c_{S}(E)S(H(E)-S) - e_{S}(E)S \]

Where \(c(E)\) is the function for colonization rate and \(e(E)\) for the extinction rate. Both are species-specific functions of the abiotic environment such that \(H\) and \(S\) also depend on \(E\), the abiotic environmental conditions. A specialist persists over the landscape in a dynamic equilibrium between habitat availability, colonization, and extinction if its occupancy \(S\) is larger than zero:

\[ S(E^{*}) = H(E^{*}) - \frac{e_{S}(E^{*})}{c_{S}(E^{*})} \]

and the distribution limit is defined by \(S(E^{*}) = 0\), such that it is located where \(H(E^{*}) = \frac{e_{H}(E^{*})}{c_{H}(E^{*})}\). Distribution limits of a habitat specialist are therefore determined by its intrinsic response to the abiotic environment (the ratio \(\frac{e_{S}(E)}{c_{S}(E)}\)), in conjunction with the habitat response to the abiotic environment (\(\frac{e_{H}(E)}{c_{H}(E)}\)).

Graphical representation of range limits

We provide three examples below illustrating how metapopulation theory can reveal some of the complexities of distribution dynamics under a changing climate.

A specialist’s persistence and therefore occupancy is jointly affected by abiotic conditions and habitat availability (occupancy) such that it can be represented graphically with \(\frac{e}{c}\) and \(H(E)\) curves. For simplicity, let’s assume linear relationships. Distribution limit occurs at the position along the environmental gradient where the habitat occupancy curve (green line) crosses the extinction to colonization ratio. An example is illustrated in Figure 1. For a landscape composed of suitable habitat patches, the habitat occupancy is 1 and does not vary with environmental conditions. The specialist’s intrinsic response is less favourable with increasing environmental conditions. Its occupancy for given environmental conditions is defined by the difference between habitat availability and the extinction to colonization ratio curves (\(S^{*} = H-\frac{e}{c}\)). The effect of environmental conditions on its occupancy can be graphically represented at \(E_{0}^{*}\) and \(E_{1}^{*}\). The difference between habitat availability and the extinction to colonization ratio curves (\(S^{*}\); shown by the arrows in Figure 1) is reduced with increasing environmental conditions, illustrating a decrease in the specialist’s occupancy and persistence (\(S(E_{0}^{*}) > S(E_{1}^{*})\)).

Interaction of the specialist and of its habitat’s response can cause indirect distribution dynamics

In a bottom-up system such as predator-prey or a habitat specialist, the response to environmental change does not only depend on the focal species but also on the response of the associated one. The covariation in the response to the environment between the two levels is therefore of critical importance. For instance, the net effect of less favourable environmental conditions to a specialist could be detrimental, without effect, or favourable depending on the effect of the environment to its habitat (Figure 2). Figure 2 A illustrates that specialist occupancy decrease (S) is amplified as environmental conditions harm simultaneously the specialist and its habitat. Conversely, stable specialist occupancy is caused by an equivalent increase of habitat availability or as one level benefits as much as the other suffers (Figure 2 B). An increase in specialist occupancy despite less favourable environmental conditions may occur if one level benefits more than the other suffers (Figure 2 C). Thus, the interaction between levels may have indirect (and counterintuitive) effects on specialist response.

Habitat mismatch affects species distribution shifts

Range limits of a habitat specialist is jointly affected by abiotic conditions and the availability (occupancy) of its habitat. Range shift in response to environmental changes is therefore not only determined by its intrinsic response to the environment, but also by the response of the habitat. As a result, a mismatch between the species response to the environment and its realized distribution may arise, in particular when different trophic levels are not responding at the same rate to environmental change. An example is illustrated in Figure 3. The distribution may shift in the geographic space, for instance toward the north, but it should stay the same in the environmental space if both levels respond similarly (Figure 3, dark shaded area). That said, if a delay or any other factor prevents the habitat from tracking the new environmental conditions, then the habitat curve will shift (Figure 3, green dashed line), and so will the distribution limit (light shaded area). Such mismatch could either benefit or harm the specialist distribution; in this example, the specialist expands to less favourable environmental conditions. The response of the habitat to changing abiotic conditions does influence the specialist distribution, both in extent and in the position of its distribution limits in both the environmental and geographical space.

Metapopulation dynamics may precipitate species decline

The projection of range shifts with correlative SDMs assumes an instantaneous response to environmental change. An implicit assumption is also that a reduction in habitat occupancy translates into an equivalent reduction in the specialist’s range, leading to extinction (Thomas et al. 2004). Metapopulation dynamics may, however, precipitate the decline of a species before the complete disappearance of suitable conditions. Consider a landscape where abiotic conditions are spatially heterogeneous, such as temperature in a mountainous area. The progressive change in this environment, like climate warming, will have two effects on the distribution of suitable patches: the first direct consequence is a reduction in habitat occupancy \(H(E)\), and indirectly follows the increase of the extinction rate with the shrinking of suitable patches. Some favourable patches may also disappear, thereby reducing the landscape connectivity. A non-linear decline of occupancy therefore arises from a linear change in environmental conditions as the ratio \(\frac{e(E)}{c(E)}\) within the specialist’s persistence function increases (Figure 4). This metapopulation effect may not be important at first while suitable habitat is abundant and patches are large, but increases as habitat occupancy decreases, supporting an acceleration of metapopulation prevalence loss to a constant abiotic environmental shift (Hanski and Ovaskainen 2000; Ovaskainen and Hanski 2002).

Spatially explicit landscapes

Analytical tools from metapopulation theory can be used to interpret range limits in spatially explicit heterogeneous landscapes. Metapopulation capacity can be evaluated for realistic landscapes where patch coordinates and size are considered. Metapopulation capacity is measured as the first eigenvalue of the landscape matrix \(M\), where elements \(m_{ij} = exp(-\alpha d_{ij})A_{i}A_{j}\) for \(j \neq i\) and \(m_{ii} = 0\) (Hanski and Ovaskainen 2000). \(\frac{1}{\alpha}\) describes the average dispersal distance, \(d_{ij}\) is the distance between patch \(i\) and \(j\), and \(A_{i}\) is the area of patch \(i\) (refer to Hanski and Ovaskainen (2000) for the full description). Metapopulation capacity is a measure of a species’ ability to maintain itself regionally as a function of connectivity and local extinctions. It provides the means to evaluate conditions for persistence given the spatial arrangement of patches and their size.

Climate change can profoundly alter landscapes as experienced by species; not only does it influence the amount of suitable habitats, but also the capacity of species to persist when colonization and extinction prevail. Consider a mountainous landscape inhabited by a high elevation habitat specialist. The landscape is marked by a steep elevational gradient in temperature where warm temperatures at low elevations exceed the species’ tolerance. The landscape would therefore be divided between suitable cold habitats on mountain tops and unsuitable warmer habitats at the bottom. The topography will not only determine the total surface of suitable conditions, but also the frequency distribution of patch sizes and of distances among mountain tops. As a result, it will influence the connectivity of the landscape and the distribution of patch specific extinction rates.

A schematic example is provided in Figure 5, inspired by the case study that will follow in the next section. Fixing a lower climatic range limit in a hypothetical mountainous landscape, we find nine suitable habitat patches of various sizes, distributed at various distances one from another (Figure 5, left panel). Habitat patches here represent high elevation mountain tops. The warming of climatic conditions causes an elevational shift of lower range limits resulting in the contraction of habitat patches. An equal contraction between patches produces important changes to the landscape’s structure (Figure 5, right panel). The number of patches declines to six for a 63% reduction of total habitat area. Patches become generally smaller from contraction and fragmentation, and the smallest patches go extinct. Further, not only smaller patches are assumed to support smaller population sizes, have superior extinction risks, and produce fewer colonizers (Hanski and Ovaskainen 2000), but the loss and the fragmentation of patches alter species dispersal ability through the loss of connectivity (Huang, Pimm, and Giri 2019). As a result, the metapopulation capacity declines by 82%.

The decrease in metapopulation capacity surpasses that of habitat amount, adding a spatial structure perspective to the assumptions made by correlative approaches. The overall effect of climate warming is not only to modify patch areas, but to change species’ ability to colonize and occupy these patches.

Case Study: Bicknell’s Thrush in North-Eastern America

We illustrate the concepts presented in the previous section with a case study of the Bicknell’s Thrush (Catharus bicknelli), a threatened bird species in Canada (COSEWIC 2009). Bicknell’s Thrush is the smallest Nordic thrush within the Catharus genus and is visually similar to the Grey-cheeked Thrush (Catharus minimus). It migrates in Northeastern America from its wintering grounds in the Greater Antilles and feeds on invertebrates and small fruits (Townsend et al. 2020). Populations are small and were reported to be declining in Canada (COSEWIC 2009). The dispersal of Bicknell’s Thrush is not known with certainty, although it has been suggested that adults nest near the site of previous successful nesting while few yearlings are observed to come back to their site of birth (Rimmer et al. 2001; Collins 2007; Studds et al. 2012). The Bicknell’s Thrush is known to be associated with very dense balsam fir (Abies Balsamea) forests, mostly at high elevations, resulting in a fragmented and highly restricted range (COSEWIC 2009, @cadieux_spatially_2019). This habitat may be ephemeral, as natural disturbances, forestry and stand succession could lead to local extinctions. Furthermore, its distribution in mountainous areas is highly contingent on climate elevation gradients. Climate change could therefore pose a major threat to the persistence of this species as favourable climatic conditions within isolated habitat patches could shrink rapidly (Rodenhouse et al. 2008). Unfavourable abiotic conditions are predicted to increase at the edges of mountaintop fir forest patches with the warming of climate and the limited response capacity of boreal tree species (Talluto et al. 2017; Vissault et al. 2020).

In the following section, we project the changes to the Bicknell’s Thrush breeding range in response to climate forcing using a standard correlative approach. We then leverage the projections using the concepts developed above to analyze the total amount of favourable habitat, the distribution of patch areas, their connectivity, and the metapopulation capacity. Finally, we compare Bicknell’s Thrush favourable landscapes under climate-only change and climate-induced forest change scenarios to illustrate arising climate-habitat mismatch. Thereby, we wish to reveal the joint effects of these two components of Bicknell’s Thrush’s distribution and demonstrate their importance on distribution dynamics.

Methods

Studied region

The Bicknell’s Thrush breeding range was projected for the region where the majority of the Canadian occurrences are identified (COSEWIC 2009; Townsend et al. 2020). Populations are primarily found in the province of Québec, specifically in the Appalachian Mountains in the southeast and the Laurentians Mountains north of the St. Lawrence River. The landscape is composed of boreal, mixed and temperate forests, with their distributions mainly driven by climatic latitudinal and elevational gradients. Mean annual temperature ranges from -4.0 to 7.5 °C in this region, but the Bicknell’s Thrush occupies locations with a more restricted range because of its preference for high-elevation areas. Annual precipitation ranges from 730 to 950 mm.

Data

Distribution data consisted of 6,079 observations of nesting behaviour sampled from 1994 to 2020 and was provided by the le Regroupement QuébecOiseaux (SOS-POP 2021). It contains observations from various sources, including scientific surveys and citizen science. The region of interest was rasterized on a grid of 250 x 250 m cells, where an observation within a cell was defined as a presence and the other cells were left empty. By gridding the region of interest, we considered the locations where one or more observations were made as a single presence, accounting for any potential effects of temporal and spatial pseudo-replication resulting, for example, from multiple sightings of the one individual in the same location.

Temperature, precipitation, elevation, and balsam fir biomass were used to model occurrences. This selection of variables was motivated by expert knowledge as best reflecting Bicknell’s Thrush preference for high elevation and fir dominated habitats (COSEWIC 2009; Townsend et al. 2020). Mean annual temperature and total annual precipitation were interpolated from climate station records for the 1981-2010 period to produce a time series of annual means (McKenney et al. 2013). Data from a georeferenced 10 km climate grid (McKenney et al. 2013) were projected to each 250 m grid cell centroid and adjusted for differences in latitude, longitude and elevation with spatial regression using BioSIM v11 (R’egni‘ere and St-Amant 2007; R’egni‘ere et al. 2017). BioSIM is capable of interpolating climate parameters at specific locations given that digital elevation mapping, which is used as an independent variable in the model, is provided. Forest composition in individual grid cells was obtained from LANDIS-II biomass outputs at simulation time = 0 (see below) which was initialized using provincial ecoforestry provincial maps and temporary forest inventory plots (see Boulanger and Pascual Puigdevall 2021). Absolute fir biomass was considered along with relative biomass to describe Bicknell’s Thrush preference for dense fir stands (Cadieux et al. 2019). Elevation data was obtained using the elevatr R package, then was rasterized at a 250 m resolution (Hollister et al. 2021).

Breeding range model

We estimated the number of observations per cell of the Bicknell’s Thrush using downweighted Poisson regression (Renner et al. 2015); a point process model for presence only data where locations of presences and of quadrature points (spatially random data points necessary to estimate the species distribution) are modelled as a function of environmental variables. In a downweighted Poisson regression, large weights are assigned to quadrature points and small weights to observations such that presence location points comprise a very small portion of the data used to estimate the model. The effect is similar to applying a spatial scaling so that the response is modelled as the number of observations per cell.

We modelled observation records as a function of climate, elevation, and forest composition with 250m resolution as

\[ \text{Presence Points} = \alpha + \beta_1(\text{temperature}) + \beta_2(\text{temperature}^{2}) \] \[ + \beta_3(\text{precipitation}) + \beta_4(\text{elevation}) + \beta_5(\text{firBiomass}) + \beta_6(\text{firRelativeBiomass}) \] \[ + \beta_7(\text{firBiomass} \times \text{firRelativeBiomass}) + \varepsilon \]

where \(\varepsilon \sim Poisson(\lambda)\)). Temperature was considered quadratically to describe both warm and cold limits. Other variables are taught to describe broad preferences and were therefore considered as linear relationships (COSEWIC 2009; Townsend et al. 2020). Absolute fir biomass was also considered in interaction with relative biomass to describe both stand development and composition. We randomly positioned quadrature points to cover most environmental variability and to maximize the accuracy of the likelihood estimation (Renner et al. 2015). We used the fitted model to predict the number of observations per cell that we then converted into the Bicknell’s Thrush breeding range. The breeding range consists of all cells with a predicted density of observation superior to 1 individual per \(km^2\) (i.e., 0.00625 observations per cell).

We assessed model predictive performance using the area under the receiver operating characteristic curve (AUC, Guisan and Thuiller 2005). AUC is essentially a diagnostic tool to measure the quality of prediction of a model. A perfect prediction yields an AUC of 1 while a random prediction yields an AUC of 0.5 (the calculation of the AUC was performed with the auc function of the R package pROC, Robin et al. 2011).

Scenarios

We projected the Bicknell’s thrush breeding range at a 250 m resolution for two scenarios to contrast the impacts of climate with forest composition dynamics over the 2020-2100 time period. We used the Bicknell’s Thrush model along with calibration conditions for the breeding range projection of 2020. We then used climate and forest composition scenarios for the 2040, 2070, and 2100 projections.

The Bicknell’s Thrush breeding range distribution was first projected over time under intermediate climate change conditions using the RCP 4.5 climate forcing scenario (van Vuuren et al. 2011), while keeping forest composition and elevation constant. Future temperature and precipitation projections for 2021-2040, 2041-2070 and 2071-2100 periods were obtained for the RCP 4.5 scenario from the Canadian Earth System Model version 2 (CanESM2). Such anthropogenic climate forcing is increasingly considered as one of the most likely scenarios given current and pledged global climate policies (Hausfather and Peters 2020). Projections were first downscaled to a 10 km resolution using the ANUSPLIN method, and then the BioSIM v11 model was used to interpolate them to a 250 m resolution (R’egni‘ere and St-Amant 2007; McKenney et al. 2011). BioSIM was used to project daily maximum and minimum temperatures (°C), precipitation (mm) by matching georeferenced sources of weather data (in this case the CanESM2 projections over the 10 km Australian National University Spline grid; Hutchinson et al. 2009) to 15,000 random spatially georeferenced points over Quebec, adjusting the weather data for differences in latitude, longitude, and elevation between the source of weather data and each random point using spatial regressions. Universal kriging using elevation as a drifting variable was then used to interpolate climate variables to the 250m grid. As BioSIM stochastically generate future daily weather time series using 30-yrs future climate normals, we averaged results from 30 BioSIM simulations to compute future climate variables that were assigned to the last year of the projection period (e.g., 2021-2040 period became 2040).

Second, we projected Bicknell’s Thrush breeding range over time by only considering climate-induced changes in forest composition (hereafter forest change) under RCP 4.5, i.e., keeping climate variables and elevation constant in the model. Projections of forest composition for the commercial forests of Québec in 2040, 2070, and 2100 were obtained from Boulanger and Pascual Puigdevall (2021) which were produced using the LANDIS-II forest landscape model (FLM, Scheller et al. 2007). LANDIS-II is a spatially-explicit, raster-based FLM that accounts for stand (e.g., interspecies competition, mortality, establishment) and landscape-level processes (e.g., disturbances, seed dispersal, and forest succession). In Boulanger and Pascual Puigdevall (2021), simulations were run at a 10-year time step from the 2020 biomass initial conditions up to 2150 under the RCP 4.5 climate scenario. In these simulations, climate-induced changes in stand dynamics as well as in wildfires were considered. Business-as-usual harvesting as well as spruce budworm outbreaks were also simulated. More details about model parameterization, calibration and results can be found in Boulanger and Pascual Puigdevall (2021).

Analyses

We assessed the impacts of climate-only change and climate-induced forest change on Bicknell’s Thrush persistence by contrasting different aspects of landscape structure from the original and forecasted landscapes. Analyses were run for the southern part of the Québec Province (\(410,080 km^2\)). Breeding range may change with respect to habitat occupancy (here, fir-stand occupancy), the spatial structure of suitable patches, or the species’ ability to occupy available suitable patches. Isolating the effect of the different elements helps to identify the drivers and their respective importance on distribution dynamics. We decomposed the landscape spatial structure into three complementary elements: the number of patches, the patch areas, and the inter-patch distances.

We further compared temporal trends in habitat amount (sensu Fahrig 2013) and persistence using metapopulation capacity (Hanski 2001). We contrasted habitat amount, metapopulation capacity without dispersal constraints, and metapopulation capacity with strong dispersal constraints to reveal different aspects of metapopulation response. Habitat amount alone determines occupancy in the absence of metapopulation dynamics (i.e., the expectation from correlative SDMs); contrasting it with metapopulation capacity under long-distance dispersal reveals the effect of a reduction in patch area on extinction; metapopulation capacity under short dispersal distance reveals the combined effects of reduction in patch area and change in landscape connectivity. Without sufficient knowledge of the Bicknell’s Thrush dispersal kernel, we therefore compared metapopulation capacity for extreme scenarios of dispersal within the range of plausible kernels. We thus evaluated metapopulation capacity for high dispersal limitations (average dispersal distance of 1 km) and for long average dispersal distance (average dispersal distance of 500 km).

Results: Connectivity in addition to habitat amount define realized range

The model had high performance and accurate breeding range prediction with an AUC of 0.95. Proportional fir biomass (slope ± standard error, \(\beta_6 = 3.39 \pm 0.46\)) and mean annual temperature (\(\beta_1 = 1.56 \pm 0.27\)) are best predictors of the breeding range. Furthermore, the quadratic temperature term is significantly negative (\(\beta_2 = -0.28 \pm 0.025\)) such that the model estimates maximum occupancy at 2.7 Celsius (mean annual temperature). Total annual precipitation (\(\beta_3 = -0.0064 \pm 0.00024\)) and elevation (\(\beta_4 = 0.018 \pm 0.00029\)) also have significant effects on occupancy. Fir biomass was not a significant predictor (\(\beta_5 = 0.0082 \pm 0.0081\)) but its interactions with fir relative abundance (\(\beta_7 = -0.048 \pm 0.012\)) and proportional fir biomass were such that stands of dense fir forest are associated with greater occupancy. The model shows a decrease in Bicknell’s thrush predicted occupancy at low elevations of the southern edge and of the northern edge of its distribution area (Figure 6).

Climate and habitat mismatch

Our model projected varying effects of climate change on Bicknell’s Thrush breeding range within the study region (Figure 6). The magnitude of change differed between climate-only and climate-induced forest change scenarios. Shifts at the range edges were more pronounced than within the range under the climate-only scenario, with contraction at the southern edge and expansion at the northern edge. Under the climate-only scenario, extensive expansion was projected as soon as 2040 at high elevation (>600 m) and rapidly warming (up to 3 °C between 2020 and 2040) regions. Multiple northward patches became momentarily suitable with climate warming at moderate elevation areas (500 to 600 m) because of the narrow range of suitable climatic conditions at these lower elevations. Important contraction was projected at the southern range edge with high elevation mountain tops insufficient to cope with temperature increase. As opposed, changes in forest composition are limited due to the slow demography and the limited dispersal of trees (Vissault et al. 2020). As a result, the projected changes to the breeding range under the forest change scenario were much more limited (Figure 6).

Changes in the spatial structure

Projections show that climate and forest changes have major consequences on the spatial structure of suitable patches (Figure 7). The number of patches within the breeding range in the climate-only scenario supports the initial observation of range expansion followed by a rapid contraction with a peak in number of patches in 2040, while the climate-induced forest change scenario shows a decline in number of patches (Figure 7, left panel). Median patch area for both scenarios varied between 0.125 and 0.312 \(km^2\) (minimum and maximum patch area = 0.0625 and 7805 \(km^2\) respectively) and indicates a skewed distribution with a dominance of small patches and few very large ones (Figure 7, centre panel). On the other hand, the median inter-patch distance varied between 218 and 280 km (minimum and maximum inter-patch distance = 0.25 and 809 km respectively) and shows a more balanced distribution with the landscape composed of distanced groups of regionally close patches (Figure 7, right panel). Although the distribution of patch areas in the climate-only scenario appears to remain constant through time, important decreases in the interpatch distances indicate the loss of small, isolated patches, the addition of geographically close patches, and the fragmentation of large patches (Figure 7, centre and right panels). Despite the apparent stability of the breeding range under the climate-induced forest change scenario, important changes in its spatial structure were observed. We observed a rapid decline in the number of patches and, in contrast to changes under the climate-only scenario, the median patch area constantly increased between 2020 and 2100, and the inter-patch distance marginally increased (Figure 7, centre and right panels). Results indicate that close patches became connected to form fewer, but larger patches in addition to the loss of small, isolated patches.

Persistence

We observed an initial increase of 64% (11,743 to 19,344 \(km^2\)) in habitat amount under the climate-only scenario (total change of +9% between 2020 and 2100; Figure 8, full blue line) while habitat amount remained almost stable with only a slight initial decrease of 11% (11,742 to 10,416 \(km^2\)) under the climate-induced forest change scenario (total change of -15% between 2020 and 2100; Figure 8, full orange line). Changes in Bicknell’s Thrush metapopulation capacity approximated those in habitat amount under long average dispersal distance (approximating mean field assumption, Figure 8). However, we observed important divergences in the Bicknell’s Thrush metapopulation capacity from habitat amount when dispersal was restricted (Figure 8). That is, metapopulation persistence accounting for patch size alone (long-distance dispersal) was closely approximated by habitat amount but differed when accounting for both patch size and connectivity (limited dispersal) when changes in the spatial structure of the breeding range were not explained by habitat amount alone.

Perspectives

Using theory and a case study, we show that the climate-induced changes in distribution are likely to be impacted by bottom-up interactions, demography, and landscape structure. We first derived three observations from metapopulation theory. i) A specialist’s range is impacted by changes in habitat occupancy and a habitat-abiotic mismatch affects the range limits of the specialist. ii) The interplay between habitat shrinking and connectivity loss is likely to yield precipitated range contraction and could potentially lead to extinction. iii) The direction and amplitude of the specialist’s response to environmental change vary with the degree of environmental response correlation between trophic levels. We projected the suitable environmental conditions for a well-known bird species whose distribution is jointly affected by climate and vegetation and we analyzed its spatial structure. We showed that climate-induced changes to the distribution of suitable climatic conditions differed from that of its biotic habitat. Furthermore, both the amount of habitat and the spatial structure distribution of the favourable abiotic and biotic conditions are predicted to be impacted by climate change. Thus, we expect the persistence of this species under climate change to be fundamentally affected by metapopulation dynamics. We show that the metapopulation approach complements the understanding of distribution changes by correlative SDMs. The metapopulation dynamics are fundamental to account for changes in distributions’ spatial structure and contribute to accurately capturing climate-induced change in species distribution.

Applications of the metapopulation approach

Many studies have investigated distribution change using metapopulation theory (Fordham et al. 2013; Schnell et al. 2013; Talluto et al. 2017; Huang, Pimm, and Giri 2019; Vissault et al. 2020), but few have considered the complexity arising from biotic interactions and dispersal in context of rapid environmental change. Some aspects have, however, been explored, starting with the development of the theoretical basis for metapopulation dynamics on heterogeneous landscapes. Spatially realistic metapopulation theory has allowed modelling of distribution dynamics in species living in fragmented landscapes (Hanski 1998, 1999a, 2001). The coupling of spatially explicit metapopulation models with dynamic climate change represents a significant conceptual advancement toward realistic projections (Anderson et al. 2009). Our analysis reveals distribution dynamics that previous methods fail to capture, demonstrating the importance of integrating dynamic processes. A simulation study of the Iberian lynx distribution was the first study to consider the interplay of climate change and trophic interactions using a metapopulation approach (Fordham et al. 2013). It showed that these factors could be explicitly considered together, exhibiting distribution dynamics of greater complexity and realism. Moreover, the use of the metapopulation approach has made possible the study of non-equilibrium distributions by the scaling of local processes at the entire distribution (Talluto et al. 2017). Recently, the approach was extended to non-equilibrium dynamics of range shift in response to climate change, opening the way for the study of nonlinear dynamics of migration (Vissault et al. 2020). The metapopulation framework that we propose here builds on these previous developments to advance toward simultaneously projecting changes in demography and dispersal in response to climate change and the multi-species effects of biotic interactions on the distribution of species.

The use of the metapopulation theory to inform conservation goes as far back as 1985 (Shaffer 1985) for species with patchy population structures and has since been adapted to account for specific spatial and population dynamics (Hanski and Simberloff 1997; Fordham et al. 2013; Huang, Pimm, and Giri 2019). In response to exploitation pressure from the logging companies and an extinction risk increasing rapidly, a spatially explicit metapopulation model was used to define the amount of pristine forest needed to assure the survival of the northern spotted owl (Strix occidentalis caurina) in the Northwestern United States (Shaffer 1985; Lamberson et al. 1993). More recently, the incidence function model has been used to study large-scale population dynamics in the Glanville fritillary (Melitaea cinxia) whose distribution has shrunk in Europe to become highly fragmented (Hanski, Kuussaari, and Nieminen 1994; Hanski 1994). The application of these models to case studies demonstrates the value of the metapopulation approach in describing the distribution dynamics of species while being strongly rooted in theory and simple enough to be parameterized using available ecological data (Hanski 1999b).

Metapopulation theory and models effect today how conservation priorities are defined at a variety of scales. The conservation of ecological corridors is the current focus of important initiatives worldwide including, but not limited to, Corridor Appalachien, Nature Conservancy Canada, Yellowstone to Yukon Conservation Initiative and Western Wildway Network Priority Corridor Project, while habitat fragmentation is a criterion of threat for the IUCN Red List (IUCN 2021). Metapopulation theory predicts the scaling of extinction risk with increasing habitat isolation, something other non-spatially explicit approaches do not consider. We further show that a species’ ability to access suitable habitat is a determining factor of its persistence. Equally, assisted colonization and habitat restoration are brought forward as means to support species persistence by increasing respectively colonization rates and habitat occupancy (Ricciardi and Simberloff 2009; Willis et al. 2009; Fordham et al. 2013). Ultimately, metapopulation theory’s main contribution to current conservation initiatives has been to raise attention on the effect of spatial structure of the landscape and dispersal on species persistence.

Metapopulation dynamics

We have shown using a metapopulation approach that a change in the occupancy of a habitat along an abiotic environmental gradient may impact the distribution of higher levels, such as predators or, here, habitat specialists. Therefore, a mismatch between the distribution of the habitat and of the favourable abiotic conditions may affect the position of the specialist’s range edge along an environmental gradient. This is the result of local increases or decreases in colonization and extinction rates from changes in habitat occupancy. Indeed, we observed the Bicknell’s Thrush breeding range projection from climate-induced forest change to remain stable despite important climate change. Less contraction than expected from climate-only projections were observed at the warm edge of southern local habitat patches, indicating the establishment of a mismatch. The high elevation coniferous patches persisted into warmer abiotic conditions, increasing fir occupancy under abiotic conditions where it was previously rare or absent. Furthermore, we observed no range expansion of the specialist where the climate-only scenario predicts northern expansion, revealing a decrease in habitat occupancy for climatic conditions where it was previously available. This observation is likely the result of prolonged persistence (i.e., extinction debt) of the Bicknell’s Thrush where it is already observed despite less favourable abiotic conditions, and the reduction of occupancy in favourable abiotic conditions where it is initially observed (i.e., colonization credit). As a result, non-equilibrium dynamics in Bicknell’s Thrush distribution change are predicted to be an important source of complexity. Forested habitat-abiotic, or resource-abiotic mismatch in response to environmental change is to be expected in natural systems from limitations in dispersal ability and demography (Svenning et al. 2014). Conversely, habitats that shift faster than abiotic conditions may instead decrease specialist persistence in its current range and favour environmental, but not geographical range stability. It is clear that non-equilibrium dynamics in species distributions are key elements of complexity. Hence, predictions are likely to be biased without proper models to account for it.

Correlative SDMs predict direct response of species’ range to habitat amount variations such that a decrease in habitat amount causes an equivalent contraction of the species’ range. However, we have shown that a metapopulation framework offers complementary information to extract from habitat projections. The contraction of a species’ range may be accelerated (or slowed) by metapopulation dynamics. Here, the effect of landscape connectivity interacts with habitat occupancy to generate dynamics of greater complexity. We observed changes in the Bicknell’s Thrush distribution projections in both habitat amount and in spatial structure of habitat patches. Landscape connectivity was affected by newly suitable habitat patches, the extinction of the smallest habitat patches, the fragmentation of the larger ones, and the dispersal distance. In concordance with our intuition, changes in Bicknell’s Thrush persistence were affected by metapopulation dynamics. Persistence could not be explained by changes in habitat amount alone contrasting with the assumption made by correlative SDMs (Figure 8). Furthermore, our results support Hanski (2015) in that connectivity is fundamental to species regional distribution, abundance, and biodiversity in opposition to the habitat amount hypothesis (Fahrig 2013). That is because the species’ ability to use all available habitat is affected by dispersal, which habitat amount alone does not represent.

More favourable abiotic conditions can have unexpected negative impacts on specialists if their habitats are negatively affected. We described this phenomenon as the effect of environmental response correlation between trophic levels (see Key concepts section). It is a concept unique to process-based approaches that cannot be observed directly using a correlative SDM approach as it originates from the joint effects of species-specific environmental performance and of biotic interactions. Although we have not been able to measure it directly with the Bicknell’s Thrush case, we observed an important contrast between its response to climate-only change and to climate-induced forest change: the habitat amount increased in the first scenario and declined in the second. We showed that regionally more favourable climatic conditions to the Bicknell’s Thrush may have, even if only temporarily due to colonization or extinction lags, the opposite effect on its habitat. Therefore, the resulting distribution dynamics from the interplay between trophic levels are complex to predict. Counterintuitive dynamics can arise from species’ environmental correlation. Indeed, the Bicknell’s thrush example illustrates the necessity of documenting the response between trophic levels to a rapidly changing environment as they can produce non-equilibrium dynamics when considered together. It is when the lower trophic level affects the specialist’s colonization and extinction rates asymmetrically that non-equilibrium distribution dynamics are observed. Because metapopulation models can incorporate such dynamics on specialists’ population dynamics, the resulting projections may be of greater realism.

Limitations of the current approach

Metapopulation models require few parameters making them relatively easy to parameterize. Such models have been calibrated for mammals and trees (Fordham et al. 2013; Talluto et al. 2017; Vissault et al. 2020) and can also be for birds although the dispersal component may be challenging to evaluate (Van Houtan et al. 2007; Studds et al. 2012). Even in the absence of a calibrated model, the metapopulation approach offers tools to interpret projections outputs from correlative SDMs. We showed that different aspects of the landscape’s structure could easily be described and studied. An integrated interpretation of distribution changes can be gained from scenarios of dispersal and extinction. Such scenarios can then be used to evaluate species persistence.

Several other factors could also impact the system’s response to climate warming. The model described here is best suited for habitat specialists whose presence depends on the prior establishment of another species that they do not impact, but it could also be generalized to other types of interactions (see Gravel et al. 2011 for an example of a very general model). The concepts developed in this study are more general than the specialist-habitat context in which they are presented and can apply to any bottom-up system. Positive and negative effects of the specialist on its habitat could influence the system’s response to climate change differently. For example, habitat (i.e., resource) removal by the specialist may reduce competition of habitat types and decrease response lag, accelerating the specialist’s decline at the scale of the landscape (Vissault et al. 2020). Prolonged occupancy of the habitat by the specialist may, on the other hand, increase habitat mismatch and support source-sink dynamics. In addition to biotic interactions, metapopulation dynamics at the landscape level could be affected by the interaction of climate change and natural disturbances. For instance, wildfires and insect outbreak regimes are expected to be strongly altered under climate change (Boulanger and Pascual Puigdevall 2021), and associated biodiversity (see Tremblay et al. (2018) for a case study). Both are important drivers of forest dynamics, and our results show that modification in habitat distribution is associated with the specialist response.

We hope that biodiversity actors benefit from more accurate, yet accessible methods to estimate distribution changes. Correlative SDMs are most often used to project distribution changes, but metapopulation models allow a more accurate estimation of colonization and extinction rates with a multispecies perspective. Our estimation of the Bicknell’s Thrush range projected that the biotic interactions will favour the species’ persistence where it already occurs, but will limit its progression further north where firs are not as abundant despite increases in climate suitability. The resulting effect is likely to be the regional contraction of the Bicknell’s Thrush range despite more favourable climatic conditions. Our study highlights the importance of demography, dispersal and biotic interactions on distribution change to rapid environmental change and the importance of spatial structure on the interpretation of projections.

References