Incorporating intraspecific variation into species distribution models improves distribution predictions, but cannot predict species traits for a wide‐spread plant species

60 –––––––––––––––––––––––––––––––––––––––– © 2019 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: Christine Meynard Editor-in-Chief: Miguel Araújo Accepted 11 September 2019 43: 60–74, 2020 doi: 10.1111/ecog.04630 doi: 10.1111/ecog.04630 43 60–74


43: 60-74, 2020
doi: 10.1111/ecog.04630 The most common approach to predicting how species ranges and ecological functions will shift with climate change is to construct correlative species distribution models (SDMs). These models use a species' climatic distribution to determine currently suitable areas for the species and project its potential distribution under future climate scenarios. A core, rarely tested, assumption of SDMs is that all populations will respond equivalently to climate. Few studies have examined this assumption, and those that have rarely dissect the reasons for intraspecific differences. Focusing on the arctic-alpine cushion plant Silene acaulis, we compared predictive accuracy from SDMs constructed using the species' full global distribution with composite predictions from separate SDMs constructed using subpopulations defined either by genetic or habitat differences. This is one of the first studies to compare multiple ways of constructing intraspecific-level SDMs with a species-level SDM. We also examine the contested relationship between relative probability of occurrence and species performance or ecological function, testing if SDM output can predict individual performance (plant size) and biotic interactions (facilitation). We found that both genetic-and habitat-informed SDMs are considerably more accurate than a specieslevel SDM, and that the genetic model substantially differs from and outperforms the habitat model. While SDMs have been used to infer population performance and possibly even biotic interactions, in our system these relationships were extremely weak. Our results indicate that individual subpopulations may respond differently to climate, although we discuss and explore several alternative explanations for the superior performance of intraspecific-level SDMs. We emphasize the need to carefully examine how to best define intraspecific-level SDMs as well as how potential genetic, environmental, or sampling variation within species ranges can critically affect SDM predictions. We urge caution in inferring population performance or biotic interactions from SDM predictions, as these often-assumed relationships are not supported in our study.
Keywords: climate change, genetic variation, habitat suitability, local adaptation, species distribution modeling, traits

Introduction
Discerning how and where populations will respond to climate change is now a central topic in ecological research, with great interest in applying this knowledge to inform conservation and management decisions in order to mitigate species extinction risks. A common approach is to determine how the potential distribution of a given species will shift in space with climatic changes using correlative Species Distribution Models (SDMs; Pacifici et al. 2015). Such SDMs correlate a species' occurrences to current climate in order to predict the species' relative probability of occurrence (P occ ) in space and time (reviewed by Wiens et al. 2009). Assuming that species track the modeled environmental conditions, this method allows ecologists to draw conclusions on how species' distributions will shift in the future (Elith et al. 2011, Hughes et al. 2012. Given the accessibility of global species occurrence records (gbif.org 2018), high resolution climate data (WorldClim et al. 2017, CHELSA, Karger et al. 2017, and user-friendly software (e.g. MaxEnt software packages, Phillips et al. 2006), SDMs are widely utilized to predict species' range shifts across the globe (Merow et al. 2013, Pacifici et al. 2015.
Despite their ubiquitous use, however, there has been increasing criticism of SDMs regarding their over-simplification of the factors that limit species distributions Peterson 2012, Early andSax 2014). SDMs that use predominately climatic factors to predict a species' distribution make the key assumptions that 1) the species' range is in equilibrium with its climatic niche (Veloz et al. 2012), 2) climate is indeed the main distribution driver (Araújo and Peterson 2012), 3) the climate niche is static over timescales relevant to predictions and 4) all populations respond identically to climate, such that the climate niche for the species is also that for individual populations (Wiens et al. 2009). Even given the long history of work that shows strong evidence for local adaptation to climate conditions in many plants and animals (Mayr 1956, Aitken et al. 2008, Pelini et al. 2009, Fournier-Level et al. 2011, Ruegg et al. 2018, it is poorly understood how differences in local population responses to climate may affect SDM results (but see Hällfors et al. 2016, Schwalm et al. 2016, Theodoridis et al. 2018) and thus how important this last assumption may be. More research on this is especially needed, as recent work has shown that predictions of range shifts using species-wide SDMs underestimate intraspecific genetic diversity loss (Balint et al. 2011, Alsos et al. 2012. The most common approach to constructing climate-based SDMs is to use all available data on a species' occurrences to predict its distribution, implicitly assuming the same climatic responses across populations Peterson 2012, Merow et al. 2013). Out of the thousands of SDM studies published, we could find only 30 previous studies (Supplementary material Appendix 1 Table A1 for search methods and results) that account for intraspecific differences in climate responses by separately modeling smaller units (henceforth, 'subpopulations') of a species range, which generally grouped according to presumed genetic differences, differing climate histories, or geographic regions (Table 1). Combining these intraspecific-level SDMs yields predictions of P occ over the same geographical extent as a species-level SDM, while predicting P occ for subpopulations according to their corresponding climate distributions. While broadly defined subpopulations may not, in fact, capture intraspecific variation in climate responses, testing for such effects by constructing intraspecific-level SDMs could be an important improvement, allowing for potential differences in climatic response (Hällfors et al. 2016). However, only three studies have examined multiple approaches to defining subpopulations (Table 1, Supplementary material Appendix 1 Table A1) and thus it is unclear to what extent these differences might influence SDM predictions. For the 30 studies we found that have used intraspecific-level SDMs, we summarize in columns (from left to right) 1) approaches to identifying subpopulations (Genetic = genetic groups and/or taxonomic lineages, Habitat = climatic and/or geographic groups, Other = discrete phenotypes), 2) approaches to validating the intraspecific-level SDMs (Cal. data = validation using the calibration dataset, Ind. data = validation using independent data on distribution, Niche div. = tests for niche divergence controlling for background environments), 3) methods used to compare alternative SDMs (Global = species-level model pooling all occurrences, Mult. subpop. = study includes multiple intraspecific-level SDMs, based on alternative groupings of occurrences, Subpop. > global = for studies that compared species-level and intraspecific-level SDMs, how many found that the intraspecific-level model performed best?) and 4) interpretation of mechanisms for intraspecific-level SDM differences (LA = demonstrating genetically-based local adaptation to climate, Sampling = reflecting differences in sampling intensity among subpopulations, Habitat avail. = reflecting differences in habitat availability among subpopulations). Note that we focus on studies using traditional SDM approaches calibrated with occurrence data within a species' native range and do not include related approaches that incorporate intraspecific structure through modeling ecosystem types or data from transplant experiments (Benito-Garzón et al. 2011, Gray et al. 2011, Hamann and Aitken 2013 Improvements in predictions between intraspecific-level SDMs and the corresponding species-level SDM have generally been interpreted as indicating local adaptation to climate (Table 1, Supplementary material Appendix 1 Table A1) or more broadly, differences in climate responses, be they adaptive or not. However, even in the absence of any genetically-based niche divergence, intraspecific-level SDMs could produce different predictions simply due to over-fitting, better representation of under-sampled climates, or environmental differences (including biotic interactions) among the defined subpopulations leading to different inferred climate responses. Regardless of the mechanisms involved, the extent and importance of differences between intraspecific-level and species-level SDM predictions have only been considered in very few of the intraspecific-level SDM studies we identified ( While SDMs only formally predict P occ , their outputs have been assumed to correlate to population performance, such as population persistence (Araújo and Williams 2000), functional traits (Thuiller et al. 2009), and abundance (Weber et al. 2017). However, predicting species performance with SDM output is controversial and recent studies have disagreed over the extent to which SDM output can accurately predict aspects of population performance. While some researchers have demonstrated strong links between SDM predictions and abundance (VanDerWal et al. 2009, Van Couwenberghe et al. 2013, Lee-Yaw et al. 2016, recent metaanalyses show that this relationship is stronger in vertebrates than in plants (Weber et al. 2017) or even that this relationship hardly exists at all (Dallas and Hastings 2018). Other studies further question the link to demographic rates (Thuiller et al. 2014, Csergő et al. 2017. Even the existing evidence for using distance to environmental, not geographic, centers to predict population performance (Martínez-Meyer et al. 2012) and genetic diversity (Lira-Noriega and Manthey 2014) has been recently contested (Pironon et al. 2017, Santini et al. 2018) as a possible oversimplification of the biogeographic drivers on populations (Dallas et al. 2017). Thus, the extent to which SDM output can be used to infer the distribution of population performance or other traits still needs closer examination. This is especially important for species interactions, which can be influential drivers of species range limits (Early and Keith 2018). Given that positive species interactions, such as facilitation and mutualism, can dramatically broaden species' ranges (Afkhami et al. 2014), it is particularly relevant to understand if they can be described by SDM output.
In this study, we test the ability of SDMs to 1) predict a species' P occ , contrasting global versus intraspecific-level SDMs, and 2) predict local population performance as well as positive species interactions. To address the first question, we examine how intraspecific-level and species-level SDMs differ for a circumboreal alpine-arctic plant, using broad genetic and habitat (i.e. biome) differences to construct intraspecific-level SDMs. While there is likely finer-scale differentiation within our broadly defined subpopulations, grouping occurrence data by broad patterns of genetic differentiation (Pearman et al. 2010, D'Amen et al. 2013, Serra-Varela et al. 2017) and climate or habitat differences (Sork et al. 2010, Hällfors et al. 2016, Hu et al. 2017) is a common approach for intraspecific-level SDMs. Despite the fact that genetic versus habitat-based approaches to defining subpopulations can potentially yield differing predictions, no previous study has compared these approaches (but see e.g. Marcer et al. 2016 for comparison of genetic and trait-based models). To address the second question, we test the predictive value of SDM-derived P occ for other aspects of population performance and ecological function.
We test these questions for the facilitative arctic-alpine cushion plant Silene acaulis (Caryophyllaceae). Silene acaulis is a long-lived gynodioecious perennial, and its cushion-like growth form and deep taproot are thought to be adaptations to harsh arctic-alpine conditions (Griggs 1956, Billings 1974. Individual cushions slowly grow radially outwards and are known to live 300 yr or longer (Morris and Doak 1998). Individual performance can be measured by cushion size, as larger cushions 1) produce disproportionally more fruits than smaller ones (Chardon et al. 2019a) and 2) grow faster, survive better, or both (Morris and Doak 1998). Silene acaulis is an ideal species for this work, as there is evidence for local adaptation to climate (Peterson et al. 2018) as well as genetic structure (Gussarova et al. 2015), and global trait data are available Morris 2010, Cavieres et al. 2013) to comprehensively analyze whether SDM output can predict traits. Its wide distribution (Hultén and Fries 1986) across the Northern hemisphere makes it optimal for SDMs (Pacifici et al. 2015) and particularly appropriate for this study, as greater intraspecific variation may exist relative to more narrowly distributed species.
Silene acaulis' cushion growth form makes it an important facilitator of other arctic-alpine species (Cavieres et al. 2016). The facilitative effects of cushion plants generally increase along elevational gradients, as they provide the necessary microhabitat for beneficiary species to establish at high elevations characterized by increased abiotic stress (Callaway et al. 2002). These facilitative interactions, however, can break down at extremely high levels of abiotic stress (Michalet et al. 2006, reviewed in Liancourt et al. 2017. Facilitative plants are not only important in structuring plant communities around the globe (Cavieres et al. 2016), but can also buffer responses to rapid climatic changes in alpine and arctic regions (Anthelme et al. 2014). As S. acaulis is projected to lose over half of its climatically suitable habitat in the next 20 yr (Ferrarini et al. 2019), it is particularly critical to understand how this important facilitative arctic-alpine species will respond to forecasted climatic changes.
We hypothesize that the genetic-and habitat-based intraspecific-level SDMs will provide more accurate distribution predictions than the species-level SDM and that the two intraspecific-level models will yield very similar results to each other (hypothesis 1). This would suggest that broad scale genetic and habitat differences capture at least some variation in local climate responses, which would be reflected in the intraspecific-level SDMs. We also expect that any differences seen in model predictions will not be explained simply by differences in climate across subpopulations. Second, we hypothesize that if SDM P occ captures the potential for high population performance (Araújo and Williams 2000), S. acaulis individual plant sizes will be larger in areas of higher P occ (hypothesis 2), as larger plants grow faster, survive better, or both (Morris and Doak 1998). Third, as facilitative interactions tend to be higher in climatically stressful areas (Callaway et al. 2002), we expect that high facilitative interaction strength will correspond with low predicted P occ for S. acaulis (hypothesis 3).

Climate data
We used four bioclimatic variables from the CHELSA dataset in the timeframe 1979 to 2013 (during which most of our species data is available) and at a 30 arc-sec (~1 km 2 ) resolution (Karger et al. 2017). These four variables were recently used in a S. acaulis SDM study (Pironon et al. 2015) and have been shown to be particularly important predictors in SDMs (Bradie and Leung 2017): maximum temperature of the warmest month, temperature seasonality (i.e. difference between annual mean minimum and maximum), precipitation of the wettest month, and precipitation seasonality (Supplementary material Appendix 1 Fig. A1). Using shapefile boundaries of North America, Europe and Russia (thematicmapping.org 2018), we cropped these bioclimatic variables to encompass the broad geographic regions that define S. acaulis' global distribution. To account for the distinct climates over the large land-locked bodies of water found within the species' range (e.g. Canada), we also removed the climate data of large lakes (≥ 50 km 2 ) and reservoirs (≥ 0.5 km 3 ; WWF Global Lakes and Wetlands Database).

Species occurrences
We combined geographic occurrences from two existing data sets on S. acaulis traits (see 'Species traits' below; Morris 2010, Cavieres et al. 2013), occurrences from a S. acaulis genetic study (Supplementary material Appendix 1 in Gussarova et al. 2015), and S. acaulis occurrence records from digital databases. We downloaded all 'Silene acaulis' (and listed subspecies) digital occurrence records from the databases BIEN (biendata.org 2018), GBIF (gbif.org 2018), and BioTIME (BioTIME Consortium 2018, Dornelas et al. 2018). To match the resolution and timeframe of the occurrence data to the bioclimatic data, we performed several operations. First, we filtered all data at 1 km geographic position accuracy or better and at 1979 data collection year or later, where these metadata were available. Second, we removed any exact latitude and longitude duplicate occurrences. Third, to reduce erroneous occurrences, we filtered all data to retain only biomes that contain alpine or tundra terrain within S. acaulis' geographic distribution ('Tundra', 'Temperate Conifer Forests', 'Temperate Broadleaf and Mixed Forests', 'Boreal Forests/Taiga'; Ecoregions 2017). As this filter deleted the only occurrence record from eastern Russia, we added it back in because of its rare verification of existence in this region (Gussarova et al. 2015) and our use of this record to determine geographic delineations of genetic groups. Fourth, we manually checked isolated southern or lower elevation occurrences in the USA and mainland Europe (GoogleEarth Pro 2009) and removed six occurrences in terrain where S. acaulis does not naturally grow (i.e. in forests). We recognize that this occurrence dataset does not fully represent the range of S. acaulis, as is the case in many occurrence records (Meyer et al. 2016). In particular, occurrences in the Canadian and Russian Arctic range of the species are markedly sparse, however both arctic and alpine regions are well represented in the occurrence data ( Fig. 1).

Species traits
We obtained S. acaulis trait data from a global cushion plant study on facilitative interactions (Butterfield et al. 2013, Cavieres et al. 2013, Lortie 2018) and a long-term demographic study Morris 2010, Peterson et al. 2018, Doak et al. unpubl.). These data span S. acaulis' geographic distribution, with a total of 50 sites over 8 countries in both North America and Europe (Fig. 1a, Supplementary material Appendix 1 Table A2). These data include individual plant cushion size (measured as elliptical area; n = 5890 plants), a good plant performance indicator because larger plants 1) grow faster, survive longer, or both (Morris and Doak 1998) and 2) produce disproportionally more fruits (Chardon et al. 2019a). Although these data were obtained over the span of multiple years, it has been shown across a range of sites and years that there is very little size variation between years due to S. acaulis' slow growth rate (Morris and Doak 1998; Supplementary material Appendix 1 Fig. A2). The data also include percent cover and richness of beneficiary species growing within individual cushion plants (n = 1674 plants). The strength of plant-plant facilitative interactions is commonly measured as beneficiary species percent cover and richness (Cavieres et al. 2016), from which we calculated a beneficiary species Shannon diversity index for each S. acaulis individual (vegan package; Oksanen et al. 2018). While cushion plant size can influence facilitative interactions (Chardon et al. 2018), our data show only moderate correlation (correlation of size and beneficiary species percent cover = 0.46; size and richness = 0.39; size and diversity = 0.33).
We subset the plant size data to 1) account for methodological differences between the two datasets and 2) focus on larger plants in order to best capture variation in population performance. As plant size data was recorded through either targeted sampling of larger individuals (Cavieres et al. 2013) or sampling all individuals in a population (Doak and Morris 2010; for comparison see Supplementary material Appendix 1 Fig. A8 in Chardon et al. 2018), we first retained only cushion sizes of plants above the 65th percentile overall from the latter dataset, a cutoff that best aligned the two size distributions (Supplementary material Appendix 1 Fig. A3a-b). As we specifically aimed to test if SDM output can predict cushion plant sizes (see 'Model performance' below), we then subset all data to only include plant sizes above the 40th percentile (Supplementary material Appendix 1 Fig. A3c). This provided the best fit out of a set of cutoffs tested, and using other cutoffs does not change the qualitative patterns in the results. Larger plant sizes are more meaningful population performance indicators than the full plant size distribution, as larger plants correlate better with environmental variables and produce disproportionately more fruits Doak 1998, Chardon et al. 2019a).

Species distribution models
To correct for some of the sampling bias present in the occurrence records, which are far denser in Europe than in either North America or Russia, we subsampled all records by keeping only one occurrence per 30 arc-sec cell ('gridSample' function in dismo package;  to match the resolution of the bioclimatic data (total n = 4107 occurrences; Fig. 1a). Although this does not correct for unsampled areas, it is a standard bias correction approach (Fourcade et al. 2014, Guisan et al. 2017. We then split the occurrences into four genetic groups identified by STRUCTURE analyses of multilocus AFLP markers (335 markers for 106 populations) by Gussarova et al. (Fig. 1b;corresponding to Fig. 4 in Gussarova et al. 2015). We also split the occurrences into six habitat biome groups in the Nearctic or Palearctic realms ( Fig. 1c; shapefile boundaries from Ecoregions 2017) representing broad habitat and climatic differences. While there is considerable correspondence between the habitat and genetic groupings, they are not identical, and also differ in the total number of subpopulations recognized, likely reflecting the fact that genetic groups capture post glacial history as well as current habitat effects. We likely do not capture the full extent of genetic variation within these defined subpopulations, and emphasize that in this study we aim to assess the potential for broad intraspecific differences to influence SDM predictions. We used Maximum Entropy Species Distribution Modeling (MaxEnt ver. 3.4.1; Phillips et al. 2018) to model S. acaulis' current distribution using 1) all occurrences together (species-level SDM), and separately for 2) occurrences within each genetic group (genetic intraspecific-level SDMs) and 3) occurrences within each habitat group (habitat intraspecific-level SDMs). We calibrated and projected individual SDMs only in the polygon corresponding to that subpopulation. We chose MaxEnt to create our SDMs, as it is a common and well-performing algorithm for presenceonly data (Phillips et al. 2006, Elith et al. 2010, Merow et al. 2013. We elected to only model current distribution (following Hällfors et al. 2016), as this allows the most appropriate evaluation of which SDM type (species-level, genetic intraspecific-level, or habitat intraspecific-level) can best predict S. acaulis' recorded geographic distribution and population performance.
We employed 10-fold cross-validation MaxEnt runs for each individual SDM with a jackknife test of variable importance and response curves for environmental variables. To binarize the resulting P occ of each run, we selected the maximum test sensitivity plus specificity threshold in MaxEnt, a commonly used and well-performing suitability threshold that maximizes the sum of sensitivity and specificity (Liu et al. 2005). We used these thresholds to create 'presence-background' maps, which show cells as either present (i.e. P occ above threshold) or as background (i.e. P occ below threshold), for each SDM type. This allowed us to compare the predicted binary 'presence-background' maps among SDM types. We constructed these maps retaining only those cells above the threshold in more than five of the replicates per individual SDM (following Hällfors et al. 2016). We then mosaicked maps across subpopulations to generate the final presence-background and also P occ (as indicated by cloglog output; Phillips et al. 2017) maps across the entire species distribution area for species-level, genetic intraspecific-level, and habitat intraspecific-level models. Given that output values are only relative to the modeled region and are dependent on occurrence density and sampling design, we recognize that comparing these values across models can be difficult (Merow et al. 2013).

Model performance
To test which SDM best predicts S. acaulis' recorded distribution (hypothesis 1), we evaluated the predicted binary presence-background maps with two types of validation data. First, we employed the standard approach of using the recorded S. acaulis occurrences used to calibrate the SDMs to calculate sensitivity (proportion of correctly identified presences) for the species-level, genetic intraspecific-level, and habitat intraspecific-level SDMs. However, since the occurrence data available is sparse in areas where S. acaulis commonly occurs, such as central Alaska, the Canadian tundra and Russia, we also compared model performance to an independently-derived distribution map. We used an existing global S. acaulis distribution map (digitized terrestrial locations from map 791 in Hultén and Fries (1986) to calculate standard performance metrics of sensitivity, specificity (proportion of correctly identified background points), and True Skills Statistic (TSS = sensitivity + specificity − 1; Allouche et al. 2006) for each of the three SDM types. TSS is particularly useful in comparing model accuracy (Allouche et al. 2006, Shabani et al. 2016, whereas the commonly employed area under the receiver operating curve (AUC) has been increasingly criticized (Lobo et al. 2008, Jiménez-Valverde 2012, Shabani et al. 2016. Using an existing distribution map (Hultén and Fries 1986) allowed us to validate our models with data independent from those used to calibrate our SDMs, a validation approach that has yet to be employed in intraspecific-level SDM studies (Table 1; Peterson et al. 2019). Due to the longlived nature and slow growth rates of S. acaulis, this 30-yr old distribution map still reflects the habitats and climates relevant for the species. Furthermore, this map has been shown to be useful in other SDM work on plants in the Northern Hemisphere (Alsos et al. 2012). We refined this large scale and low precision map to include only biomes where S. acaulis is most likely to be found, thereby employing the same criterion we used to filter S. acaulis occurrences (see 'Species occurrences' above). This kind of filter has been shown to greatly increase the accuracy of where a species is likely to be found, illustrating its applicability to improve species distribution maps (Ocampo-Peñuela et al. 2016).
We dissected differences between SDM predictions in three steps. First, we examined P occ correlations between the species-level and intraspecific-level SDMs for each distinct subpopulation to see where in the species' range SDM type influenced predictions. Second, to examine if climate differences between subpopulations cause SDM dissimilarities, we compared how predicted 1) individual subpopulation climate spaces (i.e. climate in cells predicted as present) and 2) subpopulation regional climate conditions (i.e. all cells from occurrence and background points) differ between the three SDM types. We focused on the two climate variables identified as most important by MaxEnt's analysis of variable contribution and jackknife test of variable importance. Third, we calculated Warren's I (function 'modOverlap' in package fuzzySim; Barbosa 2015), a statistic based on Schoener's D and Hellinger distance (Warren et al. 2008), to quantify niche similarity. We computed this statistic to compare the climate niches predicted by the species-level SDM with the intraspecific-level SDMs for 1) the species entire range and 2) individual genetic and habitat regions. As the habitat intraspecific-level SDM projects to fewer cells than the other two SDM types (Supplementary material Appendix 1 Fig. A5), we used only those cells when computing Warren's I between the habitat SDM and the other two SDM types.
To test if P occ can predict S. acaulis population performance (hypothesis 2), we fit linear mixed models (LMMs with function 'lmer' in package lme4; Bates et al. 2015) on cushion plant size using linear and quadratic P occ from each SDM type as predictor variables and a random effect of site. We calculated additional model summary outputs with the packages lmerTest (Kuznetsova et al. 2017) and MuMIn (Bartón 2018). We log-transformed size to meet LMM assumptions of residual distribution.
In order to test if facilitative interactions between beneficiary species and S. acaulis can be predicted by P occ (hypothesis 3), we fit LMMs on beneficiary species percent cover, richness and diversity using linear and quadratic P occ from each SDM type as predictor variables and a random effect of site. We log-transformed percent cover [log(cover + 1)] to meet LMM assumptions of residual distribution, a transformation that was not necessary for richness or diversity.

Results
We found that the genetic intraspecific-level SDM predicts the highest proportion of true presences (sensitivity, or cell overlap with recorded distribution) when compared against both the S. acaulis occurrences used to construct the models and the Hultén and Fries (1986) distribution map (hypothesis 1; Supplementary material Appendix 1 Table  A3a). The genetic SDM also has the highest TSS value relative to the habitat intraspecific-level and the species-level SDMs, with the species-level SDM performing worst (Supplementary material Appendix 1 Table A3a). The three SDM types yield quite different presence-background (Fig. 2) as well as relative probability of occurrence (P occ ) predictions (Supplementary material Appendix 1 Fig. A4-A5). In contrast, the proportion of predicted background points (specificity) was similarly well predicted by all models (Supplementary material Appendix 1 Table A3a). Niche similarity is also high between SDM types, with Warren's I ranging between 0.83 and 0.92 (Supplementary material Appendix 1 Fig. A4).
Examining prediction differences by genetic subpopulation illustrates that the genetic intraspecific-level SDM outperforms the species-level SDM in all subpopulations except the Atlantic (Supplementary material Appendix 1  Table A3b). For the European and American genetic groups, the largest mismatches occur where the species-level SDM predicts high P occ while the genetic SDM predicts low P occ (Supplementary material Appendix 1 Fig. A6). Predicted P occ between the species-level and habitat intraspecific-level SDMs are generally more similar in the Palearctic realm (Europe; Supplementary material Appendix 1 Fig. A7). The species-level SDM generally overpredicts P occ in the Palearctic realm and underpredicts in the Nearctic realm (North America), broadly corresponding to better performance by the species-level SDM in the Palearctic realm and better performance by the habitat SDM in the Nearctic realm (Supplementary material Appendix 1 Table A3c).
Out of the four climate variables we used to construct our SDMs, maximum temperature of the warmest month (average percent variable contribution to MaxEnt models for genetic SDM: 55%; habitat SDM: 59%) and temperature seasonality (34%; 27%) are the two most important environmental variables across the four and six separate genetically-based and habitat-based SDMs, respectively. These two variables are also most important for the species-level SDM (temperature: 33%; temperature seasonality: 67%). Jackknife tests of variable importance in both training and testing gains for each separate SDM also identified the variable with the highest contribution as being most important in all but the Nearctic Tundra biome.
Similarities in P occ predictions between the species-level SDM and each of the intraspecific-level SDMs correspond to similarities in predicted climate niches corresponding to the predicted presences (as defined by the maximum test sensitivity plus specificity threshold in MaxEnt). The climate space for S. acaulis in the Atlantic genetic group, where P occ predictions between the genetic intraspecific-level and global SDM are most similar, shows the largest predicted climate niche overlap compared to the other genetic groups (Fig. 3). This is supported by Warren's I for the Atlantic genetic group, which is 0.90 between the two SDM types. The Beringian genetic group also shows high niche similarity (Warren's I = 0.95), whereas the European (Warren's I = 0.76) and American (Warren's I = 0.75) show lower niche similarity between SDM types. When compared to the species-level SDM, the habitat SDM also shows the largest difference in predicted climate spaces where similarity between P occ predictions is low, most notably in the Nearctic Tundra and Conifer Forests biomes (Fig. 4). Niche similarity between SDM types, on the other hand, is higher in the Nearctic realm (Tundra = 0.92, Mixed Forests = 0.92, Conifer Forests = 0.94) than in the Palearctic ream (Tundra = 0.85, Boreal Forests = 0.81, Mixed Forests = 0.83). Furthermore, the range of available climate conditions used to construct each of the intraspecific-level SDMs overlap (Supplementary material Appendix 1 Fig. A8a, c) and the climate spaces of predicted presences is narrower and more overlapping (Supplementary material Appendix 1 Fig. A8b, d). This illustrates that differences between the intraspecific-level and species-level SDMs may not be due primarily to distinct climates among the genetic or habitat subpopulations.
SDM output only poorly predicts S. acaulis performance (hypothesis 2) and strength of facilitative interactions (hypothesis 3), and this very weak relationship (marginal R 2 : 0.01-0.15) is not improved by fitting intraspecific-level SDM P occ values. P occ values from the species-level SDM best predict both S. acaulis cushion plant size and beneficiary species percent cover, with a peak at median to high P occ (Supplementary  material Appendix 1 Table A4, Fig. 5a-b). The other measures of facilitative interaction strength, beneficiary species richness and diversity, cannot be significantly predicted by SDM output (Supplementary material Appendix 1 Table A4, Fig. 5c-d).

Discussion
We critically evaluated the performance of three approaches to model species distributions: a traditional species-level SDM using a species-wide climate niche, and intraspecific-level models based on either genetic groups or climatically-distinct habitat types. We found that the intraspecific-level SDMs decisively outperformed the species-level SDM in predicting the distribution of S. acaulis, consistent with results from the previous few studies that have made this comparison (Table 1, Supplementary material Appendix 1 Table A1). In particular, sensitivity is conclusively higher in our intraspecific-level SDMs compared to the species-level models, as also recently found by Lecocq et al. (2019). While to date few studies have included intraspecific differences in SDM models (Table 1, Supplementary material Appendix 1 Table  A1; see also Schurr et al. 2012, Ehrlén and Morris 2015, Pironon et al. 2018, the improved accuracy of intraspecificlevel SDM distribution predictions illustrates this as a promising approach. As we found that niche similarity is high between SDM types, we emphasize that multiple evaluation metrics, as well as a close examination of predicted climate spaces, are needed to assess SDM performance and prediction differences. Using a detailed global trait dataset for the species, we also found that SDM output poorly predict S. acaulis cushion size, a measure of population performance, and facilitative interaction strength. Support for this prediction is not improved with intraspecific-level SDMs. We observed more variability in cushion size with increasing P occ , which could explain the lack of a clear relationship. Such a pattern between population performance and habitat suitability has, in fact, been described in previous work (Hengeveld 1990, Brown 1995). The traits we tested loosely follow an elevational pattern, with largest plant size and strongest facilitative interactions found at mid-elevations (Supplementary material Appendix 1 Fig.   A9a-b), whereas P occ increases with elevation (Supplementary material Appendix 1 Fig. A9c). However, model fit is neither improved by adding elevation as a fixed effect to our LMMs nor by substituting climate variables for P occ as a predictor variable (Supplementary material Appendix 1 Table A5, Fig.  A10). While the ability to predict species' traits with P occ has been examined before (Thuiller et al. 2009), recent work has shown that such results need to be interpreted cautiously, especially when considering species abundances (Dallas andHastings 2018, Santini et al. 2018) and demographic rates (Thuiller et al. 2014, Csergő et al. 2017, Pironon et al. 2018. Although biotic interactions have been successfully modeled on a geographic scale (Araújo and Rozenfeld 2014) and SDM predictions can improve when incorporating facilitative interactions (Filazzola et al. 2018), our results indicate that predicting biotic interactions from SDM output values needs to be approached with caution.
We did not find that model performance simply increased with greater subdivisions of the data, as the genetic SDM with four groups outperformed the habitat SDM based on six groups and calibrated on an overall smaller range (Fig. 2). This suggests that, at least for S. acaulis, subpopulations based Figure 3. Species-level and genetic intraspecific-level SDMs predict different climate niches. (a-d) Differences in the predicted climate niche between species-level (shown in grey) and the genetic (shown in color) SDM types varies by genetic group. Only 1% of data per genetic group are shown to improve clarity, and only the predicted presence (determined by equal testing sensitivity plus specificity threshold) cells are plotted. 69 on habitat types or geographic regions may not best capture intraspecific differences in responses to climate. Our results imply that how a species is divided into subpopulations is critical to SDM inference and accuracy. Given that only three previous intraspecific-level SDM studies have compared multiple approaches to delineating subpopulations (Marcer et al. 2016), this needs to be examined in greater detail. It is especially surprising that none of these studies used independent validation data with which to judge predictive quality. We show that this can readily be done with a distribution map independent of the digital occurrences used to calibrate our models. This evaluation approach provided results in agreement with the traditional approach of using the species' occurrences used to calibrate the model. While model sensitivity is lower when evaluating model performance against a large scale and low precision map (Supplementary material  Appendix 1 Table A3), this approach ranks the models' performance in the same order. Further, even coarse distribution maps may better capture portions of a species' range that, for whatever reason, may be underrepresented in occurrence datasets, particularly for widespread and common species. Indeed, this is what we see for S. acaulis, with limited occurrence data available in large portions of the species' range (e.g. central Alaska, Canada, Russia). We therefore discourage the common practice of only validating SDM performance with the occurrence data used to construct the models in the first place, especially when that occurrence data does not represent the full extent of the species' range.
Notably, while 73% of past intraspecific-level SDM studies attribute increased intraspecific-level SDM performance to local adaptation, most do not report results that allow for an assessment of the importance of these differences or their likely causes (Table 1, Supplementary material Appendix 1  Table A1), even though several other mechanisms could also cause these improvements. First, if there are strong climate differences in separately modeled regions, intraspecific-level SDMs may fit correspondingly different climate spaces (Meynard et al. 2017). In our study, subpopulations' predicted climate niches were substantially narrower and more overlapping than their regional climate spaces (Supplementary material Appendix 1 Fig. A8), suggesting that differences in P occ are not simply due to sharp distinctions between regional climates. Second, differences in sample size (i.e. recorded occurrences) between different regions may mean that a species-level model may perform poorly for subpopulations with lower sampling intensity due to swamping of the model fit by data from better sampled regions (Pearman et al. 2010, Hällfors et al. 2016). Our study partially supports this explanation, as we see that unequal sampling intensity across regions ( Fig. 1) corresponds to differences between species-level and habitat, but not genetic, intraspecific-level predictions. Given that the majority of species likely have biased occurrence data (Meyer et al. 2016), intraspecific-level SDMs may be useful as a way to control for bias in model fits, even when there are not local differences in climate suitability.
We emphasize that SDMs themselves are not capable of fully dissecting these different mechanisms, but examination of the calibration data and model predictions can help suggest their possible importance. In particular, dissimilarity in SDM predictions can only indicate the potential for local adaptation and resulting population-level climate response, which would need to be confirmed with direct experimental work. The next steps are to then explicitly incorporate local adaptation into predictions of range shifts with climate change (for review see Peterson et al. 2019) according to well identified subpopulations, or by adopting a mechanistic approach Figure 5. Species traits are poorly predicted by and show no trend with SDM relative probability of occurrence (P occ ). (a) Species-level SDM P occ values are the only ones that significantly predict Silene acaulis cushion area, with a peak at higher P occ . There is no significant relationship between SDM output values of any SDM type and beneficiary species percent cover (b), richness (c), or diversity (d). Colors are all as in (a), and all LMMs are fit with a quadratic P occ term. Shown are the log-transformed data used to fit the models in (a) and (b). Note that sample size for the trait data in (a) is larger, and therefore a wider range of P occ values are observed. (Angert et al. 2011). In the case of S. acaulis, climate manipulation experiments have found local adaptation to temperature between populations corresponding to the Beringian and American genetic groups (Peterson et al. 2018). Silene acaulis might also respond strongly to other climatic drivers than the ones we examined, although recent studies have identified temperature to be an important climate variable for the species (Pironon et al. 2015, Ferrarini et al. 2018.
Given that we found large inconsistencies between SDM types, we emphasize that, when possible, subpopulations should be modeled separately for more accurate predictions and that the choice of how to define subpopulations needs to be well-justified. Our results illustrate the necessity of examining potential intraspecific variation in responses to climate, which, if present, violates a foundational assumption of SDMs built using a species' full climatic niche. Species traits or performance can differ with the various local climate conditions found within its range (Emery et al. 2015, Amburgey et al. 2018 and predictions for locally modeled populations often do not match those from specieslevel models (Hällfors et al. 2016, Schwalm et al. 2016. Practitioners using SDM outputs for conservation planning should be particularly wary of predictions generated from single SDMs using large scale distribution data, and aim to compare outputs from multiple SDM types. Particular care should also be taken with relating SDM output to species traits, as this relationship does not hold in many systems and, at least for S. acaulis, is not improved with better performing intraspecific-level SDMs.