Influence of climate, soil, and land cover on plant species distribution in the European Alps

Although the importance of edaphic factors and habitat structure for plant growth and survival is known, both are often neglected in favor of climatic drivers when investigating the spatial patterns of plant species and diversity. Yet, especially in mountain ecosystems with complex topography, missing edaphic and habitat components may be detrimental for a sound understanding of biodiversity distribution. Here, we compare the relative importance of climate, soil and land cover variables when predicting the distributions of 2,616 vascular plant species in the European Alps, representing approximately two-thirds of all European flora. Using presence-only data, we built point-process models (PPMs) to relate species observations to different combinations of covariates. We evaluated the PPMs through block crossvalidations and assessed the independent contributions of climate, soil, and land cover covariates to predict plant species distributions using an innovative predictive partitioning approach. We found climate to be the most influential driver of spatial patterns in plant species with a relative influence of ~58.5% across all species, with decreasing importance from low to high elevations. Soil (~20.1%) and land cover (~21.4%), overall, were less influential than climate, but increased in importance along the elevation gradient. Furthermore, land cover showed strong local effects in lowlands, while the contribution of soil stabilized at mid-elevations. The decreasing influence of climate with elevation is explained by increasing endemism, and the fact that climate becomes more homogeneous as habitat diversity declines at higher altitudes. In contrast, soil predictors were found to follow the opposite trend. Additionally, at low elevations, human-mediated land cover effects appear to reduce the importance of climate predictors. We conclude that soil and land cover are, like climate, principal drivers of plant species distribution in the European Alps. While disentangling their effects remains a challenge, future studies can benefit markedly by including soil and land cover effects when predicting species distributions.

Abstract. Although the importance of edaphic factors and habitat structure for plant growth and survival is known, both are often neglected in favor of climatic drivers when investigating the spatial patterns of plant species and diversity. Yet, especially in mountain ecosystems with complex topography, missing edaphic and habitat components may be detrimental for a sound understanding of biodiversity distribution. Here, we compare the relative importance of climate, soil and land cover variables when predicting the distributions of 2,616 vascular plant species in the European Alps, representing approximately two-thirds of all European flora. Using presence-only data, we built point-process models (PPMs) to relate species observations to different combinations of covariates. We evaluated the PPMs through block crossvalidations and assessed the independent contributions of climate, soil, and land cover covariates to predict plant species distributions using an innovative predictive partitioning approach. We found climate to be the most influential driver of spatial patterns in plant species with a relative influence of~58.5% across all species, with decreasing importance from low to high elevations. Soil (~20.1%) and land cover (~21.4%), overall, were less influential than climate, but increased in importance along the elevation gradient. Furthermore, land cover showed strong local effects in lowlands, while the contribution of soil stabilized at mid-elevations. The decreasing influence of climate with elevation is explained by increasing endemism, and the fact that climate becomes more homogeneous as habitat diversity declines at higher altitudes. In contrast, soil predictors were found to follow the opposite trend. Additionally, at low elevations, human-mediated land cover effects appear to reduce the importance of climate predictors. We conclude that soil and land cover are, like climate, principal drivers of plant species distribution in the European Alps. While disentangling their effects remains a challenge, future studies can benefit markedly by including soil and land cover effects when predicting species distributions.

INTRODUCTION
Understanding how environmental factors influence plant species distributions and diversity has been a longstanding challenge in ecology (de Candolle 1820, MacArthur 1965, von Humboldt and Bonpland 2010. In particular, it is recognized that climate has direct effects on plants' growth and physiology (Woodward 1987(Woodward , K€ orner 2003, which explains the imprint of past climate changes on extinctions, distributions, and speciation (Comes and Kadereit 1998, Mittelbach et al. 2007, Smycka et al. 2017. Given these direct effects, and the ever-increasing availability of high-resolution climatic data, climate has always been the preferred environmental predictor set to explain spatial patterns of species and diversity (Zimmermann et al. 2009, Mod et al. 2016. Nevertheless, other environmental features also influence plant species distributions.
A species' ecological niche is usually defined as a multi-dimensional space, which not only considers climate, but also other covariates such as soil or habitat structure (Hutchinson 1957). Soil indeed plays an important ecophysiological role for plants through the provisioning of nutrients and water, and its small-scale influence on plant distribution, in terms of variations in chemical and physical soil properties and bedrock types, is well documented (K€ orner 2003, Coudun and G egout 2005, Coudun et al. 2006, Piedallu et al. 2011, Bertrand et al. 2012, Dubuis et al. 2013, Scherrer and Guisan 2019. Land cover is also considered relevant in models aiming at predicting species distributions, as it adds realistic small-scale information on habitat fragmentation and human influence, which are not represented in other more commonly used bioclimatic variables (Von Holle and Motzkin 2007, Bradley and Fleishman 2008, Hopkins and Whiles 2011, Wilson et al. 2013, Cord et al. 2014, Tasser et al. 2017). Yet, soil and land cover are often neglected in research that evaluates spatial patterns of species and diversity (Thuiller et al. 2014, Mod et al. 2016, Scherrer and Guisan 2019. Soil has likely been included less often because relevant information is rather difficult to obtain at high spatial resolution. Land cover is often considered an indirect effect, and thus less prioritized in analyses of species distributions, or neglected due to its usual high correlation with climate ). However, land cover is decisive to realistically depict the impact of global changes on ecosystems, especially when predicting at fine-grained spatial scales. Therefore, combining climate, soil and land cover in plant distribution studies should provide new insights and an improved conceptual understanding of plant species distributions.
Several studies have investigated the importance of climate and soil on plant distributions (Lobo et al. 2001, Coudun et al. 2006, Hageer et al. 2017, Buri et al. 2020, plant abundance (Walthert and Meier 2017) and vegetation growth (Chakraborty et al. 2019). However, such studies have so far mainly focused on trees, analyzed small numbers of species or ignored land cover effects. In addition, these studies have not analyzed how the importance of predictors varies across space, although such information allows for a better understanding of spatial constraints on species distributions. To summarize, while climate, soil and land cover are recognized as important drivers, their respective influence on plant species distributions, and the spatial variation thereof, remain largely unexplored.
Mountain ecosystems represent natural laboratories for exploring such environmental effects on species distribution since they offer complex environmental conditions (K€ orner 2007). Abiotic gradients are steep and include extreme conditions; soil properties show considerable small-scale variations Guisan 2019, Buri et al. 2020) and land cover is particularly structured by the complex patterns of different human land uses and land-abandonment histories (Lackner and Psenner 2006, Becker et al. 2007, Gehrig-Fasel et al. 2007). Climate change is expected to have profound impacts on mountain ecosystems (Pauli et al. 2012, Rogora et al. 2018, which are centers of plant diversity and endemism (Smycka et al. 2017). Unravelling the combined role played by climate, soil and land cover in shaping plant species distributions in mountain ecosystems is, therefore, of principal importance for sound assessments of global change impacts.
Here, we explored the importance of climate, soil, and land cover in predicting the distributions of 2,616 vascular plant species across the entire European Alps (Fig. 1), encompassing two-thirds of the vascular plant species inventoried in this mountain range (~4,000 species). We investigated the relative influence of the three drivers along elevation gradients and explored plants' modeled responses along main climate, soil, and land cover gradients for 145 specific species grouped by altitudinal preferences. All influencing factors were obtained from point-process models (PPMs) by applying an innovative "predictive partitioning analysis." Based on these analyses, we tested the following hypothesis: (1) We hypothesize climate to be the dominant driver of plant distributions, although possibly with a decreasing influence with increasing elevation as climate becomes more extreme at high altitudes, and less spatially variable due to reduced surface area.
(2) We hypothesize soil to be the second most important driver with increasing importance with altitude. This expectation is based on the view that bedrock types are more variable and have a more direct effect on soils at high elevation. Soils range from very raw to well-developed, and are more directly dependent on topographic exposure and bedrock type at high elevation, leading to widely differing soil chemistries.
(3) Finally, we hypothesize land cover to be the least influential driver, with decreasing importance with elevation mainly due to a progressive reduction in habitat fragmentation and human impact, resulting in a lower modifying effect from human activities on species patterns relative to the effects from climate and soils.

Study area and species observations
This study was conducted across the European Alps, encompassing an enlarged version of the official Alpine Convention perimeter (Permanent Secretariat of the Alpine Convention 2009; Fig. 1). The perimeter was extended by considering Switzerland entirely, and by adding two French departments, i.e., Ain and Bouchesdu-Rhône. The total area covers 294,994 km 2 and encompasses a large range of environmental conditions and landscape types, including continental, oceanic, and mediterranean climates, and elevations ranging from À10 to 4,810 m above sea level.
The data on vascular plant observations used in this study was compiled from more than 210 individual sources (for a complete list, see Appendix S1: Table S1), the most important being: the National Data and Information Centre on the Swiss Flora (InfoFlora;~48%), the French National Alpine Botanical Conservatory (CBNA;~19%), the French National Mediterranean Botanical Conservatory (CBNMED;~5%), and information from the Global Biodiversity Information Facility not originating from any of the other sources considered (GBIF, available online;~2%). 4 To merge all data sets, we unified the species taxonomy using Flora Alpina (Aeschimann et al. 2004) as a backbone and checked for synonymy across different data sources (see Appendix S1: Text S1). GBIF provides a huge amount of geo-referenced species observations. Yet, many observations have considerably imprecise coordinates, duplicated records, or misleading coordinates from raster cell centroids (Gaiji et al. 2013). Strict filtering was thus applied to the compiled data to remove imprecise observations (for more details on the procedure, see Appendix S1: Text S1).
The pooled set of observations was filtered by means of two criteria: (1) for each species, observations were reduced to one observation within each 1 9 1 km grid cell to reduce pseudoreplication of species presences and (2) only species with at least 60 filtered observations across the study area were kept. Our final observational data set comprised 2,636,017 unique entries for 2,616 species (Fig. 1). For each, we extracted additional information from Flora Alpina about elevation preferences (colline, montane, subalpine, alpine, and nival species), vegetation types (herbaceous, shrub, and tree species), and geological preferences (calcareous, siliceous, and mixed species; see Appendix S1: Table S2, S3 for numerical summaries).

Environmental data
The candidate list of environmental predictors used for analysis and modeling of species distributions consisted of 27 climate, 25 soil, and 8 land cover variables. All predictors were kept at or aggregated (i.e., geographical mean) to a 1-km spatial resolution and projected to the standard Lambert Azimuthal Equal Area projection for Europe (EPSG:3035).
Climate variables.-Climate variables (means, covering the time period 1979-2013) were extracted from the Climatologies at High resolution for the Earth's Land Surface Areas (CHELSA) portal (Karger et al. 2017; available online). 5 From this data portal, we downloaded the 19 bioclimatic predictors as well as layers related to snow cover and plant phenology (for a detailed list, see Appendix S1: Table S4). Based on six original bioclimatic variables, three additional layers were generated to better represent seasonal variations, namely, summer to winter quarter differences in average daily mean temperatures (bio10Àbio11), summer to winter quarter differences in average daily precipitation sums (bio18Àbio19), and annual range in monthly precipitation sums (bio13Àbio14).  Martin 1987, Duscher et al. 2015. From ESDAC, we included only predictors of ecological relevance for plant species with full spatial coverage, namely, the 3D Soil Hydraulic Database of Europe (T oth et al. 2017), the European Soil Database derived data (Hiederer 2013a, b) and the Topsoil Organic Content for Europe (Jones et al. 2005). SoilGrids variables were chosen based on the same criteria. In total, SoilGrids and ESDAC provided a list of soil predictors describing different soil properties, including pH, organic carbon content, depth available to roots/bedrock, soil texture, water content and bulk density (for a detailed list, see Appendix S1: Table S4). Where information was available for topsoil and subsoil, we calculated profile averages before spatially aggregating the data to a 1-km resolution. Additional bedrock type layers were generated based on the two soil lithology polygon maps GLiM and IHME1500. The layers were constructed in three steps for each source. In the first step, we classified soil lithology classes into three basic categories, namely, calcareous, siliceous, and mixed (see Appendix S1: Table S5, S6). In the second step, the newly classified map was converted to a 100-m grid for each category. Finally, we aggregated the 100-m grids to a 1-km resolution by calculating proportions of the three basic categories per 1-km pixel.
Land cover variables.-Predictors relevant for habitat structure at the landscape scale were taken from the CORINE Land Cover 2012 map available at a 100-m spatial resolution (Bossard et al. 2000; data available online). 8 Land cover classification is available at three levels of detail. We generally kept the least detailed level 1 but used level 2 details to split the level 1 category "forest and seminatural areas" into "forests," "grasslands," and "open spaces with little or no vegetation." This led to a total of seven classes: "built-up areas," "agricultural areas," "forests," "grasslands," "open spaces with little or no vegetation," "wetlands," and "water" (see Appendix S1: Table S7, S8). Similarly to soil variables, we aggregated the 100-m grid to a 1-km resolution by calculating proportions of each of the seven classes per 1-km pixel. Additionally, we added a remote-sensingbased 1-km raster of vegetation productivity (MODIS Enhanced Vegetation Index [EVI] 2000[EVI] -2017 to the land cover variables (data available online. 9 For each of the 18 yr, we extracted the layers of the summer months June, July, August, and September, and averaged them across all layers to generate an 18-yr mean EVI estimate for summer.

Species distribution modeling
Variable selection.-For each of the three sets of environmental variables, a separate principal component analysis (PCA) was computed among all variables to obtain PCA axes that summarize the main trends per predictor class climate, soil, and land cover (for a description of the PCA axes and correlations with original predictors, see Appendix S1: Figs. S1-S3, Text S2). We kept six, nine, and six PCA axes, respectively, for climate, soil, and land cover for further analyses in order to cumulatively explain >85% of the total variance per predictor class, as well as to keep a balanced number of axes within each group. Furthermore, to avoid collinearity among predictors, we ascertained in the follow-up analyses that the predictor inter-correlation was |r| < 0.7 among the selected principal components of any category (Appendix S1: Fig. S4). Fewer axes for each predictor class could have been chosen, although we would have missed the representativeness of environmental conditions that some axes might bring despite their smaller variations as a result.
Model framework.-For each species, we generated 10 random sets of predictors by sampling (without replacement) our 21 PCA axes to account for model replicates (Fig. 2), and balanced environmental representativeness. Each set consisted of two PCA axes per environmental predictor class (climate, soil, and land cover) resulting in six PCA axes per set, while ascertaining that (1) each set had a unique combination of axes and (2) every PCA axis per set was considered at least once. Such a sampling procedure of PCA axes was implemented to avoid the arbitrary choice of specific covariates that might neither be relevant nor ensure equal representativeness of all three driver types. Although the direct effects of climate, soil and land cover on plant distribution cannot be measured as their exact influences are unknown for most plants, the selected PCA axes stand for proxy effects of all abiotic information available. Therefore, our sampling procedure of PCA axes as covariates assured both the inclusion of this proxy information, as well as a balanced representation of all PCA axes among and between the three driver types. For each predictor set and species, seven models were calibrated, namely, (1) FULL (climate, soil, land cover; six PCA axes), (2) ClimSoil (climate, soil; four PCA axes), (3) ClimLand (climate, land cover; four PCA axes), (4) SoilLand (soil, land cover; four PCA axes), (5) Clim (climate alone: two PCA axes), (6) Soil (soil alone: two PCA axes), (7) Land (land cover alone: two PCA axes; Fig. 2).
Our observational data set originated from a range of different sources that often lack sampling design. We therefore assumed to have bias in our presence observation sampling. This problem is known to be an issue FIG. 2. Summary of the modeling framework. For each of the 2,616 species, we calibrated seven different models representing varying combinations of two PCA predictors out of the three drivers (climate, soil, land cover). Each of these seven models was repeated 10 times by varying combinations of predictors from the three predictor sets. Model performance was assessed by fivefold environmental block cross-validation. In total 183,120 models were calibrated. PPM, process-point model; CLIM, climate. (Edwards et al. 2006, Albert et al. 2010, especially when using presence-only data (Stolar andNielsen 2015, El-Gabbas andDormann 2018). In order to correct for sampling bias, every model run was implemented with two potential observer bias covariates (following Warton et al. 2013): "Distance to Roads" and "Distance to Cities" (built with OpenStreetMap and the Global Roads Open Access Data; available online). 10,11 Both covariates had an intercorrelation of |r| < 0.7 with the principal components of climate, soil, and land cover (see Appendix S1: Fig. S4).
Model calibration.-We employed point-process models (PPMs) for our analyses; an algorithm specifically designed for species presence-only data, that models density of observations as a log-linear function of the environmental covariates (Renner et al. 2015). The algorithm, via a Poisson-GLM process, estimates the species' model response in environmental space based on density values, that are here approximates by aggregation of repeated 1-km observations along the environmental gradients. No background points or pseudo-absences required selection; rather, we used "quadrature points" to maximize the model log likelihood (Renner 2013, Renner et al. 2015. Each of the seven models was calibrated independently as a "down-weighted Poisson regression" (DWPR, following Renner et al. 2015) using a generalized linear model (McCullagh 1984) with second-order polynomials, elastic net regularization (Zou and Hastie 2005) and fivefold environmental block cross-validation (BCV: see Fig. 2 and Model evaluation). In order to have a consistent setup, we set the number of quadrature points to a high value (290,316), by extracting them on a regular 1-km mesh across the study area, rather than identifying individual minima with diagnostic tools (Renner andWarton 2013, Renner et al. 2015).
Down-weighted Poisson regression with elastic net was executed using the R package glmnet (Friedman et al. 2010). Elastic Net represents a type of regularization and variable selection mixing lasso and ridge regression approaches (Zou and Hastie 2005). It penalizes non-relevant environmental predictors that might lead to overfitting by shrinking their effects or removing them completely. Elastic net regularization requires two parameters: alpha (a), which sets the balance between lasso and ridge, and lambda (k), which sets the penalty coefficient level. For each BCV fold, we set a to 0.5, and determined the optimal k by testing 100 different values (see cv.glmnet function) and selecting the one for which model fit performed best under a new fivefold cross-validation. Elastic net regularization was only applied to PCA predictors, ensuring that the two observer bias covariates were not down-weighted or removed during model calibration.
Model evaluation.-Fivefold environmental BCV generates preliminary blocks of similar environmental conditions either spatially (Roberts et al. 2017), or according to observational patterns (Brun et al. 2019), to partition calibration (fourfold) and evaluation folds (onefold). For each model, (1) an independent PCA was computed on its environmental (PCA) predictors to summarize environmental conditions via the two first principal components and (2) species observations were then partitioned into 10 blocks of 5 environmentally stratified folds (2 blocks per fold), using the two principal components, observation coordinates and partitioning around medoids (PAM) clustering (following Brun et al. 2019). This allowed the number of observations to be environmentally balanced within each independent fold.
The BCV was adapted such that the partitioning was applied to presences only. Model performance was then evaluated for each fold by projecting the model calibrated on the remaining four folds to the study area, and by estimating the Boyce Index (BI) spatially. The Boyce index is considered one of the most appropriate metrics to assess the performance of presence-only models (Hirzel et al. 2006, Guisan et al. 2017. The index ranges from À1 to +1, with positive values indicating projections in accordance with the observations from the validation data, values around zero indicating random agreement, and negative values indicating systematic errors (Hirzel et al. 2006, Thuiller et al. 2018. Partitioning sub-model contributions.-We used BIs from our models to run a predictive partitioning analysis (based on Borcard et al. 1992), which allowed the independent single contribution of each driver to be separated from the total contribution of the BI FULL model predictive accuracy. Our approach differed from traditional variation partitioning approaches (Heikkinen et al. 2004) in that we used measures of model evaluation (i.e., BI) rather than model fit (e.g., R 2 ). For each species, we thus computed the partitioned influence of climate, soil, and land cover in predicting species' spatial distributions, by subtracting from the Boyce index of the FULL model, each of the other three models' indices (i.e., BI CLIMp =BI FULL À BI SOIL-LAND ; BI SOILp = BI FULL À BI CLIM-LAND ; for BI LANDp = BI FULL À BI CLIM-SOIL ). As a result, the absolute difference between the FULL model and the model lacking a specific driver represents the independent contribution of this driver to the prediction accuracy of the FULL model. We quantified both absolute and relative (to the full model) influences in our analyses. The relative influence is represented by the percentage of contribution of each driver effect to the sum of all effects (i.e., BI ALL = BI CLIMp + BI SOILp + BI LANDp , R CLIMp = BI CLIMp /BI ALL 9 100, R SOILp =BI SOILp /BI ALL 9 100, R LANDp = BI LANDp /BI ALL 9 100).
Model projections.-All model variants were projected across the entire European Alps at 1 9 1 km, and corrected for observer bias by setting the two bias covariates to 0 following the recommendations of Warton et al. (2013). Furthermore, to assess the species response along environmental gradients in more detail, we predicted the FULL model of a subset of 145 species (selected to balance strict elevation preferences; see Appendix S1: Table S9) along the first PCA axis of each of the three drivers while holding the other predictors constant at their mean/median (following Elith et al. 2005, Zurell et al. 2012. To that end, we considered FULL models to be those that included at least one first PCA axis and had a BI FULL >0.4. Each of the three axes/drivers was then represented by the median species response with species being grouped by elevation preferences. Cross-species analyses.-The following analyses were also performed for models with BI FULL >0.4. The relationships between species' elevation distribution and relative or absolute influences of prediction accuracy (of climate, soil, and land cover) were analyzed using second-order polynomial linear regressions. For each species, the elevation values used in this analysis were extracted by overlaying species observations with the Digital Elevation Model over Europe from the GMES RDA project (EU-DEM,). 12 We kept the 95th percentile of elevation values as a reference, so that the partitioning of species along the elevation gradient is less influenced by biased elevation sampling and better represents their ecological cold limit.
The spatially resolved relative influence of each predictor category was calculated with the averaged FULL model projections. The procedure was undertaken for all 2,616 species, and for the subset of 145 species, in order to prevent potential patterns arising from the study's sampling bias. For each species, the retained models (BI FULL > 0.4) among the 50 model versions per species (10 replicates 9 5 BCV) were projected to the study area. The average projected probability per species (among model versions) was then multiplied by the species' BI CLIMp , BI SOILp , and BI LANDp , generating for each species three grids representing the absolute influence of climate, soil, and land cover. For each driver, the species' layers were then averaged, resulting in three grids that represent the overall absolute influence of climate, soil, and land cover on plant species distributions across the European Alps. For each grid, the relative influence was then obtained by dividing each pixel value by the sum across the three grids, scaled to range from 0 to 255, and combined with an RGB composite.
Having disentangled the results of the relative influence along elevation bands (Fig. 4a), we found that although the contribution of climate remained highest across all elevation bands, it declined with elevation (decreasing from~62.2.0% to~41.0%). In contrast, the effect of soil increased progressively with elevation from 17.0% to~26.4% in the alpine belt and stabilized at higher elevations. The effect of land cover first decreased slightly from colline (~20.8%) to subalpine (~18.0%) and then strongly increased toward nival elevations (~34.0%). Analyzing the ungrouped relative effects on species along the 95th percentile of their elevation distributions (Fig. 4b) revealed similar trends in relationships, as seen in Fig. 4a with R 2 adj CLIMp~5 .24%, R 2 adj SOILp~2 .78%, R 2 adj LANDp~3 .38% (for more information on all plots and statistics, see Appendix S1: Fig. S5, Table S10).
The absolute contributions of climate, soil, and land cover predictors on plant species distributions along elevation gradients (Appendix S1: Fig. S6, Table S10) followed the same order, but the apparent hump-shaped effects appear to be influenced by higher FULL model predictions at low and high elevation (Appendix S1: Fig. S7). Moreover, the general trend of grouped and ungrouped relative effects for all 2,616 species was supported by the same trends for herbaceous, shrub, calcareous, siliceous, and mixed species (Appendix S1: Figs. S8, S9). Nevertheless, no significant relative influences were found for tree species, as this vegetation type only occurs at low to mid-elevation. Overall clearer relationships were found for species types with a higher number of species, i.e., herbaceous, calcareous, and mixed species.
The map of average contributions of the three drivers to the 145 species pool showed similar trends, as seen in Fig. 4, except in the southwestern European Alps (Fig. 5a-c). Climate mostly affected prediction accuracies at low elevation (outer rim and the non-center valleys of the study region), while land cover had the largest influence at high elevations (center of the study region). However, land cover was also important to some degree at low elevation showing, locally, the greatest influence in urban/built-up areas, as well as additional effects strongly mixed with soil in the center valleys of the European Alps (in yellow). The greatest effects of soil showed 12 https://www.eea.europa.eu/data-and-maps/data/eu-dem Article e01433; page 6 YOHANN CHAUVIER ET AL. Vol. 91,No. 2 strong mixed effects with climate overall. They were generally concentrated at mid-elevation, but a specific pattern is evident in the southwestern Alps. Finally, the map of average contributions to the 2,616 species pool (Appendix S1: Fig. S10) generally showed similar results. Nevertheless, climate displayed here less spatially relative influence than soil overall, although it remained higher in Switzerland and the southwestern Alps. The climate PCA axis 1 was found to be positively correlated with all predictors related to temperature and FIG. 4. Relationships between relative climate (blue), soil (red), and land cover (green) effects on prediction accuracy and species elevation in meters. The influence is represented (a) per elevation belt and (b) per species plotted at their 95th elevation percentile. In panel a, for each elevation band, three Wilcoxon paired-sample tests were performed to test whether the effects of R CLIMp , R SOILp , and R LANDp differ significantly among species. Obtained P values were adjusted for multiple comparisons applying a "false discovery rate" correction. Panel b shows similar trends to panel a and summarizes second-order polynomial linear relationships between the relative influences and species' 95th elevation percentiles. Statistics are detailed in Appendix S1: Fig. S5, Table S10. Colored shadings show 95% prediction interval (*P < 0.05, **P < 0.01, ***P < 0.001). negatively correlated with predictors related to snow and precipitation. The soil PCA axis 1 was positively correlated with predictors related to organic carbon, and negatively correlated with predictors related to soil pH, soil depth, and water content. The land cover PCA axis 1 was more positively correlated with predictors related to human land cover (urban/built-up areas and agriculture) and open spaces, and more negatively correlated with predictors related to forests and productivity (Appendix S1: Figs. S1-S3 for exhaustive correlations between all predictors and axes). Overall, analyses of the 145 species responses along the first PCA axis of the three drivers showed that plant species responded according to their altitudinal distribution. However, while this was complete for climate (Fig. 6a), this was less so for soil (Fig. 6b) and land cover (Fig. 6c).
(1) Responses of colline, montane, subalpine, alpine, and nival species followed successively decreasing values of the climate PCA axis 1.
(2) Responses of the same species along the soil PCA axis 1 showed that mainly colline and montane species responded toward positive PCA values.
(3) Responses along the land cover PCA axis 1 showed that nival and alpine species did not respond markedly at all, while colline, montane, and subalpine species showed clearer responses toward positive PCA values.

DISCUSSION
Our analyses demonstrate that other drivers besides climate need to be considered when modeling species distributions, especially in topographically heterogeneous regions. We not only measured the relative influence of each driver in the European Alps, but also how their respective importance varied with elevation. Although climate exerts the strongest influence across all elevation belts, its influence on species distributions was greater at low altitude and declined toward high elevations (Figs. 4, 5). This trend is most likely explained by an association among three different aspects. (1) As altitude increases, the diversity of climatic conditions gradually declines due to decreasing surface area (or number of habitats) at higher elevations, and to the extensive pattern of climate (Appendix S1: Fig. S11). Unlike soil, climate drives plant species occurrence through large-scale variability. Therefore, as elevation increases, the decrease in surface area shrinks climate variability and, thereby, climate influence. (2) Climate becomes more extreme toward higher elevations, increasingly limiting plant growth and distributions (Maestre et al. 2009, Dvorsk y et al. 2017, and progressively selecting cold-adapted species, especially endemics (Pauli et al. 2003). Coldadapted species are therefore less influenced by climatic conditions than by competition, which prevent them from occupying lower elevations. Biotic interactions are indeed generally known to have a high influence on plant species distribution (Meier et al. 2010), with competition more strongly limiting plant distribution toward lower elevations (i.e., in less stressful abiotic environments) and facilitation expanding the distribution toward higher and more stressful environments (Choler et al. 2001, Callaway et al. 2002, Maestre et al. 2009, Louthan et al. 2015. (3) Other non-climatic predictors (i.e., soil and land cover), which provide primarily smallscale information, gain importance in explaining the distribution of higher elevation plants (Figs. 4,5). The higher the elevation, the more geological barriers occur, preventing species dispersal to neighboring mountain ranges, and promoting endemism specific to soil, bedrock types, and habitats (McGraw and Levin-Madrono 1998, Brown 2001, Vetaas and Grytnes 2002, Pauli et al. 2003. In contrast to climate, our study revealed an increasing influence of soil properties with elevation (Fig. 4), along with increasing soil spatial variability (Appendix S1: Fig. S12). This trend reflects the fact that soil depth becomes shallower, while soil chemistry becomes more variable and more directly influenced by bedrock type and by the more complex topography compared to the lowlands (K€ orner 2003, Becker et al. 2007, Bales et al. 2011, Leuschner and Ellenberg 2017). Yet, this increasing effect remains fairly weak and stabilizes at mid-altitude (Fig. 4), which may reflect a lack of welldeveloped soils (i.e., high levels of organic carbon and humus) at the alpine and nival belts (Leuschner and Ellenberg 2017). Furthermore, the distinct soil influence pattern found in the southwestern European Alps (Fig. 5a) seems to originate from very heterogeneous soil conditions occurring within this region, as shown by the RGB composite map of soil PCA axes (see Appendix S1: Fig. S13). Still, we believe that more accurate and informative spatial explicit soil data sets are still needed to improve ecological analyses and modeling (Mod et al. 2016, Scherrer andGuisan 2019). We expect future studies to confirm our results or even reveal increased effects compared to climate, as soil maps with improved spatial and thematic resolution become available.
Land cover generally shows an increasing influence toward higher elevations (Fig. 4, 5), which appears to originate from a growing presence of a few dominating land cover classes that are not related directly to human activities (open spaces and grasslands; see Appendix S1: Fig. S14). Although this influence could be explained by an increase in land cover classes toward higher elevation, the modified level 1 land cover classification used in this study shows a clear opposing trend (Appendix S1: Fig. S14) that rejects such an assumption. A moderate local influence of land cover at colline and montane elevations was also visible (Fig. 4). Interestingly, while the influence of climate decreased with elevation, this trend was not as linear as we expected (Fig. 7a). Instead, the influence of climate reached a plateau at the colline and montane belts (Fig. 4a, b), possibly originating from an increasing human footprint on species redistributions at lower elevations (Nogu es-Bravo et al. 2008, Marini et al. 2009). This would cause climate to exert a lesser effect than expected on current distributions at the cost of increasing effects from man-made disturbances. COR-INE land cover classes that are indicative of artificial environments are indeed more abundant at colline elevation and fairly present at montane elevation (built-up and agriculture; see Appendix S1: Fig. S14), and our results show that artificial land cover clearly affects species distributions in the lowlands, especially urban areas (Fig. 5a vs. Fig. 5b). Nevertheless, our two classes of artificial land cover do not seem to fully capture this anthropogenic influence. Land cover influence at colline FIG. 6. Model responses of 145 species along PCA axis 1 of (a) climate, (b) soil, and (c) land cover. Each line in the main panel represents the median response of plant species of the same altitude preference. For more information on the description of PCA axes, see Appendix S1: Figs. S1-S3. For more information on the 145 species, see Appendix S1: Table S9. and montane elevations is only moderately higher than at the subalpine level, and much lower than higher elevations (Fig. 4a, b). It appears that using a classification of artificial land cover, as heterogeneous and complex as that of genuine CORINE land cover (see Appendix S1: Table S8), would have better captured human influence that is known to strongly govern species distributions (Lippitt et al. 2008, Schwoertzig et al. 2016. We therefore hypothesize that there are two opposing land cover effects along the elevation gradient. An anthropogenic effect, mainly at low elevation, and a natural one at high elevation that structures species' occurrences conditionally to climate. Consequently, and based on the relative influence of soil (see small panel in Fig. 7b based on Fig. 4b), four hypothetical climate influence trends may be distinguished in response to the prevalence of natural and anthropogenic land cover (Fig. 7b). Based on these trends, the influence of anthropogenic land cover gradually decreases toward higher elevations (Nogu es-Bravo et al. 2008), whereas natural land cover grows in importance. Although these inferences remain visual and conceptual, they pave the way for future studies to explicitly disentangle these two land cover aspects. It should be noted that, as an alternative approach, we could have computed four PCAs to obtain axes for Climate, Soil, Natural Land Cover, and Anthropogenic Land Cover, and run the overall analyses with these four sets of variables. While this might have permitted us to further disentangle land cover influence, such an approach would have required computation time and resources beyond what was available.
The assessment of influential drivers on single species distributions remains a key question in ecology since this is not known for many species. Regarding climate PCA axis 1 (Fig. 6a), species' responses follow the known gradient of ecological niches along elevation and confirm that species do not perform well outside their preferred climate niche. Findings on land cover PCA axis 1 (Fig. 6c) show that low elevation species (colline and montane) are positively associated with human-related land cover and low productivity sites. For soil PCA axis 1 (Fig. 6b), the responses of all species groups are more difficult to interpret, except for colline and montane species, which seem to respond strongly to highly acidic soils and organic carbon. This may be explained by mixed effects with land cover found at low elevations (in Fig. 5, yellow patterns in the center valleys), which could be associated with agricultural areas. Unfortunately, the increasing soil effect with elevation found in this study cannot be entirely reflected with only 145 species. It is possible that this effect is due to the relatively low resolution of our original soil layers (1 9 1 km), which mostly originate from sources that use model interpolation (ESDAC and SoilGrids). These sources likely lack the fine-scale variability of edaphic properties (Mod et al. 2016, Scherrer andGuisan 2019) that is known to be even higher in topographically heterogeneous regions. Furthermore, although our PCA methodology (i.e., 21 PCA axes with replicated sampling) effectively summarizes the full range of environmental conditions, the approach has its limits. First, overall, the interpretation of the 145 species model responses to climate, soil, and land cover PCA axes 1 remains very complex. Identifying the contributions of the original predictors (i.e., 27 climate, 25 soil, and 8 land cover variables) to species' responses is only feasible via correlation analyses and strong inferences. Second, the full extent of species' responses cannot be explained by PCA axis 1 alone, and interpretation may prove even more complex when analyzing species' responses on the remaining PCA axes.
Other possible limitations of the study are related to (1) potential coordinate imprecisions regarding species observations and (2) the obvious sampling bias toward Switzerland and the South-Western Alps (Fig. 1) that might misrepresent the spatially relative influence shown in Fig. 5. Regarding (1), species observations were compiled from 210 individual sources. Although such a high number of sources certainly allowed the large scope of our study, this might have additionally introduced imprecise or misplaced species observations. Regarding (2), observational bias was taken into consideration by FIG. 7. Conceptual framework for the relative influence of climate on plant species distributions with elevation. The figure is built around the hypothesized decreasing relative effect of (a) climate alone (with and without human impact) and (b) climate alone + full land cover (Full LC) with different LC effects singled out along elevation. Full LC draws the exact inverse trend of soil-relative influence (see inset based on Fig. 4b).
adding distance to roads and cities as two independent bias covariates in the model calibration. While the implementation of such covariates has resulted in general improvements in previous studies , applying this type of correction may only be partly efficient in mountain environments. Indeed, mountain accessibility is not only defined by roads and cities, which are present mostly at lower elevations (Yackulic et al. 2013), but also according to climate gradients and topography. Considering that two-thirds of our species occur at the colline and montane belts (Appendix S1: Table S2, S3), the use of such bias covariates was a logical choice, although certainly not entirely adapted to the study's environment. Future improvements might consider more representative accessibility covariates such as the raw sampling effort of the study area, as well as independent testing data sets (Stolar and Nielsen 2015), to improve estimates of relative influence and evaluation performances, respectively.
Interpreting the independent effects of climate, soil, and land cover on plant species distributions remains a complex task. These three predictor groups are known to interact in their effects on species distributions, and we have indeed shown that they are partly interdependent on each other (Fig. 3b). Climate and land cover both influence soil attributes, especially soil moisture and nutrients (Post and Kwon 2000, Teuling et al. 2007, Seneviratne et al. 2010, Macdonald et al. 2012, whereas land cover is intrinsically linked to climate , Yang et al. 2017. Climate is therefore the main driver responsible for structuring natural habitats and their specific properties such as soil. This complicates the interpretations and analysis when attempting to disentangle the three drivers (Pearson et al. 2004, Jetz et al. 2007). These potential interactions could have been explicitly addressed with other multivariate methods such as structural equation modeling (SEM). Implementing such a method would have potentially allowed the large unexplained shared effects (latent variables) between climate, soil, and land cover (measured variables) to be disentangled and interpreted (total unexplained BI SHARED ¼ BI FULL À ðBI CLIMp þ BI SOILp þ BI LANDp Þ = 0.43).
While our findings focus on the large system of the European Alps, questions remain regarding the generalization of similar influence patterns to other mountain systems or continents. In particular, analog ecological patterns are often demonstrated along both elevation and latitudinal gradients, such as decreasing plant biodiversity or increasing climate extremes (Stevens 1989, Weir and Schluter 2007, Schemske et al. 2009). Consequently, we may well ask whether the elevational trends we observed also exist from the equator to the poles. In principle, similar influence gradients might be found on the condition that soil heterogeneity and natural land cover increase toward higher latitudes. Given the increasing availability of plant observations from online collaborative platforms, and global high-resolution climate, soil, land cover, and remote sensing data (Drusch et al. 2012, Fick and Hijmans 2017, Karger et al. 2017, Kobayashi et al. 2017, we expect future studies to investigate such global trends of influence.