Global patterns and drivers of alpine plant species richness

Aim: Alpine ecosystems differ in area, macroenvironment and biogeographical history across the Earth, but the relationship between these factors and plant species richness is still unexplored. Here, we assess the global patterns of plant species richness in alpine ecosystems and their association with environmental, geographical and historical factors at regional and community scales. Location: Global. Time period: Data collected between 1923 and 2019. Major taxa studied: Vascular plants. Methods: We used a dataset representative of global alpine vegetation, consisting of 8,928 plots sampled within 26 ecoregions and six biogeographical realms, to estimate regional richness using sample- based rarefaction and extrapolation. Then, we evaluated latitudinal patterns of regional and community richness with generalized additive models. Using environmental, geographical and historical predictors from global raster layers, we modelled regional and community richness in a mixed- effect modelling framework. Results: The latitudinal pattern of regional richness peaked around the equator and at mid- latitudes, in response to current and past alpine area, isolation and the variation in soil pH among regions. At the community level, species richness peaked at mid-latitudes of the Northern Hemisphere, despite a considerable within- region variation. Community richness was related to macroclimate and historical predictors, with strong effects of other spatially structured factors. Main conclusions: In contrast to the well- known latitudinal diversity gradient, the alpine plant species richness of some temperate regions in Eurasia was comparable to that of hyperdiverse tropical ecosystems, such as the páramo. The species richness of these putative hotspot regions is explained mainly by the extent of alpine area and their glacial history, whereas community richness depends on local environmental factors. Our results highlight hotspots of species richness at mid- latitudes, indicating that the diversity of alpine plants is linked to regional idiosyncrasies and to the historical prevalence of alpine ecosystems, rather than current macroclimatic gradients.


| INTRODUC TI ON
More than 200 years after the attempt by Alexander von Humboldt to formulate a unified theory of the natural world, understanding the global patterns of diversity remains one of the greatest challenges in biogeography and macroecology (Brummitt et al., 2020;Keil & Chase, 2019;Kier et al., 2005;Kreft & Jetz, 2007;Kreft et al., 2008;Weigelt et al., 2016). In particular, mountains have been revealed as centres of biodiversity, with a disproportionately high species richness in comparison to their corresponding lowland regions (Antonelli et al., 2018;Muellner-Riehl et al., 2019; Rahbek, Borregaard, Antonelli, et al., 2019). Along the elevational gradient of mountains, the compression of life zones brings different biomes into proximity, with the alpine belt representing the outpost for plant life above the climatic tree line. Alpine ecosystems, governed by low-temperature regimes, cover c. 3% of land outside Antarctica and are distributed across all continents and latitudes Testolin et al., 2020). Despite ongoing efforts to monitor changes in the biota of mountain summits in the face of climate change Pauli et al., 2012;Steinbauer et al., 2018), we still lack a picture of the global patterns of plant diversity in alpine habitats, let alone an understanding of its major drivers.
According to the general latitudinal diversity gradient, biodiversity is expected to peak at the equator (Hillebrand, 2004).
Among the possible explanations for this pattern (Lomolino et al., 2017), latitude is normally interpreted as a proxy for climatic conditions and available metabolic energy, which might have an effect on speciation rates (Wang et al., 2009). Whether this general rule also applies to alpine ecosystems, however, is still a matter of debate. Regardless of their latitude, alpine ecosystems are determined by low-temperature conditions, hence low energy input. Therefore, lowland and alpine thermal conditions from polar to equatorial latitudes are increasingly decoupled from one another (Testolin et al., 2020). Besides having a lower energy input compared with the lowlands, alpine ecosystems are also highly heterogeneous in their topoclimates (Quinn, 2008), which might weaken the correlation between latitude and primary productivity (Testolin et al., 2020). For these reasons, plant diversity in alpine areas might decouple from major climatic gradients.
Alpine areas are also isolated from each other, forming fragmented systems of "sky islands" surrounded by lowland environments that limit species dispersal (McCormack et al., 2009).
Following the ecological principle of the species-area relationship (Lomolino, 2000b) and its application to the theory of island biogeography (MacArthur & Wilson, 1967), the extent of alpine habitats and their isolation could have affected rates of colonization, speciation and extinction of plants (Heaney, 2000;Steinbauer et al., 2016).
These processes might have resulted in biodiversity patterns linked to the historical and current abundance of alpine habitats at the global scale. Although it has been reported widely that the biogeographical history of mountains has shaped diversity patterns of cold-adapted plant species in alpine regions Harris, 2007;McGlone et al., 2001;Sklenář et al., 2014), a major unresolved question is the extent to which the interplay of ecological drivers and historical contingencies dictates the patterns of alpine plant diversity at the global level (Nagy & Grabherr, 2009). The significance of these drivers might shift from global to local spatial scales and can reveal new patterns and relationships that are not evident at regional scales at which alpine plant diversity patterns have been studied so far (Jiménez-Alfaro et al., 2014;Lenoir et al., 2010;Moser et al., 2005;Vonlanthen et al., 2006).
Here, we compiled a dataset of 8,928 vegetation plots with 5,325 vascular plant species sampled by botanical experts in alpine ecosystems over the past 100 years and representative of global alpine vegetation. By analysing the data at both regional and community levels, we investigate: (a) the global latitudinal patterns of alpine plant species richness; and (b) the relative influence of environmental, geographical and historical factors in driving such patterns. We also evaluate how those patterns and drivers change between regional and community levels and how they relate to hotspots of alpine plant diversity recognized at the global scale.

| Study system and data collection
We considered as zonal alpine vegetation any plant community dominated by graminoids, forbs and dwarf shrubs above the climatic tree line (Körner, 2003). In addition to strictly zonal habitats, snow-patch vascular plant communities and vegetation on rocks and screes are also found ubiquitously in the alpine belt and were included in our study. We did not consider vegetation from polar climates owing to the absence of elevational tree lines and their distinct environments (Quinn, 2008;Walter & Box, 1976). Therefore, the alpine vegetation included in the present study corresponds to the "mid-latitude alpine tundra" and the "tropical alpine biome" groups as defined by Quinn  (Supporting Information Table S1). Datasets from different sources were standardized by identifying a minimum common set of plot attributes, including size, elevation and geographical coordinates. When the geographical coordinates were missing for small, clearly delimited areas, we estimated plot locations from maps (i.e., Mount Jaya; Hope et al., 1976) or by randomly assigning the coordinates of raster cells with the same elevation (±10 m) as the plots in that area (i.e., Mount Wilhelm and Drakensberg; Brand et al., 2015;Wade & McVean, 1969), using the SRTM-3 digital elevation model at 30 m resolution (Farr et al., 2007;NASA & JPL, 2013). Species cover values with discrete scales were transformed to the mean value of the corresponding percentage interval. Species names were harmonized using the Taxonomic Name Resolution Service (Boyle et al., 2013) online tool (https://tnrs.biend ata.org/) with default settings, updating the names to the most recent nomenclature and merging subspecies and varieties to the species level by summing their respective cover values.
The initial dataset, consisting of 10,408 plots, was filtered further by removing plots with tree species or incomplete taxonomic identification. When taxa identified to the genus level or higher taxonomic rank represented ≥ 10% of the plot vegetation cover, the corresponding plot was discarded; otherwise, we removed those taxa from the plot record (3,086 plots from which at least one taxon was removed; median number of taxa removed = 1). Each plot was then assigned to a region based on its location. Regions were defined based on the approximate extent of ecoregions (Olson et al., 2001), which represent an ecologically meaningful framework for identifying distinct geographical units at the global scale. Given that the names of some ecoregions did not reflect the presence of an alpine vegetation belt, we renamed these regions after the main mountain ranges where the plots were located, following Körner et al. (2017) (Supporting   Information Table S2). For the analyses, we retained only regions with ≥ 60 plots and removed extremely small or large plots (<.25 or >400 m 2 ). To filter out compositional outliers, we performed a detrended correspondence analysis (DCA) on each regional dataset, excluding those plots whose score on the first axis (DCA1) was larger or smaller than 10 times the width of the interquartile range from the median. After removing the outliers, the gradient length of DCA1 ranged from 3.6 to 9.9 standard deviation units of species turnover within different regions (Supporting Information Table S3), indicating different, yet high, degrees of regional beta diversity. Finally, to assess the representativeness of our dataset, we compared the climatic space of the plots against the climatic envelope of global alpine areas (Testolin et al., 2020;Supporting Information Figure S1). The

| Diversity measures
Given that the number of samples differed considerably among regions, we estimated regional species richness using samplebased rarefaction and extrapolation (Chao et al., 2014) with the software R (R Core Team, 2020) and the package iNEXT (Hsieh et al., 2016). This technique allows a statistically sound comparison of diversity across groups with different sample sizes through the construction of sampling curves for species richness. These curves can be rarefied (interpolated) to smaller sample sizes or extrapolated (predicted) to larger sample sizes (Chao et al., 2014;Hsieh et al., 2016). Here, we estimated the regional richness for a unique sample size of 180 plots, corresponding to approximately three times the smallest regional sample (Figure 1b,c). We chose 180 plots as a trade-off between the loss of data in intensively sampled regions versus the inclusion of all regions in the analyses. As such, these estimates should not be interpreted as representing the total regional species pools, but rather as comparable estimates of regional richness. Given that our global dataset comprised plots of different sizes, we evaluated the effect of plot size on the species richness estimates. To do this, we compared the same estimates using three subsets of different plot sizes (small, <10 m 2 ; medium, ≥ 10 and <100 m 2 ; and large, ≥ 100 m 2 ). For those F I G U R E 1 Overview of the alpine vegetation dataset and regional species richness. (a) Spatial distribution of alpine vegetation plots. (b) Number of plots collected in this study (N) and estimated species richness (S est ) for a comparable number of 180 plots in 26 alpine regions and six biogeographical realms. (c) Rarefaction curves of species richness for each region. Dashed lines indicate extrapolated values beyond the available number of plots. Continuous lines indicate that regional estimates were interpolated from larger sample sizes. The shaded areas represent the 95% bootstrap confidence intervals [Colour figure can be viewed at wileyonlinelibrary.com] regions where ≥ 60 plots of each of the three different sizes were available (Alborz Mountains, Central and Eastern Alps, Colombian and Ecuadorian Andes, Eastern African Mountains, South Central Rocky Mountains and Western Carpathians), we compared richness estimates obtained from the different subsets. We found that, regardless of the subset used, the relative differences among regions were largely preserved, especially for those datasets comprising large numbers of plots (e.g., Central and Eastern Alps and Western Carpathians). An exception was the region of the Colombian and Ecuadorian Andes, where regional richness estimates were highly dependent on plot size (Supporting Information Figure S2). However, large and small plots both resulted in lower species richness estimates compared with medium-sized plots, suggesting that the differences are driven by different vegetation types being sampled with differently sized plots.
For each plot, we calculated community richness as the total number of species. We evaluated latitudinal patterns of regional and community richness using generalized additive models (GAMs), with a smoothing term for latitude. Our alpine regions were characterized by very different extents, and plot size varied widely. To account for different regional extents and plot sizes in the evaluation of species richness patterns, we also fitted GAMs on the residuals from ordinary least square regressions of ln(regional richness) as a function of ln(current local alpine area) and ln(community richness) as a function of ln(plot size). The procedure for calculating local alpine area is described in section 2.4.

| Environmental predictors
To analyse the drivers of species richness, we retrieved a set of environmental variables linked to plant diversity in the alpine belt from online sources. We calculated several climatic predictors at the plot level using digital sources at c. 1 km resolution. We used data from CHELSA (Karger et al., 2017) within the time frame of the growing season, defined as days with mean temperature > .9°C (Paulsen & Körner, 2014). Given that daily temperature data were not available, we estimated the growing season using monthly averages, including the months with a mean temperature > .9°C. Although this might have resulted in a sharper delimitation of season lengths, it probably had little effect on our global analyses. We included the mean temperature, precipitation, growing degree days and mean potential evapotranspiration of the growing season, all of which have been reported to have positive effects on photosynthetic activity and species richness in alpine areas (Körner, 2003;Moser et al., 2005;Nagy & Grabherr, 2009). Growing degree days (i.e., the sum of monthly temperatures > .9°C multiplied by the total number of days in such months) were calculated using the "growingDegDays" function of the R package envirem (Title & Bemmels, 2018). The mean potential evapotranspiration of the growing season was estimated with the "hargreaves" function of the R package SPEI (Beguería & Vicente-Serrano, 2017), using maximum and minimum monthly values of temperature and monthly precipitation. The monthly values of potential evapotranspiration obtained were then averaged across months with a mean temperature > .9°C.
Together with climate, soil pH is known to be a significant driver of species richness in the alpine belt (Vonlanthen et al., 2006) and is a good surrogate for the dominant bedrock, effectively differentiating calcareous and siliceous substrates Lenoir et al., 2010). We derived estimates of soil pH from the SoilGrids database, averaging the values estimated at 5 and 15 cm depths (Hengl et al., 2017). When the pH value was missing for a plot location (45 plots), we assigned the value of the closest pixel to the plot. Despite the limitations posed by the use of global datasets to estimate fine-scale soil properties (Hengl et al., 2017), the obtained values covered a wide span of soil pH variation in alpine environments (4.40-8.35) and were, therefore, useful to distinguish dominant bedrock types.
Regional values of the predictors computed at the plot level were then estimated as the average of all vegetation plots within a region. For climatic predictors and soil pH, we also calculated the standard deviation of the predictor to test for the effect of environmental heterogeneity.

| Geographical and historical predictors
In addition to environmental variables, large-scale geographical factors, such as area and isolation, are known to influence the current diversity of island-like ecosystems (MacArthur & Wilson, 1967;Whittaker et al., 2008Whittaker et al., , 2017, as do their historical changes caused by climatic fluctuations Weigelt et al., 2016). We delimited alpine area as the portion of land with a mean temperature of the growing season between 3.5 and 6.4°C or with a length of the growing season between 1 and 3 months (Paulsen & Körner, 2014). We did this both for current climatic conditions and considering climate during the Last Glacial Maximum (LGM) (Supporting Information Figure S3).
Alpine areas were calculated at two scales, reflecting the spatial extents of each regional sample (local area) and the total continuous alpine area extending beyond the samples (total area). We defined the local area as the extent of the alpine area contained within the convex hull formed by all plots of each region. In some cases, the coarse resolution of the climatic datasets used to estimate the alpine areas failed to detect any alpine patch within the hulls. Therefore, we applied a 5 km buffer around each hull to include at least some alpine area patches for all regions. The total alpine area for each region was estimated as the continuous extent of all alpine patches intersected by the hulls, reflecting the total extent of alpine habitats available to species dispersal (Supporting Information Figure S4). Calculating alpine areas at two scales also allowed us to estimate the completeness of the regional samples by calculating the percentage of the local alpine area encompassed by the samples over the total alpine area (Supporting Information Table S3). Given that the local and total ln-transformed areas were highly correlated (Pearson's r > .8), only the former was retained in the subsequent analyses.
In addition, we estimated the current and LGM isolation as the minimum distance from the centroid of each alpine region to the boundary of the nearest alpine area ≥ 1,000 km 2 . We set a threshold of 1,000 km 2 to exclude smaller alpine patches that could still be part of the same alpine region, that is, islands of the same archipelago . If an alpine region had a total area ≥ 1,000 km 2 , isolation was set to zero. Current and LGM alpine areas and isolation were ln-transformed.
Given that past climatic changes could have affected current diversity patterns (Graham et al., 2014), we also calculated the velocity of climate change since the LGM as a measure of regional climatic instability (Loarie et al., 2009;Sandel et al., 2011), using the "gVoCC" function of the VoCC package (Molinos et al., 2019) with current and LGM mean annual temperatures. The latter was calculated as the average of the two PMIP3 climatic datasets derived using the CCSM4 and MIROC-ESM climate models (Sandel et al., 2011;Weigelt et al., 2013Weigelt et al., , 2016. Finally, we included biogeographical realms (Keil & Chase, 2019) as a proxy for differences in evolutionary history. Owing to the lack of regional data, we did not account for differences in the geological history of mountains. However, we acknowledge that this could influence speciation and partly explain species richness (Whittaker et al., 2008).

| Statistical analyses
We tested the influence of environment, geography and history on estimated regional richness by fitting individual Poisson generalized linear mixed-effects models (GLMMs) to each predictor with the "glmer" function of the R package lme4 (Bates et al., 2015). Initially, we tested univariate relationships to select a set of significant variables to be used in subsequent multivariate modelling. To account for uncertainties in regional richness estimates, we weighted the observations by the inverse of their 95% confidence interval width. We also scaled the predictors by subtracting their mean and dividing by their standard deviation across the regions to ensure model convergence.
To control for overdispersion and reduce the risk of type I errors, we added an observation-level random effect (OLRE) to the models, that is, a unique level of a random effect for each data point that models the extra-Poisson variation present in the data (Harrison, 2014). The ratios between the sum of squared Pearson residuals and the residual degrees of freedom of the fitted models with OLRE indicated no additional overdispersion. Next, we analysed the correlations among significant predictors with the Pearson correlation coefficient. We found that some of our regional variables were strongly correlated with one another (Supporting Information Figure S5), limiting our ability to distinguish partial contributions. However, we built alternative multivariate models by retaining only the significant, uncorrelated predictors. Finally, we checked for the presence of spatial autocorrelation in model residuals with the Moran's I test implemented in the "testSpatialAutocorrelation" function of the DHARMa package (Hartig, 2020) and found none (Supporting Information Table S4). We also fitted a null (intercept-only) model to compare the goodness of the fits. Models were compared using the Akaike Information Criterion corrected for small sample sizes (AICc) in addition to marginal and conditional R 2 (mR 2 and cR 2 ), calculated with the "r.squaredGLMM" function of the MuMIn package (Barton, 2019). Given that the only random effect in the models was an OLRE, mR 2 = cR 2 .
We modelled community richness by fitting Poisson GLMMs including the environmental, geographical and historical predictors. Growing degree days and precipitation of the growing season were highly correlated with temperature and evapotranspiration of the growing season, respectively (Pearson's r > .6). Likewise, area and isolation-related variables were highly correlated with one another. Thus, to avoid multicollinearity issues, we retained temperature and evapotranspiration of the growing season, in addition to current and LGM local areas, and excluded the other variables. We also accounted for different plot sizes by adding their ln-transformed values to the model and controlled for overdispersion by adding an OLRE. Given that the plot-level predictors were derived from digital datasets at 1 km resolution, we randomly selected one plot for each .01° cell (c. 1 km). We repeated the process 999 times and obtained as many random subsets of 2,534 plots, that is, one plot for each .01° cell. Before modelling, all predictor variables were scaled to ensure model convergence. We then fitted the GLMMs to each of the 999 subsets. Given that the Moran's I tests highlighted strong spatial autocorrelation of the models' residuals, we re-fitted the models including random intercepts for .05 (≈5 km) and .1 (≈10 km) degree cells, which largely resolved the issue (Moran's I ≈ 0; p > .05 for 75% of model fits). We also tested for regional effects by fitting another model to the 999 subsets, with an additional random intercept for regions. Finally, we averaged the fixed effect coefficients of the resulting models using weights based on their AICc with the "model. avg" function of the MuMIn package. The two resulting averaged models (with and without the random intercept for region) were compared using mean AICc, mR 2 and cR 2 , obtained by calculating the weighted average of the respective indices for the 999 fits.

| RE SULTS
According to sample-based rarefaction and extrapolation of regional richness (estimated for 180 plots), the richest alpine regions in this study were the Colombian and Ecuadorian Andes (Neotropics; 543  (Figure 3a). The same patterns emerged even when different regional extents and plot sizes were accounted for, with an additional peak of regional richness at F I G U R E 2 Latitudinal patterns and drivers of estimated regional species richness. (a) Regional plant species richness estimated for 180 plots (S est ). The scatterplot on the right represents the latitudinal trend. The three horizontal grey lines on the map and the scatterplot represent the equator and the tropics. The black line represents a generalized additive model (GAM) fit. (b) Single-predictor models of regional species richness. The dots represent the regional plant species richness estimated for 180 plots (S est ). The error bars represent the 95% bootstrap confidence intervals of the richness estimates. Black lines represent the individual Poisson generalized linear mixed-effects model (GLMM) fits. The grey bands are the 95% bootstrap confidence intervals. Marginal R 2 (mR 2 ) and model significance are reported. Significance codes: <.001 (***); <.01 (**); <.05 (*). The numbers of the main coldspot and hotspot regions are reported according to  Figure S6a,b).
At the community scale, species richness was positively related to the evapotranspiration of the growing season (β = .13 ± .05; p < .001), velocity of climate change (β = .10 ± .05; p < .001) and LGM alpine area (β = .19 ± .05; p < .001), whereas it was negatively related to soil pH (β = −.12 ± .05; p < .001). Nearctic plots were generally poorer in species than plots in other realms (β = −.44 ± .22; p < .001; Figure 3b; Supporting Information Table S5). Overall, the fixed effects explained 22% of the variance, and the inclusion of the random effects controlling for the spatial aggregation of plots at 5 and 10 km increased the explained variance to 58%. The inclusion of regions as an additional random effect increased the total explained variance further to 65% and left as significant fixed effects the mean temperature of the growing season (β = .04 ± .03; p < .05) and soil  Table S5).

| Regional patterns and drivers
Our results, based on a representative sample of global alpine vegetation, showed a latitudinal pattern of plant species richness, with peaks at mid-latitudes and around the equator. The highest estimate of regional richness was detected in the Colombian and Ecuadorian Andes (Neotropics). This region is home to the páramo ecosystem, a centre of plant diversity within the tropical biodiversity hotspot known to host the richest alpine flora in the world (Madriñán et al., 2013;Myers et al., 2000). At higher latitudes, we also found that the Pamir and Altai Mountains (Eastern Palaearctic) exhibited regional richness comparable to the páramos, representing hotspots of alpine plant diversity outside the tropics. This is consistent with previous studies that highlighted the high plant diversity of the Altai, Pamir, in addition to other Central Asian mountain systems (Agakhanjanz & Breckle, 1995;Brummitt et al., 2020;Kier et al., 2005;Körner, 1995;Xing & Ree, 2017). Nevertheless, these putative mid-latitude alpine hotspots are generally excluded from global centres of biodiversity, despite their importance as refugia for cold-adapted plant species . When accounting for the extent of the local alpine area, the Drakensberg (Afrotropics) and the Australian Alps (Australasia) emerged as centres of alpine plant richness of the Southern Hemisphere. Indeed, the high-elevation plateau of the Drakensberg has been widely recognized as a continental hotspot of botanical diversity (Brand et al., 2019;Carbutt, 2019), and the Australian Alps have been listed among the main national areas of plant species richness (Bell et al., 2018;Crisp et al., 2001). Other Africa (Anthelme & Dangles, 2012) and is an active volcano, with eruptions that limit development of vegetation (Nagy & Grabherr, 2009).
Moreover, the Northern Scandes were completely glaciated during the Pleistocene glacial maxima and located far from the Southern European glacial refugia further south (Lenoir et al., 2010).
Our models showed that the regional richness of alpine ecosystems is mostly independent of macroclimatic gradients. An analogous decoupling pattern has also been reported for the global diversity of grasses as a response to biogeographical history and the adaptation of certain lineages to cold and arid environments (Visser et al., 2014). Indeed, global alpine areas are climatically constrained toward low-temperature conditions (Körner, 2003;Nagy & Grabherr, 2009;Paulsen & Körner, 2014;Testolin et al., 2020).
Thus, although alpine plants respond to changes in temperature and light because of topography, large-scale richness patterns of alpine vegetation seem to be largely independent of energy gradients that determine species diversity at lower elevations (Hillebrand, 2004).
Moreover, global alpine areas are subjected to different amounts of precipitation and are differentiated along a gradient of humidity (Körner, 2003;Nagy & Grabherr, 2009;Testolin et al., 2020).
Although our dataset encompasses a large portion of the variation in water availability of global alpine areas (Supporting Information Figure S1), the effect of precipitation on regional richness was not significant. This suggests that the association of water availability with plant species richness might be restricted to local scales and especially to arid regions, where precipitation is the main factor limiting plant growth (Palpurina et al., 2017).
Contrarily to macroclimate, we found a positive effect of the extent of current alpine area and a negative effect of isolation. The importance of these factors is consistent with the predictions of the theory of island biogeography (MacArthur & Wilson, 1967), which posits that larger, less isolated islands are characterized by lower extinction rates and greater chances of being colonized by new species. Nevertheless, the historical legacy of the extent and isolation of alpine areas during the LGM also left a strong imprint on regional richness patterns that is independent of their current geographical characteristics. The extent of alpine areas during the LGM was the second strongest predictor of regional richness and, together with the current area, explained almost 80% of the variance. This is consistent with recent refinements of the theory of island biogeography that incorporate the effect of Late Quaternary climate oscillations on oceanic islands Weigelt et al., 2016).
Pleistocene glacial-interglacial cycles acted like a "historical sieve" (Körner, 1995) on alpine plant diversity. During glacial periods, downslope shifts of the alpine belt resulted in increased surface area and connectivity of tropical alpine archipelagos, in addition to colonization and diversification processes in mid-latitude mountain ranges that favoured in situ speciation (Flantua et al., 2020). The high species richness found in the Andes is probably the result of multiple contingencies related to South American tropical diversity and strong past connectivity of these mountains , which are not co-occurring in any other tropical region. Moreover, in Central Asian mountains, the emergence of habitat corridors during glacial periods resulted in extensive long-distance dispersal, with the consequent admixture of previously isolated floras (Agakhanjanz & Breckle, 1995;Agakhanyantz & Lopatin, 1978). Indeed, the Pamir Mountains are a continental hub for floristic migrations (the Pamir Knot) that connects south-central Asian ranges to the northern Siberian mountains (Agakhanjanz & Breckle, 1995), whereas the Altai Mountains connect diversity between Euro-Siberian and Central Asian floristic regions (Chytrý et al., 2012). Our results also show a positive relationship between regional richness and soil pH variability (a surrogate for bedrock heterogeneity) largely driven by the Pamir and Altai Mountains.
This finding confirms the effect of habitat heterogeneity on species richness inherent to larger areas (Lomolino, 2000a) through the occurrence of more diverse bedrock types (Moser et al., 2005).

| Community patterns and drivers
At the community scale, the latitudinal pattern of species richness was less pronounced than at the regional scale, with a single peak at mid-latitudes of the Northern Hemisphere, but a wide range of values within all regions. While controlling for plot size, we found a positive effect of evapotranspiration of the growing season and a negative effect of soil pH on community richness. The former is consistent with the species-energy hypothesis, which states that more productive communities (i.e., where higher temperatures and solar radiation support greater photosynthetic rates) are also richer in species (Wright, 1983). The latter could be explained by the absence of strongly acidic soils (pH < 4) in our dataset. Furthermore, soils with high pH values can be linked to reduced nutrient availability in harsh conditions and the confounding effect of reduced precipitation (Chytrý et al., 2007;Palpurina et al., 2017), explaining the lower species richness in our dataset. Despite the underlying causes of these effects, our results are in line with the role of energy-driven processes and bedrock mineralogy as determinants of vascular plant species richness in alpine communities (Lenoir et al., 2010;Moser et al., 2005;Vonlanthen et al., 2006). Nevertheless, evolutionary and historical factors might also affect current patterns of community richness (Ricklefs & He, 2016). The Nearctic realm exhibited lower community richness than any other biogeographical realm, possibly owing to limited evolutionary radiation of the North American temperate flora (Qian & Ricklefs, 2000). In addition, the velocity of climate change and the LGM extent of alpine area were both positively related to community richness, indicating that the greater availability of alpine habitats in the past influences plant species richness at the community scale (Pärtel & Zobel, 1999).
Large-scale environmental factors, however, explained only a limited proportion of the variation in community richness compared with regional and sub-regional effects, suggesting that dispersal-related processes and other spatially structured factors strongly influence local richness patterns (Dormann et al., 2007). In alpine landscapes, these effects are regulated by elevational and meso-topographical gradients that affect microclimatic conditions and local plant diversity (Bruun et al., 2006;Jiménez-Alfaro et al., 2014;Scherrer & Körner, 2011). A weak influence of global macroclimatic gradients on local communities has also been detected for functional diversity across plant formations (Bruelheide et al., 2018), but it had not been tested before on a single ecosystem. Our results, therefore, support the dominant role of within-region effects linked to postglacial spatial configuration and historical contingencies, rather than macroclimatic factors, when explaining the global variation of alpine local communities.

| Data constraints and assumptions
Despite including several alpine regions across all continents and latitudes, our dataset lacked information about some outstanding centres of alpine plant diversity, such as the Himalayas and Hengduan Mountains (Ding et al., 2020;Favre et al., 2015;Muellner-Riehl et al., 2019;Xing & Ree, 2017) or the Caucasus (Agakhanjanz & Breckle, 1995;Körner, 1995), owing to the limited availability of community data from these areas. Therefore, our results related to the latitudinal patterns of regional and community richness could be refined further by the future addition of data from currently missing regions. Nevertheless, our aim was not to present a complete census of global alpine regions, but to assess their richness patterns and the corresponding drivers using a representative sample. In this respect, the collection of georeferenced vegetation plots presented here encompasses all continents and a wide range of latitudes; it represents regions with markedly different biogeographical history and vegetation types growing on different substrates, and covers a large portion of the climatic envelope of global alpine areas (Supporting Information Figure S1). Despite the lack of some remarkable alpine regions, our dataset allowed us to highlight the presence of extratropical alpine diversity centres and the importance of historical factors in shaping the current alpine plant richness patterns.
We also note that the use of heterogeneous surveys from different collectors might create issues related to different sampling effort among regions. We controlled for sampling effort in two ways. First, we used rarefaction and extrapolation techniques that assumed that the spatial distribution of plots in each region was representative of the regional diversity. Although this assumption is difficult to prove without additional data, we note that our regional samples were selected to capture the local heterogeneity of vegetation types and covered a wide range of elevations in all regions (Supporting Information Table S3), thus increasing the probability that our regional richness estimates were correlated with regional species pools. Second, we explicitly quantified the proportion of the alpine area sampled in each region, thus allowing inter-regional comparisons even when the samples covered very different extents or only a small fraction of the total available alpine area (e.g., Ladakh Range, Pamir Mountains or Southern Cordillera Occidental Peru).
In addition, we note that taxon concepts might not be applied consistently across all datasets, that is, they are the result of "lumping" and "splitting" of taxa delimitations that change with time and from place to place (Rouhan & Gaudeul, 2014;Wiser, 2016).
Although a harmonized species nomenclature cannot account fully for this taxonomic bias (Wiser, 2016), it still represents the most effective tool to address taxonomic inflation in macroecological studies (Isaac et al., 2004). Indeed, by correcting misspelt names and merging synonyms, we assume that the main sources of error relevant to the estimation of species richness in different regions were removed, and remaining issues about potential pitfalls in the geographical distribution of species (Boyle et al., 2013) are not pertinent to the present study.

| Conclusions
Overall, we found that the latitudinal distribution of plant species richness in alpine ecosystems is decoupled from the general latitudinal diversity gradient and that it relates to regional idiosyncrasies, rather than macroclimatic gradients. Although our results are conclusive enough to support that current and historical effects of area, isolation and environmental heterogeneity exert an overarching influence on vascular plant richness in global alpine ecosystems, we are still far from understanding the processes behind such effects.
Future alpine research should, therefore, consider local information about soil biotic and abiotic composition, topographical features and microclimatic variation at the regional scale. Additionally, further efforts should be oriented toward the collection of plant community data from underrepresented regions. Indeed, this work is the starting point for defining global hotspots of alpine plant diversity, and further investigations including patterns of endemism, functional variation and phylogenetic diversity are still needed. This type of information, together with dynamic regional diversity models accounting for spatio-temporal connectivity, will provide a better understanding of the patterns we have found here and a tool for the effective conservation of alpine biodiversity in response to climate change.