Origins of global mountain plant biodiversity: Testing the ‘mountain‐geobiodiversity hypothesis’

Our objective is to analyse global‐scale patterns of mountain biodiversity and the driving forces leading to the observed patterns. More specifically, we test the ‘mountain geobiodiversity hypothesis’ (MGH) which is based on the assumption that it is not mountain‐uplift alone which drives the evolution of mountain biodiversity, but rather the combination of geodiversity evolution and Neogene and Pleistocene climate changes. We address the following questions: (a) Do areas of high geodiversity and high biodiversity in mountains overlap, that is can mountain geodiversity predict mountain biodiversity? (b) What is the role of Pleistocene climate change in shaping mountain biodiversity? (c) Did diversification rate shifts occur predominantly with the onset of more pronounced climate fluctuations in the late Neogene and Pleistocene fostering a ‘species pump’ effect, as predicted by the MGH?


| INTRODUC TI ON
Global comparisons of the biodiversity of mountains have fascinated natural scientists for a long time. Famously, the 1802 expedition by Alexander von Humboldt (*1769 to mount Chimborazo in the Cordillera Occidental range of the Andes led him to develop the first map of biodiversity (Humboldt, 1807;Wulf, 2016), which summarizes the vegetation and abiotic parameters of eleven elevational zones found on the Chimborazo in comparison with other mountains he had previously visited (e.g. Alps, Pyrenees, on Tenerife). Humboldt appreciated vegetation from the perspective of climate and locality rather than taxonomic categories, a completely new approach at his time, that still has repercussions on our understanding of ecosystems today (Wulf, 2016). He was also aware that the biotic and abiotic world would have possibly shared a common temporal history (Humboldt, 1807). Humboldt regarded nature as being composed of a 'net of life', which he considered should be investigated in an integrative manner from all disciplines rather than solely through the eyes of a botanist, a geologist or a zoologist. Recently, the term 'geobiodiversity' or 'geo-biodiversity' has marked the increased interest in the union of biological and geological sciences (e.g. Pătru-Stupariu et al., 2017;Mosbrugger, Favre, Muellner-Riehl, Päckert, & Mulch, 2018; compare also Fan et al., 2013;Štursa, 2013;Kumarasami, Ramkumar, & Jyotsana, 2013;Izakovičová, Miklós, & Oszlányi, 2014; first mention probably by Lee et al., 2004; in more remote context probably by Gressel, 2001; previously also referred to as 'ecodiversity', for example Naveh, 1994;Barthlott, Lauer, & Placke, 1996), although these terms have thus far not been formally defined in the literature. Building on Humboldt´s idea of exploring the relationship between geodiversity and biodiversity in mountains in particular, our study aims at providing a global comparison of mountain geodiversity and vascular plant diversity, looking both at short-term (ecological) as well as long-term (evolutionary) processes.
Studies of the correlation between patterns of plant diversity and contemporary abiotic conditions at geographical meso-to macroscales have focused predominantly on the environmental variables relating to two main hypotheses, namely the 'energy-water' hypothesis and the 'habitat heterogeneity' hypothesis (Chen, Jiang, Ouyang, Xu, & Xiao, 2011;Currie, 1991;Gaston, 2000;Hawkins et al., 2003;Jimenez, Distler, & Jorgensen, 2009;Kreft, Jetz, Mutke, & Barthlott, 2010;Turner, Lennon, & Lawrenson, 1988;Wright, 1983). For example, a study of vascular plants of the Tibeto-Himalayan region found habitat heterogeneity, water and energy variables to be important environmental determinants of vascular plant richness patterns (Mao et al., 2013). However, as correlations between plant richness and environmental variables vary with life form and taxonomic scale, the explanatory power of these variables also changes on different spatial scales due to the proportion of life forms and the asymmetric effects of these drivers (Mao et al., 2013). Irrespective of any focus on a particular group of organisms, studies have shown that high levels of abiotic features ('geodiversity') may result in and thus predict high levels of biodiversity (e.g. Bailey, Boyd, Hjort, Lavers, & Field, 2017;Braun, Mutke, Reder, & Barthlott, 2002;Najwer, Borysiak, Gudowicz, Mazurek, & Zwolínski, 2016;Parks & Mulligan, 2010;Zwolínski, Najwer, & Giardino, 2018). Geodiversity encompasses 'the natural range (diversity) of geological (rocks, minerals, fossils), geomorphological (land form, processes) and soil features' (Gray, 2004, p.8) as well as the atmosphere/climate and hydrosphere (e.g. Parks & Mulligan, 2010). The fact that on a global and regional scale geodiversity may explain 50%-70% of the biodiversity distribution (compare Antonelli et al., 2018 for a study on mountain terrestrial tetrapods) led to the formulation of the 'geobiodiversity' concept (see below). In addition to contemporary conditions, deep-time historical processes, such as geological events (e.g. mountain building) or long-term climatic fluctuations that potentially affect the rates of in situ speciation and/or extinction of lineages, need to be considered (Ricklefs, Latham, & Qian, 1999), which may act in parallel on species richness patterns, but do not necessarily need to do so (Whittaker, Willis, & Field, 2001).
According to the 'geobiodiversity concept', biodiversity research requires a holistic or systemic approach that links (a) biodiversity, (b) Main conclusions: Our analyses point towards an important role of historical factors on mountain plant species richness. Mountain systems characterized by small elevational ranges and strong modifications of temperature profiles appear to harbour fewer radiations, and fewer species. In contrast, mountain systems with the largest elevational ranges and stronger overlap between today´s and LGM temperature profiles are also those where most plant radiations and highest species numbers were identified.

K E Y W O R D S
biodiversity, climatic fluctuations, geobiodiversity, geodiversity, gradient, heterogeneity, Humboldt, mountains, Pleistocene, vascular plants geodiversity (in the broad sense) and (c) their historical evolution. The inclusion of historical processes is crucial because past geological events, climatic fluctuations, biotic radiations and dispersal/migration events, in particular during the Neogene and Pleistocene, had a major impact on the present-day distribution of biodiversity. Based on this concept, Mosbrugger et al. (2018) developed the 'mountain-geobiodiversity hypothesis' (MGH) to explain the high levels of biodiversity found in the Tibeto-Himalayan region. The MGH proposes that three boundary conditions are required to maximize the impact of mountain formation and surface uplift on regional biodiversity patterns and are key for the origination of montane biodiversity. These are (a) the presence of lowland, montane and alpine zones (full elevational zonation), (b) climatic fluctuations for a 'species pump' effect and (c) high-relief terrain with environmental gradients. The MGH thus considers all major factors responsible for high biodiversity in a given mountain region (Muellner-Riehl, 2019). First, the presence of steep ecological gradients would have allowed local lineages to adapt to a high variety of niches along the elevational gradient and a range of differently pre-adapted lineages to immigrate into mountains. Second, climatic fluctuations would have fostered a 'species pump' effect (here used to describe repeated cycles of connectivity and isolation as drivers of divergence; for original use of the term, with refugia acting as 'species pumps' during climatic fluctuations, see Haffer, 1997), with high-relief terrain providing geographical barriers for allopatric speciation. Third, high-relief terrain, together with the presence of a full elevational zonation, would have provided refugia and short migration distances into favourable habitats during climatic changes, fostering connectivity and reducing the danger of extinction (Muellner-Riehl, 2019). Inherently, the MGH also aims at explaining the frequently observed time-lag between the onset of mountain uplift (beyond sea level) and biotic diversification/radiations (Mosbrugger et al., 2018). It remains yet untested whether the postulations of the MGH will hold true for other mountain regions than the Tibeto-Himalayan region as well. In this case, a higher geological diversity in species-rich mountains would be expected, with species radiations dating back in time no earlier than the Pleistocene (ca. 2,6 Mya; rapid and profound climatic fluctuations), or end of Miocene (ca. 7 Mya; pronounced wet-dry fluctuations), as found for the Tibeto-Himalayan region´s biota (Mosbrugger et al., 2018).
Here, we present an approach to evaluate the MGH on a global scale, incorporating plant species diversity, geodiversity, as well as climatic and evolutionary history. We thus complement and expand the recent approach by Antonelli et al. (2018), who investigated global determinants of mountain (tetrapod) biodiversity, by explicitly considering the temporal dimension. Specifically, we ask the following questions: 1. Do areas of high geodiversity and high biodiversity in mountains overlap, that is can mountain geodiversity predict mountain biodiversity? 2. What is the role of Pleistocene climate change in shaping mountain biodiversity?
3. Did diversification rate shifts of mountain plants occur predominantly with the onset of more pronounced climate fluctuations (wet-dry and/or warm/cold) in the late Neogene and Pleistocene (from ca. 7 Mya) fostering a 'species pump' effect, as predicted by the 'mountain-geobiodiversity hypothesis'?

| Biodiversity and geodiversity/climate data
First, we compiled vascular plant species inventory data from mountain systems across different floristic realms (including equatorial, subtropical, temperate and polar zones; see Figure 1 and Table S1.1).
For each region, we extrapolated the number of vascular plant species per 10,000 km 2 using the power model of the species-area relationship (Arrhenius, 1921) with a biome-specific constant (z) that describes the extent to which species richness varies with area size following Kier et al. (2005). If a region overlapped with more than one biome, species richness was calculated as a weighted average according to the relative area of overlap with each biome (e.g. for the Tibeto-Himlayan region, separate species numbers for the Qinghai-Tibetan Plateau, the Himalayas and the Hengduan Mountains were not available, and therefore combined). Second, for each study area, we extracted a number of environmental variables. We used the GTOPO30 digital elevation model (U.S. Geological Survey, 1996) to calculate the elevational range and calculated the average net primary productivity (NPP) from MODIS data (MOD17A3; Running, Mu, & Zhao, 2011). In addition, we quantified geodiversity by computing a geodiversity (GD) index, using globally harmonized geological, soil, hydrological and topographical datasets in grid cells of 10 × 10 km.
A geological dataset derived from the Global Lithological Map database (Hartmann & Moosdorf, 2012) harmonized into sixteen lithological classes was used to compute a lithological index, based on the number of the lithological formations in each grid cell. We derived a soil index for each grid cell based on the number of soil types from the SoilGrids repository (Hengl et al., 2017) which is based on the World Reference Base for Soils system (IUSS Working Group WRB, 2015). SoilGrids provides interpolated gridded soil information based on machine learning techniques, and is trained with more than 150,000 soil profile descriptions. For the hydrological index we used the global lakes (Messager, Lehner, Grill, Nedeva, & Schmitt, 2016) and rivers (Lehner, Verdin, & Jarvis, 2008) and calculated lake area and river length per grid cell. The slope index was based on the Multiple-Error-Removed Improved-Terrain Digital Elevation Model elevation database (Yamazaki et al., 2017). From this we computed, reclassified and summed the slope range and the standard deviation of the slope for each grid cell. All indices were reclassified into five equal interval classes and combined into a GD index (at 10 × 10 km resolution), which could theoretically vary from 1 to 20. However, the realized GD index only ranged from 2 to 14. The minimum realized GD value of 2 was because each cell contains at least a slope index value and a soil index or a hydrological index value. The lithology dataset contained no data in some of the grid cells because ice, glaciers and water were excluded as they are non-lithology units. The maximum realized GD value of 14 indicates that combinations of very high sub-indices (≥15) do not occur. For each mountain polygon, we used all GD index values (at 10 × 10 km resolution) and applied Zonal Statistics in ArcGIS Pro to obtain the mean for each mountain range.
Finally, to assess the impact of past climatic fluctuations on present-day species diversity, we calculated the region-specific change in mean annual temperature, a strong predictor of species diversity in mountain systems (Antonelli et al., 2018), between the Last Glacial Maximum (LGM, c. 21,000 years ago) and today. We used the temperature anomaly (ΔT) of the LGM relative to the present as a rough proxy for Quaternary climatic fluctuations (the last 2.6 million years), as it covers almost the full Quaternary temperature range. For each region, we extracted mean annual temperature both for the present-day and the LGM (calculated as the average of three circulation models: CCSM, MIROC & MPI-ESM) using the Worldclim dataset (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) at a 2.5 arc-minute resolution. The mean and mode of temperature differences for each mountain range were calculated using the raster (Hijmans, 2018) and rgeos (Bivand & Rundel, 2018) packages in R (v. 3.5.1, R Development Core Team, 2016). To visualize the changes in available environmental space in each mountain system, we plotted the relative frequency distributions of present-day mean annual temperatures during the LGM. We accounted for the fact that present-day temperature regimes might have been available for plant species elsewhere (i.e. in neighbouring regions/valleys) during climatic fluctuations by constructing a 200 km buffer around each study area.
We used generalized linear models (GLMs) with Poisson error distribution and log link to relate plant species richness to contemporary and historical environmental variables (NPP, geodiversity, elevation and temperature changes since the LGM). We evaluated all single and multiple-parameter GLMs, using the Akaike information criterion (AIC, Burnham & Anderson, 2002) to identify the best-fit model. Prior to the analyses, all parameters were standardized to allow for a comparison of the relative importance of different predictors. All analyses were performed using the R statistical software package (v. 3.5.1, R Development Core Team, 2016).

| Timing of mountain plant radiations across the globe
We gathered studies on time-calibrated plant radiations for different mountain systems of the world (a list of the data sources is found in Appendix S1). We distinguished between studies that performed diversification rate analyses (Table S1.3 in Appendix S1) and those that did not (Table S1.4 in Appendix S1). For the former, we identified whether shifts in diversification rates occurred, and for all studies, we determined with which mountain system radiations were predominantly associated. We ignored species-poor clades resulting from (usually recent) dispersal to other mountain ranges. Thus, our scoring may not correspond to the complete current distribution of each taxon.
Potential drivers (e.g. uplift, climate modification or variation, polyploidization, etc.) or their precursors (biome shifts, key opportunities, pre-adaptations, etc.) of diversification were then attributed to each radiation as suggested in each study by the authors. Some of these may function both as precursors or drivers depending on the relative timing of their evolution with regard to radiation. This is, for example the case for key innovations (a trait allowing a species to conquer a new adaptive zone; Hodges & Arnold, 1995;Simões et al., 2016), which may evolve repeatedly during a radiation (thus would not be an attribute of all taxa involved). Such a key innovation was named as 'convergent', and was considered as a driver rather than a precursor, which would be characterized by a single evolution of a trait and its persistence throughout the radiation. As an adaptive trait, a key innovation may allow a taxon to benefit from key opportunities (opening new environmental F I G U R E 1 Global digital elevation model (GTOPO30) and mountain systems included for the evaluation of vascular plant species diversity and geodiversity/climate. Colours indicate the main biome for each mountain system: tropical (green), subtropical (orange), or temperate (dark blue) spaces fostered by geological or climatic modifications; Hughes & Eastwood, 2006;Uribe-Convers & Tank, 2015;Winkworth, 2004), the interaction between key innovations (intrinsic/trait) and opportunities (extrinsic/environment) potentially leading to radiations.
These drivers or precursors were then scored as either formally 'tested' or 'suggested', the former implying the use of some targeted analysis. Typically, drivers and precursors were deemed 'tested' if higher diversification rates were associated to them (e.g. key innovations), whereas they were considered as 'suggested' when only a temporal correlation was found (e.g. timing of diversification vs. timing of mountain uplift or climatic fluctuations). Biogeographical analyses highlighting dispersal into new regions or biomes (e.g. Páramo) and associated with consecutive increased diversification allowed a key opportunity to be deemed as 'tested'. Moreover, biome or vegetation shifts (during or prior to radiation) were considered as 'suggested' unless ancestral character reconstructions were provided and shifts in biome or vegetation were associated with shifts in diversification rates. In the absence of diversification rate analyses (Table S1.3), drivers or precursors were classified based upon the evidence and analyses presented. For example, niche shifts were scored as 'tested' if the differences in niche dimensions were presented within a phylogenetic framework. Table 1 provides an overview of a number of studies suggesting or testing the influence of drivers and precursors of diversification for different mountainous areas. Climatic drivers or precursors may include climatic variations (Pleistocene) and modifications over a large time span (e.g. aridification, Miocene cooling). Geological drivers or precursors include mountain uplift and sometimes ruggedness. Biotic drivers and precursors may refer to key innovations and opportunities, pre-adaptations, pollinator shifts, hybridization, biome and niche shifts, and shifts in vegetation distribution.

| Global determinants of mountain plant species diversity
We compiled species diversity estimates for 16 mountain systems around the world, including all continents except Australia and Antarctica, revealing a variation in plant species richness of more than one order of magnitude. The highest plant species diversity was recorded in the Tibeto-Himalayan region, the Andes and the Albertine Rift mountain systems, whereas the lowest species numbers were found in the temperate mountain ranges of Europe (Table S1.1). It should be noted that the Tibeto-Himalayan region as delineated here comprises the Qinghai-Tibet Plateau as well as parts of the Himalayas and Hengduan Mountains, thus combining both temperate and subtropical elements. We used 2-D area (i.e. area of the polygon) for the estimation of species diversity, which could potentially underestimate the true area depending on the slope of the mountain range (Wang, Chen, & Chen, 2017).
However, given that we exclusively compare mountain regions, we consider this effect negligible.
The best-fit GLM combined three variables (NPP, Elevation, Temperature anomaly) to explain plant species richness in mountain ranges (Table S1.2, AIC = 7037.1, R 2 MF = 0.626). Importantly, the temperature differences between the LGM and today (ΔT) had a stronger effect than either of the two contemporary environmental predictors (NPP, Elevation; Table 2). The GD index and Elevational range were strongly correlated (r = .70) and given that Elevational range showed better performance in single predictor models (Table S1.2), GD was not considered in the multivariate regression models. This was further corroborated by the variance inflation factor in a GLM combining all variables (Elevational range: 6.605; GD: 7.407). One important difference might be that the harmonization and categorization of the environmental variables for the GD index may be associated with some loss of information (in contrast to the elevation data).   Colombia, Figure 2c).

| Timing of mountain plant radiations across the globe
We found 47 studies investigating plant radiations in mountains, including 30 that tackled diversification dynamics, the majority of which detected shifts in diversification rates (Table 1) In a very few cases, as in the family Ericaceae (Table S1.

| Geodiversity alone is not sufficient to predict mountain plant diversity-the role of climatic fluctuations
We F I G U R E 2 Relative frequency distribution of available present-day temperatures during the Last Glacial Maximum (LGM) for selected mountain systems, indicating the difference between past and present climate. Present-day mean annual temperatures were extracted for each study area (right side of each plot) and compared to the distribution during the LGM (left side) including a 200 km buffer region to account for the availability of suitable temperatures in neighbouring lowlands (plots for all other mountain systems can be found in Appendix S2) (a) (b) (c) Kissling et al., 2012), despite evidence for the influence of historical contingency from both experimental and comparative evolutionary biology (Blount, Lenski, & Losos, 2018). Assuming niche conservatism (Losos, 2008)

| Interaction between climatic changes and mountain topography modifies the effect of geodiversity on biodiversity
The comparison of past and present-day climate (temperature) suggests that the interaction between climatic changes and mountain topography may further modify the effect of geodiversity on biodiversity (e.g. the 'species-pump' effect; see also Muellner-Riehl, 2019).
The topography of mountains strongly affects habitat availability in different elevational belts (Elsen & Tingley, 2015), which in turn im-  Schwery et al. (2015) are not shown. Whenever only age ranges were given, we calculated mean age between maximum and minimum of this age range. Climatic changes are represented by deep-sea oxygen-isotope records, modified from Zachos, Dickens, and Zeebe (2008). THR, Tibeto-Himalayan region for radiations, because vertical displacement would lead to less variation in connectivity (Flantua et al., 2019). In contrast, strong fluctuations (such as between the LIG and LGM) would decrease the probability of persistence of suitable habitats within the mountain range, leading to a stronger filtering effect (extinction) ultimately characterized by fewer radiations detectable today. In addition, such a filtering effect should be stronger for mountain systems with less extensive elevational range. Indeed, mountain systems characterized by small elevational ranges and strong modifications of temperature profiles (Figure 2 and Figure S2.1) appear to harbour fewer radiations (e.g. European mountains, North American mountains; Figure 3). In contrast, mountain systems with the largest elevational ranges and stronger overlap between today´s and LGM temperature profiles are also those where most radiations were identified (e.g.

Neogene and Pleistocene fostering a 'species pump' effect, as predicted by the 'mountain-geobiodiversity hypothesis'
Our compilation of studies on plant radiations in mountain systems around the world revealed that shifts to higher diversification rates (or starts of radiations) show a concentration from the late Miocene towards the Pleistocene (Figure 3; Tables S1.3 and S1.4). Globally, temperatures have been falling since the Eocene, with climatic cycles becoming more pronounced from the late Neogene (from 7 to 8 million years) onwards, and intensifying in the course of the Pleistocene. processes. In many cases, mountain uplift appears as the 'go-to' interpretation, particularly in the absence of diversification rate analyses.
For example the diversification of Rhodiola (Zhang, Meng, Allen, Wen, & Rao, 2014) was suggested to be driven by geological processes in the Tibeto-Himalayan region, whereas its peculiar habitat (incl. high alpine screes and rock crevices) and morphology (incl. traits related to Crassulacean metabolism) also suggest a potential role of intrinsic processes. In contrast, studies stemming from the Andes tend to highlight the role of climatic variations, intrinsic drivers (e.g. key innovations) and precursors (e.g. key opportunity, pre-adaptations; e.g. in Lupinus, Drummond, Eastwood, Miotto, & Hughes, 2012) more often than geological processes (but see review by Luebert & Weigend, 2014). This difference may find its source in mountains system-specific discrepancies within the scientific community.
Many types of biotic factors have been suggested as drivers or precursors of radiations in mountain systems. Some are strongly connected to morphology or physiology. For example, key innovations may include iteroparity for Lupinus in mountains of the New World (Drummond et al., 2012), berry-like fruits in Tripterospermum (Matuszak, Favre, Schnitzler, & Muellner-Riehl, 2016) and the Andean Campanulaceae (Lagomarsino et al., 2016), or the cushion life form typical of some Saxifraga (Ebersbach et al., 2017) and Androsace (Roquet, Boucher, Thuiller, & Lavergne, 2013) Knox & Palmer, 1995). Similarly, key innovations may help a taxon to access a newly formed habitat, such as the wood structure of Alchemilla in the Afroalpine (Gehrke, Kandziora, & Pirie, 2016). Some studies also suggest that further biotic drivers may contribute to diversification, such as hybridization   Najwer et al., 2016).

| Challenges encountered in mountain-
The

| CON CLUS IONS
About 200 years ago, Alexander von Humboldt for the first time comprehensively synthesized the observations and ideas about the interplay of biotic and abiotic nature at his time, supported by countless researchers over almost two decades, in his monumental 'Kosmos' (Humboldt, 1845(Humboldt, -1850. His expedition to mount Chimborazo, and his journeys to other mountain regions around the world, led him to develop the first map of biodiversity, the famous 'Naturgemälde' (Humboldt, 1807), which was of fundamental importance to his later work and general insights into the 'net of life' (Wulf, 2016). Following his approach of exploring the impact of geodiversity on biodiversity on mountains in particular, our study constitutes an up-to-date global comparison of mountain geodiversity and vascular plant diversity, incorporating ecological, evolutionary and geological processes, taking advantage of methodological advances and a wealth of the most recent case studies on organismic evolution.
Our study highlights that the complexity of plant evolution in mountain systems is not only determined by contemporary abiotic and biotic factors, but also influenced by historical factors, which need to be better integrated. In addition, the role of mountain topography in affecting species distributions and the dynamics of speciation and extinction under climate change need to be considered. Ultimately, only a standardized approach, both on the geodiversity as well as the biodiversity side, will lead to a comparable global-scale application, and a similar recognition of the importance of geodiversity to nature and the human population (incl. ecosystem services; compare Reynard & Brilha, 2018) as is nowadays recognized for biodiversity.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data on which this paper is based are in part submitted to Dryad (GIS Shapefiles, Geodiversity summary statistics, MODIS MOD17 Net Primary Production (NPP) data, R script; https://doi. org/10.5061/dryad.v1kh40j), and in part listed in the Appendix S1 (Tables S1.1., S1.3. and S1.4). The datasets derived are found in the Appendix S1 (Table S1.2), as well as in Tables 1 and 2.