Temporal variability is key to modelling the climatic niche

Niche‐based species distribution models (SDMs) have become a ubiquitous tool in ecology and biogeography. These models relate species occurrences with the environmental conditions found at these sites. Climatic variables are the most commonly used environmental data and are usually included in SDMs as averages of a reference period (30–50 years). In this study, we analyse the impact of including inter‐annual climatic variability on the estimation of species niches and predicted distributions when assessing plant demographic response to extreme climatic episodes.

Most of these studies generally use correlative SDMs that statistically correlate species' presences with the environmental conditions of sites where species occur (Franklin, 2010;Guisan et al., 2017;Peterson et al., 2011). These models are deeply rooted in Hutchinson's environmental niche concept (Guisan & Zimmermann, 2000), which implies that species' observed ranges are the geographical translation of species' environmental requirements in the N-dimensional space (i.e. the fundamental niche, sensu Hutchinson, 1957), discarding the areas where competitors impede the species presence (Soberón & Nakamura, 2009) and areas out of reach by dispersal (Barve et al., 2011;Pulliam, 2000). Several limitations of this method have been recognized, linked to both the underpinning theory and the commonly used methodologies. Some of the most criticized limitations include the equilibrium assumption (Guisan & Thuiller, 2005), which neglect any lag between environmental changes and species distributions (Blonder et al., 2015;Svenning & Skov, 2004); the niche conservatism assumption (Pearman et al., 2008), since these models implicitly consider that the niche remains constant when projected in time and space ; the lack of explicit inclusion of biotic interactions (Wisz et al., 2013); or the frequent use of exclusively climatic predictors (Franklin, 2010;Mod et al., 2016;Thuiller, 2013), partly encouraged by their widespread availability in contrast to the relative lack of fine-grain resolution of other important and landscape-heterogeneous environmental variables (e.g. soils for plant distributions, Mod et al., 2016). Finally, SDMs usually include climatic variables as monthly or annual averages of reference periods of 30 to 50 years (Fick & Hijmans, 2017;Hijmans et al., 2005;Karger et al., 2017), in most cases, irrespective of whether modelled species' lifespan is significantly longer or shorter than these periods (Roubicek et al., 2010). This disregard for temporal variability of the climatic predictors at sites where occurrences have been collected, has been suggested to have a big impact on the accuracy of SDMs (Bateman et al., 2012;Zimmermann et al., 2009).
Characterizing the environmental requirements using exclusively averages over reference periods could lead to underestimations of the real size of species' niche in the N-dimensional space. Not accounting for climate variability implies that, for example, two species withstanding the same average temperature for a given period have the same thermal niche, even if one of them lives under a constant climate every year, while the other suffers wide fluctuations in temperatures between years. Accordingly, the impact of including inter-annual variability in niche modelling may not be very important for species inhabiting areas with little temporal climatic variability (Hannah et al., 2014;Lembrechts et al., 2019). However, for species withstanding wide inter-annual climatic variability, differences between the niche characterized from climatic averages and the niche characterized including inter-annual variability is likely to be substantial (Hannah et al., 2014;Lembrechts et al., 2019). Similarly, species with small geographic distributions would be expected to show greater relative underestimation of their niche size when ignoring inter-annual variability than species with larger geographic distributions, where spatial variability could compensate for temporal variability (Hannah et al.,2014).
Underestimations of niche size caused by the non-inclusion of the inter-annual climate variability could have implications on different facets of ecology and biogeography. Since environmental suitability and population dynamics are theoretically related (Pulliam, 2000), niche breadth underestimation could contribute, among other ecological factors, to decouple climatic suitability from demographic processes (Csergő et al., 2017;Thuiller et al., 2014). For instance, environmental suitability estimated from average climates could predict population absence or decline under some environmental conditions where the population intrinsic growth rate is positive during some years which ensures the persistence of the population (Dullinger et al., 2012;Eriksson, 1996). Also, climatic variability could be particularly relevant when predicting biodiversity changes under future climate change scenarios in which the recurrence of climatic extremes is increasing (Coumou & Rahmstorf, 2012;IPCC Working Group1, 2014). Depending on the magnitude of extreme climatic events, the position of populations in the environmental space could move away from niche centres and may even fall outside of the species niche (i.e. become sink populations), reducing populations' capacity to survive if these conditions persist. But if niches are wider than traditionally considered, for a given change in climatic conditions the probability of populations to be pushed out of their niche space would be reduced.
Despite these substantial implications, studies that accounted for inter-annual climatic variability are generally rare (Hannah et al., 2014;Zimmermann et al., 2009). Some mechanistic modelling approaches (Dormann et al., 2012;Kearney & Porter, 2009;Kearney, Matzelle, & Helmuth, 2012) may implicitly consider temporal variability in species distribution, by characterizing the physiological relationship between species responses and environmental variables.

K E Y W O R D S
climatic averages, climatic suitability, demography, inter-annual variability, niche modelling, species niche In turn, range dynamics or dynamic occupancy models (Kéry, Guillera-Arroita, & Lahoz-Monfort,2013;Royle & Kéry, 2007;Schurr et al., 2012) may include occurrences information from yearly re-surveys. However, these approaches are not possible for the vast majority of species, due to the lack of experimentally-derived physiological response curves (Araújo et al., 2013) or the absence of temporal systematic data collection . Other studies have explicitly considered inter-annual climate to assess the relevance of temporal climatic variability on SDMs (Bateman et al., 2016;Bateman et al., 2012;Germain & Lutz, 2020;Margalef-Marrase, Pérez-Navarro, & Lloret, 2020;Zimmermann et al., 2009), but none of them has still provided an explicit comparison of the direct impact of including inter-annual climatic variability on niche size, and its impact when assessing the relationship between environmental niches and demographic responses.
Here, we attempted to fill these gaps by: (a) developing a novel procedure to characterize species niches based on long-term yearly climatic resolution and long-term climatic averages in a common environmental space; (b) assessing whether this inclusion of inter-annual variability improves the correspondence between niche-based environmental suitability modelled from averaged climate and population-level processes, illustrated by decay responses to an extreme drought event of Pinus halepensis' populations in southeast Spain; and (c), quantifying the relative change in niche size when considering inter-annual climatic variability versus averaged climate in 42 species with different distribution extents, corresponding to contrasting climatic ranges across the Mediterranean basin.

| MATERIAL AND ME THODS
We quantified species' climatic niches by projecting geographic occurrences into a two-or three-dimensional environmental space and then estimating kernel densities (Blonder et al., 2014;Broennimann et al., 2012). We calculated an "inter-annual variability-based" niche by considering for each species' occurrence the climatic values for every single year of a reference period (e.g. 1979-2013). In contrast, we obtained the "average-based" niche by considering only the average climate of every species' occurrence for the same period. (PNOA for year 2016) the extreme drought event, as the percentage of cover of dead and highly affected trees (dry canopy surface higher than 50%). These decay percentage values were then contrasted in the field for 14 plots of around 6,300 m 2 (total sampled surface of 87,964.59 m 2 ,  This two-dimensional environmental space was used to represent both the average Pinus niche and the Pinus inter-annual niche by projecting the scores of the climatic values of Pinus occurrences into this environmental space. We then used kernel density functions to determine the Pinus density for each cell of the two-dimensional environmental space, by applying a Gaussian kernel, selecting optimal bandwidth by cross-validation (Duong & Hazelton, 2005) and removing values below the 0.05 percentile. Finally, climatic suitability was estimated by dividing each density value in the niche space by the maximum density value within the niche, obtaining values ranging between 0 and 1.

| Climatic suitability and demographic responses in
We then extracted the climatic suitability of Pinus populations affected by drought for the extreme year 2013-2014, for which we had population decay data. Climatic conditions for this year were obtained from monthly precipitation and maximum, minimum and mean temperature records from between 68 and 114 weather stations of the Spanish Meteorological Agency (AEMET). This climatic dataset was translated into 1 km 2 resolution-maps following Ninyerola, Pons & Roure (2000) and using latitude, longitude and elevation as explanatory variables for climate. We then applied the "biovars" function (dismo package; Hijmans et al., 2016) to convert them into the final bioclimatic variables format. Finally, we selected the same 12 bioclimatic variables as used for niche characterization and translated these values into the two or three-dimension environmental space.

| Suitability-decay analyses
We applied GLMs with populations decay (binary or continuous) as response variable, and populations' climatic suitability during the extreme year (estimated from inter-annual or average niches) as explanatory variables. We produced four alternative GLMs: binary decay versus climatic suitability from average climatic niche, binary decay versus climatic suitability from inter-annual climatic niche, continuous decay versus climatic suitability from average climatic niche, continuous decay versus climatic suitability from interannual climatic niche. Each model used binomial error distribution and logit link. These models allowed us to compare the accuracy of the relationship between population decay and climatic suitability depending on niche characterization approach (inter-annual or average). These models were applied both for two-dimensional and three-dimensional niches.
In order to assess whether differences in accuracy between models were biased by the possible non-random locations of the observed populations in the environmental space, we compared the predictive capacity of both models (average niche based and inter-annual niche based) in 30,000 simulated decay datasets obtained by randomly creating subsets from the continuous decay dataset. In particular, differences in climatic suitability are expected to be maximal in the part of the environmental space not shared by average and inter-annual niches ( Figure S3), since average niche predicts a suitability of 0 while inter-annual niche predicts positive suitability values. So, for each randomly simulated dataset, we estimated the percentages of populations located in the non-shared space (ranging from 0% to 100%). Simulated dataset sizes ranged between 118 and 263 populations. For each simulated decay dataset, we also calculated the explanatory capacity of average and inter-annual niche-based climatic suitability, by applying GLM with continuous population decay as response variable and average or inter-annual climatic suitability as explanatory variables (see above explanation).
Finally, we tested the differences of models explanatory capacity by applying GLM with R 2 as response variable, niche model (average or inter-annual climate based), percentage of population located in the non-shared space between two niches, the interaction of these two variables, and population size for each simulated dataset as explanatory variables. This last was included in order to correct the possible low performance of dataset with lower number of populations.

| Influence of range size on change in average versus inter-annual niche
We characterized the climatic niche of 42 Mediterranean plant species with range sizes varying from the whole Mediterranean basin to very small ranges for endemic species (i.e. Mediterranean, occidental Mediterranean, Iberoafrican and Iberian Peninsula SE endemic species, Table S2), using the average and the inter-annual approaches described in the previous section. In addition, these species were selected for living under particularly variable climate between years. Species' occurrence data were collected from GBIF (GBIF, 2019, http://www.gbif.org) and from the herbarium of Jardí Botànic de Barcelona (https://museu cienc ies.cat/es/area-cient ifica/ depar tamen tos-cient ifico s/jardi n-botan ico/). Species' occurrence records were filtered to remove taxonomic inconsistencies (unappropriated species synonyms) and to reduce possible sampling bias and spatial autocorrelation. To do this, we selected each species filtering distance applying the ecospat.mantel.correlogram function (ecospat package; Di Cola et al., 2017) and selecting the minimal distance at which spatial autocorrelation is not significantly different than 0. Then, we retained only one occurrence per pixel of side equivalent to the selected filtering distance. For those species where the statistically selected filtering distance was lower than 1, we kept 1 km of pixel side, in order to maintain the spatial concordance between occurrence and climatic datasets. Species' final datasets ranged from 59 to 6,032 observations, with a total record of occurrences of 78,978.
Average and inter-annual climatic datasets were collected from CHELSA database (Karger et al., 2017) for the period 1979-2013.
We selected the same 12 bioclimatic variables as selected above for P. halepensis. We obtained the inter-annual climatic dataset by extracting the climatic values of every occurrence for each year of the period 1979-2013 (that is, each species occurrence × 35 years), and the average climatic dataset by estimating the mean climate of the 35-year period for every species occurrence. In this case, we calculated the common environmental space by using a PCA built using the inter-annual climatic data from all the occurrences of all 42 analysed species (PCA-occ sensu Broennimann et al., 2012) (Table S6) since we wanted to compare different species niches among them.
We estimated average and inter-annual niches by projecting the climatic values associated to each species occurrences into the first two PCA axes defining the climatic space, and using kernel density functions to determine species' density for each cell of the two-dimensional environmental space as described above, and removing density values below the 0.05 percentile.
We then calculated the ratio between the niche size (area) based on average and inter-annual approaches (niche area ratio = inter-annual niche area/average niche area). In order to test whether the relevance of including inter-annual variability varies depending on species average niche size, we applied a linear model (lm) with niches' area ratio as response variable and average-based niche area logarithmically transformed and species distribution range as response variables. Finally, we tested the linearity of the relationship between inter-annual niche size and average niche size using a log-log model (i.e. lm with log of inter-annual niche area as response variable and log of average niche area as explanatory variable) (Smith, 1980). We also analysed possible changes in niche size/background size ratio (i.e. the proportion that a species niche occupies of the available environmental space) when including inter-annual variability, by testing the linearity of the relationship between average niche size/average background size ratio and inter-annual niche size/background niche size ratio using a log-log model (see Supplementary Material niche/ background ratio explanation).

| Climatic suitability and demographic responses in Pinus halepensis
The first two PCA axes explained 60% of the variability of the 12 bioclimatic variables ( Figure S2). The inter-annual niche of P. halepensis was 42% larger than the average niche. During the extreme climatic year, 93.3% of unaffected forests and 63.3% of highly affected forests were inside the P. halepensis inter-annual niche (Figure 2a).
These values decreased to 52.8% for unaffected forests and 21.2% for highly affected ones when niche was calculated with average climate ( Figure 2b).
Models relating species decay with species climatic suitability during the extreme year show that populations suitability significantly explained species decay irrespective of the decay dataset (binary or continuous) or the type of niche (average or inter-annual) ( Figure 3, Table S2). However, models with inter-annual measures of climatic suitability better fitted with decay records, had considerably lower AIC, and slightly higher R 2 (particularly for continuous dataset, see Figure 3).
The explanatory capacity of models that include suitability from inter-annual climatic suitability was always higher than that of models including average climatic suitability, independently of the proportion of populations located in the non-shared area defined by the average-based and inter-annual variability-based niches ( Figure S3). Nevertheless, difference in explanatory capacity between the two type of models increased as the percentage of populations located in the non-shared area increased ( Figure S4 and Table S3).
In case of three-dimensional species niche, the explained variability of the environmental space increased up to 76.15%. Climatic suitability derived from three-dimensional niche was highly correlated to two-dimensional niche-based climatic suitability (Pearson correlation > 0,8 for average climatic suitability and > 0,95 for inter-annual climatic suitability, see Table S4), leading to similar results when explaining Pinus populations decay after extreme drought ( Figure S6 and Table S5).

| Influence of range size on change in average versus inter-annual niche
The two first PCA axes explained 60.8% of the variability of the 12 climatic variables describing the distribution of the set of 42 Mediterranean species ( Figure S7). Linear models across species with different distribution ranges showed a significant negative  effect of the species average niche size logarithmically transformed on inter-annual/average niche size ratio, implying that species with smaller average-based niche area (corresponding to species with more restricted distribution ranges), increased more their niche area when considering inter-annual climatic suitability than species with larger average-based niche area (i.e. wider distribution ranges) ( Figure 4 and Table S7). Accordingly, per distribution categories the endemic species were the ones which increased more their niche size when considering inter-annual variability. (Table S7). Log-log model results relating inter-annual variability-based niche size as a function of average-based niche size showed a regression slope significantly lower than 1 and higher than 0 (Table S8), which is associated to a positive but saturated relationship between average niche area and inter-annual niche area (Smith, 1980); so the higher values of average-based niche area, the lower increase of inter-annual niche size ( Figure S8). Similarly, log-log model relating the average-based niche/background ratio and inter-annual variability-based niche/ background ratio, also showed a regression slope lower than 1 and higher than 0 (Table S9), illustrating a saturated relationship rather than linear between the two ratios. Therefore, the higher the ratio of the study area occupied by species niche when including average climate, the lower the increase in the proportion that the inter-annual variability-based niches occupy within the inter-annual background region ( Figure S10 and S11).

| Including inter-annual variability in niche characterization
This study highlights some space. In addition to niche analyses in the environmental space, interannual variability could also be incorporated more generally into species distribution modelling approaches that allow for a hierarchical data structure (Kuznetsova et al., 2017;Wang & Maintainer, 2016) or by independently considering yearly climatic conditions at each site.
Finally, including climatic variability when building SDMs is particularly promising in climate change studies since it will help to better determine future species distribution under a climatic context of increasing variability ).

| Climatic suitability and demographic responses
Our results emphasize the importance of climatic variability to explain the relationship between demographic processes and climatic suitability, specifically under extreme climatic events. We found that P. halepensis' climatic suitability estimated from niches characterized by inter-annual climatic variability better explained species decay under extreme events in comparison to climatic suitability estimated from niches that only considered climatic average (Figure 3). It is theoretically assumed that if species' climatic niches are properly represented, there should be a correspondence between population demographic processes (such as growth or mortality rates) and niche parameters (Csergő et al., 2017;Pulliam, 2000;Thuiller et al., 2014). Multiples studies have tried to assess this relationship (Csergő et al., 2017;Thuiller et al., 2014), particularly under extreme events (Lloret & Kitzberger, 2018;Margalef-Marrase et al., 2020;Pérez Navarro et al., 2018), but this correspondence does not always emerge (van der Maaten et al., 2017;Thuiller et al., 2014).
Very often these niche-demography relationships, particularly in case of parameters related to growth, are affected by many local aspects-such as favourable microhabitats or biotic interactions-that impede identifying a clear relationship between niche and population performance (Csergő et al., 2017;Pearson & Dawson, 2003).
Nevertheless, neglecting temporal variability in this kind of studies could also introduce a substantial error when predicting population performance in fluctuating environments (Niehaus et al., 2012).
Although suitability derived from average-based niche may be robust enough to explain dramatic demographic responses under extreme climatic episodes, as we showed in our results (see also Lloret & Kitzberger, 2018;Pérez Navarro et al., 2018;Sapes et al., 2017), the better the characterization of the niche variability, the better the species fitness' predictions.
In addition, the inclusion of inter-annual variability when predicting species responses was especially relevant for populations located close to the niche margins where such populations can be included or excluded from the niche depending on yearly variation in climate (i.e. in the non-overlapping space between inter-annual and average-based niches) ( Figure S3). This area of the species niche corresponds to unsuitable climatic conditions according to average-based niche but suitable climatic conditions when climatic variability is included in niche characterization. The average-based niche does not show a good predictive capacity for demographic responses in this non-shared area (between the two types of niches), by considering this area unsuitable for species presence.
Thus, if the study sample was composed only by populations from this non-shared area, the only model that could predict populations response correctly would be the one including inter-annual climatic variations ( Figure S4). This mismatch also implies that populations located in this area, which in many cases (but not necessarily) also correspond to populations living at geographic range margins, could be wrongfully located out of the niche space when using average-based species niche.

| Influence of range size on change in average versus inter-annual niche
These above-mentioned implications of including inter-year variability could be even more important for species with narrow distribution ranges, which implies smaller average-based niches compared with those with broader distributions. Our results showed that, in the Mediterranean basin, species with smaller average-based niche areas increased relatively more their niche area when considering inter-annual climatic variability (Figure 4, Figure S8, S11 and Tables S8   and S9), probably because in species with narrow distribution ranges, spatial climatic variability does not compensate for temporal variability. Therefore, endemic and rare species tend to increase their niche size more when including yearly climatic variability, particularly when inhabiting highly fluctuating climates.

| Other implications of inter-annual variability versus average-based climatic models
This study highlights limitations, barely considered to date, of using climatic averaged datasets of 30-50 years' periods (Fick & Hijmans, 2017;Hijmans et al., 2005;Karger et al., 2017) when modelling species' niches and distributions (Hannah et al., 2014;Zimmermann et al., 2009). Not accounting for the whole temporal climatic resolution at sites of species occurrences when modelling climatic niche could result in errors affecting the characterization of species geographical distribution and losses of predictive power.
Underestimations of then niche envelope could lead to underestimations of species' environmentally suitable areas in geographic space, and such errors could also be propagated into related management plans (e.g. selection of favourable areas for protecting species or assisted migrations; Guisan et al., 2013). Importantly, these possible niche underestimations could vary depending on species distribution ranges, being particularly relevant for rare species. In addition to the error derived from the non-inclusion of climatic variability, using systematically 30-50 years averaged climatic periods could be especially pernicious when characterizing the climatic niche of short-lived species, as in pests (Jaime et al., 2019)-since these emerge explosively as a consequence of specific climatic conditions that appear in particular years of the period-, or in case of non-sessile organismswhich can buffer their environment variability by movement, thus only perceiving the climate of the locations during specific favourable time periods (Lembrechts et al., 2019)-. In these cases, it is particularly important to know the specific moment of species occurrence and use the corresponding temporal resolution (Bateman et al., 2016).
Estimates of invasion risk or extinction rate could also be affected by niche characterization errors. In this case, niche size underestimation could lead to underestimate the environmentally suitable habitat of invasive species in new regions, or to overestimate extinction risk (Feeley & Silman, 2011) as populations and species could be able to survive under changing climates or extreme events more than expected with average-based niches.

ACK N OWLED G EM ENTS
We are grateful to Oficina de Impulso Socioeconómico del Medio

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13207