The methodological framework for simulating the dispersal of bryophytes under changing climate conditions is presented in Fig. 4. A grid of pixel-specific environmental conditions and dispersal kernels, combining information on species dispersal traits, local wind conditions, as well as landscape features affecting dispersal by wind, is generated and used as input in simulations of species dispersal in the landscape under changing climate conditions.

Data sampling

The European bryophyte flora includes 1817 native or naturalized species⁴¹. Because information on bryophyte species distribution is scarce and very heterogeneous, challenging the application of climatic suitability models⁴², we selected 10 species based upon their representativeness for each of the four main biogeographic elements (i.e., groups of species sharing similar distribution patterns), namely the Arctic-Alpine, Atlantic, Mediterranean, and wide-temperate elements (Supplementary Table 2). For each of these species, we downloaded data from the Global Biodiversity Information Facility (https://www.gbif.org). We excluded data collected before 1960, which represented, on average, 41 ± 12% of the data available, for two reasons. First, old records often lack sufficiently precise location information. Second, we wanted to avoid a potential mismatch between old observations and current climate conditions used for modeling. To complete these data and generate a dataset across the entire range of each species in Europe, we specifically performed a thorough literature review to document their occurrence from more than 600 sources. Only points that were separated by at least 0.1° from each other were subsequently retained for modeling (“ecospat.occ.desaggregation” function in Ecospat 3.1⁴³) to avoid sampling bias and reduce the risk of spatial autocorrelation. Altogether, the number of observations available for each species ranged between 55 and 34,035 (database available from Figshare, https://doi.org/10.6084/m9.figshare.8289650).

Average spore diameter was recorded for each species from Zanatta et al.⁴⁴ and references therein. Species unknown to produce sporophytes were assigned a spore size of 150 µm to take dispersal through larger asexual propagules into account. Spore settling velocities V_t and release height (0.03, 1 and 10 m, which roughly correspond to habitat preferences for ground-dwelling, saxicolous, and epiphytic species, respectively) were determined for each species (Supplementary Table 2) following Zanatta et al.⁴⁴.

Nineteen bioclimatic variables, averaged over the period from 1970 to 2000, were retrieved from WorldClim 1.4 at a resolution of 30 arc-seconds⁴⁵. Although snow is an important driver of species distributions in Arctic regions⁴⁶, the lack of sufficiently detailed information on snow precipitation across Europe prevented us from implementing this variable.

Given the spatial grain of our study, the hypothesis that some species will persist in small microhabitats, where temperatures can be cooler and humidity higher than in the surrounding environment, cannot be rejected. Data at finer scales for both present and future conditions would therefore be desirable⁴⁷. Recently developed methods to generate fine-grained climatic data taking into account microclimatic effects modulated by microtopographic variation in the terrain, vegetation cover and ground properties using energy balance equations cannot, however, yet be implemented across large spatial scales⁴⁸.

For future climate conditions, a wide range of GCMs have been described and their variation represents the largest source of uncertainty in future range prediction studies⁴⁹. No criterion exists to evaluate GCMs, whose performance may vary among regions and variables⁵⁰. Due to computational constrains associated with our migration simulations (see below), we followed Didersky et al.⁵¹. and selected two GCMs that reflected the highest and lowest levels of predicted changes due to climate change for two angiosperm species in Europe⁵⁰, namely MPI-ESM-LR⁵² and HadGem2-ES⁵³. For each GCM, we analyzed two climate change scenarios. These scenarios are expressed by the representative concentration pathways (RCPs), using values comparing the level of radiative forcing between the preindustrial era and 2100. The moderate scenario RCP4.5 assumes 650 ppm CO₂ and 1.0–2.6 °C increase by 2100, and refers to AR4 guideline scenario B1 of IPCC AR4 guidelines. The pessimistic scenario RCP8.5 assumes 1350 ppm CO₂ and 2.6–4.8 °C increase by 2100, and refers to A1F1 scenario of IPCC AR4 guidelines⁵⁴. Climatic data for each GCM and each RCP were averaged for each of the four time periods considered, i.e., 2010–2020, 2020–2030, 2030–2040 and 2040–2050.

Monthly average and daily maximum wind speeds measured at 10 m as well as predicted wind speeds for the same ten-year time periods between 2010 and 2050, were computed from EURO-CORDEX (https://euro-cordex.net). Canopy height data were obtained from the global scale mapping of canopy height and biomass at a 1-km spatial resolution⁵⁵. Wind speed and canopy height were sampled for each pixel and each time-slice to generate kernel maps through time (see below).

Deriving climatic suitability maps

The correlation among the 19 bioclimatic variables was computed from 50,000 random points. To avoid multicollinearity, five bioclimatic variables with a Pearson correlation value of R < 0.7 (as recommended in ref. ³⁰) were selected. These variables were: mean of monthly temperature range, temperature seasonality, mean temperature of warmest quarter, precipitation of wettest month and precipitation of warmest quarter. Since the geographic background should not only reflect the extant, but also the potentially occupied range in the past⁵⁶, and since, in bryophytes, models built from large geographic backgrounds are recommended⁵⁷, pseudo-absence points were sampled from a random selection from all points within the studied area excluding available presence points across Europe.

To account for model uncertainty, we generated ensemble models⁵⁸ using generalized linear model (GLM)⁵⁹ and boosted regression trees (GBM)⁶⁰ with the package biomod2 3.3–7⁶¹. Following Barbet-Massin et al.⁶², 10,000 pseudo-absence points were sampled for GLM and then down-weighted to give them same overall prevalence as presences. For GBM, we sampled a number of pseudo-absence points identical to the number of datapoints. For GLM, the default parameter set (selection procedure via AIC, quadratic model, interaction level = 0, interaction level between variables considered, logit function) was used. For GBM, 5000 trees were included, the maximum depth of each tree was set to 5, the fraction of the training set observations randomly selected to propose the next tree in the expansion was set to 0.8. All other parameters were set to default (Bernouilli distribution, minimum number of observations in the terminal nodes of the trees = 5, shrinkage = 0.001, Number of cross-validation folds = 3). Ten replicates were run and, for each run, 80% of the data was used to calibrate the models, whereas the remaining 20% was kept aside to evaluate the performance of the model using the AUC and TSS metrics. We generated a consensus model of the 10 replicates for each of the GLM and GBM models, wherein each individual model contributed proportionally to its goodness-of-fit statistics. Finally, we computed the suitability at each pixel based on the average of the two GLM and GBM consensus models. Because, despite our thorough literature survey to document species distributions, the number of actually sampled points is a dramatic under-estimation of the actual number of occupied pixels by the species across the study area, all pixels identified as climatically suitable by binarized climatic suitability model projections were employed as initial distribution points for migrations during the first time slice. This could lead to an overestimation of the number of source pixels and raises the issue that, like in any hybrid correlative-mechanistic model, datapoints are employed both for inferring the niche and initiating dispersal simulations, whereas datapoints are themselves the result of a dispersal process⁶³. If, due to dispersal constraints, a species is absent from climatically suitable conditions, climatic suitability models may therefore underestimate species niche range⁶³. Although bryophytes are extremely good dispersers, so that, unlike in some angiosperm species¹⁸, there is no mismatch between the predicted and observed northern limit of distribution³⁹, the present analyses suggest that there is a time-lag of more than a century before newly suitable areas are fully colonized. Nevertheless, our datapoints were sampled across the entire European range of the species, encompassing the full range of climatic conditions that they can experience, so that the potential failure to incorporate localities where the species had not time to disperse yet would not affect the boundaries of our global niche estimate.

The ensemble model was then projected onto future climatic layers using two GCMs and two RCPs per GCM (see above). A key issue with modeling responses to climate change is that we do not fully understand how models made under current conditions will transfer to future conditions. Models developed using too many predictors may run the risk of overfitting to local conditions, restricting the predictive power of the model^64,65. Tests of transferability across taxa and geographic locations have, however, failed to demonstrate consistent patterns, and a general approach to developing transferable models remains elusive⁶⁶. Here, we compared the ROC and maxTSS values computed from the test sets (20% of the data) to those observed at the level of the entire dataset, assuming that these statistics at the level of the entire dataset would substantially drop at the level of the test sets in case of severe issues of overfitting.

The continuous suitability index was transformed into a binary presence/absence model, using maximum TSS to reclassify values.

Dispersal simulations under changing climate conditions

The MigClim model

Dispersal simulations under changing climate conditions were performed with a modified version of Migclim⁶⁷ adapted for wind dispersal. To simulate dispersal under climate change, MigClim requires information on species dispersal capacities, a map of species initial distribution, a map of present climate conditions, and maps of future climate conditions at p intervals that divide the period between time present and the end of the simulations, set by the user, into p climatic periods. In MigClim, source pixels are represented by actually occupied pixels and target pixels are pixels that newly become climatically suitable under climate change. Dispersal simulations are performed from source pixels into target pixels as follows (see Fig. 2 in Engler et al.⁶⁸):

1. 1.

   For each target pixels, all the potential source pixels located within a user-defined range are identified.

2. 2.

   The probability that the target pixel is colonized from all the potential source pixels is computed through the probability Pcol (see below). Optionally, long-distance can be added to the simulation, with a user-defined range and probability.

3. 3.

   These steps are repeated nDisp times, with nDisp typically set to 1 year, until the end of the first climatic period.

4. 4.

   At the end of each of the p climatic periods, pixels that are no longer suitable due to changes in environmental conditions have their values reset to zero, and climatic suitability is updated to reflect environmental change, potentially resulting in a series of newly suitable target pixels.

To define Pcol, MigClim implements a dispersal kernel, which is a vector indicating the probability of dispersal P(x) as a function of the distance from the source. Since dispersal from a source pixel could take place in any direction, MigClim implements a coefficient of diffusion called Surface_j, which corresponds to the number of pixels belonging to a same distance class from the source, to compute the probability that a diaspore from a source pixel ends-up in a target pixel and not in any other pixel located at the same distance range:

To account for the number of diaspores produced by a source pixel j, MigClim implements a parameter called Successful Seeds, which accounts for the number of seeds produced by a source pixel j and allows for turning individual dispersal event probabilities into species spread rates.

Finally, P(pixel_j) values are computed at increasing distances from the source and combined to generate a global probability of colonization Pcol from n potential source pixels:

where P_mat(i) is a probability that is function of the time as the source pixel I became occupied and represents the increase in reproductive potential of source pixel i over time.

Implementing the Wald model in MigClim for simulating dispersal by wind

We developed a new version of Migclim, MigClim 1.7⁶⁹, designed for wind dispersal. While a single kernel was employed across the landscape until the end of the simulations in the previous implementation of MigClim, we employed a wind-dispersal kernel that was sampled for each pixel individually to account for variations in wind conditions and was modified at the same time as the p climatic change intervals to take future wind conditions into account.

We employed the Wald model⁷⁰ to infer dispersal kernels. The WALD model was initially developed⁷⁰ and largely used for wind-dispersed seeds^34,35, so that its use for smaller particles could be questioned. Bryophyte spore-trapping experiments in fact revealed that the tail of the dispersal kernel is, beyond hundreds of meters, not distance-dependent, suggesting that, once a spore is airborne, it could disperse over hundreds to thousands of kilometers, regardless of the distance from the source⁷¹. Spatial genetic structures consistently show, however, significant isolation-by-distance patterns for all distance classes, evidencing that realized colonization rates are distant-dependent⁷² and justifying the implementation of a mechanistic model such as WALD. Furthermore, the WALD model assumes that (i) the slippage velocity between the particles and surrounding air is zero, leading to an infinite drag coefficient, so that the particles and surrounding air parcels are tightly coupled, and that (ii) the diaspore terminal velocity is reached instantly after release. These conditions are precisely met in small particles, which (i) are characterized by low Reynolds numbers, and hence, high drag coefficients, and (ii) almost readily reach terminal velocity after release. The WALD model has thus also been applied to small particles such as pollen grains and spores^73,74.

The Wald model⁷⁰ defines the probability P(x) of colonization at distance x from the source depending on intrinsic (e.g., settling velocity, height of release) and extrinsic (e.g., wind speed) parameters, across the distance range between the source and target pixels, as follows:

With \(\mu ^\prime = \frac{{H\bar U}}{{{\rm{Vt}}}},\lambda ^\prime = \left( {\frac{H}{\sigma }} \right)^2\;{\mathrm{and}}\,\sigma ^2 = 2{{K}}h\frac{{\sigma _w}}{{\bar U}}\)

where x is the distance from the source, Ū is the horizontal mean wind speed at the height of seed release, H is the release height, h accounts for canopy height, Vt is the diaspore terminal velocity, K is von Karman’s constant (0.4), and σ_w is a turbulence parameter corresponding to the standard deviation of the vertical wind velocity.

Starting from the centroid of a source pixel, we finally integrate the Wald model over the shortest and largest distances between the source and target pixels to obtain the probability of colonization of the latter.

Parameter estimation

We derived the turbulence parameter σ_w from wind speed data and canopy height⁵⁵. σ_w = 1.25 u*, where u* is the wind-induced friction velocity depending on canopy height. Since wind speed is typically measured over short vegetation (h_s, set at 0.3 m), we first inferred σ_w above taller vegetation of variable height h from the wind-induced friction velocity measured above short vegetation, \(u^{\ast}_{s}\). Hypothesizing that, at the top of the atmospheric surface layer (~200 m), the mean velocity is not affected by the texture of the ground vegetation,

Following Nathan et al.³⁵, \(u^{\ast}_{s}\) was estimated using von Karman’s formula from the measured wind speed \(\bar U_s\):

where K is von Karman’s constant (0.4), w is the height, at which the wind was measured (here 10 m), and Z0_s = 0.1 h_s.

The friction velocity u* for taller vegetation of height h was then derived using Eq. (5).

To derive the mean wind speed at the height of release H, we implemented either the logarithmic or exponential wind profile⁷⁵. When the height of release H is roughly higher than the canopy height h, the logarithmic wind profile describes the decline in horizontal wind speed with decreasing height above the surface, due to the surface resistance, as:

with Z0 = 0.1 h and d = 0.7 h.

In contrast, when the height of release H is below the canopy, we implemented the exponential wind profile:

with the mean wind speed at canopy height \(\overline {U_h}\) derived from Eq. (6), and α derived from Gualtieri and Secci⁷⁶ as α = 0.24 + 0.096Z₀ + 0.016log²Z₀, where Z₀ = 0.1 h

The parameters Ū, h and σ_w are sampled for each pixel and each time-slice (10 years intervals) to generate kernel maps through time.

We determined “Successful Seed” empirically following the calibration method of Engler and Guisan⁶⁷. Although “Successful Seed” was determined once on the basis of a single empirical study⁷¹ and kept constant across species, this study reported observed colonization rates at distances of hundreds of meters from the source colony, giving us a unique opportunity to make the link between our deterministic models and actual observations, increasing the realism of our approach. P_mat was set to 1.

Finally, in addition to short-distance dispersal events with a probability defined by the kernel described above, any pixel located at >10 km from a potential source could be colonized by LDD. The maximum LDD distance was set to unlimited based on phylogeographic evidence³⁹. Following Robledo-Arnuncio et al.³¹, we employed the results of previous Approximate Bayesian Computation methods for LDD inference from genetic structure data in bryophytes^39,77 to define the range of LDD probability values, set to 0, 10⁻⁴, 10⁻³, 10⁻² and 10⁻¹.

Migclim simulations

We modeled the dispersal of a species under a climate change scenario over a period of 40 years, from 2010 to 2050. Starting with an initial distribution for the year 2010, the climatic suitability of cells was updated every 10 years to reflect the projected changes in climatic conditions under the considered climate change scenario. Since our simulations run over 40 years, we need four different climatic suitability maps. The wind layers were updated at the same 10 years intervals as the climatic data to produce series of spatially and temporally explicit kernel maps. We assume that our species disperse once a year, and hence, our simulations performed a total of 40 dispersal steps between 2010 and 2050. For each 10 years climatic period, pixels were identified as potentially suitable based on the binarized climatic suitability model projections. While climatic suitability thus drove colonization probability, a recent study raised the intriguing idea that spread rates at the migration front increase as climatic suitability decreases as a response to the need to seek for more suitable habitats⁷⁸. In bryophytes, however, such a mechanism would be unlikely as inadequate resources and investment in environmental stress defence typically result in shifts from sexual to asexual reproduction⁷⁹.

For each species, we ran a sensitivity analysis by testing the impact of variation of the free parameters described above: two values of horizontal windspeed Ū (monthly average and daily maximum), three values of spore release height Z₀ (0.03, 1 and 10 m), and four values of LDD probabilities (see above). For each parameter combination, 30 MigClim replicates were performed.

We computed the ratio between the predicted loss of suitable area (fraction of initially suitable cells that became unsuitable by 2050) and the simulated effective colonization rate (fraction of newly suitable cells by 2050 that were effectively colonized) using two extreme values of the LDD probability range, that is, 0 and 0.1.

To determine the time-lag of the colonization of newly suitable habitats, the analyses were run for 500 years, keeping the environmental parameters at their 2050 values.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.