The mossy north: an inverse latitudinal diversity gradient in European bryophytes

It remains hotly debated whether latitudinal diversity gradients are common across taxonomic groups and whether a single mechanism can explain such gradients. Investigating species richness (SR) patterns of European land plants, we determine whether SR increases with decreasing latitude, as predicted by theory, and whether the assembly mechanisms differ among taxonomic groups. SR increases towards the south in spermatophytes, but towards the north in ferns and bryophytes. SR patterns in spermatophytes are consistent with their patterns of beta diversity, with high levels of nestedness and turnover in the north and in the south, respectively, indicating species exclusion towards the north and increased opportunities for speciation in the south. Liverworts exhibit the highest levels of nestedness, suggesting that they represent the most sensitive group to the impact of past climate change. Nevertheless, although the extent of liverwort species turnover in the south is substantially and significantly lower than in spermatophytes, liverworts share with the latter a higher nestedness in the north and a higher turn-over in the south, in contrast to mosses and ferns. The extent to which the similarity in the patterns displayed by spermatophytes and liverworts reflects a similar assembly mechanism remains, however, to be demonstrated.


Table of contents
Atlas of Flora Europeae provides a good representation of vascular plant diversity patterns in western Europe vascular plants, but not in eastern Europe 5-6 . Therefore, eastern MGRS pixels were discarded in following analyses.

Environmental predictors
As environmental predictors we used the 35 macroclimatic variables of CliMond (https://www.climond.org/) 7 , as well as monthly and annual potential evapotranspiration (http://www.cgiar-csi.org, see ref. 8). To avoid multicollinearity, we run a Pearson correlation analysis eliminating one of the variables in each pair with a correlation value higher than 0.8. The final set of six variables used to run the models were 'mean diurnal temperature range', 'temperature seasonality', 'precipitation seasonality', 'mean moisture index of warmest quarter', 'mean moisture index of wettest quarter' and 'annual potential evapotranspiration'.

Background selection and statistical modelling
For each species, we generated 10 sets of pseudo-absences equalling the number of presences and sampled from pixels not adjacent to reported occurrences, that were later used to produce an ensemble model 7 using three different techniques: generalized linear models (GLM) 10 , Maxent 11 , and Random Forests (RF) 12 , as implemented in the R 13 package BIOMOD 2.0 14 . The performance of the models was assessed by randomly splitting 10 times the data into a 70% dataset to generate the models and a 30% dataset to estimate their predictive accuracy. After elimination of all models with an AUC<0.8, we generated for each species an ensemble model, consisting in a weighted mean of the models predictions, where the contribution of each individual technique was proportional to its predictive accuracy.

Evaluation statistics and binarization of species' potential distribution
Because different ensemble models could generate different models 15 , we generated two ensemble models per species, including either (1) only models with an AUC>0.8 (ROC consensus model), or (2) models with a true skill statistic (TSS) > 0.7 (TSS consensus model). The contribution of each model to the final ensemble model was proportional to their goodness-of-fit statistics.
If stacking of the binary models reduces the over-prediction 2 , the selection of an appropriate threshold still can reduce error rates in both individual and ensemble SDMs 16 . To take into account the effects of the selected threshold used to generate binary models on species' over-prediction, we produce three different binary models For the TSS consensus model, we also generated three richness models: (1) optimizing the TSS statistics; (2) applying a maximum of 5 % of omission error; and (3) applying a maximum of 10 % of omission error.
We stacked (i.e., summed) the different binarised SDMs predictions so that we obtained 6 stacked species distribution models (S-SDMs) per species group (bryophytes and vascular plants), which represent the potential species richness for both groups (Figs. S1 and S2). To compare differences between models, we calculated the Pearson correlation coefficients by pairs of S-SDMs (Tables S1 & S2).
As correlations were very high, we discarded possible differences between S-SDMs, so in the following analyses we employed only the S-SDMs obtained by optimizing the ROC statistics for all the subsequent analysis.

Over-fit testing for bryophytes
By definition, rare species distributions usually have few observations, and are prone to over-fitting when modelled. To ensure that models were not over-fitted for large set of predictors, we employed the 'ensembling of bivariate models' approach presented by Lomba et al. 18 . It involves generating bivariate models (using only two independent variables every time) and averaging them with a weighted ensemble approach. We have run this approach for bryophytes and climatic variables to compare the results obtained and discard over-fitting problems. As in the original S-SDMs we generated 300 models for every species; in order to be comparable we have run 480 bivariate models per species (4 pseudo-absence sets, 4 replicates, 10 groups of pairs of independent variables, and 3 modelling techniques). We generated a consensus (only models with an AUC>0.8) and a binary model per species (optimizing the ROC statistic). The final richness model generated was compared with the original richness model to discard over-fitting. The Pearson correlation coefficient between both approaches was 0.99, hence we consider that over-fitting problems can be disregarded in the modelling approach followed.    For the first comparison we produced a dataset of observed species richness for 45 UTM squares (100x100 km) scattered across Europe and intensively sampled according to a previous review of the literature ( Figure S3, see Appendix S5 for the complete list of references).
Then we generated a set of variables describing extant macroclimatic conditions, environmental heterogeneity, spatial patterns and historical factors across Europe (see Table 1 of the main manuscript). We selected the variables that best describe patterns of bryophyte species richness across the 45 UTM squares through a multimodel inference approach 5 approach. We employed the MuMin R package 6  consequently projected onto Europe to predict bryophyte species richness across the continent ( Figure S3).
On the other hand, we downloaded all the collections available for bryophytes in the Global Biodiversity Information Facility (GBIF) database across Europe. We counted the number of collections per UTM squares (100x100 km) to generate a map of sampling effort for bryophytes in Europe ( Figure S3).
Finally we calculated the Pearson correlation coefficients between the three maps generated on the 45 pixels with observed richness values from bibliography. Those values were highly correlated with the potential richness obtained by S-SDMs and MEM, and very low correlated with sampling effort. In conclusion, the potential bryophyte richness (S-SDMs) map is not biased by sampling effort, and it is a good tool to study biodiversity patterns for bryophytes in Europe.
We also compare the S-SDMs obtained for vascular plants and observed richness patterns. We generated a map of observed richness patterns ( Figure S4) for vascular plants with all the data available in the study 9 .

Supplementary Methods 4: Comparison of potential richness maps
Potential richness (stacked species distribution models, S-SDMs) maps were employed to compare the spatial patterns of species richness among mosses, liverworts, ferns and spermatophytes. The first step was two couple pixel by pixel richness patterns analyses.
First, we calculated and mapped local Lee's L bivariate spatial association 1 using our own implementation of this statistics with the R language. This function is now included in the 'spdep' package 2 . Lee's L statistics is a measure of the pixel-to-pixel correlation between two spatial variables, corrected for the spatial autocorrelation of To complement this analysis, we mapped the residuals of a Procrustes 4 analysis with the "vegan" R package. Procrustes rotation allows the comparison of the ordinations of two data matrices through an algorithm that minimises the sum of squared distances between corresponding points belonging to the two matrices.
Procrustean comparison was applied to two canonical correspondence analyses (CCAs) 5 performed with the ade4 R-package 6 and conducted using the matrices of spermatophytes, mosses, liverworts and ferns species probabilities in each 100 km pixel against a set of environmental variables (Table 1)  Inertia of spermatophytes, ferns, mosses, and liverworts composition explained by environmental variables was 65.05%, 68.67%, 67.44%, and 66.13%, respectively. The CCA gets the most important gradients of the independent variables to explain the richness patters for each group (Fig. S6), and the Procrustes analysis (Fig. S7) compare both CCA to see if they are the same. High differences between gradients (high values of residuals) indicate that environmental variables have different influence for the taxonomical groups compared. For deriving a value of c, we assume that total species presences are the same in expected and observed SR sets, i.e. ΣES i =ΣOS i . Hence, For the estimation of the value of z, SR was calculated for groups of contiguous pixels ranging in area from ∼10 4 km 2 (1 pixel) to ∼64*10 4 km 2 (64 pixels), completely nested within them and covering almost all Europe. Pixels covering sea surfaces were excluded from this analysis. The species-area relationship was calculated via regression of SR on area (Dengler 2009 Figure S8d: Mean and SD annual precipitation (mm/m 2 ) by latitudinal band.