Frontiers of Biogeography

Climatic of Sphagnum species Abstract Peat mosses (genus Sphagnum ) dominate most Northern mires and show distinct distributional limits in Europe despite having efficient dispersal and few dispersal barriers. This pattern indicates that Sphagnum species distributions are strongly linked to climate. Sphagnum -dominated mires have been the largest terrestrial carbon sinks in Europe over the last few millennia. Understanding the climatic drivers of Sphagnum species distributions is important for predicting the future functionality of peatlands. We used MaxEnt, with biologically relevant climatic variables, to model and clarify the current distributions of 45 Sphagnum species in Europe. We used a dataset of 238 316 records from across Europe (30° to 90° N, -30° to 63° E; Sahara to the Arctic, Azores to Ural mountains). We used annual degree-days, annual water balance and their monthly standard deviations (i.e. seasonality) as climatic predictors over a range of spatial resolutions (from 10 to 200 km pixel size). With these climatic predictors, we produced reasonably accurate projections of the distribution of 45 species (overall AUC >0.8). Large pixels (100 and 200 km) resulted in loss of detail, but smaller pixels (10-50 km) performed well across fit measurements. Projected distributions at the 50 × 50 km resolution showed the largest resemblance to published distribution maps. Suitable climate for many Sphagnum species was associated with the northern, western and mountainous parts of Europe. We found that annual water balance was an important indicator of Sphagnum presence. Limits in relation to annual water balance were the same as reported by bioclimatic peatland models from North America. Most Sphagnum species were limited to annual degree-days between -5000 °C y -1 and 5000 °C y -1 . Seasonality in both climate variables separated species, with degree-day seasonality having a stronger influence than water balance seasonality. High degree-day seasonality as a consequence of cold temperature sets a northern distribution limit to some species. The results suggest that the future of Sphagnum diversity in Europe is most strongly dependent on changes in water availability and in seasonal temperature variation. degree-days, GBIF data, peat mosses, peatlands, precipitation, species distribution models


Introduction
Understanding the distribution of peat mosses (Sphagnum species) is important as they are the main drivers behind the carbon accumulation in northern peatlands over millennia (Clymo and Hayward 1982). In Europe there is a wide diversity of mire types circumscribed by climatic zones. Shallowpeated tundra mires occur in northern arctic and alpine areas. Boreal bogs, fens and huge aapa mire complexes cover large parts of the landscape from Fennoscandia eastwards. More locally in the boreal and temperate zones, fens associated with groundwater discharge can host a high species richness of Sphagnum (but with less dominance). Maritime bogs occur in coastal areas in the northwest, particularly in Norway, Scotland and Ireland. Further south in Europe, mires become small and scattered and quite variable. In the Mediterranean and semi-arid areas with summer drought, peatlands are confined to depressions with standing water all year. While Sphagnum species are most prominent in peatlands (see Joosten et al. 2017 for a detailed European overview), they are also common in oceanic wet heaths and in swamp forests. In Europe there are 58 species (Hodgetts 2015, Séneca andSöderström 2009). Even though Sphagnum species have effective long-distance dispersal (Sundberg 2010(Sundberg , 2013, they vary in their distribution in Europe from south to north and from oceanic to continental areas. Experiments have indicated that northern species of Sphagnum tolerate cold autumns and winters better than species with more southern distribution (Campbell and Rydin 2019). There is also a large variation in drought resistance and growth response to temperature among Sphagnum species (Breeuwer et al. 2008, Rydin andJeglum 2013). Hence, we expect that the biogeography of Sphagnum can be predicted by physiologically relevant aspects of the climate (cf. Woodward and Williams 1987).
Recent distribution models for bryophytes include both country-wide and continental-wide scales (Mateo et al. 2016, Mateo et al. 2013, Oke and Hager 2017, Sergio et al. 2007). Oke and Hager (2017) modelled Sphagnum distribution to predict and project peatland distribution over North America. They showed that Sphagnum and peatland presence have a climatic limit where the difference between precipitation and potential evapotranspiration is -500 mm y -1 . Popov (2016Popov ( , 2018) modelled species distribution of two sub-genera of Sphagnum in European Russia and Finland and found that most Sphagnum species were strongly positively associated with high humidity or precipitation and negatively with high summer temperature. This pattern is expected, as Sphagnum lacks the ability to regulate water loss and depends on water being freely available Climatic drivers of Sphagnum species distributions 4 (Rydin and Jeglum 2013). The fact that species composition differs between oceanic and continental areas (Daniels and Eddy 1990, Gignac and Vitt 1990, Flatberg 2013 indicates that the biologically relevant features of the climate are not only annual temperature and water balance, but also their seasonal variation. However, in focusing on north-east Europe, Popov (2016Popov ( , 2018 did not cover the main area of Sphagnum diversity, namely the highly oceanic Atlantic coast. Species distribution modelling (SDM) uses species observations and a set of environmental predictors to produce environmental suitability maps (Guillera-Arroita 2017). SDM is a useful tool to assess ecological preferences and tolerances of species to better understand their biogeography.
The usefulness of a model for a particular question is dependent not only on the model inputs, but also on the extent of the investigated area, model resolution and the assessment of model fit (Guillera-Arroita 2017, Leroy et al. 2018, Seo et al. 2009, Yackulic et al. 2013. In general, continental-scale models rely on databases with presence-only species data (e.g., Mateo et al. 2013), and maximum entropy modelling (MaxEnt) is widely used (Elith et al. 2011). Within the MaxEnt framework, there are a large set of tools to handle variation in survey intensity and quality across countries (Leroy et al. 2018). Although these advances have improved the use of mixed data sources, there are still challenges with SDM in general. For example, inappropriate variable selection, combined with incautious interpretation of evaluation statistics, is known to produce false patterns (Fourcade et al. 2018). Biologically relevant predictors at an appropriate scale are likely better than a large set of automatically produced variables (e.g., the bioclimatic variables from WorldClim; Hijmans et al. 2005). Furthermore, the resolution (pixel size) is known to affect the results and quality of species distribution and species richness models. Low-resolution projections can over-estimate species ranges (Seo et al. 2009), and species richness predictions are scaledependent (Pineda and Lobo 2012). The effects of resolution are particularly relevant when modelling multiple species with a variety of range sizes; such as the Sphagnum genus.
We used MaxEnt with relevant climatic predictors (annual water balance and degree-days and their seasonal variability) to model the distributions of the Sphagnum species of Europe. Our aims were: (1) to test the relationship between climatic variables and Sphagnum species distributions, (2) to produce climate suitability maps for each Sphagnum species and (3) to map the association between climatic variables and modelled Sphagnum species richness. We performed an array of models at different spatial resolutions to explore model sensitivity, in terms of predictive performance, to changes in pixel size.

Geographical extent
We used an extent of 30 to 90° latitude and -30 to 63° longitude, of which geographical Europe comprises the main part. This region is bounded by geographical barriers that limit Sphagnum distributions (Kyrkjeeide et al. 2016), viz oceans (Arctic and Atlantic) and deserts (northern Africa, the Arabian Peninsula and western and central Asia). Sphagnum dispersal is limited by the Gulf stream to the west and the Ural mountains to the east. South to north the region climatically spans the arid to polar division of the Köppen climate classification (Peel et al. 2007). The western parts are oceanic and strongly influenced by the Gulf stream. There is consequently a large temperature gradient from north to south, which is tempered by a west-east continentality gradient. The interaction of the two affects the seasonality of temperature and precipitation.

Records data and taxonomy
Sphagnum species record data were downloaded from the database of the moss flora of Russia 1 (Ivanov et al. 2017), Biocase 2 , GBIF (Occdownload Gbif.Org 2017) and the Finnish Biodiversity Information Facility 3 . We used only data from 1980 onwards. This is a period of extended and extensive bryological and mire related recording in several countries. For example, the first bryophyte atlas in Britain and Ireland (Hill et al. 1991), the Wetland Inventory in Sweden (Gunnarsson et al. 2014) and the Inventory of Swiss Mire Landscapes of Particular Beauty and National Importance (Hammer et al. 2009). For records without geographical coordinates, we looked for coordinate data in the locality information and used the coordinates of the centre of the locality. Where we only had the name of the locality, we georeferenced the locality to within 5 km using online resources (Geo 4 , Google maps 5 and Yandex 6 ). These databases were further supplemented with published records from the literature (Appendix S1).
Nomenclature followed Hodgetts (2015) with the exception of S. beothuk (Andrus) which had previously not been recognised in Europe. Synonymous taxa were found using The Plant List (2010). All doubtful taxa were removed. There are three cases in which the taxonomy has shifted somewhat over time.
(1) Sphagnum beothuk was recently split from S. fuscum in Europe and has only been recorded in Britain, Ireland, Norway and Sweden; therefore, a few records listed as S. fuscum in our data from other regions may in fact be S. beothuk. Hill (2017)

Climate data
We used four climatic variables: annual and seasonal degree-days and water balance. These variables were chosen based on Sphagnum ecology. To calculate the climatic variables for the geographical extent, we downloaded raster layers of mean monthly maximum and minimum temperature and precipitation from Worldclim (v. 1.4) (Hijmans et al. 2005) at a resolution of 30 arc seconds (~1 km 2 at the equator) and monthly means of potential evapotranspiration (PET) from Trabucco and Zomer (2009). These layers were re-projected to an Albers Equal Area grid with centre point 60° N; 16.5° E. The re-projected variables were used to calculate degree-days and water balance. Degree-days are normally based on daily data and a set base temperature. We did not have access to data with such high temporal resolution and instead define degree-days operationally as the mean of monthly max and min temperature multiplied by the number of days in the month.
This measure combines temperature with time, indicating which areas are warmer or colder for longer or shorter durations. We did not use the base temperature, as its inclusion is related to phenological signals (e.g., to initiate growth in the spring) which is not relevant for Sphagnum.
Water balance was calculated as precipitation -PET for each month. The monthly values of both variables were summed to give annual degree-days and annual water balance, respectively. To represent the seasonal variability, we used the standard deviation of the monthly data as Degree-7 Day Seasonality and Water Balance Seasonality. The climate data were then aggregated to pixels of 10 × 10, 25 × 25, 50 × 50, 100 × 100 and 200 × 200 km by calculating the mean values for the aggregated pixels.

MaxEnt models
Species distribution models were constructed using MaxEnt v. 3.4.1 (Phillips et al. 2017) in the dismo v. 1.1-4 R package (Hijmans et al. 2017). MaxEnt is an algorithm that minimises the distance between the probability densities of environmental predictors taken from species occurrences and background points (Elith et al. 2011). This minimal solution is the relative probability that a grid cell is contained within a collection of presence points and is equivalent to the relative occurrence rate of a species (Merow et al. 2013). It has been shown to be robust to variation in sample size and, in comparison with other modelling frameworks, has one of the best predictive powers for presence/background data (Wisz et al. 2008). We used MaxEnt's raw output for all calculations following Merow et al. (2013).
Compared to classical statistical approaches (standard GLMs) MaxEnt can better handle collinearity and find the best set of parameters. This is particularly true when predictive accuracy is in focus (De Marco and Nóbrega 2018). Correlations of the climate layers in our study at each resolution were generally weak or moderate (r < 0.6), except for annual degree-days vs seasonal water balance (200 × 200 km pixels, r = 0.78; 10 × 10 km pixels, r = 0.74). Taken together, collinearity is unlikely a major issue in our analyses.
Despite its strengths, uncritical use of MaxEnt has been shown to be problematic. Underlying record data bias and overfitting of environmental predictors as well as their spatial resolution have all been shown to affect SDM quality in ways that may not be detected quantitatively (Beck et al. 2014, Merow et al. 2014, Pineda and Lobo 2012. We used a number of techniques to counteract such effects, which are described below.

Model fitting procedure
We converted the cleaned data records to spatial grids with the same projection and resolutions as the climate rasters. For each species, we retained 10% of the presence points per resolution as a test data set. We split the species into a high presence group (≥30 records) and a low presence group (20 Climatic drivers of Sphagnum species distributions 8 to 29). No model was produced for species with fewer than 20 presences since sample sizes below 20 have been shown to perform poorly (Soultan and Soufi 2017). For each species in each group, we partitioned the data in two different ways appropriate for the sample size, producing two models. For each species in the high presence group, we partitioned the presence data into 10 randomly assigned groups (k-folds), and we also partitioned these data into four spatially explicit blocks. Spatially explicit partitioning has been shown to improve model projections for biased data (Radosavljevic andAnderson 2014, Roberts et al. 2016). Presences were assigned to each block using the ENMeval v. 0.3.0 R package (Muscarella et al. 2014). For species in the low presence group, we used n-1 jackknifing, following Shcheglovitova and Anderson (2013), and a model with no partitioning. For species with low prevalences jackknifing has been shown to improve model performance.
As pseudo-absences we randomly selected 10 000 background pixels. When fewer than 10 000 pixels were available, we used the whole environmental background. We used the 'all features' setting of MaxEnt and set the β-regularisation to 2. Doubling the regularisation reduces model complexity and likelihood of overfitting, reducing the effects of spatial bias (Cao et al. 2013, Radosavljevic and. Quantitatively comparing performance of MaxEnt models built using background data is complex, and multiple metrics have been recommended (Mouton et al. 2010). Lawson et al. (2013) recommend using continuous evaluation metrics and characterise prevalence-dependant and independent metrics as favouring calibration and discrimination respectively. Discrimination is the ability to separate presences and absences based on model outputs, and calibration the numerical match between observed and predicted probabilities.
Here we used the area under the curve (AUC) of the receiver operator curves (ROC), log-loss and omission rate for the test and training data sets. AUC is a measure of the probability that a presence is ranked higher than an absence for a binary classifier and is prevalence independent. Log-loss is a measure of the difference between two probability distributions and is equal to the negative loglikelihood, and is prevalence dependent. Omission rate is the proportion of presences misclassified as being absences (false negatives; Merow et al. 2013). Evaluation statistics were calculated from the average values from the partitioned models. For each species at each resolution, we selected the model which had the most favourable evaluation statistics across both test and training datasets.
We used MaxEnt´s raw output for all model evaluations and further analyses (Merow et al. 2013).
A qualitative assessment of the model was done by visually comparing the rasterised outputs with published species range maps from Daniels and Eddy (1990) to check that the projected distributions were reasonable.

Species distributions
We examined the relationship between species distributions and climate at the finest resolution (10 × 10 km). The relative importance of the predictors in the final model is presented as "permutation importance", which is determined by randomly permuting the values of a predictor among the training points and calculates the decrease in training AUC (Philips and Dudik 2008). The larger decrease, the greater the influence of the predictor on the model. The output is normalized so that the sum of all predictor importance is 100 percent. We illustrated the relationship of the species distributions to the climate factors by first conducting a principal component analysis (PCA) on the rasterised raw outputs. Thereafter we fitted thin-plate GAM splines to the PCA using the ordisurf function from the vegan v. 2.5-4 R package (Oksanen et al. 2019

Projected species richness
To project species richness from model output, we converted the rasterised climate suitability maps to presence/absence maps. To determine if a species was projected to be present in a pixel, we used appropriate threshold values following recommendations of Liu et al. (2013). For the high presence group, we used the maximum sensitivity plus specificity calculated using the training presence points; for the low presence group, we used the mean predicted suitability of up to 10 000 random points when available. Projected species richness for each pixel was then calculated as the sum of presences from a raster stack where each raster is a presence/absence map of a species. There was a strong north-south gradient in annual degree-days, while degree-day seasonality had a stronger east-west gradient (Fig. 1). Annual water balance had the highest values in western Norway and the Alps. Water balance seasonality was highest in the south and east, likely as a consequence of higher rates of potential evapotranspiration in the summer.

Model evaluations
Most species had enough presences to be assigned to the high presence group (Table 1) The quantitative evaluation showed that models performed well across all species (Table 1). Both training and test AUC were generally above 0.8 for most of the models (Table 1, Fig. S1). AUC did decrease with increased pixel size and the lowest AUCs were found at the 200 × 200 km resolution.
Only S. girgensohnii and S. squarrosum had test AUC < 0.8. Test log-loss decreased with pixel size for most species, with the lowest values again at the 200 × 200 km resolution. However, training log loss was lowest at 50 × 50 km resolution. Omission error varied little across models but was slightly lower at the 25 and 50 km resolutions. At the coarsest resolutions (100 and 200 km), the projected areas lost details (Fig. 2), and the species responded less to the seasonality variables (flatter curves; Fig. 3). Given that, the results are fairly robust to change in resolution between 10 and 50 km (Appendix S2), we present climatic suitability maps for the finest resolution, 10 × 10 km.        (Fig. 4). Figure 6. Projected species richness at pixel size 10 × 10 km, calculated from thresholded mean raw outputs of the selected MaxEnt models. Threshold values were the maximum sensitivity plus specificity calculated from the presence points used to create the models.

Projected species richness
Species richness as predicted by climate suitability is strongly concentrated towards the north and west (Fig. 6). Most strongly it is centred on the Fennoscandian region, particularly Norway. To the south and east, a higher species richness is associated with mountainous regions, such as the Dinarides of the Balkans. Around the coasts of Africa, the Black Sea and the Mediterranean region, suitable climate is patchy.  Marginal response curve outputs (Fig. 8) show a similar response for most species regarding annual degree-days, with the highest probability of presence around 0 degree-days. For most species, the probability of presence in relation to annual water balance rapidly rises in the range -500 and 500 mm y -1 . Species responses vary markedly above 500 mm. Three species have high marginal responses even at low water balance values: S. angermanicum, S. annulatum and S. rubiginosum.
These species are typically northern and scattered in occurrence. In relation to degree-day seasonality, most species had ranges which fell between variabilities of 100 °C y -1 to approximately 350 °C y -1 , but there appear to be groups of species with different peaks. As the response value of one group decreases, another group of species increases. Water balance seasonality had little influence on the models (Table 2). It should be noted that this predictor was correlated with annual degree-days, which showed high contribution to the models. These predictors cannot be completely separated, but annual degree-days gives a better model fit than water balance seasonality.
Regardless of the low contribution of water balance seasonality to the model, there was a general pattern of increasing probability of presence from 0 mm y -1 to a general peak at around 40 mm y -1 .
After this point, most species probabilities remain high or, in some cases, decline slightly. Several species reach peaks above 40 mm y -1 up to approximately 100 mm y -1 .

Distribution maps of suitable climate for Sphagnum
Using only four climate variables, we produced projections that agree reasonably well with known Sphagnum species distributions (Daniels and Eddy 1990, Séneca and Söderström 2009, Masing et al. 2010). Our selected variables are known to be closely linked to the biology of plants in general and Sphagnum in particular (Franklin 1995, Gignac et al. 2000, Rydin and Jeglum 2013 Hager 2017). The use of a few, mechanistically based predictors should avoid problems of overfitting and low biological realism (Fourcade et al. 2018, Merow et al. 2014).
Structured partitioning of presences has been shown to improve model predictions in structured data (Radosavljevic and Anderson 2014). We found that the structured partitioning models (block partitioning) were selected most frequently at finer resolutions for the high presence group, likely a consequence of western Europe being more densely recorded than eastern Europe. Such structural bias is further evidenced by species associated with western, oceanic regions having better evaluation metrics for the 10 k-fold method (random partitioning).
Spatial resolution is known to affect MaxEnt predictions (Seo et al. 2009). The lowest resolutions (pixels > 50 x 50 km) resulted in loss of details and were uninformative for projection at a relevant landscape scale. They also led to overestimated ranges of the species (cf. Seo et al. 2009) compared to published distributions (Daniels and Eddy 1990, Séneca and Söderström 2009, Masing et al. 2010. AUC values remained high across all resolutions (Table 1) and for other statistics (test omission, test logloss) the differences between resolutions from 10 to 50 km were subtle. Within this range the results seem robust to choice of pixel size.
When comparing the projected species distributions with the literature (Daniels and Eddy 1990), the impression was that the 50 × 50 km resolution gave the result that was most similar to published maps across all species, and for some applications this could be a reasonable resolution. However, discrepancies between results from models solely based on climatic predictors and observed species distributions are to be expected as the models do not account for biological interactions or dispersal limitations (Pineda and Lobo 2012). The good agreement (high AUC, low omission rates and a reasonable replication of Daniels and Eddy (1990) distribution maps) in our results is aided by the fact that Sphagnum are little affected by interactions such as disease and herbivory (Rydin and Jeglum 2013), and by the high dispersal ability of most Sphagnum species (Mikulášková et al. 2015, Sundberg 2013. Hence, our results are in line with Kyrkjeeide et al. (2016), who suggested that spore-producing Sphagnum species are able to fully occupy their climate niche.

The relationship between climate and Sphagnum species distribution
Sphagnum species are generally favoured by a wet climate (Rydin and Jeglum 2013), as indicated by the importance of annual water balance as predictor for Sphagnum presence. The surprising result that at such a low value as -500 mm yr -1 the climate starts to become suitable for most species is probably an effect of the extraordinary capacity of Sphagnum to store water in extracellular spaces and in hyaline cells (Rydin and Jeglum 2013). Earlier studies on climatic requirements for peatlands in Canada found similar values (Gignac et al. 2000, Oke andHager 2017) suggesting that there is a limit for Sphagnum species and for the formation of high latitude peatlands in general.
This limit is related to the physiological constraints of Sphagnum: their photosynthesis is strongly dependent on tissue water content, and even though they can store water they are unable to control water losses (Rydin and Jeglum 2013). A suitably wet climate is achieved by high precipitation (as in oceanic regions) or by low temperature that reduces evaporation (in the north or at high altitude).
However, Sphagnum can grow in areas with drier macroclimate by microclimatic and hydrological buffering, which allows non-peatland species to occupy, for example, shaded swamp forests and fens with groundwater discharge in relatively dry climates (Joosten et al. 2017). The seasonality of water balance contributed little to our models. The weaker influence of fluctuation in water balance may be due to its unequal effects on Sphagnum presence across our study area. Water balance seasonality may be similarly high in some boreal and Mediterranean areas (Fig. 1). However, in the north the seasonality is largely an effect of variation in precipitation, whereas in the Mediterranean, it is an effect of variation in PET that leads to long dry periods that adversely affect Sphagnum physiology and ultimately its survival (Gerdol and Vicentini 2011). Species with a western distribution are associated with low degree-day seasonality (Fig. 7), as expected in an oceanic climate, but somewhat surprisingly not with low water balance seasonality. This indicates that oceanic climates are best captured by an even temperature and a high, but seasonally fluctuating precipitation.
As for temperature, our models suggest that suitable climate for Sphagnum falls within a range between -5000 and 5000 °C y -1 , but with rather narrow peaks around 0 °C y -1 . The amount of warmth in the environment on an annual basis is biologically important as there are physiological limits to growth and survival (Franklin 1995). Photosynthesis of Sphagnum begins to be impacted by temperature above 35 °C (Haraguchi and Yamada 2011), but provided that the mosses remain moist, they can survive in places with high air-temperatures owing to the cooling effects of evaporation (Rydin 1984, Van der Molen and Wijmstra 1994, Dyukarev et al. 2009), though there are exceptions (Lange 1973). High temperatures are most likely in regions where water is scarce and evaporative demand is high. Such areas are unlikely to support Sphagnum due to a lack of freely available water, not because of temperature limitations. Species with a northern distribution are associated with lower annual degree-days and also with higher degree-day seasonality (Fig. 7).
In the north, seasonal variation is caused by low temperatures during winter and late autumn (before snow cover gives protection). Comparing species occupying similar micro-habitats, northern species in general tolerate such conditions better than their southern counterparts (Campbell and Rydin 2019).
Our models suggest that reduced wetness would be the most detrimental climatic change, probably for all species of Sphagnum. It is predicted that water availability will decrease substantially in southern Europe with climate change (e.g., Ruosteenoja et al. 2018), which potentially can limit Sphagnum species southern distribution. In addition to the direct effect on growth, a drier climate can lead to increased cover of vascular plants in peatlands (Rydin and Jeglum 2013), which will further hamper Sphagnum performance (Berendse et al. 2001, Bengtsson et al. 2021

Patterns of Sphagnum species richness
Our models show that both projected and observed species richness in Sphagnum is highest in European north and west, while the number of species decreases southwards and eastwards. This pattern agrees well with studies examining Sphagnum species richness within eastern Europe (Popov 2016(Popov , 2018. Most Sphagnum species have a suitable climate in the boreal zone where temperature is not too high. It is not unlikely that this is a legacy from the Sphagnum species radiation that occurred in northern boreal ecosystems (Shaw et al. 2019), but in addition, a positive response to high annual water balance leads to increasing suitability for most species towards the Atlantic.
Despite our attempts to reduce sampling bias by collecting data not only from GBIF, but also from primary publications (see Materials and Methods), the richness maps may be affected to some extent by variation in sampling efforts. North-western Europe comprises species-rich areas that are intensively sampled (e.g., the British Isles, most of Fennoscandia and the Netherlands). But there are also areas with low predicted (and real) species richness that are intensively sampled (e.g., eastern England). Interestingly, a species richness cold-spot is predicted at the northern borders of Finland, Norway and Sweden (Fig. 6). The area does not stand out in the climatic maps (Fig. 1). It is a remote area with differing record intensities from the three countries, but field studies would be needed to say with certainty whether the low species richness is an effect of a peculiar climate or an effect of recording bias.
Whilst Sphagnum distributions are closely allied to climate variables, the genus can be present in areas decoupled from its predicted climate niche, owing to favourable local topography and hydrology (Cerrejón et al. 2020). To the south and east in Europe, projected Sphagnum species richness is high in mountain regions even though they are not rich in peatlands (Joosten et al. 2017 (2013) report 14 and 21 species for these regions, respectively. Our projections indicate Sphagnum presence but not particularly high richness (Fig. 6). We also project low climatic suitability for Sphagnum along Georgia's Black Sea coast despite this being an area of high Sphagnum productivity (Krebs et al. 2016).

Conclusions
Sphagnum species distributions can be well understood based on a few climatic variables linked to their physiology. Sphagnum presence is to a large extent governed by water availability, leading to higher species richness to the north and west in Europe, with mountainous regions providing suitable climate for several species to the south and east. Sphagnum species distributions are also related to climate seasonality, for example as a result of different tolerance to low autumn and winter temperature. Agreement between our models and bioclimatic models of Sphagnumdominated peatlands implies changes in humidity may have the strongest impact on both Sphagnum diversity and peatland function. The coldest regions of Europe are likely under-sampled and there are too few records of several newly described species, calling for directed sampling campaigns to further improve the understanding of Sphagnum biogeography.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1. Line graph of the evaluation parameters of AUC and Log-loss for all species and methods across resolution sizes. Table S1. Excel table of model evaluations for all produced models.
Appendix S1. Literature sources for supplemental records.