Diego Nieto-Lugilde

Jonathan Roger Michel

Henri Lenoir

Sylvain Abdulhak

David Aeschimann

Stefan Dullinger

Jean-Claude Gégout

Antoine Guisan

Harald Pauli

Julien Renaud

Jean-Paul Theurillat

Jonathan Lenoir

Wilfried Thuiller

Jérémie Van Es

Pascal Vittoz

Wolfgang Willner

Thomas Wohlgemuth

Niklaus E Zimmermann

Jens-Christian Svenning

Tree cover at fine and coarse spatial grains interacts with shade tolerance to shape plant species distributions across the Alps

Keywords: biotic interactions, biotic modifiers, facilitation, light competition, plant-plant interactions, shade tolerance, spatial grain, species distribution models Ecography (Cop.)

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

One important goal of ecology is to provide reliable forecasts of biotic responses to climate change [START_REF] Pearson | Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?[END_REF][START_REF] Meier | Biotic and abiotic variables show little redundancy in explaining tree species distributions[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF]. At the same time, obtaining a mechanistic understanding of the determinants of species distributions is still one of the greatest challenges [START_REF] Thuiller | A road map for integrating eco-evolutionary processes into biodiversity models[END_REF]). Climate has long been recognised as the main driver of species distributions [START_REF] Pearson | Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?[END_REF]. However, many recent studies have shown that other factors, such as dispersal limitation, remnant population dynamics, transient eco-evolutionary dynamics, environmental disturbances, edaphic conditions, and biotic interactions, may modify species-specific responses to climatic drivers [START_REF] Pounds | Biological response to climate change on a tropical mountain[END_REF], Root et al. 2003[START_REF] Parmesan | Ecological and Evolutionary Responses to Recent Climate Change[END_REF][START_REF] Lenoir | Going against the flow: potential mechanisms for unexpected downslope range shifts in a warming climate[END_REF][START_REF] Normand | Postglacial migration supplements climate in determining plant species ranges in Europe[END_REF][START_REF] Bertrand | Disregarding the edaphic dimension in species distribution models leads to the omission of crucial spatial information under climate change: the case of Quercus pubescens in France[END_REF]). Among these, biotic interactions have been recently highlighted as a key factor that should be better incorporated into both species distribution models (SDMs) and biodiversity models [START_REF] Kissling | Towards novel approaches to modelling biotic interactions in multispecies assemblages at large spatial extents[END_REF][START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF][START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF]. Although previously thought to be of highest importance at fine spatial grains (see Eltonian Noise Hypothesis in [START_REF] Soberon | Colloquium Papers: Niches and distributional areas: Concepts, methods, and assumptions[END_REF]Nakamura 2009, Cooper and[START_REF] Cooper | opinion: Habitat data resolution and the detection of species interactions[END_REF], recent studies have shown that biotic interactions can also significantly alter species distributions at coarse spatial grains and broad spatial extents [START_REF] Bullock | Geographical separation of two Ulex species at three spatial scales: does competition limit species' ranges?[END_REF][START_REF] Araújo | The importance of biotic interactions for modelling species distributions under climate change[END_REF][START_REF] Gotelli | Macroecological signals of species interactions in the Danish avifauna[END_REF][START_REF] Meier | Co-occurrence patterns of trees along macro-climatic gradients and their potential influence on the present and future distribution of Fagus sylvatica L[END_REF].

In plants, both correlative [START_REF] Pellissier | Species distribution models reveal apparent competitive and facilitative effects of a dominant species on the distribution of tundra plants[END_REF][START_REF] Meier | Biotic and abiotic variables show little redundancy in explaining tree species distributions[END_REF], 2011[START_REF] Kunstler | Effects of competition on tree radial-growth vary in importance but not in intensity along climatic gradients[END_REF], 2012[START_REF] Meineri | Modeling alpine plant distributions at the landscape scale: Do biotic interactions matter?[END_REF] and experimental studies [START_REF] Klanderud | Climate change effects on species interactions in an alpine plant community[END_REF][START_REF] Hector | Effects of Dominance and Diversity on Productivity along Ellenberg's Experimental Water Table Gradients[END_REF] have provided evidence for the importance of biotic interactions for species distributions. For example, shrub cover and abundance affect herb and dwarf shrub richness and abundance in tundra habitats [START_REF] Pajunen | Impact of shrub canopies on understorey vegetation in western Eurasian tundra[END_REF], and the co-occurrence patterns of competing trees alter the distribution patterns of European beech [START_REF] Meier | Co-occurrence patterns of trees along macro-climatic gradients and their potential influence on the present and future distribution of Fagus sylvatica L[END_REF]). However, the importance of taller species (such as trees) in driving herb and shrub species distributions across large geographic regions has never been explicitly assessed (see [START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF]) despite the well-known roles of light availability and shade tolerance in determining local plant community composition (Kobe et al. 1995, Valladares and[START_REF] Valladares | Shade Tolerance, a Key Plant Feature of Complex Nature and Consequences[END_REF]. Furthermore, due to the sessile nature of plants, direct plant-plant interactions are commonly thought to mainly occur within short distances, and there are thus few studies of the impact of plant-plant interactions at coarser spatial grains [START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF], but see [START_REF] Bullock | Geographical separation of two Ulex species at three spatial scales: does competition limit species' ranges?[END_REF].

Trees have been referred to as biotic modifiers [START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF]) because of their abilities to modify microclimatic conditions [START_REF] Pinto | Landscape attributes drive complex spatial microclimate configuration of Brazilian Atlantic forest fragments[END_REF]) and soil properties [START_REF] Thomas | Ecology of Woodlands and Forests: Description, Dynamics and Diversity[END_REF]. In areas with high diurnal and seasonal climatic variation, dense tree cover dampens the microclimatic variability in the understory compared with open areas (von [START_REF] Arx | Spatio-temporal effects of forest canopy on understory microclimate in a long-term experiment in Switzerland[END_REF]. Additionally, trees add organic matter to the soil, often increasing nutrient availability [START_REF] Augusto | Impact of several common tree species of European temperate forests on soil fertility[END_REF]. However, the most important effect of trees is that they reduce light availability in the understory. As some species are adapted for photosynthesis under low-light conditions (Kobe et al. 1995, Valladares and[START_REF] Valladares | Shade Tolerance, a Key Plant Feature of Complex Nature and Consequences[END_REF], such species may benefit from a dense tree cover, whereas other (shade-intolerant) species are excluded. Hence, trees may have varying impacts on co-occurring herb and shrub species, ranging from competition to facilitation, depending on the shade tolerance of the interacting species. Despite their sessile nature, the effects of trees may propagate to coarser spatial grains (e.g., 1 km 2 ) via plant population dynamics, by affecting colonisation probabilities and local extinction risks of herbs and shrubs. For instance, an open gap in a forested landscape could be too small to support a viable population of a light-demanding species [START_REF] Canham | Different responses to gaps among shade-tolerant trees species[END_REF] or too distant from source populations and/or too short-lived to be colonised. Similar processes may also occur for shade-tolerant species and small patches of dense tree cover in open landscapes [START_REF] Jamoneau | Drivers of plant species assemblages in forest patches among contrasted dynamic agricultural landscapes[END_REF]. These examples suggest the fine-grain occurrence of understory plants may be affected not only by tree cover at the same grain but also by landscape-grain tree cover through metapopulation-like dynamics.

Changes in the role and importance of biotic interactions along abiotic gradients have been discussed for decades, and we use this perspective as an additional way to study the effects of tree cover on herb and shrub distributions [START_REF] Paine | Intertidal community structure[END_REF][START_REF] Kaufman | Diversity of New World Mammals: Universality of the Latitudinal Gradients of Species and Bauplans[END_REF][START_REF] Maestre | Refining the stress-gradient hypothesis for competition and facilitation in plant communities[END_REF]. A frequent approach is to study their effects along the elevational ranges of the species [START_REF] Normand | Importance of abiotic stress as a range-limit determinant for European plants: insights from species responses to climatic gradients[END_REF][START_REF] Kunstler | Effects of competition on tree radial-growth vary in importance but not in intensity along climatic gradients[END_REF][START_REF] Roux | Biotic interactions affect the elevational ranges of high-latitude plant species[END_REF]. The asymmetricabiotic-stress-limitation hypothesis (AASLH) proposes that abiotic environmental and biotic interaction drivers constrain the upper and lower elevational limits, respectively [START_REF] Paine | Intertidal community structure[END_REF][START_REF] Normand | Importance of abiotic stress as a range-limit determinant for European plants: insights from species responses to climatic gradients[END_REF][START_REF] Ettinger | Climate determines upper, but not lower, altitudinal range limits of Pacific Northwest conifers[END_REF], assuming a stronger effect of competition in habitats not physiologically limited by abiotic conditions (i.e., the lower elevational limit). However, biotic interactions may ameliorate the limiting physiological conditions, thereby extending the upper elevational limits of certain species [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF]. The stress gradient hypothesis (SGH) predicts varying biotic interactions as environmental stress decreases, with an increase in competition and a decrease in facilitation [START_REF] Bertness | Positive interactions in communities[END_REF]Callaway 1994, Maestre et al. 2009), supporting the role of biotic interactions in influencing both upper and lower elevational limits of species ranges [START_REF] Defossez | Do interactions between plant and soil biota change with elevation? A study on Fagus sylvatica[END_REF].

Here, we (1) incorporate tree cover into SDMs (cf. [START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF] to measure its importance as a driver relative to climate and soil for fine-grain herb and shrub species distributions across a large region, the Alps, and (2) assess whether tree cover effects are affected by species-specific shade tolerances. Then, we (3) investigate whether tree cover effects occur only at fine spatial grains (plot grain: ~10 m) or if there are also detectable coarser grain (landscape grain: ~1 km) tree cover effects on fine-grain herb and shrub species distributions. Finally, we (4) evaluate whether tree cover affects herb and shrub elevational range limits, and if it does, whether the patterns are consistent with the AASLH or the SGH.

MATERIALS & METHODS

Study area

The study area covers the entire Alps (43°29'49" -48°53'00"N; 4°46'28" -17°04'00"E), as delineated by the Alpine Convention Boundary (Fig. 1; http://www.alpconv.org). Elevation ranges from sea level up to 4,810 m a.s.l., and mean elevation is approximately 1,037 m a.s.l. Mean annual temperature and precipitation range from -11.6 to 16.6 °C and from 515 to 2,883 mm yr -1 , respectively (WorldClim database; [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]. Soil conditions are diverse, with large areas of calcareous as well as siliceous bedrocks, and soil pH ranging from 3 to 8 (European Soil Portal; http://eusoils.jrc.ec.europa.eu/). The main vegetation types are forest, agricultural land, scrubland, grassland, and scattered vegetation on rocks, with forests covering 40.4 % of the study area (European Forest Data Center; http://forest.jrc.ec.europa.eu).

Data and study species

We compiled vegetation plot data from two different sources: (1) the Alps Vegetation Database (n = 35,735 plots from 1900 to 2009; [START_REF] Lenoir | The Alps Vegetation Database -a geo-referenced community-level archive of all terrestrial plants occurring in the Alps[END_REF]) and (2) a survey of the French National Alpine Botanical Conservatory at Gap (n = 1,782 plots from 1980 to 2009; CBNA). This combined database covers a large elevational gradient from 0 up to 3,412 m a.s.l., distributed across the entire Alps. We used Flora Alpina (Aeschimann et al. 2004) as the taxonomic reference when integrating the databases. Note that each vegetation plot is an exhaustive list of all vascular plant species co-occurring within a given location at a given point in time and thus provides reliable information on both presence and absence data for a given species.

To avoid mismatches between plot locations and the raster-based predictor variables (see below), we selected those plots with a spatial reference of at least 500-m accuracy. To reduce pseudo-replication, we randomly selected one plot among those with identical spatial coordinates (i.e., plots from time-series and spatially nested plots). To estimate plot-grain tree cover, we selected plots having abundance-dominance indices and excluded plots having presence-absence data only. The final dataset comprised 6,935 plots surveyed after 1980, with plot sizes ranging from 1 to 500 m 2 (see Appendix 1 in the supporting information for an assessment of the overall representativeness of the final dataset for all predictor variables relative to the original dataset). As our objective was to study the impact of tree cover on the distribution of lower-stature plant species, the smallest plots (1 m 2 ) were also relevant because they were the most common plot size for herbs and shrubs in open habitats. Consequently, we removed trees and tall shrub species-i.e., life forms coded 'p' and 'n-p' in Flora Indicativa [START_REF] Landolt | Flora indicativa[END_REF])-as well as 634 herb and shrub species not listed in Flora Indicativa. Finally, to limit noise in subsequent modelling analyses and to lessen the likelihood of overfitting due to the limited number of occurrences, we selected 960 herb and low-stature shrub species with more than 20 occurrences within the 6,935 selected plots. The total number of occurrences per modelled species ranged from 20 to 3,451, with an average of 264 occurrences.

We completed the database with descriptor and indicator values for most of the taxa from the Flora Indicativa database [START_REF] Landolt | Flora indicativa[END_REF]. [START_REF] Landolt | Flora indicativa[END_REF] ranked most of the terrestrial plant taxa in the Alps according to their optimal occurrence along key environmental gradients using an ordinal scale ranging from 1 to 5. For the purpose of this study, we used one plant trait descriptor for life forms (LF) and six plant indicator values for air temperature (T), continentality (K), soil moisture (M), soil reaction [soil pH] (R), nutrient availability (N), and light (L). LF describes the Raunkiaer life forms for plants and was used to select herbs, shrubs and trees. We classified species as shade tolerant (L = 1-2), semishade tolerant (L=3), or shade-intolerant (L = 4-5) and, according to their optimal elevational belt, as colline (T = 4-5), montane (T = 3-3.5), subalpine (T = 1.5-2.5), or alpine and nival (T=1).

Predictor variables

Because our objective was to investigate the effect of tree cover in SDMs on lower-stature plants at two different spatial grains (plot and landscape), we calculated two sets of predictor variables for climate, soil, and tree cover.

Plot grain-For five of the six indicator values (V: T, K, M, R, N) in a given plot, we computed the mean indicator value (V PMI ) by averaging the values of all species recorded in the respective plot. Note that the averaged values were not weighted by abundance (see Diekmann 2003 for a discussion of the use of weighted and unweighted averaged values). These climate and soil factors (T PMI , K PMI , M PMI , R PMI , and N PMI ) are well known to affect plant species distributions (e.g., [START_REF] Skov | Potential impact of climatic change on the distribution of forest herbs in Europe[END_REF]Svenning 2004, Bertrand et al. 2012) and were used as predictor variables for the plot-grain SDMs (~10-m grain). To avoid circularity when fitting SDMs for a given focal species, we computed V PMI after removing the indicator values of the focal species from all plant assemblages in which it occurred. Because the V PMI s were calculated based on the species observed in each plot, these values represent the vegetation at the time when each plant assemblage was recorded. T PMI characterises the mean temperature during the growth period, whereas K PMI characterises the air humidity, daily and annual variation in temperature, and minimum temperature. M PMI reflects the moisture in the soil during the growth period of the plants and was included in the climate models to provide information regarding water availability. R PMI and N PMI represent the pH and the amount of nitrogen in the soil, respectively.

To calculate tree cover at the plot grain (TC plot ), we used field information provided by the vegetation surveys. We calculated the cover percentage for all species taller than 2 m (recorded as vegetation layers 1, 2, 3, and 4 in the database) in each plot by transforming the original abundance-dominance indices to approximate percentage values and then summing the cover percentages for all the species in these four vegetation layers, as implemented in Turboveg [START_REF] Hennekens | TURBOVEG, a comprehensive data base management system for vegetation data[END_REF].

Landscape grain-To study the effect of coarse-grain tree cover on plot-level species occupancy, we prepared a similar set of predictor variables, but at the landscape grain (see Appendix 2 for similarities and differences between landscape-and plot-grain variables).

Here, we used freely available data from various sources that are frequently used in SDMs, enabling their use in a broader set of studies where environmental information at finer spatial grains is not available. Using GRASS 6.4.2RC1 for Windows (GRASS Development Team 2011), we created a 30″ resolution grid in the geographic coordinate system and datum WGS84 (i.e., approximately 925 m × 625 m or 0.58-km 2 cells) that encompassed the spatial extent of the final dataset. Then, we computed climate, soil, and tree cover variables at this spatial resolution.

We used monthly mean climatic data for current conditions (~1950-2000) from the WorldClim database [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]) at 30″ resolution to calculate three climatic variables: growing degree days (GDD), absolute minimum temperature (AMT) and water balance (WBAL). The WorldClim data were estimated from long-term monitoring and subsequent spatial interpolation [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF], offering a rough approximation of these metrics. To represent soil conditions, we downloaded two variables from the European Soil Portal -European Commission - [START_REF] Jones | Estimating organic carbon in the soils of Europe for policy support[END_REF][START_REF] Panagos | European Soil Data Centre: Response to European policy support and public data requirements[END_REF]): the organic content in the topsoil (OCT) and soil reaction (pH). OCT was used as a surrogate for nutrient availability. Further details in the calculation of these variables are provided in Appendix 2.

The landscape-grain tree cover data were derived from the "Forest Map 2000" provided by the European Forest Data Centre -European Commission (http://efdac.jrc.ec.europa.eu/). This raster layer contains information on forest presence/absence at 25-m resolution and is in ETRS89 Lambert Azimuthal Equal Area projection. We re-projected the original data using the nearest neighbour method at the original resolution (~25 m). Subsequently, we resampled to 30″ resolution, computing the surface ratio (ranging from 0% for treeless areas up to 100% for completely forested areas) covered by forest in each grid cell. We refer hereafter to this variable as the landscape-grain tree cover (TC land ).

Tree cover and species distribution

To assess the importance of tree cover on species distributions relative to climate and soil, we fit species distribution models (SDMs) using four different combinations of predictors for all selected species in the final dataset: climate, soil and tree cover (CST); climate and soil (CS); climate and tree cover (CT); and climate alone (C). First, we ran models with predictor variables calculated at the plot grain. Then, we investigated whether the effect of tree cover was also observable at a coarser grain by fitting the same combinations of predictors with the independent variables calculated at the landscape grain. In a post-hoc analysis, we included TC land in the CST models at the plot grain to evaluate their combined effect. All the predictor variables included in each model had low correlations (see Table A1), with the exceptions of K PMI vs. M PMI (Spearman correlation, r s = -0.767), GDD vs. AMT (r s = 0.921), GDD vs. WBAL (r s = -0.898), and AMT vs. WBAL (r s = -0.779). Most importantly, TC plot and TC land , which constitute the variables of greatest interest here, showed low correlations with soil and climate variables.

We used Generalised Linear Models (GLMs) with a binomial family and a logistic link function in which all the predictor variables were included as linear and quadratic terms. However, to avoid interpreting the results solely based on a single algorithm, we performed supplementary analyses using Maxent [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF]) and Boosted Regression Trees (BRT; [START_REF] Ridgeway | The state of Boosting[END_REF][START_REF] Elith | A working guide to boosted regression trees[END_REF]) (see Appendix 3 for more details on both methods). For each modelled species, the response variable is a binary variable (0/1) of presence (1) and absence (0) data. Given the high quality of the data, the absence data are reliable. For Maxent, which was designed to compare presences with a background sample, we used all selected plots (including those where the modelled species was present) as background. Our objective with this approach was to obtain comparable results between modelling approaches while minimising the effect of sampling bias [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data[END_REF]). The results from Maxent and BRT were consistent with those obtained with GLMs (see Appendix 3). Previous studies have also used data accumulated over long periods to infer and study biotic interactions [START_REF] Araújo | The importance of biotic interactions for modelling species distributions under climate change[END_REF][START_REF] Hof | How biotic interactions may alter future predictions of species distributions: future threats to the persistence of the arctic fox in Fennoscandia[END_REF][START_REF] Kunstler | Competitive interactions between forest trees are driven by species' trait hierarchy, not phylogenetic or functional similarity: implications for forest community assembly[END_REF][START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF], despite the complication that the strength and direction of the interactions could change over time (e.g., due to changes in environmental conditions). To assess whether interannual variation could affect the results, we performed a sensitivity analysis with Logistic Mixed Models with survey year as a random effect. The results were highly consistent with those based on GLMs, indicating that such temporal effects did not bias the results and the reasoning (see Appendix 4).

To study the performance of models, we randomly selected 75% of the presence-absence data from the final dataset for each of the 960 modelled species to calibrate the models (training dataset) and withheld the remaining 25% to evaluate them (testing dataset). Although the testing datasets were not used to fit the models, they might not be spatially independent from the training datasets. This spatial dependency could bias model evaluation, overestimating the performance of the models; thus, caution should be taken in using the models to make projections, but it is less likely that this bias would affect the comparisons of the models. We used two threshold-independent measures, the nonparametric Area Under the Curve of the Receiver Operating Characteristics plot (AUC) and the parametric point-biserial correlation coefficient (COR), and one threshold-dependent measure, the true skill statistic (TSS). AUC has been widely used but is known to be sensitive to prevalence, whereas COR and TSS have been suggested to be unbiased by prevalence issues [START_REF] Allouche | Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS)[END_REF][START_REF] Liu | Measuring and comparing the accuracy of species distribution models with presenceabsence data[END_REF]. To calculate TSS and to obtain presence-absence predictions from the models, we used the species-model-specific threshold that minimised the difference between sensitivity (proportion of presences correctly predicted) and specificity (proportion of absences correctly predicted) in absolute terms for the training dataset. By using this threshold, our objective was not to assess the ability of the models to predict potential areas of suitability but rather to assess their ability to predict areas the species actually occupies.

Based on this set of 960 modelled species, we first evaluated the significance of incorporating tree cover into the SDMs, in addition to climate and soil. Because of the high number of samples, we used pairwise comparisons of model performance among the four candidate models (C, CS, CT, and CST) with species as the sampling unit (Blach-Overgaard et al. 2010). Boxplots were drawn with boxes showing the interquartile range of the difference in model performance between two candidate models of a given species for each of the six possible combinations of pairwise comparison (C:CS, C:CT, C:CST, CS:CT, CS:CST, and CT:CST) and for each of the three measures of model performances (AUC, COR, and TSS). The notches inside the boxes indicate the 95% confidence interval around the median of the difference between the two candidate models. If the notches of a box do not overlap the zero horizontal line, there is strong evidence (95% confidence) that the median difference in model performance between the two candidate models differs from zero. Additionally, the statistical significance of these differences was also tested using Student's paired t-test, with the sample units (species) assumed to be independent.

Finally, we assessed the relative importance of the variables for the herb and shrub distributions by using hierarchical partitioning. For each species, we ran all possible models with different combinations of the six predictors, modelled as linear terms. Then, the independent contribution of each variable was calculated and partitioned based on the entire set of models.

Tree cover and community composition

If light-mediated plant-plant interactions drive herb and shrub distributions, then the proportion of shade-tolerant species within the community should increase with tree cover; i.e., tree cover should negatively correlate with the mean indicator values of the plot for light (L PMI ). Due to the nested nature of TC plot within TC land , the two variables are correlated, and their effects on herb and shrub distributions could be confounded (r s = 0.387). We disentangled their independent effects by regressing L PMI against TC plot and/or TC land , using partial regressions to assess the effect of each tree cover variable while accounting for the other.

Tree cover and the elevational range limits

To estimate the effect of tree cover, TC plot , on the elevational range limits of lower-stature plants, we compared two extreme scenarios (i.e., low vs. high tree cover) using either real observations or simulations. For the real observations, we first selected two subsets of vegetation plots among all the plots in the dataset: one subset of plots with low TC plot values (≤50 %; n = 1,869) and another subset of plots with high TC plot values (>50 %; n = 4,158). Because the plots were not evenly distributed across the elevational gradient, we performed stratified random sampling to obtain two subsets of plots that were uniformly and equally distributed between 0 and 2200 m a.s.l., selecting 50 plots at each 50-m elevation interval whenever possible (n = 1,238 for each subset; mean TC plot = 15.7% and 81.5% for the low and high subsets, respectively; see Appendix 5). Then, for each species within each subset, we ranked its occurrences along the elevational gradient and selected the 1 st and 99 th percentiles as the lower and upper elevational limits. To reduce uncertainties due to the stratified random sampling, we replicated this computation 10 times and then computed the mean values. Finally, for each species, we calculated the elevational differences between the low and high tree-cover areas at both its lower and upper elevational limits.

For the simulations, we used the most complex model (CST) at the plot grain to predict the elevational limits for each species under the low and high tree-cover scenarios. All predictor variables were set to their original values except for TC plot , which was set to 10% and 90% to simulate low and high tree-cover conditions, respectively. Then, we used the predicted occurrences of a given species (translated to presence-absence predictions by applying the species-model-specific thresholds) and applied the same approach as that applied with the empirical data to calculate the elevational differences between the low and high tree-cover scenarios at both the lower and upper elevational limits. We also compared these results among the different shade-tolerance groups of species using Student's t-test and assuming the sample units (species) to be independent.

All analyses were performed using R 2.14.0 (R Development Core Team 2011) and the dismo, gbm, hier.part, and ncf packages [START_REF] Hijmans | Species distribution modeling[END_REF][START_REF] Oksanen | Community Ecology Package[END_REF][START_REF] Ridgeway | gbm: Generalized Boosted Regression Models Root TL, et al. Fingerprints of global warming on wild animals and plants[END_REF][START_REF] Walsh | hier.part: Hierarchical Partitioning[END_REF].

RESULTS

Fine-and coarse-grain tree cover effects for plot-level species distribution

Species distribution models with predictor variables calculated at the plot grain for the 960 herb and shrub species generally had higher performance than models with predictors at the landscape grain (Table 1). For both of these, the climate-only models (C) performed well, but including soil or tree cover (CS and CT models) improved the overall performance, except for TSS (Table 1 and Fig. 2). Furthermore, including tree cover in addition to climate and soil (CST models) resulted in an additional improvement in performance over the CS models (Table 1 and Fig. 2; paired t (959) = 5.273 and p<0.001 and paired t (959) = 4.359 and p<0.001 for AUC at the plot and landscape grain, respectively). Importantly, using TC land with TC plot as predictors with the climatic (C) and soil (S) variables improved the performance according to COR (paired t (959) = 2.947; p<0.005).

Focusing on the CST models, which showed the best model performances at both the plot and landscape grains, we found different patterns in the relative importance of each predictor variable (Fig. 3). At the plot grain, T PMI was by far the most important predictor variable, followed by TC plot and M PMI . Hence, TC plot was more important than some climatic (K PMI and M PMI ) and soil variables (N PMI and R PMI ). At the landscape grain, AMT was the most important predictor variable, followed by WBAL and GDD. Tree cover at the landscape grain (TC land ) was less important than at the plot grain but still more relevant than OCT (Fig. 3). The response curves of different species along the tree-cover gradient, as estimated by the CST models, differed among the colline, montane and subalpine groups of herb and shrub species (Fig. 4). As expected, the shade-tolerant species showed an increasing probability of occurrence with higher values of tree cover, whereas shadeintolerant species showed a decreasing trend with higher tree cover, with the semi-shadetolerant species having an almost flat response curve.

Fine-and coarse-grain tree cover effects on plot-level community composition

Both TC plot and TC land had negative effects on L PMI (Fig. 5) and together explained a large percentage (58.5%) of the lower-stature plant species composition in terms of their shade tolerance. Partial regressions showed that TC plot and TC land had significant independent effects (Fig. 5), with a stronger unique effect of TC plot (46.4% explained variance) than TC land (9.3%).

Tree cover modifications of elevational range limits

Testing the effect of changes in tree cover at the plot grain on the elevational range limits of lower-stature plants, we found contrasting patterns depending on species shade tolerance (Fig. 6). Based on the empirical data, the elevational range of shade-intolerant species was narrower under high tree cover conditions due to negative differences at the upper limits for the colline and montane species groups and due to both negative and positive differences at the upper and lower limits, respectively, for the subalpine species group. In contrast, the elevational range of shade-tolerant species was greater under high tree cover due to negative differences at the lower limits for the colline species group and due to both negative and positive differences at the lower and upper limits, respectively, for the montane and subalpine species groups. Semi-shade tolerant species showed weaker patterns with no clear trend. The outputs from the simulations were largely consistent with the results based on observations, although with a tendency to predict weaker tree-cover effects than were empirically observed (Fig. 6).

DISCUSSION

This study provides the first broad-scale quantitative assessment of the role of tree cover in the distributions of multiple herb and shrub species and their elevational range limits at two contrasting spatial grains (coarse vs. fine). Notably, we found that tree cover was correlated with distribution and community patterns not only at the vegetation plot level but also at a much coarser grain, near 1-km. These results are thus in line with recent findings regarding the effects of biotic interactions at large spatial scales [START_REF] Gotelli | Macroecological signals of species interactions in the Danish avifauna[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF], exemplifying how local interactions could be observable at coarser grains (but see [START_REF] Cooper | opinion: Habitat data resolution and the detection of species interactions[END_REF]. Additionally, the results suggest that tree cover could lead to important modifications in the elevational range limits of plant species based on their shade tolerances, suggesting an important role of biotic interactions such as competition and facilitation.

Plant-plant interactions are important at coarse spatial grains

Tree cover improved the overall performance of the models for herb and shrub species distributions (Fig. 2). The pairwise comparisons between the candidate models suggest a secondary but important role of tree cover relative to climate and soil variables. The analysis of variable contributions (Fig. 3) for the most complex models (CST) showed that tree cover is more important than soil for models at the plot grain. At this grain, tree cover played a stronger role than at the coarse landscape grain (Fig. 3), likely reflecting the direct effect of the nearest neighbouring trees (Fig. 4) on herb and shrub distributions.

When the predictor variables were measured at the landscape grain, climate remained the main driver, but soil became more important than tree cover (Fig. 3). The changes in variable importance across spatial grains could reflect the differences between the two sets of predictors (i.e., their time frame, origin, and the particular variable represented). Regardless of this grain-dependent change, tree cover still improved model performance at the coarser grain (Fig. 2). Additionally, including TC land with TC plot in the CST models provided small but significant improvements in model performance. Interestingly, curves along the tree cover gradient at the landscape grain were similar to those at the plot grain, although more flattened (Fig. 4), suggesting that shade-tolerant species might be favoured not only by dense tree cover in their direct vicinity but also by being in a forested landscape, and vice versa for shade-intolerant species. These tree cover effects extended up to 1 km for most species and even further in some cases (see Appendix 6 for a supplementary analysis).

The community-level analysis confirmed the structuring role of tree cover in terms of shade tolerance. Indeed, the significant negative effects of tree cover at both plot and landscape grain on the community-based index of light conditions were consistent with the results of the SDMs (Fig. 5). Furthermore, partial regression analyses showed (1) that tree cover in the direct vicinity is more important than tree cover in the surrounding areas and (2) that their effects act independently of each other, with only a small proportion of shared variation (3.2%). These patterns suggest that tree cover could affect plot-level species occurrences not only via direct local processes [START_REF] Gravel | Shade tolerance, canopy gaps and mechanisms of coexistence of forest trees[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF]) but also via landscape-scale effects, e.g., metapopulation dynamics [START_REF] Eriksson | Regional Dynamics of Plants: A Review of Evidence for Remnant, Source-Sink and Metapopulations[END_REF][START_REF] Ehrlén | Long-term assessment of seed limitation in plants: results from an 11-year experiment[END_REF]). In accordance with [START_REF] Bullock | Geographical separation of two Ulex species at three spatial scales: does competition limit species' ranges?[END_REF], who showed that competition between related plant species could have effects at coarse spatial grains, the results contradict the generalised idea that biotic interactions primarily matter at local and fine spatial grains [START_REF] Pearson | Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?[END_REF]. Our findings highlight the importance of considering biotic interactions (e.g., plant-plant, plant-pollinator or plant-grazer interactions) in large-scale studies of plant distributions and diversity patterns [START_REF] Greve | Continental-scale variability in browser diversity is a major driver of diversity patterns in acacias across Africa[END_REF], including predictive models projecting future climate and land use changes [START_REF] Araújo | The importance of biotic interactions for modelling species distributions under climate change[END_REF]Luoto 2007, Thuiller et al. 2013).

Plant-plant interactions relate to shade tolerance

Species response curves along the tree cover gradient (Fig. 4) clearly differed in shade tolerance. Notably, shade-tolerant species were positively correlated with tree cover. This pattern suggests that these species benefit from the altered abiotic conditions that a dense tree cover offers, resulting in a wider elevation range than under low tree cover, in otherwise unsuitable conditions (Figs. 4 & 6). However, shade-intolerant species were negatively correlated with dense tree cover at both elevational range limits and thus had a narrower elevation range than at low tree cover despite the regionally suitable climatic conditions (Figs. 4 & 6).

For shade-intolerant species, tree cover showed an asymmetric pattern regarding elevational limits. For subalpine shade-intolerant species, an increase in tree cover increased the lower elevational range limits, whereas it did not modify the upper limits (Fig. 6). Such a trend would be consistent with the asymmetric-abiotic-stress-limitation hypothesis (AASLH), suggesting that tree cover could act as the limiting factor, through light competition, that sets the lower elevational range limits of shade-intolerant plant species distributions [START_REF] Normand | Importance of abiotic stress as a range-limit determinant for European plants: insights from species responses to climatic gradients[END_REF]. For colline and montane shade-intolerant species, an increase in tree cover decreased the upper elevational range limits (Fig. 6). Such a trend suggests the negative effects of biotic interactions through competition towards the upper elevational range limits of plant species distributions, which is inconsistent with both the stress gradient hypothesis (SGH) and AASLH, since SGH propose positive effects through facilitation [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF] and AASLH purely abiotic constraints at the upper elevational range limit [START_REF] Kikvidze | Climatic drivers of plant-plant interactions and diversity in alpine communities[END_REF][START_REF] Defossez | Do interactions between plant and soil biota change with elevation? A study on Fagus sylvatica[END_REF][START_REF] Pottier | The accuracy of plant assemblage prediction from species distribution models varies along environmental gradients[END_REF].

For shade-tolerant species in the montane and subalpine vegetation belts, an increase in tree cover increased the upper and decreased the lower elevational range limits (Fig. 6), most likely reflecting facilitation processes. Although these results would support the importance of biotic interactions in determining both limits of species ranges [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF], they are a deviation from the SGH, which predicts an increase in competition and a decrease in facilitation as abiotic stress decreases [START_REF] Maestre | Refining the stress-gradient hypothesis for competition and facilitation in plant communities[END_REF]. Such a pattern (deviating from the SGH) towards the lower elevational range limit most likely reflects the fact that most lowland areas are human-altered open habitats (see Appendix 7).

Implications of the use of biotic interactions in species distribution models

The addition of tree cover generally improved the model predictions of herb and shrub species distributions compared with models including only climate and soil variables, especially for the threshold-independent measures (Table 1). These results are thus in line with the increasingly recognised importance of including biotic factors as important drivers in SDMs [START_REF] Gilman | A framework for community interactions under climate change[END_REF][START_REF] Mcmahon | Improving assessment and modelling of climate change impacts on global terrestrial biodiversity[END_REF][START_REF] Kissling | Towards novel approaches to modelling biotic interactions in multispecies assemblages at large spatial extents[END_REF][START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF][START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF][START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF]). However, several factors have prevented their extensive inclusion in SDMs so far. On the one hand, knowledge of interacting species and the nature of the interaction are required but not always available. To manage this problem, exploratory analyses on co-occurrences have been performed [START_REF] Pellissier | Species distribution models reveal apparent competitive and facilitative effects of a dominant species on the distribution of tundra plants[END_REF][START_REF] Kunstler | Effects of competition on tree radial-growth vary in importance but not in intensity along climatic gradients[END_REF], 2012), or interactions have been inferred between congeners or species with similar environmental requirements [START_REF] Meineri | Modeling alpine plant distributions at the landscape scale: Do biotic interactions matter?[END_REF]. However, these approaches should be carefully implemented because they can lead to misinterpretations if some important variables or species were not included in the analyses, if the supposed interactions do not apply, or if the interacting species do not have overlapping distributional ranges [START_REF] Ovaskainen | Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions[END_REF][START_REF] Meineri | Modeling alpine plant distributions at the landscape scale: Do biotic interactions matter?[END_REF]). On the other hand, biotic interactions can involve multiple species, which would reduce the impact of any abiotic variable included in these models [START_REF] Kissling | Towards novel approaches to modelling biotic interactions in multispecies assemblages at large spatial extents[END_REF]. As an alternative, we used a synthetic index that summarises the information regarding biotic interactions through one or few variables, representing the effect itself (as a modulator in the sense of [START_REF] Linder | Biotic modifiers, environmental modulation and species distribution models[END_REF]) instead of the numerous species that generate it [START_REF] Boulangeat | Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances[END_REF].

Our results also have implications for the use of SDMs to predict plant responses under climate change scenarios and biodiversity conservation. Such responses are increasingly recognised as not being as straightforward [START_REF] Dullinger | Extinction debt of high-mountain plants under twenty-first-century climate change[END_REF] as initially thought [START_REF] Lenoir | Climate-related range shifts -a global multidimensional synthesis and new research directions[END_REF], for example simple poleward and upward range shifts [START_REF] Parmesan | Ecological and Evolutionary Responses to Recent Climate Change[END_REF][START_REF] Wilson | An elevational shift in butterfly species richness and composition accompanying recent climate change[END_REF][START_REF] Lenoir | A Significant Upward Shift in Plant Species Optimum Elevation During the 20th Century[END_REF]. Indeed, the results of the present study suggest that biotic factors can be important additional drivers that influence species-specific range shifts. The long life cycles of trees, however, implies long time periods to generate such vegetation dynamics (Lenoir et al. 2008, Svenning and[START_REF] Svenning | Disequilibrium vegetation dynamics under future climate change[END_REF], and therefore, dynamic models accounting for transient population dynamics should ideally be used to predict the future of plant distribution and diversity [START_REF] Araújo | The importance of biotic interactions for modelling species distributions under climate change[END_REF][START_REF] Meier | Co-occurrence patterns of trees along macro-climatic gradients and their potential influence on the present and future distribution of Fagus sylvatica L[END_REF], 2012[START_REF] Dullinger | Extinction debt of high-mountain plants under twenty-first-century climate change[END_REF][START_REF] Hof | How biotic interactions may alter future predictions of species distributions: future threats to the persistence of the arctic fox in Fennoscandia[END_REF][START_REF] Thuiller | A road map for integrating eco-evolutionary processes into biodiversity models[END_REF]. Relative variable importance (I) for 960 herb and shrub species as independent variance explained by each variable (as linear predictors) in a hierarchical partitioning with all the predictor variables as linear terms. Boxplots show the per cent contributions across all species. Diamonds represent mean values. Variables are ordered by decreasing relative importance. Variables are coloured by scale, with the plot-grain variables in light grey: air temperature (T PMI ), continentality (K PMI ), soil moisture (M PMI ), nutrients (N PMI ), soil reaction (R PMI ) and tree cover at the plot-grain (TC plot ); and the landscape-grain variables in dark grey: growing degree days (GDD), absolute minimum temperature (AMT), water balance (WBAL), organic content of the top soil (OCT), soil reaction (pH) and tree cover at the landscape grain (TC land ). Mean response curves of all 960 species to the tree cover gradient using CST models at both plot and landscape grain. The black line represents the mean response curve, and the light grey area shows the 99% confidence intervals. To compute these response curves, all other predictor variables were set to their mean values. The total number of herb and shrub species in each shade tolerance group (L1-L2 for shade; L3 for semi-shade; and L4-L5 for light) is at the bottom of the graph. Pure effect of plot-grain tree cover (TC plot ) on the plot mean indicator value for light (L PMI ) after accounting for the effect of landscape-grain tree cover (TC land ) on L PMI and the pure effect of TC land on L PMI after accounting for the effect of TC plot on L PMI . All regression analyses were performed with standardised values. Hexagons represent point densities in a grey scale. Differences in the elevational range limits of different species between low and high tree cover conditions for 717 herb and shrub species: 193 from the colline belt, 302 from the montane belt and 222 from the subalpine belt. Squares indicate observed differences based on empirical data between high (> 50 %) and low tree cover (≤ 50 %), subtracting the latter from the former. Triangles show the same differences based on simulated data obtained from the predictions of the CST models at the plot grain. The simulation was run with the predictor variables for each plot at their calibration values, except for tree cover, which was set at 90% and 10% to simulate high and low tree cover, respectively. Bars indicate the 99% confidence intervals for the mean values calculated with one-sample Student's t-tests.

Horizontal lines indicate no difference in the elevational distribution of the species between the high and low tree cover conditions. Median, minimum and maximum values for the Area Under the Curve of the Receiver Operating Characteristic plot (AUC), correlation coefficient (COR), and true skill statistic (TSS) for 960 herb and shrub species modelled with different sets of predictor variables (C: climate; S: soil; and T: tree cover). Evaluation was performed on the test datasets with 25% of the total plots.

Figure 1 .

 1 Figure 1.Spatial distribution of the 6,935 studied vegetation plots across the Alps. The background colour indicates decreasing tree cover from green (100%) to light grey (0%) at the landscape grain (30").

Figure 2 .

 2 Figure 2. Pairwise comparisons of model performances for 960 herb and shrub species modelled with different sets of predictor variables (C: climate; S: soil; and T: tree cover) among four candidate models (C, CS, CT, and CST) and based on three different measures: Area Under the Curve of the Receiver Operating Characteristic plot (AUC), correlation coefficient (COR), and true skill statistic (TSS). Boxplots show the difference in model performance between two candidate models of a given species for each of the six possible combinations of pairwise comparison (C:CS, C:CT, C:CST, CS:CT, CS:CST, and CT:CST) and for each of the three measures of model performances (AUC, COR, and TSS). Notches inside the boxes indicate the 95% confidence interval around the median of the difference between the two candidate models. The horizontal line represents no difference between the candidate models, and the two candidate models differ in performance when there is no overlap between the notches of a given box and this line. Note the change in scale for AUC.

  Figure 3.

  Figure 4.

  Figure 5.

  Figure 6.

Table 1

 1

Nieto-Lugilde et al. Page 22 Ecography (Cop.). Author manuscript; available in PMC 2015 August 17. Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts Europe PMC Funders Author Manuscripts Nieto-Lugilde et al. Page 23

Ecography (Cop.).Author manuscript; available in PMC 2015 August 17.

ACKNOWLEDGEMENTS

We acknowledge two anonymous reviewers and the subject editor (Robert P. Anderson) for their insightful comments and suggestions on earlier versions of the manuscript. Universidad de Granada (Spain) funded D.N.L. (Perfeccionamiento de doctores Fellowship). J.-C.S. was supported by the European Research Council (ERC-2012-StG-310886-HISTFUNC). W.T. received support from the European Research Council under the European Community's Seven Framework Programme FP7/2007-2013 Grant Agreement no. 281422 (TEEMBIO). We additionally consider this study a contribution to the Center for Informatics Research on Complexity in Ecology (CIRCE), funded by Aarhus University and the Aarhus University Research Foundation under the AU IDEAS program.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.