Microhabitat analyses support relationships between niche breadth and range size when spatial autocorrelation is strong

724 –––––––––––––––––––––––––––––––––––––––– © 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: Dan Warren Editor-in-Chief: Miguel Araújo Accepted 6 January 2020 43: 724–734, 2020 doi: 10.1111/ecog.04798 43 724–734


Introduction
The idea that species with large tolerance to environmental variables are more widespread than specialized species is one of the key hypotheses in macroecology, and has been corroborated by a large body of evidence (Brown 1995, Slatyer et al. 2013, Kambach et al. 2019. A comprehensive meta-analysis observed consistent, positive relationships between niche breadth and range size across a broad range of organisms and using multiple measures of niche, suggesting that this relationship can be a general ecological pattern (Slatyer et al. 2013).
However, recent studies have highlighted that correlations between niche breadth and range size can emerge just as a by-product of strong spatial structure of environmental variables (Boucher-Lalonde and Currie 2016, Di Marco and Santini 2016, Moore et al. 2018b, Saupe et al. 2019. For instance, Moore et al. (2018b) found significant relationships between range size of South African plants, and the breadth of their climatic niche. Nevertheless, they also noted that similar relationships can be produced by spatial null-models in which random layers, having a spatial structure similar to climatic variables, were used instead of true climate. The possibility that spatial autocorrelation positively biases the niche breadth-range size (NB-RS) relationship casted the doubt that some of the observed NB-RS patterns could be artefactual (Moore et al. 2018b). Unfortunately, spatial autocorrelation is one of the most pervasive features of natural systems, and many parameters measured in ecological studies (e.g. species distribution, climate, distribution of resources, soil features…) generally are characterized by strong autocorrelation (Wagner and Fortin 2005). This makes it challenging to test the robustness of the NB-RS relationships.
The NB-RS relationship has often been evaluated by measuring the species' bioclimatic niche on the basis of broad scale bioclimatic variables (hereafter: macro-climatic; Slatyer et al. 2013, Moore et al. 2018b. This is probably related to both the growing interest on bioclimatic (Grinnelian) niche during the last years (Soberón and Nakamura 2009), and to the easy availability of these data over broad spatial extents, particularly through geographic information systems (Karger et al. 2017). However, macroclimatic data have some limitations. First, species often interact with their environment at scales finer than most of available macroclimatic data, thus macroclimatic layers only provide an imperfect representation of the conditions actually perceived by study organisms (Potter et al. 2013, Bennie et al. 2014, Scheffers et al. 2014. Second, spatial autocorrelation of macroclimatic data is generally very high, thus assessing whether NB-RS relationships are genuine, or are just a by-product of autocorrelation, can be particularly difficult (Di Marco andSantini 2016, Moore et al. 2018b). Micro-habitat analyses provide a complementary approach to measure species niche. When performed at the scale at which individuals interact with their environment, microhabitat analyses provide a measure of conditions actually experienced by individuals, and can take into account processes such as microhabitat selection and species thermoregulation, which allow individuals to maintain conditions within their eco-physiological tolerance limits in landscapes characterized by strong heterogeneity (Scheffers et al. 2014, Ficetola et al. 2018, Moore et al. 2018a, Rubalcaba et al. 2019a. Furthermore, microhabitat variables often are characterized by strong spatial heterogeneity, as nearby sites can have very dissimilar microhabitats. The reduced spatial structure of microclimate may be a key feature to assess the goodness of NB-RS relationships (Moore et al. 2018a). However, we are not aware of studies assessing whether microhabitat analyses can help distinguishing between true NB-RS relationships and the artefactual effects of spatial autocorrelation.
In this study, we analyzed patterns of niche breadth across the European terrestrial salamanders [genus Hydromantes, subgenera Atylodes and Speleomantes; see Wake (2013) for details on nomenclature], comparing patterns obtained at the macroclimatic and at the microhabitat level. During summer, these salamanders move to underground spaces, where microhabitat features can show strong variation over just a few meters (Ficetola et al. 2018, Mammola and Leroy 2018, Sánchez-Fernández et al. 2018, Supplementary material Appendix 1 Fig. A1). Salamanders select the microhabitats with abiotic conditions (e.g. temperature, humidity, light) within their ecophysiological limits, thus analyses of the microhabitat within sectors occupied by salamanders provide a good representation of their niche requirements at fine scale (Lunghi et al. 2016, Ficetola et al. 2018. The aims of this study were 1) to test the validity of the NB-RS hypothesis at different levels, using macroclimatic and microclimatic niche; 2) to test whether a positive NB-RS relationship can be generated by spatial autocorrelation of environmental variables. We predict that this pattern would be strongest when niche is measured using the environmental variables with highest autocorrelation (Moore et al. 2018b). Finally, 3) to assess whether fine-scale niche analyses can help to ascertain whether NB-RS relationships are genuine, or emerge just as a by-product of spatial autocorrelation (Moore et al. 2018b), since microclimatic variables are expected to have weaker spatial autocorrelation.

Study system
European terrestrial salamanders (genus Hydromantes) are a clade comprising eight allopatric species of Plethodontid salamanders living in Italy and France; most of them are endemic to small areas, although range size shows a variation over three orders of magnitude across species (50-18 000 km 2 ; Fig. 1). These salamanders are active at the surface during cool and wet seasons (from autumn to spring). However, during summer external conditions become too harsh (dry and hot) for outdoor activity and salamanders move to belowground environments, where they select sectors with conditions of temperature, humidity and light intensity within their preferred range (Ficetola et al. 2018). All species are associated with sectors characterized by low light intensity, cold temperature and high humidity, still species have different realized niches, with some species living in darker and more humid sectors, and other being more tolerant to dryness and light (Ficetola et al. 2018(Ficetola et al. , 2019. In this study, we used distribution data covering the whole range of all the species, and cave surveys performed during summer, to measure niche breadth at both the macroclimatic and microhabitat scales.

Macroclimatic data
To analyse the macroclimatic niche, we used an updated version of the database of Ficetola et al. (2018), which gathered high-resolution distribution records covering the whole range of all the European Hydromantes species. For each presence locality, we extracted four bioclimatic parameters: mean annual temperature, temperature seasonality, summed precipitation and precipitation seasonality, obtained from the 30 arc-second resolution Chelsa-Clim database, which provides accurate estimates of precipitation in landscapes with high topographic complexity (Karger et al. 2017). Previous studies have shown that these variables allow an excellent representation of salamander niche, and that variation in these parameters can strongly affect fitness-related traits (Lunghi et al. 2018). We only considered a single locality per each 30 arc-second cell; localities within the narrow hybrid zone between H. ambrosii and H. italicus (Ficetola et al. 2019) were excluded from the study. Temperature seasonality and summed precipitation were log-transformed before analyses to improve normality.

Microhabitat data
To measure species niche at the microhabitat scale, we surveyed caves covering the whole range of all the European Hydromantes species (Fig. 1). We performed surveys in early summer (mostly June-July 2011-2018), when the detection of salamanders is the easiest. Seasonal variation of underground microclimate is very weak (Lunghi et al. 2015, Mammola 2019; previous analyses showed that these salamanders select similar microhabitat conditions through the year, and that niche estimates obtained from data collected in summer are similar to estimates from the other seasons (Lunghi et al. 2015). Furthermore, these salamanders generally are at equilibrium with their environment for water and temperature and, in the field, temperature differences between air and body temperature generally are < 0.5°C (Spotila 1972, Lunghi et al. 2016, suggesting that air conditions are excellent proxies of operative conditions of individuals (Sunday et al. 2014, Lunghi et al. 2016). Surveys were always performed during the central hours of sunny and dry days. We subdivided caves in 3-m longitudinal intervals (hereafter: sectors). The size of sectors roughly corresponds to the home range size of individuals during their underground activity (Lanza et al. 2006), therefore observations are unlikely to represent transient individuals. Surveyed sectors covered the whole cave, or all the sectors until the first empty one after the last salamander. Within each sector, we performed visual encounter surveys to detect active salamanders, and recorded four abiotic variables affecting their distribution. Relative air humidity (%; accuracy: 0.1%) and air temperature (°C; accuracy: 0.1°C) were measured using EM882 and Lafayette TDP92 multi-function devices. Furthermore, we measured minimum and maximum incident light (illuminance, measured in lux; accuracy 0.01 lux) in the portions of the sector receiving more and less light, respectively, using the EM882. A few caves were surveyed multiple times; for these caves we considered only one survey, selecting that in which the largest number of sectors was assessed. See Ficetola et al. (2018) for additional details on sampling. Incident light was log-transformed before analyses to improve normality.

Statistical analyses
Environmental variables were strongly collinear (Supplementary  material Appendix 1 Table A1), thus we used an approach based on principal component analysis to measure  (Boulangeat et al. 2012, Kambach et al. 2019. For macroclimatic data, we reduced the dimensionality of the four bioclimatic variables using principal component analysis (PCA). The first two PCA axes explained 90.9% of the overall climatic variation across presence localities. The first PCA axis was positively related to mean temperature and precipitation seasonality, and negatively related to annual precipitation and temperature seasonality; the second PCA axis was positively related to mean temperature and temperature seasonality, and negatively related to annual precipitation and precipitation seasonality (Supplementary material Appendix 1 Table A2a). Following Boulangeat et al. (2012) we calculated niche breadth at the bioclimatic level as the area of the minimum convex polygon (MCP) encapsulating the presence points in which each species occurred across the climatic space. In order to control for differences in sample size among species, and to down-weight the influence of outlier plots, we calculated niche breadth as the average MCP of 500 random draws of five presence points in which each species occurred (averaged area of the polygons; Boulangeat et al. 2012, Kambach et al. 2019. As the number of five presences is somewhat arbitrary, we repeated analyses using 10 presence points, obtaining very similar results. Previous studies and preliminary analyses on our data have shown that results are highly robust and consistent, regardless to the number of presence points (5, 10 or 15 points) (Boulangeat et al. 2012). We used the same approach to calculate niche breadth on the basis of microhabitat parameters, by performing a PCA on the four microhabitat variables (variation explained by the first two PCA components: 81.6%). In this case, the first PCA axis was positively related to temperature and light, and negatively related to humidity; the second PCA axis was positively related to temperature, and negatively related light and humidity (Supplementary material Appendix 1 Table A2b).
The range area of each species was calculated on the basis of the range polygons of species; polygons were generated from the presence records of each species using α-hulls (Ficetola et al. 2014). α-hulls are particularly suited to describe species ranges, which can have non-convex structures (Burgman and Fox 2003). Subsequently, we assessed the relationship between niche breadth and range size. Related species may share similar biological traits because they inherit them from their ancestors, thus data cannot be analysed by tools assuming that each species is statistically independent (Freckleton et al. 2011). Therefore, relationships between niche breadth and range size were assessed using generalised phylogenetic least squares (PGLS), a method taking into account phylogenetic non-independence that accurately estimates parameters with low type 1 error (Freckleton et al. 2011). To integrate phylogeny into analyses, we used the calibrated tree published by Carranza et al. (2008). Both niche breadth and range size were log-transformed to improve normality. In PGLS, we estimated regression parameters and phylogenetic signal simultaneously, using maximum likelihood to assess phylogenetic signal (λ). PGLS was built using the caper package in R 3.5 (Orme et al. 2018).

Assessing the potential effects of spatial autocorrelation: null-models
Environmental data are often characterized by strong spatial structure, and this can lead to spurious relationships between range size and niche width. We therefore used null-model tests to assess the robustness of relationships between niche width and range size. At the macroclimatic level, we used the spatial null-model developed by Moore et al. (2018b) to generate randomized climatic layers with the same mean, variance and spatial structure of the real macroclimatic variables. We modelled the expected value of each macroclimatic variable as a continuous Gaussian field, which was defined by an intercept (fixed effect) equal to the average value of the variable, and using Matérn correlation to describe its spatial autocorrelation. The parameters of the Matérn process were estimated using the integrated nested Laplace approximation, as implemented in the R-INLA package (Bivand et al. 2015, Rue et al. 2017. After fitting one model per macroclimatic variable, we used them to simulate 250 realizations of each environmental layer (pseudo-macroclimate); see the Supplementary material of Moore et al. (2018b) for additional details and complete scripts. We then used the aforementioned PCA-based approach to estimate the niche breadth of each species on the basis of the randomly generated (pseudo-macroclimate) layers. To assess whether the observed relationship is significantly greater than the values observed under the null-model, we compared the observed standardized coefficient of the PGLS with the ones obtained under the spatial null model (Moore et al. 2018b). The observed relationship was considered significant if falling outside the 95% CI generated through null-models.
The Moore et al. (2018b) spatial null-model has been developed to generate raster layers with specific spatial structure, and can be successfully applied to macroclimatic layers. However, microhabitat variables have features that do not allow the direct application of the Moore et al. (2018b) approach. First, each cave has one single set of associated coordinates (corresponding to cave entrance), still multiple sectors are measured within each cave. Second, the spatial null-model is optimized for the analysis of grid layers with regularly spaced cells, while surveyed caves are irregularly spaced. We thus modified the Moore et al. (2018b) nullmodel to make it applicable to microhabitat data and to generate pseudo-microclimatic data with the same features and spatial structure of the true data. First, we built an empty grid (resolution: 30 arc-seconds, i.e. the same resolution of macroclimatic data) covering the whole range of the study species. For each grid cell where at least one surveyed cave occurs, we extracted the minimum and the maximum value recorded at the four microhabitat variables across all the sectors surveyed within the grid cell. Therefore, we generated eight grid layers: max and min temperature; max and min humidity; max and min values observed for both maximum and minimum incident light. When multiple caves were sampled within the same cell, we considered the sectors of all the caves. We then applied the Moore et al. (2018b) spatial null-model to these eight grid layers, to generate random microhabitat variables retaining the same spatial structure of the observed microhabitat variables; the result was a set of eight simulated micro-climatic layers (i.e. min and max values of temperature, of humidity…; Supplementary material Appendix 1 Fig. A2). For each microhabitat variable, the values observed at cave i under the spatial null-model were then drawn from a uniform distribution with min and max values constrained between the min and max values of that variable observed under the null-models in the cell where cave i is located. The distribution of min and max values of microhabitat variables showed relatively high autocorrelation, therefore constraining the simulated microhabitat variables between the min and max values observed in nature allowed to generate variables characterized by autocorrelation levels comparable to the ones of true variables (Supplementary material Appendix 1 Fig. A2). In order to evaluate whether this spatial nullmodel produces data with a spatial structure similar to the true data, for each macro-and microhabitat variable we compared the correlogram of true data with the one of variables generated by the null-model, using the EcoGenetics package in R (Roser et al. 2017).
Environmental variables showed a high degree of multicollinearity (Supplementary material Appendix 1 Table A1), and a null-model that does not match the structure of the real data could lead to suspect conclusions. We therefore used the condition index to measure the collinearity of both real and simulated environmental variables; values of condition index > 30 indicate the occurrence of collinearity (Dormann et al. 2013). For microhabitat, real variables showed a higher collinearity than the simulated ones (Supplementary material Appendix 1 Fig. A3). To confirm that our conclusions were not biased by collinearity, we generated a second set of simulated variables, mimicking the collinearity of true variables. In the real microhabitat variables, temperature was strongly correlated to humidity, the two light measures were related to each other, while light and temperature were uncorrelated (Supplementary material Appendix 1 Table A1b). Therefore, the simulated collinear dataset included: simulated temperature, simulated maximum light, and two simulated variables collinear to temperature and maximum light. The collinearity between simulated variables was controlled using decay parameters of 8.1 and 2.3 for temperature and light, respectively, as this ensured obtaining a collinearity structure matching the one of real data (see Supplementary material Appendix 1 Fig. A4 for additional details).
Differences in spatial autocorrelation among microclimatic variables is an additional potential issue, as they could affect the multi-dimensional nature of the niche, potentially influencing PCA results. For microclimatic variables, we therefore repeated the NB-RS analyses using univariate tests. Univariate niche breadth was measured as the range of each variable, after excluding the 5% most extreme values to take into account differences in sample size among species. We then used PGLS to assess univariate NB-RS relationships, and compared the observed NB-RS relationships with those obtained under the spatial null-models.

Macroclimatic niche versus range size
Macroclimatic analyses of niche were based on 573 salamander detections in 303 different grid cells (Supplementary material Appendix 1 Table A3, A4). At the macroclimatic scale, niche breadth varied over one order of magnitude across species. We found a strong, positive relationship between niche breadth and range area (PGLS: t 6 = 2.81, p = 0.03, R 2 0.57); the species with the smallest range (H. sarrabusensis) also showed the narrowest niche, while the species with broadest range (H. italicus) showed the broadest macroclimatic niche (Fig. 2a).
The spatial autocorrelation of macroclimatic variables was very high (Moran's I always > 0.9; Fig. 3a). The pseudo-climatic data, generated using the spatial null-model, showed very similar, strong autocorrelation (Fig. 3b). Simulated macroclimatic variables also showed collinearity values similar to those observed for the true variables (Supplementary material Appendix 1 Fig. A3a). We observed a positive relationship between simulated macroclimatic data and range area. Using simulated macroclimatic variables, the median effect size of the NB-RS relationship was 3.433, and 100% of simulated effect sizes were positives (95% simulation envelopes: 0.63-9.06; Fig. 4a). The relationship between macroclimatic niche and range area observed with real data was not significantly different from that produced by null-models (Fig. 4a).

Microhabitat niche versus range size
In microhabitat surveys, we surveyed 353 caves and detected salamanders in 575 out of the 1669 surveyed cave sectors (Supplementary material Appendix 1 Table A3, A4); niche breadth varied over one order of magnitude across species (Fig. 2b). We observed a strong, significant relationship between niche breadth and range size also at the microclimatic scale (PGLS: t 6 = 3.04, p = 0.02, R 2 0.61). Also in this case, the species with the smallest range (H. sarrabusensis) showed the narrowest niche, while the two species with broadest range (H. strinatii and H. italicus) showed the broadest microclimatic niche (Fig. 2b).
Although the spatial autocorrelation of all the microhabitat variables was significantly higher than zero (for all variables, p < 0.005), Moran's I values were much lower than those observed for macroclimatic parameters (range: 0.03/0.59; Fig. 3a). The pseudo-microhabitat data, generated using the spatial null-model, showed autocorrelation similar to the one of environmental variables (Fig. 3b). The relationships between simulated microhabitat and range area was generally positive; using random microhabitat variables the median effect size was 0.93 (95% simulation envelopes: −1.37/3.01), and 79% of simulated effect sizes were positive (Fig. 4b). Nevertheless, the true relationship between microhabitat niche and range area was significantly stronger than the one obtained through null-models (98% of simulated effect sizes were lower than the observed relationship; Fig. 4b). The collinearity of simulated microclimatic variables was lower than the one observed in the true variables (Supplementary material Appendix 1 Fig. A3b). However, differences in collinearity did not affect our conclusions, as results were nearly identical when we used the null-model with the same collinearity of real data (Supplementary material Appendix 1 Fig. A4, A5).
Results of univariate tests, analysing separately the four microclimatic variables, were generally weaker than the results of the multivariate analysis. The NB-RS relationship was positive for all the microclimatic parameters, but was significant for temperature only (Supplementary material Appendix 1 Table A4). The observed univariate relationships were generally more positive than the ones of spatial nullmodels, but were weaker than the multivariate one, and laid within the 95% simulation envelopes (Supplementary material Appendix 1 Fig. A6).

Discussion
The positive relationship between range size and niche breadth has been proposed as a general pattern, with important implications for both ecological theory and management (Slatyer et al. 2013). However, recent analyses showed that strong autocorrelation can inflate the likelihood of generating positive NB-RS relationships, thus challenging the generality of this pattern (Moore et al. 2018b). Our study provides a unique combination of data describing species distribution at multiple scales, including exhaustive macroclimatic and distribution data, and detailed (3 m resolution) information on the microhabitat actually selected by individuals. This allows comparing NB-RS relationships obtained using environmental data with different resolution and autocorrelation levels. On the one hand, we confirm that strong autocorrelation of macroclimatic variables makes it difficult disentangling true NB-RS relationships from the mere by-product of autocorrelation. On the other hand, multivariate microhabitat analyses can be an important step to solve this conundrum, given their weaker autocorrelation and the direct links between microhabitat features and the ecophysiology of individuals (Rubalcaba et al. 2019b).
In European terrestrial salamanders, the relationship between range size and niche width was evident. Despite the limited number of species, this taxon shows a remarkable variation in range size: the species with the smallest distribution (Hydromantes sarrabusensis) is among of the European vertebrates with smallest range (just 50 km 2 ), while the range of H. italicus is 400-times larger, a feature that makes this species relatively widespread within amphibians (~60% of world's amphibians have smaller ranges; IUCN 2019). Differences in range size closely match the strong variation in niche breadth, which spanned over one order of magnitude across species (Fig. 2). The narrow niche of localized species was not a by-product of small sample size, given that niche breadth was measured using a permutation procedure that allows taking into account differences in sample size (Boulangeat et al. 2012, Kambach et al. 2019. Two main hypotheses can thus explain the broader niche of large ranged species. First, generalist species are able to exploit a broader range of environmental conditions, and can thus occupy large ranges (Slatyer et al. 2013). Second, strong spatial autocorrelation of environmental data could artificially lead to higher values of niche breadth in species with large ranges (inflated niche breadth-range size relationship; Moore et al. 2018b).
The macroclimatic analyses do not allow teasing apart these two hypotheses, given that randomly generated climatic layers produced relationships similar, and sometimes stronger, to the ones obtained with the true variables (Fig. 4a). However, the support of these contrasting hypotheses can be better evaluated by combining macroclimatic and microhabitat data. Both macroclimatic and microhabitat niches are described using abiotic scenopoetic parameters, and both include variables representing thermal energy and water availability which, along with light, are the main determinant of salamander activity (Spotila 1972, Buckley et al. 2012, Fisher-Reid et al. 2012, Ficetola et al. 2018). On the one hand, we found a strong positive relationship between niche breadth measures obtained using the macroand microhabitat approaches (Fig. 2c), in accordance with Figure 3. Spatial autocorrelation (Moran's I) of (a) macroclimatic and microhabitat variables used to measure niche breadth, and (b) variables randomly generated using the spatial null-model. Error bars are 95% confidence intervals. macroecological predictions on the processes determining species' distributions (Barnagaud et al. 2012). On the other hand, the autocorrelation of microclimatic data is moderate to weak (Fig. 3a), thus in this case it is less likely that positive NB-RS relationships emerge just as a by-product of spatial structure. Besides the weak autocorrelation, microhabitat analyses also have the advantage of providing direct information on the conditions actually experienced by individuals (Lunghi et al. 2016), and can be a good proxy of operative eco-physiological conditions, which can be used to obtain measures of fundamental niches (Kearney and Porter 2009).
The fact that analyses performed using macro-and microclimatic variables produced similar NB-RS relationships supports the hypothesis of a genuine effect of niche breadth on range size (Slatyer et al. 2013). The macro-and microhabitat approaches describe non-identical aspects of the niche. For instance, close to the surface differences between air temperature and mean annual temperature can be strong (Lunghi et al. 2015), thus transferring the results at these two levels is challenging (Ficetola et al. 2018, Mammola 2019. Nevertheless, some studies found similarity between niche parameters measured at fine (body temperature, microclimate) and broad spatial scales (Kozak andWiens 2007, Fisher-Reid et al. 2013). Therefore, it is possible that species with large tolerance for temperature are active in caves with a variety of temperature conditions, and can also have ranges spanning across broad climatic gradients.
Using null-models to test ecological hypotheses can be challenging, as inappropriate null-distributions may lead to excessively conservative or liberal outcomes (Ulrich et al. 2009, Zurell et al. 2010. The generation of spatial nullmodels is particularly complex for microhabitat, because of the difficulty of modelling its fine-scale variation using regular grid cells. For instance, despite simulated microhabitat matched well the overall pattern of spatial autocorrelation (Fig. 3), some realizations showed diverging autocorrelation values, and this could in principle affect the conclusions of analyses. Nevertheless, our results remained highly consistent when we used alternative models with different assumptions (Supplementary material Appendix 1 Fig. A5), supporting the robustness of our conclusions. Results of univariate tests were similar among microhabitat variables, as they all showed positive NB-RS relationships. However, patterns at the four microhabitat variables were non-identical, with temperature and humidity showing larger autocorrelation than light (Fig. 3). In turn, temperature was the only variable with a significant univariate relationship, suggesting that autocorrelation may play a role in driving NB-RS relationships also for microhabitat parameters, and that also the outcome of microhabitat analyses should be considered with caution. None of the univariate relationship showed effect size significantly higher than the null-models. This is apparently in contrast with the multivariate test (Fig. 4), which showed a higher effect size than the univariate ones. The multivariate NB-RS relationship was not driven by one single microhabitat variable, but instead represented the joint effect of the four variables (Supplementary material Appendix 1 Table A6). In fact, measures of niche breadth obtained from single microhabitat variables are not necessarily correlated. For instance, H. sarrabusensis was the species with the smallest geographic range, and also exhibited the narrowest niche for light intensity, while another small-ranged species (H. supramontis) showed the narrowest niche for temperature (Supplementary material  Appendix 1 Table A3, A4). The multivariate measure of niche breadth was well correlated to all the univariate measurements (Supplementary material Appendix 1 Table A6). The different univariate measurements represent tolerance to different limiting factors (temperature, light, water availability), each of which can affect species fitness in different ways (Peterman and Semlitsch 2014, Lunghi et al. 2018, Ficetola et al. 2019; see below for discussion). The multivariate measurement jointly considers the most important drivers of salamander distribution and helps in better representing the multivariate nature of the niche (Soberón and Nakamura 2009). In fact, if a species has a great tolerance to dry environments, but only a narrow thermal niche, it will probably be able to exploit narrower ranges than species with good tolerance to both variables. This can explain why the multivariate test shows the strongest effect size, and the fact that it is significant even though none of the univariate relationships is stronger than the null-models. Terrestrial salamanders underwent a complex pattern of range dynamics since the end of glaciations, and their present-day distribution is determined by the interplay between geological features and climatic fluctuations during the Quaternary (Chiari et al. 2012, Cimmaruta et al. 2015. Tolerance to the different microhabitat parameters can provide different advantages and can impact both the fitness and dispersal of individuals. In turn, these effects can have consequences on species ranges (Peterman and Semlitsch 2014). For instance, tolerance to temperature and/or humidity can increase salamanders outdoor activity during warm and dry seasons, thus enhancing the probability of long-term dispersal between mountain massifs, and also allows species persistence in arid areas of the range (Peterman and Semlitsch 2014). Conversely, tolerance to light can increase activity during daytime or nearby the surface, improving foraging performance or facilitating dispersal among subterranean habitats (Lunghi et al. 2018, Ficetola et al. 2019. Therefore, multivariate analysis of niche can better describe the multifaceted processes that can enable the more tolerant species to exploit broad geographical ranges (Welsh and Droege 2001).
Plethodontid salamanders are a peculiar group, with small ranges and strongly dependent on belowground ecosystems, still our conclusions can be relevant for a very broad range of organisms. Despite being rarely the focus of macroecological studies, belowground environments host a very large number of species, that have a key role in the functioning of ecosystems (Beck et al. 2012, Bardgett andvan der Putten 2014). Furthermore, belowground organisms are not the only ones strongly dependent on microhabitats. For instance, microclimatic conditions can be more important than macroclimate for many ectothermic vertebrates and arthropods (Scheffers et al. 2014, Sunday et al. 2014, thus our conclusions can be relevant for several other taxa. Previous work suggested weak correspondence between macroclimatic and microhabitat niche in terrestrial salamanders (Ficetola et al. 2018), while in the present study we found excellent match between niche breadth measured at the macro-and microhabitat scales (Fig. 2c). However, Ficetola et al. (2018) focused on the average conditions occupied by populations (niche position), while here we considered niche breadth. Position and breadth are different properties of the niche, and species with similar niche breadth can have different positions, or viceversa (Rannap et al. 2009, Gómez-Rodríguez et al. 2015. Obtaining accurate measures of niche position on the basis of macroclimatic data may be particularly difficult for species performing active selection of specific microhabitats, because of the strong difference between broad-scale and fine-scale measurements of climatic parameters (Potter et al. 2013, De Frenne et al. 2019. Conversely, the match between analyses at the two scales could be better for niche breadth, given that a broad tolerance can allow species to exploit both a wide range of microhabitats, and a range encompassing multiple climatic regimes.
A growing number of studies are showing that niche breadth can impact multiple population and species parameters (Gómez-Rodríguez et al. 2015, Saupe et al. 2015, Qiao et al. 2016, Ralston et al. 2016, Rolland and Salamin 2016. However, it is also evident that a complex structure of environmental data, such as autocorrelation, makes it challenging testing the hypotheses on relationships between niche breadth and species distributions. The combination of macroclimatic and microhabitat data helps exploring niche variation at multiple scales, thus providing a unique opportunity to understand how environmental conditions shapes the distribution of species, and allows a better identification of the complex factors that actually determine the range of species.

Data availability statement
Species distribution and microhabitat data used in this study are available as Supplementary material Appendix 1 Table A4. To avoid illegal poaching on protected species, we degraded the quality of coordinates reported in the tables (Lunghi et al. 2019). Macroclimatic data are available from Karger et al. (2017).