Trait‐based projections of climate change effects on global biome distributions

Climate change will likely modify the global distribution of biomes, but the magnitude of change is debated. Here, we followed a trait‐based, statistical approach to model the influence of climate change on the global distribution of biomes.


| INTRODUC TI ON
Climate change is affecting ecosystems globally. Effects on terrestrial plants include changes in physiology, species composition and biomass production (Boone et al., 2018;Esquivel-Muelbert et al., 2019;Wadgymar et al., 2018), which in turn may translate into changes in vegetation and, ultimately, the global distribution of biomes. For example, studies have been reporting a poleward shift of temperate and boreal forest at high latitudes (Bjorkman et al., 2018;Myers-Smith & Hik, 2018), and the encroachment of tropical forest into tropical grasslands (Esquivel-Muelbert et al., 2019;Stevens et al., 2017). Such changes in the distribution of biomes may alter local biotic interactions, changing the structure of ecological communities and the provision of ecosystem services (e.g. Coops et al., 2018;Lavorel and Grigulis, 2012).
Process-based dynamic global vegetation models (DGVMs) are typically used to forecast biome shifts due to climate change.
Generally, DGVMs describe how climate affects plant physiology (e.g. plant photosynthesis and evapotranspiration) through a hypothesized mechanism, which in turn determines the population dynamics and vegetation structure. This way, DGVMs help to understand mechanisms that may result in biome changes (Gonzalez et al., 2010;Sakschewski et al., 2015;Scholze et al., 2006;Sitch et al., 2008). Yet, due to the complex nature of ecological processes, DGVMs simplify real-world variability in plant traits into a limited number of plant functional types (Sitch et al., 2008;Yang et al., 2016). Moreover, while these mechanistic models explicitly include processes that underpin scenario specific changes (Connolly et al., 2017), the combination of many processes modelled simultaneously may result in emergent properties that may be difficult to validate.
As an alternative, correlative plant functional trait-based models may also contribute to predict biome changes due to climate change. These statistical models rely on the direct link between community mean plant traits and the local environment (Boonman et al., 2020), where an equilibrium state is assumed and processes are implicitly described via the impact of the environment on individual plant fitness (e.g. effects on growth, reproduction and survival; Keddy, 1992;Violle et al., 2007). Similar to how DGVMs use physiological processes as an intermediate step between climate and vegetation structure, this intermediate step allows us to include trait variation of different species in similar environments and trait adaptation of similar species in different environments. In addition, this intermediate step introduces a causal link between climate and biomes (e.g. the conversion of water availability to respiration via specific leaf area, and of temperature to metabolism via wood density and height), therefore reducing the risk of establishing spurious relationships (e.g. Fourcade et al., 2018;Santini et al., 2021). Assuming that these models capture the causal relationships between climate and community mean trait values, the relationships are in principle transferable to future scenarios (Dormann et al., 2012;Yang et al., 2016Yang et al., , 2019. Here, the environmental constraints on traits are used as a backbone for predicting biome changes in response to climate change (van Bodegom et al., 2014;Yang et al., 2016Yang et al., , 2019. As argued in the field of distribution modelling (Dormann et al., 2012), these two families of modelling approaches should not be seen as competing, but as complementary approaches with their own strengths and weaknesses. Although the mechanistic approach is common for large-scale vegetation predictions, the fundamentally different statistical, trait-based models may provide an independent line of evidence and may bring new prospects for understanding climate change impacts on vegetation distributions (Lavorel and Grigulis, 2012;Yang et al., 2015;Yang et al., 2016). The comparison between two fundamentally distinct modelling approaches, each with its own limitations, provides insight in the models' ability to capture the right mechanisms and make meaningful predictions (e.g. Newbold et al., 2017). While biomes and future trait distributions have already been predicted using a trait-based approach, this technique has not yet been used to predict the change in global biome distributions due to climate change (Bodegom et al., 2014;Madani et al., 2018).
In this study, we aim to predict the impact of climate change on the global distribution of biomes. We project future biome distribution changes via global community mean trait-environment relationships of specific leaf area, plant height and wood density. Specifically, we investigate how predicted changes in geographic patterns of plant traits can translate into biome contractions, expansions or shifts. By comparing the results of our study to the results of process-based models predicting global vegetation changes (e.g. Alo & Wang, 2008;Gonzalez et al., 2010;Scholze et al., 2006;Sitch et al., 2008), this study may bring new insights into the underlying causes of vegetation change (Boonman et al., 2020;Madani et al., 2018;van Bodegom et al., 2014). In addition, this may help to assess ecological consequences of climate change and may ultimately aid the allocation of large-scale conservation efforts (Laughlin, 2014).

| Modelling approach
We followed a three-step trait-based approach to model the influence of climate change on the global distribution of biomes ( Figure 1). This paragraph provides an overview of these steps, and each step is discussed in more detail in the following sections. As the first step, we predicted the current global distributions of plant community mean specific leaf area (SLA), height and wood density at a 0.5-degree resolution in the WGS84 coordinate system based on recently developed global trait-environment relationships (Boonman et al., 2020). Second, we applied a clustering model that predicts the probability of occurrence of each biome as a function of the combination of the three traits. We trained this model based on the predicted present-day global trait patterns and a map of the present biome distributions (Olson et al., 2001) and predicted current biome distributions based on the highest probability of occurrence of all biomes per grid cell. Third, we projected changes in plant community mean traits and corresponding changes in global biome distributions under two alternative future greenhouse gas emission scenarios, representing a low (RCP 2.6) and extreme (RCP 8.5) emission scenario corresponding to an increase of global mean annual temperature of +1.2°C and +3.5°C in 2070, respectively.

| Trait modelling
We established global trait-environment relationships using the plant community mean trait dataset from Boonman et al. (2020), which includes georeferenced, locally measured trait data representative of natural plant communities. As plant traits, we selected specific leaf area (SLA), plant height and wood density, which are related with plant form and growth rate (Griffin-Nolan et al., 2018), are known to vary with environmental gradients and across geographic regions and biomes (Charles-Dominique et al., 2018;Chave et al., 2009;Freschet et al., 2011;Griffin-Nolan et al., 2018) and could be modelled with sufficient predictive power (Boonman et al., 2020).
For each trait, we first calculated community mean trait values for each community in each study and then combined multiple community means in one grid cell by taking their mean, in order to limit spatial pseudo-replicates. This led to 361, 217 and 125 0.5-degree grid cell average community means for SLA, plant height and wood density, respectively. We fitted the trait-environment relationships at a resolution of 0.5 degrees using six environmental variables: minimum temperature of the coldest month, humidity index, precipitation in the driest quarter of the year, precipitation seasonality and, averaged to a depth of 30 cm, soil cation exchange capacity and soil pH ( Figure 1, step 1). We selected these variables for their ecological relevance while minimizing collinearity (for details see Boonman et al., 2020). We obtained the current bioclimatic variables (averages for 1979-2013) from CHELSA version 1.2 (Karger et al., 2017) and the soil characteristics from SoilGrids250 m (Hengl et al., 2017). The latter were resampled to match the 0.5-degree resolution of the trait and climate data. The humidity index was calculated as the mean annual precipitation divided by the mean annual potential evapotranspiration (Zomer et al., 2008).

F I G U R E 1
Modelling overview. Trait-environment relationships are established and, using an ensemble modelling approach, global plant trait distributions are predicted (step 1). These are linked to observed biome distributions with a classifying Gaussian Mixture Model (GMM) (step 2). For future predictions of biome distributions, future climate predictions are used as input for the fitted ensemble model resulting in future plant trait distributions, which, in turn, are used as input for the fitted GMM resulting in future biome distributions (step 3) We used four modelling techniques (generalized linear model, generalized additive model, random forest and boosted regression trees) and combined their predictions using an ensemble forecasting approach, calculating the current global trait distributions as the average predictions of the four models weighted by the models' predictive performance (cross-validated R 2 ; Appendix S1). Specifics on model settings can be found in Boonman et al. (2020). Further, we assessed the applicability domain of the models using a multivariate environmental similarity surface (MESS) analysis (Appendix S1). This analysis quantifies per grid cell the difference of the most extrapolated environmental predictor (i.e. the predictor with a grid cell value that is furthest outside that predictor's range of values within the dataset compared to all other included predictors) and the environmental range of that predictor covered by locations in the plant trait dataset, while considering the distribution of these data within the global environmental range (Elith et al., 2010).

| Biome modelling
To predict biomes from trait combinations, we fitted a Gaussian Mixture Model (GMM) (Witte et al., 2007). This clustering model links trait combinations to biomes by estimating the probability of occurrence of each biome based on densities of trait combinations found in a specific biome. We calibrated the GMM by overlaying the predicted current global trait distributions with an observed biome map from Olson et al. (2001) (Figure 1, step 2). For this calibration step, we excluded grid cells where at least one of the traits was predicted to have an extreme value (i.e. > the 99th quantile of the predicted values), which resulted in the removal of 3% of all grid cells. From the observed biome map, we omitted the biome "rock and ice" because soil data on ice are not available, and the biomes "lakes," "flooded grassland and savanna" and "mangroves" because water-dependent plant communities were not covered by the trait dataset (Boonman et al., 2020). We included the following biomes: 1) tundra; 2) boreal forest and taiga (hereafter boreal forest); 3) (sub)tropical moist broadleaf forest (hereafter tropical moist forest); 4) (sub)tropical dry broadleaf forest (hereafter tropical dry forest); 5) (sub)tropical coniferous forest (hereafter tropical coniferous forest); 6) temperate coniferous forest; 7) temperate broadleaf and mixed forest (hereafter temperate mixed forest); 8) (sub)tropical grassland, savanna and shrubland (hereafter tropical grassland); 9) temperate grassland, savanna and shrubland (hereafter temperate grassland); 10) montane grassland and shrubland (hereafter montane grassland); 11) Mediterranean forest, woodland and scrub (hereafter Mediterranean woodland); and 12) desert and xeric shrubland (hereafter xeric shrubland). We used untransformed trait predictions rather than the underlying community mean trait observations to calibrate the GMM to ensure a sufficiently large number of replicates per biome to properly train the model, enabling us to predict more biomes (Yang et al., 2016).
We then applied the GMM to predict the present-day distribution of biomes by assigning each 0.5-degree grid cell the biome with the maximum probability of occurrence ( Figure 1, step 2). A clustering model assumes perfect matches between the vegetation types as recorded in the original trait dataset and the observed biome map used to train the GMM. However, biomes naturally include smallscale variations, which may be defined as vegetation mosaics intrinsic to biomes, for example grassland or heathland patches in the temperate forest biome. Nevertheless, we checked for (dis)similarities between the recorded vegetation type and the observed biome on the Olson global map and found a good match between the recorded vegetation type and the biome or vegetation type common to the matching biome (Appendix S2). We assessed the performance of the GMM by using a 10-fold split-sample cross-validation (80% training; 20% testing) and calculating biome-specific True Skills Statistics (TSS) values (Allouche et al., 2006). The overall TSS value was calculated as the average value weighted by the relative extent of each biome.

| Forecasting
We extracted climate projections for 2070 (averages from 2061 to 2080) from three divergent general circulation models (GCMs), CCSM4, CNRM-CM5 and MRI-CGCM3 (Varela et al., 2015), forced with the emission scenarios RCP 2.6 and RCP 8.5 (from now on referred to as low and extreme climate change scenarios). We downloaded the climate data associated with the emission scenarios from CHELSA version 1.2 (Karger et al., 2017). We selected these scenarios to cover a wide range of vegetation changes for the considered range of possible projected changes in climate (Thuiller et al., 2019). For each combination of GCM and climate change scenario, we projected future trait distributions using the fitted ensemble model ( To draw sufficiently grounded conclusions on biome expansions, contractions and shifts in response to climate change, the extent of biomes should be large enough to have a good predictive accuracy (Table 1), which is why we focused on the seven most widespread biomes (jointly covering 88% of the global terrestrial vegetated surface and each covering more than 5% of the land surface): tundra, boreal forest, tropical moist forest, temperate mixed forest, tropical grassland, temperate grassland and xeric shrubland. Biome surface areas were calculated by grid cell size in km 2 using the "area" function of the raster package.
To account for the uncertainty in biome predictions for both the present and future, and thus relax the assumption that cells correspond to the biome with the highest probability of occurrence, we generated 100 alternative biome maps for the present as well as each of the two future scenarios. Each map was generated by resampling each grid cell using the predicted probability of occurrence of each biome per grid cell as weights. The results on biome changes are based on the comparison of each of the 100 current maps to each of the 100 alternative future biome maps per future scenario (Appendix S3). We further assessed the reliability of independently predicted future traits and their combinations by assessing the applicability domain of the future predictions (Appendix S1), comparing observed and projected among-trait correlations (Appendix S1) and checking the degree of overlap between hyper-volumes built with projected community mean trait combinations and observed plant trait values (Appendix S1).

| Trait predictions
Projected changes in trait distributions were less pronounced under the low than the extreme climate change scenario (Figure 2).
Overall, community mean SLA was predicted to decrease by up to 21% and 31% under the low and extreme climate change scenarios, respectively. Yet, it increased by up to 12% and 18% (low and extreme climate change scenario) in extreme dry and tropical wet areas (Figure 2d,g). Community mean plant height generally increased by up to 207% under the low climate change scenario and up to 388% under the extreme climate change scenario, but decreased up to 59% and 60% (low and extreme climate change scenario) in parts of Western Europe, the West Siberian Plain, the boreal shield of Canada and around the Mato Grosso Plateau (Figure 2e,h). Finally, community mean wood density overall increased up to 7% and 11% under the low and extreme climate change scenarios, respectively, but decreased up to 3% and 4% (low and extreme climate change scenario) in the most xeric areas of the world (Figure 2f,i). TA B L E 1 Confusion matrix with 12 biomes. The numbers represent the total number of grid cells averaged over the model runs predicted (row) and observed (column) to be a specific biome. Per observed biome, the shade of the cell indicates the relative difference in number of grid cells predicted to be a specific biome, where darker shades represent larger differences. The colour of each biome matches the colours in the biome maps (Figure 4). The first column of the table shows the relative surface area (km 2 ) of each biome according to the observed biome map. The second column shows the predictive accuracy of the GMM as TSS values per biome. Bold values indicate the biomes with a large enough extent and good predictive accuracy, which are considered in analyses on future biome distributions.

F I G U R E 2
Trait predictions under the current climate, and the relative change per trait for the low (RCP 2.6) and extreme (RCP 8.5) climate change scenarios, using fitted trait-environment models and future climate (Appendix S4). The first row indicates predictions for Specific Leaf Area (SLA), the second row for height, and the third row for wood density (WD). The green squares represent the grid cells with community mean trait observations. The most extreme future changes (values smaller than the 0.01 percentile and larger than the 0.99 percentile) were removed from these plots to enhance interpretation. No predictions were made for Greenland or the south of Egypt, as there are no soil predictions and precipitation seasonality is zero, respectively F I G U R E 3 Density plots showing the locations of the biomes in twodimensional trait space for each of three pairwise combinations of traits. Colours represent biomes and the centre(s) of each contour plot corresponds to the most common combination(s) of traits for that specific biome The predictive accuracy (TSS) varied greatly between the 12 terrestrial biomes (  Figure 3 or 2) biomes identified during the original trait measurements, for example Mediterranean woodland being confused with temperate grassland (Appendix S2).
Despite fine-scale differences, the weighted average TSS across all biomes was 0.51, resulting in a global pattern of predicted biomes broadly consistent with the actual distribution of biomes (Figure 4a,b and Appendix S5).

| Biome predictions
Tropical biomes were projected to flourish under climate change (Figures 4 and 5).  (Figures 4 and 6). Yet, we also predicted tropical grassland expansion at the expense of tropical moist forest, which reflects the battle of occurrence between tropical moist forests and grasslands, both having a high probability of occurrence in tropical climates (Figures 4 and 6 and Appendix S8).
Xeric shrublands were projected to expand by 6.4% (4.9%-8.0%) to 13.4% (10.8%-15.3%; low and extreme climate change scenario) in the Kalahari and Gobi regions, Australia and North America, replacing tropical and temperate grassland (Figures 4-6). This is mostly reflecting projected decreases of community mean SLA (Appendixes  4 and 6). Specifically, we projected that 51.2% (49.3%-53.1%) to 55.5% (53.6%-57.3%) of the current extent of boreal forest will shift northwards under low and extreme climate change, respectively (Figures 4 and 6). As a result, the southern border of the Eurasian boreal forest would be displaced northwards (Figure 4) (Figures 4 and 6). Alongside these shifts we found decreases in extent of 5.5% (3.6%-7.2%) to 11.9% (10.5%-13.2%) for boreal forest and 18.5% (15.8%-21.0%) to 26.0% (24.0%-28.9%) for temperate mixed forests (low and extreme change scenarios; Figure 5). Declines in temperate mixed forest were expected to occur mostly in Europe, the east of North America and China, along the southern edge of the biome (Figure 4) (Figures 4 and 6). This reflects projected decreases in community mean SLA and increases in community mean height and wood density in these areas (Appendixes S4,S6,S7).

| D ISCUSS I ON
In this study, we have built upon global models linking climate to plant traits and plant traits to biomes (e.g. Boonman et al., 2020;Madani et al., 2018;van Bodegom et al., 2014) to project biome distribution changes under a low and extreme climate change scenario.
Our results suggest that tropical biomes will flourish, while temperate biomes will shift northwards and reduce in extent, particularly tundra. Despite the limited number of traits used to characterize biomes, overall our projections of biome contractions, expansions and shifts are robust to uncertainty in biome classification. This indicates that trait-based models are a promising complementary approach to mechanistic models to assess the effect of climate change on biome distributions, especially as more data and traits become available in the future.

| Projected biome changes
In the northern regions, we projected a climate change-induced  where changes have been linked to increasing temperatures and an increased atmospheric CO 2 level (Alo & Wang, 2008). Further, although woody cover increases in the tundra regions, we expect overall reductions in the area of boreal and temperate mixed forests, which will likely be replaced by temperate grassland (Figures 5 and   6). These results are also largely concordant with those produced by DGVM studies (Alo & Wang, 2008;Gonzalez et al., 2010;Scholze et al., 2006;Sitch et al., 2008). More specifically, when comparing our results to those of Park et al. (2015), who used biomes and climate change scenarios comparable to this study, our results are highly similar: we predict a 5.5% decrease in boreal forest surface area under low climate change, which matches their predicted decrease of boreal woody plant habitats of 5.3%.
We projected an increase in forest extent in tropical regions, caused by the woody encroachment into tropical grasslands  (Scholze et al., 2006). The fact that fire is not included as a predictor in our model may partly explain our difficulty to distinguish between the tropical grassland and tropical moist forest biomes in the clustering process (Appendix S8).

| Implications
Projected changes in biome distribution may have implications for ecosystem services. For example, changes in tundra and tropical moist forest extent may alter the global carbon budget, desertification of grasslands may have important ramifications for disease control, climate regulation, food provisioning and soil erosion control (Millennium Ecosystem Assessment, 2005), and overall biome changes may disrupt key biotic interactions (e.g. Coops et al., 2018;Zhou et al., 2017). Note, however, that by relying on traits instead of species, we recede from making assumptions on species adaptability The overlap between the outer and inner circle depicts the amount of surface area of the current biome that remains the same in the future. The overhang of the inner circle shows the proportion of area that was another biome in the current situation (which can be found by tracing back the arrows) but changed into this biome in the future. Biome colours correspond to the colours in Figure 4. BoF, boreal forest; Tun, tundra; TrMF, tropical moist forest; TeBF, temperate mixed forest; TrG, tropical grassland; TeG, temperate grassland; Des, xeric shrubland range expansions and delayed local extinctions in relation to climate change, but also because information on the age of vegetation structures (e.g. young versus old forests) was not included, our results must be interpreted as predictions of the distribution of biomes at equilibrium with the future climate.
We also acknowledge that additional factors may play a role in shaping vegetation patterns. First, we did not include vegetation succession nor the effect of increased carbon dioxide concentrations on plant growth, while DGVMs do include these. Second, we did not consider biotic feedback like herbivores, which influence vegetation by grazing, browsing and trampling and may thus reduce woody encroachment in low-rainfall, open areas (Stevens et al., 2017). Third, human impact, such as the increasing human land demand for food production, will likely counteract the pro-

| Outlook
Even though the models in this study use only three traits and do not include complex mechanistic feedback, they were able to pre-  (Conradi et al., 2020).
For example, the biome "dry deciduous forest" in the Indian subcontinent would more accurately fit the "mesic deciduous savanna" classification based on vegetation structure and trait composition (Ratnam et al., 2019). This confusion may have affected our results, explaining why, for example, our models predict tropical grassland in the Indian subcontinent instead of tropical dry forest (Figure 4).
Besides improvements on the biome classification side, improving trait data may also increase the accuracy of our models' predictions. For example, the extrapolation of wood density as shown in Appendix S1 may be an additional cause for the confusion between tundra and boreal forest when training the GMM, as trait predictions for these areas with similar, extrapolated climatic conditions may have resulted in largely similar traits. Additional data for the three traits included in the present study would help to reduce uncertainty of the used trait-environment relationships and to project trait changes more accurately (Figure 4c,e,g and Appendixes S2 and S9). Additionally, we are aware of the strong correlation between the extent and the predictability of biomes, where biomes with a low extent and similar traits are confused with biomes that have similar traits but a larger extent. Considering more traits may certainly improve the discriminative ability of the biome classifying model, even for those biomes with a lower extent. Good candidates for inclusion would be traits that are known to differ between biomes (e.g. root traits; Guerrero-Ramírez et al., 2020), respond to environmental factors that are expected to change into the future (e.g. hydraulic traits; Griffin-Nolan et al., 2018) and traits that are important for ecosystem functioning (e.g. diaspore size; Díaz et al., 2016). The inclusion of traits with distinct values for specific vegetation types may even allow for more refined, regional vegetation projections.
Our modelling framework can be revisited and updated as more trait data and more accurate maps of observed biomes become available. While the prediction of biome shifts remains challenging and uncertain, our findings provide an independent, complementary line of evidence for large-scale biome changes due to climate change compared to DGVMs. Our statistical modelling approach can also highlight mismatches with process-based models warranting further investigation on hypothesized mechanisms. As an example, our correlative model cannot predict vegetation that is highly dependent on factors other than climate (e.g. fire) due to the lack of suitable predictors for future scenarios (e.g. fire intensity).
The differences found between our predictions of tropical moist forests and grasslands and the ones from DGVMs that do include fire can contribute to the debate between the existence of alternative stable states for tropical forest and savanna where some parties argue that the different biomes only exist due to the presence of fire (Lasslop et al., 2020) and others argue that other factors like soil determine the distribution of the two biomes (Lloyd & Veenendaal, 2016).
Ultimately, our correlative models may have the greatest po- where results of one may complement those of the other which contributes to understanding mechanisms and output. Our study might help to develop more accurate models, which can provide an independent approach for large-scale biodiversity assessments, opening up possibilities to map and project the provisioning of ecosystem services.

ACK N OWLED G EM ENTS
Funding was provided by the ERC project CoG SIZE 647224. A. B-L was supported by a Juan de la Cierva-Incorporación grant (IJCI-2017-31419) from the Spanish Ministry of Science, Innovation and Universities.

CO N FLI C T O F I NTE R E S T
All authors had no conflict of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13431.

DATA AVA I L A B I L I T Y S TAT E M E N T
The community mean trait data that support the findings of this study are available through the publicly accessible DANS EASY repository (https://doi.org/10.17026/ dans-xf5-qdxd). In addition, we provide the R scripts for the GMM, the simulation of biome predictions, the global biome maps and the biome probability maps (https://doi. org/10.17026/ dans-xf5-qdxd) in the same DANS EASY repository.