Area Not Geographic Isolation Mediates Biodiversity Responses of Alpine Refugia to Climate Change

Climate refugia, where local populations of species can persist through periods of unfavorable regional climate, play a key role in the maintenance of regional biodiversity during times of environmental change. However, the ability of refugia to buffer biodiversity change may be mediated by the landscape context of refugial habitats. Here, we examined how plant communities restricted to refugial sky islands of alpine tundra in the Colorado Rockies are changing in response to rapid climate change in the region (increased temperature, declining snowpack, and earlier snow melt-out) and if these biodiversity changes are mediated by the area or geographic isolation of the sky island. We resampled plant communities in 153 plots at seven sky islands distributed across the Colorado Rockies at two time points separated by 12 years (2007/2008–2019/2020) and found changes in taxonomic, phylogenetic, and functional diversity over time. Specifically, we found an increase in species richness, a trend toward increased phylogenetic diversity, a shift toward leaf traits associated with the stress-tolerant end of leaf economics spectrum (e.g., lower specific leaf area, higher leaf dry matter content), and a decrease in the functional dispersion of specific leaf area. Importantly, these changes were partially mediated by refugial area but not by geographic isolation, suggesting that dispersal from nearby areas of tundra does not play a strong role in mediating these changes, while site characteristics associated with a larger area (e.g., environmental heterogeneity, larger community size) may be relatively more important. Taken together, these results suggest that considering the landscape context (area and geographic isolation) of refugia may be critical for prioritizing the conservation of specific refugial sites that provide the most conservation value.


INTRODUCTION
Rapid anthropogenic climate change is altering temperature and precipitation patterns around the world, threatening biodiversity (IPCC, 2014;Díaz et al., 2019). While some changes in biodiversity patterns are inevitable, the impacts of climate change on biodiversity may be mitigated in part via the conservation of refugia: geographic locations or habitat types that are buffered from the most intense and immediate effects of environmental change (Tzedakis et al., 2002;Keppel et al., 2012).
Refugia have played a key role in maintaining biodiversity during past (quaternary) climatic changes (Taberlet and Cheddadi, 2002;Byrne, 2008;Carnaval et al., 2009). For example, the climatic stability, rugged terrain, and complexity of soils and microclimates of the Klamath-Siskiyou region of southwest Oregon allowed the region to serve as a refugium for ecological communities that required cool and moist conditions in the past (Whittaker, 1960;Olson et al., 2012;Copeland and Harrison, 2015). Currently, refugia have the potential to maintain biodiversity in the face of anthropogenic climate change (Ashcroft, 2010;Dobrowski, 2011;Ashcroft et al., 2012;Keppel et al., 2012Keppel et al., , 2018. While maintenance of taxonomic diversity has been the traditional focus of refugial conservation efforts, the role of refugia in preserving functional and phylogenetic diversity may be more important as functional traits and phylogenetic identity are thought to more directly control ecosystem function and stability (Díaz and Cabido, 2001;Cavender-Bares et al., 2009;Funk et al., 2017). Despite the recognized importance of climate change refugia in general, there is still a need to better identify specific sites that may serve as high-quality refugia given the economic and political constraints on conserving refugial habitats (Keppel et al., 2012(Keppel et al., , 2015Mokany et al., 2017).
Importantly, refugia may largely differ in their capacity to buffer against environmental changes and therefore promote the persistence of threatened species (Keppel et al., 2015). Within geographic locations or habitat types that are broadly considered refugia, individual sites may vary in their ability to withstand the impacts of climate change based on their landscape context (Gaston and Blackburn, 1996;Ashcroft et al., 2009;Keppel et al., 2015). Specifically, island biogeography theory states that site area and geographic isolation mediate rates of dispersal and ecological drift (MacArthur and Wilson, 1967), with subsequent impacts on taxonomic, functional, and phylogenetic diversity (Graham and Fine, 2008;Spasojevic et al., 2014;Fernández-Palacios et al., 2015;Carmona et al., 2016). Differences in area and isolation among refugial sites may also alter how those sites respond to shifts in the niche selection regime imposed by climate change (Keppel et al., 2015). For example, larger refugial sites are more likely to contain higher levels of environmental heterogeneity and thus may be able to provide optimal environmental conditions for a greater diversity of species (Currie et al., 2004), functional strategies (Spasojevic et al., 2014), and/or phylogenetic lineages (Leibold et al., 2010) even with changes in regional climate (Keppel et al., 2015). Similarly, larger sites may also be able to support larger populations or communities that are less susceptible to demographic stochasticity (Loreau and de Mazancourt, 2008). On the other hand, less isolated refugial sites may be more likely to be connected to other refugia via dispersal and thus may be able to retain species (Leibold et al., 2004), functional strategies (Schleuter et al., 2012), or phylogenetic lineages (Thorpe et al., 2008;Eldridge et al., 2018) negatively affected by changes in regional climate through source-sink dynamics. Understanding how the area and isolation of refugia mediate changes in taxonomic, functional, and phylogenetic diversity caused by climate change is therefore a key next step in prioritizing the conservation of specific refugial sites that provide the most conservation value.
Alpine tundra is an excellent system for understanding how site area and isolation may affect biodiversity change in climatic refugia. Found at the tops of mountains, alpine tundra is characterized by meadows of wildflowers and grasses (alpine by definition has no trees) and harsh environmental conditions (Bowman and Seastedt, 2001;Korner, 2003Korner, , 2004. Alpine tundra in western North America currently displays an archipelago-like distribution, with "sky islands" of suitable high-elevation habitat surrounded by a matrix of unsuitable low-elevation habitat (Pewe, 1983). These disjunct sky islands formed when the large continuous alpine communities present during the last glacial maximum began to fragment and retreat upslope in response to natural climatic warming (Pewe, 1983). Thus, alpine sky islands serve as an ex situ refugia for a formerly widely distributed paleo-habitat type, as evidenced from molecular data on a wide range of organisms (e.g., Skrede et al., 2006;Shafer et al., 2011). The upslope movement and subsequent fragmentation of tundra habitat produced individual sky islands that vary in their size and geographic isolation (McCormack et al., 2009). This variation in landscape context may serve as a critical mediator of shifts in alpine taxonomic, functional, and phylogenetic diversity in response to climate change as mountain ecosystems in western North America are currently warming faster than the global average, a trend that is expected to accelerate (Rangwala and Miller, 2012;Pepin et al., 2015). Understanding how sky island area and isolation affect biodiversity change will not only help to identify the best refugial sites for this threatened ecosystem but also potentially generalize to other systems with island-like distributions (e.g., serpentine soils).
Here, we explore if the landscape context (site area and geographic isolation) mediates changes in the biodiversity of refugial alpine sky islands in the Colorado Rocky Mountains over a 12-year period. Specifically, we explored how site area and isolation mediate changes in species richness, Shannon's diversity index, community-level phylogenetic diversity, community weighted mean (CWM) trait values, and functional trait dispersion between our two sampling periods. Based on the principles of island biogeography (MacArthur and Wilson, 1967;Brown and Kodric-Brown, 1977;Connor and McCoy, 1979;Ottaviani et al., 2020a), we predict: 1) that site area will mediate biodiversity change as larger sites have the potential to support a wider range of habitat types and species and/or larger populations; and 2) that geographic isolation will mediate biodiversity change as potential local extinction or declining populations may be offset by dispersal from nearby areas of tundra.

Field Sampling
During the summers of 2007 and 2008, we sampled community composition at seven alpine sites in the Colorado Rocky Mountains (see Spasojevic et al., 2014, for additional details). Using topographic maps, we selected sites to capture variations in size and geographic isolation while considering accessibility. All sites sampled were sky islands: alpine tundra habitat completely surrounded by less-suitable habitat for alpine populations such as rock outcrops, subalpine meadows, or coniferous forests. We sampled species composition and abundance in a series of 1-m 2 plots within each site, where the number of plots was proportional to the total site area and the arrangement of plots were standardized across major topographic gradients. In alpine tundra, topography plays a key role in shaping biodiversity patterns via snow deposition (Walker et al., 2001), and thus this sampling design enabled us to capture the full range of alpine vegetation types including snow-bed (high snow accumulation), moist meadow (intermediate snow accumulation), and fell-field (very low snow accumulation) communities. Within each plot, plant composition was estimated using visual estimations of percent cover of all species with the assistance of a 1 m × 1 m frame containing a 10 cm × 10 cm string grid. Across all sites, we sampled species composition in 153 plots, with the number per site ranging from eight at Greenhorn Mountain to 40 at Niwot Ridge (Figure 1) with an average of 23 plots per site. Plot locations were predetermined and spaced evenly along topographic gradients at each site at ∼200-m intervals. In 2019 and 2020, we revisited these sites, keeping plot number and sampling locations as consistent as possible based on GPS coordinates recorded in 2007/2008. Given the relatively low accuracy of GPS units (±3 m), we examine site-level patterns of biodiversity changes instead of plot-level changes, as we are more confident that we accurately resampled each site, if not individual plots (N = 14 in subsequent analyses, seven sky islands in each sampling period).

Spatial Variables
During site visits, the approximate alpine area was delineated on a topographic map, and the area was later calculated using ImageJ (Rasband, 2007). Geographic isolation was calculated as the area (km 2 ) of non-alpine tundra habitat within a 10-km radius from the center of each site (not including the site area) using ArcGIS 9.0. The classic theory of island biogeography by MacArthur and Wilson (1967) suggests that area more strongly influences extinction rates and geographic isolation more strongly affects colonization rates. However, advances to this theory acknowledge that geographic isolation can influence Frontiers in Ecology and Evolution | www.frontiersin.org extinction rates (the rescue effect: Brown and Kodric-Brown, 1977) and area can influence colonization rates (the passive sampling effect: Connor and McCoy, 1979). Importantly, area in the alpine tundra is often correlated with topographic heterogeneity, where larger sites have a higher probability of capturing a greater variation in topographic variation that will influence biodiversity patterns (Walker et al., 2001).

Trait Measurements
To quantify changes in CWM trait values and functional diversity metrics among sampling periods, we used a combination of trait data that have been collected from previous studies in this system (Spasojevic and Suding, 2012;Spasojevic et al., 2014), from additional measurements conducted at the Niwot Ridge LTER (one of our sites), and from measurements of any new species that were not encountered in our first sampling period. All trait measurements followed established protocols (Perez-Harguindeguy et al., 2013) and were collected on a minimum of 10 individuals for each species. Importantly, in our fist sampling period, we compared trait values among the most abundant species present in both northern and southern sites (approximately 25% of our total species pool) and found that there was little variation in mean trait values along this broad north to south geographic gradient (Spasojevic et al., 2014), suggesting that using trait mean values from this database will allow us to examine broad-scale patterns of functional diversity.
We focused on four putatively important traits: overall height, specific leaf area (SLA), leaf area, and leaf dry matter content (LDMC). Plant vegetative height, a trait often allometrically related to overall plant size (biomass, rooting depth, lateral spread) as well as to competitive interactions for light (Aan et al., 2006), was measured as length from ground level to the highest photosynthetically active tissue. We also collected a fully formed adult leaf, with no signs of damage or senescence at peak biomass. Collected leaves were hydrated overnight prior to weighting fresh mass. Leaves were then dried at 55 • C for 4 days and weighed to determine leaf dry mass. Individual leaf area was calculated from the leaf scans using ImageJ (Rasband, 2007); leaf area is associated with leaf energy and water balance, with various stressors (heat, drought, cold, and high radiation) tending to select for small leaf area (Cornelissen et al., 2003). SLA was calculated as leaf area (cm 2 ) per unit of dry leaf mass (g) and is associated with the leaf economics spectrum, which characterizes a species capacity for stress tolerance vs. resource acquisition (Westoby et al., 2002). Leaf dry matter content was calculated as the ratio of dry mass to fresh mass and is also associated with the leaf economics spectrum (Garnier et al., 2001).

Diversity Metrics
To describe differences in biodiversity among sampling periods, we first calculated species richness, Pielou's evenness, and Shannon's diversity index using the vegan package (Oksanen et al., 2019) in R version 4.0.2 (R Core Team, 2019). We then calculated two metrics of phylogenetic relatedness: mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) using the PICANTE package (Kembel et al., 2010) in R. Mean pairwise distance is a metric of relatedness that measures the mean branch distance among all species pairs in a community (Webb et al., 2002) and is more sensitive to tree wide patterns of phylogenetic clustering and evenness (Kraft et al., 2007). Mean nearest taxon distance is the mean distance separating each species in the community from its closest relative in the community (Webb et al., 2002) and is more sensitive to patterns of clustering and evenness closer to the tips of the phylogeny (Kembel, 2009). We compared these metrics to null communities by randomizing species co-occurrences 9,999 times while maintaining sample richness and species occurrence frequencies. For both of these metrics, we calculated a standard effect size (SES; as per Gurevitch et al., 1992)  For each site in each sampling period, we calculated two complementary functional diversity metrics: CWM trait values (Garnier et al., 2004) and functional dispersion (FDis; Laliberte and Legendre, 2010) for each trait. CWM trait values were calculated as the sum across all species of species' trait values weighted by their relative abundance (Garnier et al., 2004). Following Laliberte and Legendre (2010), we calculated FDis as the mean distance of each species, weighted by relative abundances, to the centroid of all species in a site. Although there are many metrics of functional diversity (reviewed in Mouchet et al., 2010;Schleuter et al., 2010), we focused here on FDis because it is independent of species richness, takes into account species abundances, and can be used for single traits or multiple traits (Laliberte and Legendre, 2010). Moreover, Ricotta and Moretti (2011) proposed a unified analytical framework that combines CWM and a close analog of FDis (Rao's Q). Functional diversity calculations were conducted using the FD package (Laliberte and Legendre, 2010) in R.

Statistical Analysis
To test whether our diversity metrics (i.e., richness, evenness, Shannon's diversity, NRI, NTD, CWMs, FDis) differed among years and if those differences were mediated by refugial area or isolation, we constructed separate generalized linear mixed models (GLMM) with each diversity metric as a response variable and with year, area, isolation, and all possible interactions of those as predictor variables, and site as a random factor-analyses were conducted in JMP version 13.0 (JMP R , Version 13.0. SAS Institute Inc., Cary, NC, 1989-2019).
Lastly, we found that FDis of SLA decreased between sampling periods (F 1,3 = 41.67, P = 0.007; Figures 4D,H) and that there as a significant year × area interaction (F 1,3 = 24.41, P = 0.02; Figure 4D) where FDis of SLA decreased more in smaller sites. Moreover, height had a significant year × area interaction (F 1,3 = 26.30, P = 0.01; Figure 4A) where FDis decreased in smaller sites and increased in larger sites between sampling periods. We did not find a significant year × isolation interaction for FD is of height (F 1,3 = 5.24, P = 0.106; Figure 4E), but did find significant year × area × isolation interaction (F 1,3 = 35.51, P = 0.01). Lastly, we found no difference among years in the FDis of leaf area (F 1,3 = 5.24, P = 0.11; Figures 4B,F)

DISCUSSION
Overall, we found that the taxonomic, phylogenetic, and functional diversity of plant communities in the alpine sky islands of the Colorado Rocky Mountains is changing, and our initial observations suggest that landscape context may mediate the ability of refugia to withstand the impacts of climate change (Gaston and Blackburn, 1996;Ashcroft et al., 2009;Keppel et al., 2015). Specifically, we found an increase in species richness, a trend toward increased phylogenetic diversity, a shift toward a more stress-tolerant functional strategy, and a decrease in the FDis of SLA over time. Importantly, we found that refugial area mediated some of these responses, where larger sites had a larger increase in richness than smaller sites, CWM height decreased in a smaller site and increased in larger sites, and FDis of SLA decreased more in smaller sites. On the other hand, none of our responses was mediated by geographic isolation, which may be due to the rarity of long-distance dispersal events (Tackenberg and Stocklin, 2008) for the long-lived, clonally reproducing species that compose alpine communities (Forbis, 2003;Rossetto and Kooyman, 2005;Herben et al., 2015). Collectively, these results suggest that considering the landscape context (area and geographic isolation) of refugia may be important for conserving refugial sites that provide the most conservation value.

Changes in Taxonomic Diversity
Our observation of increased species richness between sampling periods finds mixed support in the literature. Tundra-wide syntheses of observational studies and global change experiments have generally observed decreases in taxonomic diversity as a result of shrub expansion due to increasing temperatures and decreasing snow cover duration period (Elmendorf et al., 2012;Pearson et al., 2013), and models of climate change forecast further losses of species richness as climate change continues to accelerate (Nabe-Nielsen et al., 2017;Niittynen et al., 2020). However, research in alpine tundra specifically has found contrasting results. In non-Mediterranean European alpine systems, Steinbauer et al. (2018) found that species richness has significantly increased on 87% of 302 mountain summits since 1871 and that these increases were accelerating through time in conjunction with accelerating increases in temperature and precipitation. Steinbauer et al. (2018) attributed the increase in species richness to the movement of subalpine species into the alpine zone due to the amelioration of harsh abiotic conditions. Numerous other studies from the alpine tundra of boreal and temperate Europe corroborate these results and interpretation (Stanisci et al., 2005;Gottfried et al., 2012;Wipf et al., 2013). Studies in Mediterranean alpine tundra have found more mixed results, with some reporting decreases in species richness , and other studies showing increases largely driven by the upward movement of thermophilic generalist species (Jiménez-Alfaro et al., 2014;Evangelista et al., 2016). Fewer studies have been conducted in the alpine of the western North America, and those studies also demonstrate mixed support for species richness increases. In support of our results, at the Niwot Ridge LTER in the Colorado Rocky Mountains (one of our study sites, though we did not use longterm monitoring plots in this study), Spasojevic et al. (2013) found a significant increase in species richness over a 21-year period (sampled at irregular intervals from 1989 to 2010: 1989, 1990, 1995, 1997, 2006, 2008, and 2010) and attributed this trend to the increased establishment of subordinate alpine species as a result of increasing temperature, snow deposition, and nitrogen deposition. However, Scharnagl et al. (2019) found significant decreases in species richness using a different dataset at Niwot Ridge sampled over a 40-year period (sampled every 10 years) and attributed the change in richness to an increase in shrub cover and a decrease in average snow depth. The contrasting results may be a product of the differing sampling periods but exemplify the need for more research in North American alpine ecosystems, as these are two of the only published long-term monitoring studies focused explicitly on community composition within the Rocky Mountains (Elmendorf et al., 2012;Bjorkman et al., 2018).

Changes in Phylogenetic Diversity
In conjunction with the observed increase in species richness between our two sampling periods, we also found a marginally significant (P = 0.05) increase in phylogenetic diversity (using metrics independent of species richness). These patterns suggest that the species present within these alpine refugia are shifting toward being more widely distributed across the phylogeny, capturing greater variation in evolutionary histories potentially as a result of continued climatic warming relaxing some of the environmental filtering in the alpine. Despite the increased interest in linking evolutionary history with climate change (e.g., Botero et al., 2015;Harrison et al., 2020), few studies have considered changes in phylogenetic diversity in response to climate change within alpine tundra ecosystems. In one of the few studies to indirectly address phylogenetic diversity changes, Lesica and McCune (2004), Lesica (2014), and Lesica and Crone (2017) found that species with artic or boreal evolutionary histories have declined in the Glacier National Park (MT, United States) since 1988, and that this decline is worse for dicots than for monocots. While studies by Lesica and McCune (2004), Lesica (2014), and Lesica and Crone (2017) point to interesting trends in how coarse phylogenetic groupings may mediate responses to climate change, we were unable to find research that analyzed how community-level phylogenetic diversity in the alpine shifted in response to climate change using common metrics of phylogenetic diversity (e.g., Faith's PD, MNTD, etc.). This is an underappreciated aspect of biodiversity change that requires more attention moving forward (Gerhold et al., 2015), especially in refugia.

Changes in Functional Diversity
The shift toward stress-tolerant functional traits we observed contrasts with tundra-wide syntheses Pearson et al., 2013;Bjorkman et al., 2018;Steinbauer et al., 2018), which have generally found that increased temperatures drive shifts toward communities with greater average heights and trait values related to the resource acquisitive side of the leaf economics spectrum (high SLA and Leaf N content, low LDMC). However, in those syntheses, soil moisture plays a strong role in mediating the effects of regional temperature on functional composition and diversity. For example, the analysis of 117 tundra sites throughout the northern hemisphere by Bjorkman et al. (2018) found that higher temperatures were only correlated with increases in resource acquisitive traits in wetter tundra sites. In sites where soil moisture was limiting, temperature was correlated less positively with increasing heights and stresstolerant functional traits were favored (Low SLA and Leaf N, high LDMC). In the mountain ranges of the North American West, increased average temperatures have been coupled with an overall decline in snowpack and earlier melt-out dates (Fyfe et al., 2017;Mote et al., 2018;Rhoades et al., 2018). Taken together, these regional climatic changes may be reducing available soil moisture during the summer growing season, thus selecting for alpine species possessing more stress-tolerant functional traits. This is further reinforced by the long-term monitoring work in Glacier National Park of Lesica and McCune (2004) where shifts toward community types associated with lower soil moisture and more stress-tolerant functional strategies have been observed.

Landscape Context
In addition to these changes through time, our results suggest that a deeper consideration of landscape context is important for assessing and conserving refugia. We found that some of the variations in taxonomic richness, phylogenetic diversity, functional composition, and functional diversity, between sampling periods were mediated by site area but not by geographic isolation. Site area can mediate changes in biodiversity within refugia through several mechanisms. First, larger sites are more likely to contain a greater variety of environmental conditions (Bell et al., 1993;Hulshof and Spasojevic, 2020), which may buffer environmental change. This may be the case in alpine systems where stark topography leads to differential distribution of snow on the landscape, creating strong gradients in stress exposure and resource availability and generating a complex mosaic of distinct vegetation types that differ greatly in species composition and productivity (Bowman and Fisk, 2001;Walker et al., 2001;Bowman et al., 2003;Seastedt et al., 2004;Litaor et al., 2008). Second, larger sites may be able to support populations and communities of larger size, and as community size decreases in smaller refugia, the probability of local extinctions is predicted to increase due to demographic stochasticity (Loreau and de Mazancourt, 2008). Third, larger sites may be easier to colonize than smaller sites, potentially counteracting local extinction (Brown and Kodric-Brown, 1977).
Due to the remote locations of our sampling site, we were unable to collect data on environmental heterogeneity or community size at the scale necessary to test these alternative mechanisms. However, we found that increases in taxonomic richness and phylogenetic diversity were greater in larger sky islands and that these sky islands seemed to display increased stability in terms of functional composition and diversity. For example, CWM values for leaf traits shifted more strongly toward stress-tolerant functional strategies in smaller sky islands than they did in larger sky islands, suggesting that increased temperatures and decreased snowpack may have more intense effects on smaller sky islands. Furthermore, FDis values for SLA declined substantially for smaller sky islands while large sky islands remained relatively constant, indicating a reduction in the breadth of functional strategies that small sky islands can support.
Geographic isolation, on the other hand, did not play a strong role in mediating biodiversity change likely due to the importance of clonal reproduction for many alpine plants and physical constraints on seed dispersal. In other systems, research has demonstrated that a functional trade-off exists between a species ability to locally persist through vegetative reproduction and its capacity for sexual reproduction and seed dispersal (Rossetto and Kooyman, 2005;Herben et al., 2015). For the long-lived perennial plant species that compose the majority of alpine flora of the Rocky Mountains, clonal reproduction can be just as important as sexual reproduction for recruitment (Angevine, 1983;Eriksson, 1989), potentially obscuring the influence of seed dispersal, and thus geographic isolation, on biodiversity changes (Forbis, 2003). Furthermore, seed dispersal among sky islands is difficult due to the vast distances between areas of suitable habitat (Tackenberg and Stocklin, 2008) and the lack of long-distance dispersal mechanisms in alpine species necessary to traverse these distances (Morgan and Venn, 2017).

Limitations and Future Directions
While the trends we found are generally supported by the literature, it is important to note two potential sources of bias within our data. First, the differences in exact plot locations due to inherent inaccuracy of GPS likely led to slightly different sampling locations. In this system, where there is high betadiversity over short distances (Bowman and Fisk, 2001;Walker et al., 2001;Bowman et al., 2003;Seastedt et al., 2004;Litaor et al., 2008), some of the variations we find could be due to altered sampling locations. Because of this issue, we examined site-level patterns of biodiversity changes instead of plot-level changes as we are more confident that we accurately resampled each site, if not individual plots. Furthermore, a comparison of dominant species (those composing a cumulative sum of ∼70% of site relative abundance) for each site between sampling periods shows that the identity of dominant species is relatively stable and suggests that we sampled similar vegetation types in each time period (Supplementary Table 2). Interestingly, this table qualitatively reinforces the results of our landscape context analysis, with larger sites (Niwot, Cornwall, Colorado Mines, Mummy) showing less turnover in dominant species identity than smaller sites (Boreas, Buffalo, and Greenhorn). Second, our sampling intentionally avoided areas with shrubs to instead focus on biodiversity changes in the herbaceous community. Changes in shrub abundance are a major driver of change in the alpine (Elmendorf et al., 2012;Bjorkman et al., 2018) and may explain some of the contrasting results found with other studies that included shrubs. Specifically, our data do not capture potential shrub encroachment into the alpine and associated changes in functional diversity and composition, particularly increases in height, that might accompany this. Despite these limitations, our research suggests that the landscape context of refugia may be an important factor mediating climate change-induced shifts in biodiversity patterns.
We further highlight two important future directions for biodiversity research in climate change refugia. First, examining biodiversity patterns across scales (alpha, beta, and gamma) and linking these patterns to explicit mechanisms (environmental heterogeneity, community size) will provide a clearer picture of how landscape context mediates biodiversity change in refugia and how best to conserve refugia. Second, shifts in functional diversity patterns caused by changing climate may be highly influenced by belowground traits (Ottaviani et al., 2020b). This aspect of functional diversity deserves more attention, especially in abiotically stressful systems like alpine tundra where plants can allocate up to 80% of biomass to belowground organs (Klimešová et al., 2019).

CONCLUSION
Despite the recognized importance of refugia for conserving biodiversity in the face of climate change, there is still a need to better identify high quality refugial sites (Keppel et al., 2012(Keppel et al., , 2015(Keppel et al., , 2018. Our results suggest that landscape context has the potential to mediate the ability of refugial sites to withstand the impacts of climate change (Gaston and Blackburn, 1996;Ashcroft et al., 2009;Keppel et al., 2015). Moreover, our results also highlight that the life history of the species in a refugial landscape is important to consider. In alpine tundra, where species have long life spans (Steinger et al., 1996;Morris and Doak, 1998), clonal reproduction is common (Angevine, 1983;Eriksson, 1989), and many species lack the dispersal mechanisms necessary for longdistance movement (Tackenberg and Stocklin, 2008), we found that area played a stronger role in mediating biodiversity change than geographic isolation. In systems with short-lived annual species, geographic isolation and its influence on dispersal may be relatively more important. Future research examining the landscape context of refugia across a broader range of ecosystems and species will be critical for generalizing which sites will be most buffered from the impacts of climate change and offer the greatest potential for conserving biodiversity.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: Resampling is part of dissertation chapter that has yet to be published. We will make all data publicly available when the dissertation is published. Requests to access these datasets should be directed to MS and JH.

AUTHOR CONTRIBUTIONS
MS conducted the 2007/2008 survey and analysis. JH conducted the 2019/2020 survey, and wrote the first draft of the manuscript with input from MS. JH and MS designed the study. Both authors contributed to the article and approved the submitted version.

FUNDING
This work was supported in part by the Niwot Ridge Long-Term Ecological Research Program (NSF DEB 1457827).