Global warming will affect the maximum potential abundance of boreal plant species

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Signe Normand Editor-in-Chief: Miguel Araújo Accepted 1 February 2020 43: 1–11, 2020 doi: 10.1111/ecog.04720 43 1–11


Introduction
The biogeographical distribution of plant species is changing as a response to current global warming, and this trend is expected to continue over the next decades (Lenoir and Svenning 2015). This is especially likely in boreal and arctic environments, where changes in climate are predicted to be larger than at more equatorial latitudes (IPCC 2013). At northern latitudes, the effects of warming on vegetation have been already described in terms of upslope shifts of species ranges in mountainous regions (Felde et al. 2012), thermophilisation of communities at high altitudes (Gottfried et al. 2012) or increase of plant biomass and shifts in relative species abundances at high latitudes (Sistla et al. 2013, Vuorinen et al. 2017. For terrestrial plants, latitudinal migrations are far less reported than altitudinal shifts (Lenoir and Svenning 2015), maybe because the latter are more easily studied (Jump et al. 2009). Moreover, the response of different species to warming is highly variable, in both magnitude effect and direction (Lenoir et al. 2010, Chen et al. 2011, even within given taxonomic groups and regions (Parmesan and Hanley 2015). For instance, in the well-studied tundra biome, warming lead to an increase in the cover of shrubs, forbs and graminoids, while a decrease occurred in the cover of mosses and lichens (Walker et al. 2006, Elmendorf et al. 2012, Lang et al. 2012. Temperature sensitivity of vegetation has been widely analysed by using observed or predicted average monthly temperatures, where daily temperature data have not been available (Prevéy et al 2017), but growing degree-days may better forecast timing of the biological events (Pauling et al. 2014). In order to take advantage of best practices for conservation and ecosystem management, the effects of warming on vegetation may need to be forecasted. This is a challenging objective given the intrinsic complexity of species environmental responses and between-species positive and negative interactions, and only some studies have attempted to predict latitudinal shifts of northern vegetation under future global warming scenarios to date (but see Gignac et al. 1998, Bakkenes et al. 2002, Beauregard and de Blois 2016. Forecasting species distributions under global warming scenarios requires relying on information emerging from the current relationship between species and temperature in order to predict how projected temperature changes will affect species in the future. However, an important concern related to this approach is that temperature may not be the only factor directly limiting plant species abundance. The limiting effect of other variables, such as soil characteristics (Salemaa et al. 2008, Bertrand et al. 2012, Beauregard and de Blois 2014, light regime (Nieto-Lugilde et al. 2015), and water availability (Greenberg et al. 2015) might be more relevant in some circumstances. On the other hand, ecological facilitation may ameliorate physical stress and allow species to survive at locations that would a priori be beyond their environmental tolerances (Michalet and Pugnaire 2016). For instance, dense forest canopies buffer climatic extremes and create microrefugia for understorey plant species (Davis et al. 2019). Thus, the response of species to temperature may vary among regions with different combinations of environmental variables (Peñuelas et al. 2004), land management (García-Valdés et al. 2015, Tonteri et al. 2016) and biotic interactions (Shi et al. 2015). This includes locations where temperature may be the factor limiting abundance, as well as locations where other factors may oust temperature from limiting species abundance (Gibson-Reinemer and Rahel 2015, Beauregard and de Blois 2016). Liebig's law of the minimum predicts that only one of these environmental factors will be the active limiting constraint of species abundance at any given point in time and space (Hiddink and Kaiser 2005). However, the reality is that it is difficult to identify which factor is limiting species abundance at each point in space at present time, and it is almost impossible to predict the shift of all potential limiting variables under future global change scenarios. Therefore, if the intention is to forecast how temperature changes will affect species distribution, we need to focus on locations where temperature will actually constrain species performance.
Quantile regression models are beneficial for estimating the effect of a limiting factor when it is known that other unmeasured factors could also be the active limiting constraint at some locations (Cade et al. 1999, Cade and Noon 2003, Koenker 2005, Austin 2007, Vaz et al. 2008. For instance, it may be possible to model the limiting effect of temperature; i.e. the maximum potential abundance that a species could attain at each given temperature (Greenberg et al. 2015, Carrascal et al. 2016. Notably, this maximum potential abundance may reflect the whole potentiality of the species, either if it is determined by the fundamental niche or by facilitation allowing persistence beyond it (Bruno et al. 2003). That may be an approach to the maximum ecological response of a species to temperature, in principle not affected by the limiting effect of other environmental factors. We think this makes the method especially suited to forecast the fate of species under future scenarios of global warming for two reasons. First, as much as the model captures the isolated effect of temperature, it avoids confounding effects from other limiting environmental factors. Second, forecasting the maximum potential abundance is an implicit recognition of uncertainties on future realized distributions related to species' dispersal capacities and to reliability of land-use change scenarios. Although the technique is gaining popularity with time (Eastwood et al. 2003, Schröder et al. 2005, Vaz et al. 2008, Chessman 2012, Cozzoli et al. 2014), quantile regression is not yet widely used by ecologists, and has virtually never be used for forecasting the effects of global warming (but see Jarema et al. 2009).
The objective of this work is to forecast the effect of global warming on the maximum potential abundance of boreal forest understory plants for the upcoming decades. We rely on the hypothesis that the effect of temperature -when significant -is reflected in the maximum potential abundance that a species can attain at a certain location. Forecasts were restricted to a selected group of species with a significant response to temperature that was maintained stable when including other environmental factors in the model such as precipitation, texture and fertility characteristics of soil, stand basal area and proportion of deciduous tree species.

Study area and data
We used soil and vegetation data from locations across Finland surveyed in 1985-1986 on a systematic sampling network established by the Finnish national forest inventory (Mäkipää and Heikkinen 2003). Vegetation, tree and basic soil measurements were available from all 3009 permanent sample plots (Mäkipää and Heikkinen 2003), while more complex soil analyses were available from 453 plots (Tamminen andStarr 1990, Tamminen 2000). The 300 m 2 circular monitoring plots formed a regular network of clusters ( Fig. 1; Reinikainen et al. 2000, Mäkipää andHeikkinen 2003). In southern Finland, each cluster consisted of four plots at 400 m intervals, with one cluster per 16 × 16 km square. In the north, there were three plots per cluster at 600 m intervals, with one cluster per 24 × 32 km rectangle. In order to avoid pseudo-replication issues derived from the aggregated nature of plots, we randomly selected one plot per cluster for analysis. Moreover, we focused on forest plots on mineral soils (excluding peatlands, i.e. peat layer > 30 cm or peatland vegetation covers > 70% of forest floor, and non-forested plots, i.e. bare lands with tree growth < 1 m 3 yr −1 ), including the full range of soil fertility, silvicultural treatments, tree species and age classes. Considering these restrictions, we included 868 sample plots in the analyses.
The percentage cover of understory plants was estimated visually on 2 m 2 sampling quadrats located systematically within the plots (Fig. 1). Species plot-level abundance was estimated as the mean cover over 3-6 quadrats. If the plot included different stand compartments or land types, we used only the quadrats equivalent to the central point. We studied 25 major understory species with a minimum prevalence of 150 observations over the 868 plots. These species belonged to the plant groups of dwarf shrubs (6), graminoids (4), herbs (5), mosses (9) and lichens (1), and had a mean prevalence of 388 plots in the 868 (range 150-797; Supplementary material Appendix 1). A minimum prevalence of 150 was established to minimize error type I in quantile regression models (Cade et al. 2005).
Past temperature and precipitation data was provided by the Finnish Meteorological Inst. (FMI) at a daily resolution for a grid scale of 1 × 1 km 2 (Venäläinen et al. 2005). These values are an approach to those experienced by plant species, considering fine-scale thermal heterogeneity and differences between measured air temperatures and ground-level temperatures. We derived effective temperature sums and cumulative precipitation for each year and sample plot, and averaged them for a period of 25 yr prior to forest inventory . We used annual precipitation, since beside summer precipitation water volume of snow have importance for soil moisture (Ilvesniemi et al. 2010). The effective temperature sum (growing degree days GDD °C) was defined as the sum of the positive differences between diurnal mean temperatures and 5°C yr −1 . The GDD is a measure of heat accumulation widely used by phenological and agricultural studies to forecast state of vegetation (Pauling et al. 2014). For future temperatures, we used predictions for the period 2041-2070 derived by the FMI considering the IPCC emission scenario A1B (Jylhä et al. 2004), which is very similar to the more recent RCP 6.0 scenario (Moss et al. 2010, IPCC 2013. Other predictors included in the analyses were soil texture, soil fertility, the total basal area of the stand, and the proportion of deciduous tree species. Soil texture indicated soil water holding capacity and moisture availability. Soil texture was assessed for each vegetation quadrat using a ten-point scale ranging from bare rock to grain size < 0.002 mm, and averaged per plot as a gross measure of the physical structure of the soil. In the subset of 453 plots with more precise soil data, we confirmed that soil texture was significantly related to both clay content in soil (t-statistic = 9.610, p-value << 0.001) and soil median grain size (t-statistic = −4.981, p-value << 0.001). Soil fertility was a factor of four levels that indicated soil nutrient availability. It was defined in the field based on vegetation composition and stand characteristics, using an eightpoint scale where 1 corresponds to the highest fertility  (Hotanen et al. 2008). We merged classes 1 + 2 and 5-8 because of a low number of observations from extreme classes. In the subset of 453 plots we showed that the reclassified variable was significantly related to the carbon/ nitrogen ratio of the soil organic layer (t-statistic = 14.02, p-value << 0.001). Stand basal area was used as a proxy for light availability for the understory vegetation, and was measured from three relascope plots that included area of the primary sampling plot. The proportion of deciduous trees calculated from the basal area of all trees was used as a surrogate of biotic interactions, via canopy shading, throughfall and litter quality. All tree species were identified as a part of stand inventory, but we used two classes according to dominant species and their functional influence on understory vegetation: 1) conifer species, mainly including Scots pine Pinus sylvestris and Norway spruce Picea abies, and 2) deciduous species including silver birch Betula pendula, downy birch Betula pubescens, common aspen Populus tremula, grey alder Alnus incana, common alder Alnus glutinosa and other minority deciduous species. Correlation between predictors was moderated or low: temperature sum -precipitation (0.64) and soil fertility -basal area of tree stand (−0.38) showed the largest values, while other pairs of predictors showed absolute values below 0.28.

Forecasting maximum potential abundance (MPA) of species
The relationship between species abundance and environmental factors was analysed with quantile regression models, using the rq() function in the R package quantreg (Koenker 2019, R Core Team). At an exploratory stage, we found that most plant species displayed zero-inflated abundance data, and for most species, the abundancetemperature plot showed polygonal point clouds and heteroscedastic error distributions rather than curvilinear relationships (see example in Fig. 2). Beyond the abovementioned conceptual advantages of quantile regressions, the technique also deals well with this type of data (Cade et al. 2005, Vaz et al. 2008. Furthermore, we decided to use quadratic terms of our predictors to allow for non-linear responses. After through exploration, we disregarded the more flexible nonparametric additive version of quantile regression models because they presented large and highly heterogeneous confidence intervals. All models were defined at the upper 95% quantile (tau = 0.95), which was chosen as the best compromise solution between defining the upper response and minimizing error type I (Cade et al. 2005). Thus, the 95% quantile is an approach to the upper response at 100%, whose calculation is not computable by quantile explaining the effect of the effective temperature sum on the maximum potential abundance and on the average abundance of Dicranum polysetum, respectively. Both models include the linear and quadratic terms of effective temperature sum. Quantile regression model was developed using qr() function of quantreg R package, and generalized linear model was developed using the glm() function of stats R package (R Core Team). Blue shaded areas indicate confidence intervals at 95% for each model. The 868 sample sites included in the model are represented with black transparent points, so that darker regions indicate point superposition. The small graph is a zoom into the smaller values of abundance (0-2) in the larger graph. regression models. In order to consider the proportional nature of the response variable and to restrict predictions between 0 and 100, we followed a procedure based on Orsini and Bottai (2011) that relies on the equivariance of quantiles to monotonic transformations. We added a small number to all observed abundance values (specifically half the smallest observed non-0 abundance) and logit-transform them prior to applying rq() function. Afterwards, we back-transformed the estimated values and subtracted the 'small number'.
We developed two sets of explanatory models. Temperatureonly models related species abundance in 1985 to average effective temperature sum in the period 1961-1985 in models that did not include any other covariate. All-predictors models related species abundance in 1985 with effective temperature sum, cumulative precipitation, soil texture, soil fertility, stand basal area and proportion of deciduous tree species. We analysed the response in the two models in order to select the species for which we will attempt forecasting. Specifically, we selected those species that showed a response to temperature that was significant in both temperature-only models and all-predictors models and consistent between them (see Supplementary material Appendix 2 for more details). The consistency between the two models protects our results from the potential effect of predictor collinearity.
The predictive performance of models was assessed using the measure proposed by Li and Peng (2017), based on the check loss functions, coupled with two-fold cross validation. Training and validating datasets for cross-validation were built from plots within clusters: for each cluster we randomly assigned one plot to the training dataset and one to the validating dataset. This way we avoided spatial autocorrelation in evaluation procedures (Roberts et al. 2017). Both temperature-only and all-predictors models were fitted to the training dataset, they were used to predict the 95% quantiles associated to the predictor values of the validating set, and check losses were computed in the logit-scale between the predictions and the observed abundance values in the validating set. We repeated this evaluation procedure for 10 000 random partitions of the dataset into training and validating subsets, and calculated the mean check losses L 1,sp (temperature-only model) and L 2,sp (all-predictors model) over the plots in the validating subsets and over the 10 000 repetitions, for each species (sp). Lower L-values indicate a better predictive capacity.
We used temperature-only models to forecast the MPA of species in Finland under the IPCC A1B scenario of global warming for the 2041-2070 period. The use of temperature-only models avoided reliance on environmental variables for which there are not reliable forecasts for future (e.g. basal area of trees or tree composition, which depend on forest management). As the response of species to novel temperatures is unknown, no extrapolation was attempted to areas whose effective temperature sum is beyond the range in training models (Dormann 2007).
We calculated latitudinal shifts of the abundance distribution by identifying comparable points of the predictions for 1985 and 2040-2070. Specifically we calculated the minimum and maximum latitude of the values corresponding to 75%, 50% and 25% the maximum absolute value of MPA predicted for the species. Then, for each species, we focus on the range limits that fell within Finland for the two dates. We excluded those predicted to fall at the maximum or minimum absolute latitudes of the country, as we could not assess whether they would reach further locations if they would have the possibility. For species with predominantly southern distributions in 1985, the latitudinal shift was calculated for their northern limits (i.e. for the maximum latitude at which these 75%, 50%, 25% MPA abundance values were found). For northern species, the latitudinal shift was calculated for their southern limits (i.e. the minimum latitude at which these abundance values are found).

Significance, consistency and model predictive performance
The predicted response of species to effective temperature sum was generally consistent between temperature-only models and all-predictors models: trends were alike in 15 of the 16 species that showed significant responses to effective temperature sum in both models (Supplementary material Appendix 1, 3). Temperature-only models and all-predictors models showed similar predictive performances, especially for those species with significant and consistent responses to temperature sum (Supplementary material Appendix 1, 3, 4).

Species' response to effective temperature sum and future projections
MPA limits imposed by effective temperature sum were predicted with temperature-only models for each point in Finland in 1985 and 2040-2070, and two main types of patterns were identified (Fig. 3, Supplementary material Appendix 5). In a first group of species, MPA values were significantly higher at higher temperatures. These species reached their highest MPA in the hemiboreal and southern boreal subzones of Finland in 1985 (see Calamagrostis arundinacea, Dicranum polysetum, Maianthemum bifolium and Trientalis europaea in Supplementary material Appendix 5). In the future, their MPA will increase northwards, potentially reaching northern regions where they were rare in 1980s (see 2040-2070 predictions in Fig. 3, Supplementary material Appendix 5). Specifically, their northern range edge is predicted to move northwards over 500 km from 1985 to 2040-2070 (Table 1, Supplementary material Appendix 6). In the future climate, these species could also face changes in their abundance in southern Finland, but we cannot assess that in our data.  (Table 1).
Other species show lower abundances at higher temperature sums, indicating that they are limited by the highest temperatures (see Carex globularis, Dicranum scoparium, Hylocomium splendens, Pleurozium schreberi, Polytrichum juniperum, Solidago virgaurea, Vaccinum uliginosum, Vaccinum vitis-idaea in Supplementary material Appendix 5). According to 1985 data, the MPA of these species peaked in the middle and northern boreal sub-zones. Thus, these species are predicted to reach highest abundances at higher latitudes and to decrease under future warming scenarios. The southern range limits of these species are expected to move northwards an average of 376 km from 1985 to 2040-2070 (Table 1, Supplementary material Appendix 6).
The response of species to effective temperature sum was not related to the functional group to which the species belong: all functional groups showed a variety of patterns and thermal optimums (Fig. 4).

Species response to other environmental factors
In addition to temperature sum, understory plant species showed significant responses in MPA to tree stand composition, soil fertility and annual precipitation (Fig. 5, Supplementary material Appendix 7). Calamagrostis arundinacea, Dicranum majus, Epilobium angustifolium, Luzula pilosa, M. bifolium and T. europaea had significant positive response to high levels of soil fertility (class 1), while Calluna vulgaris, Cladina arbuscula, C. rangiferina, Empetrum nigrum and P. schreberi had highest MPA in lowest fertility class (class 4). Abundance of Vaccinium myrtillus increased and V. vitis-ideae decreased with increasing precipitation. Trientalis europaea had significant positive and P. schreberi significant negative response to increase in the proportion of deciduous trees. Higher canopy cover measured as stand basal area had positive influence on abundance of V. myrtillus, M. bifolium and a bryophyte D. majus. On the other hand, E. angustifolium, L. pilosa and Polytrichum commune showed significant negative response to increase in stand basal area.

The fate of boreal understory plants in a warmer future
We predict that locations capable of supporting the maximum abundances of fifteen boreal plant species will move an average of 6.6 km yr −1 . This means a total move of ca 460 km (range: 49-607 km) northwards from 1985 to 2041-2070 (see 75% MPA in Table 1, Supplementary material Appendix 6). We predict similar shifts for other range limits (see 25% and 50% MPA in Table 1), though the retraction of species that are currently dominating the north may be faster (mean shift of 553 km, at a pace of 7.9 km yr −1 ) than the expansion of southern species (419 km, at 6 km yr −1 ). Our results are in accordance with those of Loarie et al. (2009), as they predicted that boreal species would need to move between 0.01 and 10 km yr −1 to remain in the same temperature niche assuming there are minimal changes in other environmental conditions. Notably, our estimates fall in the upper part of the range, indicating that temperature sensitive species in Finland may experience some of the highest pressures for range shift across the boreal region. The monitoring of the future changes in the species range and abundance would provide indicator of the effects of climate change in the forests.
Our results emphasize inter-specific differences in responses to temperature and the impact of global warming on the understory vegetation of boreal forests. Species in all functional groups from dwarf shrubs, herbs and grasses to bryophytes and lichens showed significant responses to temperature, though the direction of effects varied. Calamagrostis arundinacea, M. bifolium, T. europaea and D. polysetum, showed a positive significant and consistent response to effective temperature sum, and are forecasted to potentially expand northwards. The abundance of a bryophyte species D. polysetum has increased in southern and middle boreal zone since 1950s (Mäkipää and Heikkinen 2003), and the present results suggest that the change is partly driven by Table 1. Mean latitudinal shift in northern and southern range limits of different thresholds in the maximum potential abundance (25%, 50% and 75% MPA) from 1985 to 2040-2070, and the mean of these three thresholds (Mean shift), across species. Shift are calculated when possible given data limitations (see Methods), and mean shift is the pondered mean of these three values. Thus, n indicate the number of species across which mean values are calculated. All latitudinal shifts are given in km. 8 climate warming. Based on our nation-wide vegetation data from Finland we are not able to predict their responses beyond temperature range within country. Therefore, we cannot predict their future abundance in southern Finland where novel climates are forecasted for 2040-2070. Nonetheless, since C. arundinacea, M. bifolia and T. europaeus are currently very common in the Baltic countries south from Finland (Hultén 1971), they will likely maintain or increase their abundance in southern Finland under future climates. A second group of species currently have their highest potential abundances in middle or northern boreal forests (Supplementary material Appendix 5). We predict that these species will migrate northwards, reducing their maximum potential abundance in their current areas of distribution. These predictions are supported by the decline of C. globularis observed in southern Finland from 1950 to 1995 (Reinikainen et al. 2000) and the decline of Cladina rangiferina and H. splendens in southern and middle boreal zone since 1950s' (Mäkipää and Heikkinen 2003). These changes in abundant generalist species suggests that climate warming will likely influence plant community composition and modulate the relative abundance of competing species.

Range limits
Our predictions are in terms of the maximum potential abundances of species, so that results interpretation pertains to potential changes in abundance. Realization of the predicted shift in distribution will depend on population dynamics including source-sink processes (Fordham et al. 2013), on the ability of species to cope with warming in situ by thermal tolerance (Bertrand et al. 2016) or phenological plasticity (Amano et al. 2014), on the dispersal and colonization capacity of species (Moser et al. 2011), on the impediments they will encounter to follow temperature shifts (Dullinger et al. 2015), and on the limiting effect of factors beyond temperature in future destinations. Wind-dispersed grasses and bryophytes (e.g. C. arundinacea and D. polysetum) are likely to respond more effectively to released climate limitations than short herbs with low seed production and very local dispersion (M. bifolium and T. europaeus ;Hiirisalmi 1969, Raatikainen 1990, Honnay et al. 2006. Especially in wind-dispersed species, source populations may produce abundantly sexual propagules that maintain the existence of in sink populations in less suitable habitats. On the other hand, herb, grass and fern species with high nitrogen requirements (Heikkinen and Mäkipää 2010) may not perform well in nutrient-poor destinations even if thermal environment is optimal. According to our analysis, C. arundinacea, M. bifolia and T. europaea had positive response to increasing soil fertility (Supplementary material Appendix 7), so in the future climate they are likely to colonize fertile sites in northern locations but not to thrive on poor sites. Thus, anthropogenic nitrogen deposition will likely favour the expansion of these species (Heikkinen and Mäkipää 2010). Furthermore, understory species associated with deciduous trees such as T. europaea may not easily colonise large areas dominated by conifers in the north; although we must consider that tree distribution may also change with warming. The abundance of some plant species depend strongly on light conditions and increases in clear-cut areas due to opening for light (Tonteri et al. 2016), while shade-tolerant species such as M. bifolium, V. myrtillus and D. majus, respond positively to increasing stand density (see basal area in Supplementary material Appendix 7). The interaction of these factors will determine the capacity of species to track global warming and will potentially provoke extinction debts and colonization Figure 5. Quantile regression models at 95% of the partial effect of different environmental factors on the abundance of Dicranum polysetum as an example species. Models were developed using data of 868 nation-wide plots sampled in Finland in 1985 and climate data from 1961 to 1985. Abundance of species is measured as percentage cover. The partial effect of each variable is controlled by the average effect of all other variables in all-predictors models. All variables are included with linear and quadratic terms, except for the factor fertility level, which is only included in its simple term. Confidence intervals at 95% are shown with a shaded area. Temp. sum: effective temperature sum (°C), Cumul. precip.: cumulative precipitation (mm yr −1 ), Basal area: stand basal area (ha −1 ), Soil texture: soil texture (index 1-10, from bare rock to grain size < 0.002 mm), Soil fertil.: soil fertility (index 1-4, from high to low fertility), Decid. trees: proportion of deciduous tree species (% of tree basal area). 9 credit. As the result from the interaction of all factors that affect plant species are challenging to predict (Franklin et al. 2016), our predictions on maximum potential abundances aim to give insight to species' potential range within which species realized distributions will be encountered in future climate.
We identified potential shifts of boreal understory vegetation that may be of concern for tree stand regeneration. Among species particularly problematic for forestry, C. arundinacea may expand its northern range limits and increase its abundances in northern locations. Calamagrostis arundinacea is a tall grass with high biomass production that can regrow from buds and rhizomes and produce many seeds. These traits define good dispersal capacity and may facilitate the performance of the predicted range shifts (Aubin et al. 2016). Thus, if the predicted range shift of C. arundinacea occurs, this may affect survival of the tree seedlings and stand regeneration.

Limitations of our results
We argue here that quantile regressions are capable to predict the limitation imposed by a specific factor on species abundance. Nonetheless, a prerequisite for a complete model training is the existence of locations throughout the environmental gradient where the species is actually limited by the factor of interest (Cade et al. 1999). For instance, C. rangiferina could probably reach larger abundances at northern latitudes than predicted by our models (Supplementary material Appendix 5). This and other lichen species are intensively grazed by semidomesticated reindeers in the northern half of Finland (Kumpula et al. 2000, Den Herder et al. 2003, which reduce their ground coverage well below the limits imposed by temperature. In this case, the lack of northern non-grazed plots in the survey data may prevent accurate modelling of temperature limitations to these lichen species in cold regions. Cade et al. (1999) suggested that, even if incomplete, these responses may be more consistent with the ecological theory of limiting factors than estimates by commonly used models. Nonetheless, these sources of uncertainty may require attention when these models are applied.
The focus on maximum potential abundance: a novel conceptual approach to global change biogeography Species distributions are usually defined in terms of range limits, range size, probability of occurrence or local abundance (see reviews by Ehrlén andMorris 2015, Yalcin andLeroux 2017). Here we propose the concept of 'maximum potential abundance' (MPA) as an approach that includes by definition the uncertainty derived from the effect of environmental factors not considered in the model, the dispersal limitations of species or any other contingency determining the realized abundance of a species. This approach is especially well suited for capturing the fundamental niche with respect to a specific environmental factor, in contrast to most common approaches that reflect species' realized niche. What is more, as the approach is based on biogeographical data and not on physiological experiments, it may evidence the fundamental niche of species even in the cases in which it was expanded by ecological facilitation (Bruno et al. 2003). That is, the maximum potential abundance will not estimate species realized abundance at a specific point in time and space, but the limitation that an environmental factor is exerting at that point (as defined by Liebig's law of the minimum). The realized abundance at this point will be predicted to fall between zero and the maximum potential abundance, and this range will explicitly reflect the uncertainty related to the influence of other environmental factors not accounted for Noon 2003, Greenberg et al. 2015). We believe this perspective is especially interesting to discuss the impact of specific drivers of global change under future scenarios in systems with available abundance data sampled across largescales. Howard et al. (2014) already highlighted the advantage of moving towards good quality abundance data for improving species distribution models, and Gibbons et al. (2007) stated that this type of data may not be more costly to collect than presence data in terms of time and number of observers.
Here we illustrate our proposal to use quantile regression models and maximum potential abundances in climate change biogeography in its simplest form. However, the concept may be combined with other procedures. For instance, spatialized demographic models could work on the maximum potential abundance imposed by climate at each location, in order to estimate realized abundances below that limitation. Furthermore, the idea may serve to estimate the ecological limitation imposed by environmental factors related to the habitat of a species, and to predict the impact of other global change drivers such as land use change or nitrogen deposition. In summary, we propose that quantile regression models are useful to estimate the limiting effect of any environmental factor on species numbers and have a large potential application in biogeography and global change research.

Data availability statement
A subset of the database used in the study, including all environmental variables and abundance data of two species for all 869 plots, is available from the Dryad Digital Repository: < https://doi.org/10.5061/dryad.vq83bk3pm > (Villén-Peréz et al. 2020). The whole database is archived in the vegetation database of the Natural Resources Institute Finland and available from authors on request.
Funding -This research was partially funded by the Academy of Finland (grants 140766 and 312912) and by Comunidad de Madrid in agreement with Univ. de Alcalá under the program for support of young scientists' research (CM/JIN/2019-003). During the development of this work, SV-P was supported by a contract of visiting scientist (associated with grants 140766 and 312912) and by a TALENTO Fellowship (2017-T2/AMB-6035, Comunidad de Madrid.