Climate disequilibrium of fishes along elevation and latitudinal gradients: Implications for climate tracking

Differences in realized and fundamental thermal niches can reveal how temperature constrains species distributions. Because climate‐induced range shifts assume that temperature influences distribution limits (i.e. climate equilibrium assumption), understanding the factors that determine the realized thermal niche of species is critical for understanding their climate tracking abilities.


| INTRODUC TI ON
A guiding assumption behind climate-induced range shifts is that temperature constraints on physiological processes is the primary driver of distributional limits (i.e. climate equilibrium hypothesis; Early & Sax, 2014;Freeman, 2016). Hence, warming temperatures at distribution limits would result in an expansion or contraction of a species range, as long as non-temperature factors also do not influence distribution limits. Despite numerous studies documenting climate-induced range shifts, there is large variability in the extent to which some species are shifting distributions compared with others (Freeman et al., 2018;Gibson-Reinemer & Rahel, 2015). Although physiological stress related to climate is believed to be the primary driver of high elevation and high latitude limits, especially for tropical species (Cuesta et al., 2020;Janzen, 1967), thermal limits do not always correspond to observed distribution limits (Allen-Ankins & Stoffels, 2017;Freeman, 2016;Slatyer & Schoville, 2016;Sunday et al., 2012). Thus, additional factors are likely limiting the capacity for species range shifts.
One evaluation of whether organisms will track climate change is by comparing their realized thermal niche (i.e. observed relationship between distribution limits and temperature) with their fundamental thermal niche (i.e. experimentally determined limits to temperature; Allen-Ankins & Stoffels, 2017; Sunday et al., 2012).
If a species' realized thermal niche is larger (overfilling distribution) or smaller (underfilling distribution) than their fundamental thermal niche, then factors unrelated to temperature are likely limiting the species' distribution (Slatyer & Schoville, 2016). Range shifts from climate change are unlikely for species exhibiting these "climatic mismatches", such as for species with small range sizes who underfill their distributions (Sunday et al., 2015). An alternative approach is determining if realized thermal niches match along different thermal gradients, such as those related to elevation and latitude (Kollas et al., 2014;Pellissier et al., 2013;Randin et al., 2013). The extent to which species exhibit climatic mismatches between realized and potential distributional limits can provide insight into climate tracking abilities, because thermal niche filling is frequently associated with traits of dispersal and mobility (Siefert et al., 2015;Svenning & Skov, 2004).
The ability of organisms to track and follow climatic shifts is dependent upon biological traits and the environmental gradient examined. Range shifts are limited for trait specialists who exhibit patterns of strong niche conservatism (Comte et al., 2014;MacLean & Beissinger, 2017). Sunday et al. (2015) found that highly mobile and trophic generalist species had climate-induced range shifts of greater magnitude than low mobility, trophic specialist species. In addition to biological traits, the presence of natural (steep landscape gradients; mountain ranges) or anthropogenic (dams) barriers will also restrict the ability of organisms to track temperature changes (Comte & Olden, 2017a;Gibson-Reinemer et al., 2017;LeMoine et al., 2020;Radinger et al., 2017). Organisms may be unable to shift their distributions upslope because environmental factors unrelated to temperature also change with elevation (e.g. habitat size; Buisson et al., 2008;Rahel & Nibbelink, 1999). Combining trait-based approaches and evaluations of biogeographical distribution limits is critical to inferring climate tracking abilities (Svenning & Skov, 2004).
Our study tested the climate equilibrium hypothesis for fish species along elevation and latitudinal gradients. Both freshwater and marine fishes have experienced profound range shifts from climate change (Comte et al., 2013;Pinsky et al., 2013), yet there are important differences in dispersal potential between the two groups that influences the ability of species to fill realized thermal niches.
Freshwater fishes have distributions constrained within dendritic networks where biogeographical barriers (e.g. oceans, mountains, catchment divides, waterfalls) often limit dispersal and connectivity among habitats (Carvajal-Quintero et al., 2019;Leprieur et al., 2009;Rahel, 2007), whereas marine species have few constraints on mobility and dispersal (Comte & Olden, 2017a). Thus, marine organisms should have realized thermal niches which closely align to their fundamental thermal niches, resulting in predictable range shift patterns (Sunday et al., 2012). Here, we focus on cold-edge distribution limits because cold-edge limits show stronger niche conservatism and appear to be more associated with physiological processes than warm-edge limits (Comte et al., 2014;Pellissier et al., 2013).
Using established methods for comparing fundamental and realized thermal niches across spatial gradients (e.g. Randin et al., 2013;Sunday et al., 2011Sunday et al., , 2012, our objectives were to (1) compare the fundamental and realized cold-edge limits of fish species along latitudinal gradients at a global scale and (2) compare the cold-edge elevation and latitudinal realized thermal niche limits of freshwater fishes across the Rocky Mountains-Great Plains continuum of North America. Our first hypothesis, pertaining to objective 1, was that freshwater species would underfill their cold-edge limits more than marine species due to the biogeographical barriers common in freshwater systems.
Our second hypothesis, also pertaining to objective 1, was that morphological traits related to low dispersal and non-migratory behaviour would be associated with the underfilling of cold-edge latitudinal limits. Our third hypothesis, related to objective 2, was that species with morphological traits of high dispersal would be more likely to have an elevation thermal limit that exceeds their latitudinal thermal limit (i.e. overfilling elevation), given the importance of mobility traits for occupying high elevation, small stream habitats (Bower & Winemiller, 2019;Lamouroux et al., 2002). The fourth hypothesis, pertaining to objective 2, was that more species would underfill elevation limits relative to latitudinal limits given strong environmental filtering along upstream-downstream gradients .

| Objective 1: Fundamental and realized thermal niches along latitudinal gradients
We compiled a database on the fundamental niche (i.e. cold temperature limits) of freshwater and marine fish species as estimated from thermal tolerance studies following the methods of Sunday et al. (2011Sunday et al. ( , 2012. The database was constructed from an extensive survey of the peer-reviewed literature (Web of Science and Google Scholar) combined with online, published databases of species' thermal tolerances (e.g. GlobTherm; Bennett et al., 2018;full citation list available with open access data). We only considered studies that quantified thermal minimum tolerance via one of two metrics. The first metric was lethal minimum temperature (LT min ), which is quantified as the water temperature where a predefined proportion of a population dies following experimental temperature exposure. The second metric was critical thermal minimum temperatures (CT min ), which is defined as the water temperature at which there is a loss of equilibrium for fish (e.g. mobility failure, onset of muscle spasms; Beitinger et al., 2000).
Given that cold tolerance is less studied than warm tolerance (Beitinger et al., 2000), we considered multiple study design components to maximize the number of species for which thermal minimal tolerances could be included. Thermal tolerance studies often acclimate organisms across a range of temperatures to consider thermal plasticity because the temperature an individual is acclimated to influences its observed critical or lethal limit (i.e. colder acclimation results in a lower thermal minimum; Beitinger et al., 2000). Thus, when assigning species thermal minima values, preference was given to studies that acclimated species across multiple temperatures, and the thermal minima of a given species was assigned using the lowest acclimation temperature when multiple thermal minima values existed. We excluded fish species that were considered endemic with very small latitudinal range sizes (<5°) because such species are unlikely to have poleward distributions limited by temperature (Sunday et al., 2015). We also excluded species with poleward distribution limits restricted by non-thermal barriers (e.g. edge of a continent, islands, mountain range; Sunday et al., 2012).
To determine the extent to which fish species fill out the coldedge of their latitudinal range, the associated latitude for a species' realized and fundamental thermal niche was determined based on the mean water temperature of the coldest month. Realized latitudinal limits (i.e. observed poleward limits) for marine species were based on water temperatures derived from gridded sea surface water temperatures averaged for the period 1971-2000 (https:// daac.ornl.gov/cgi-bin/dsvie wer.pl?ds_id=980). Fundamental latitudinal limits for marine species were based on gridded sea surface temperatures that corresponded to the laboratory-derived fundamental thermal minima (i.e. potential poleward limits). Water temperatures for freshwater species were derived from gridded air temperatures from the WorldClim database (http://www.world clim. com/) based on average monthly climate data for the period 1950-2000. Air temperatures were then converted to water temperatures following Punzet et al. (2012) and adjusted accordingly depending upon the climatic region in which the cold temperature limit for a species occurred.
Rather than simply averaging mean temperatures into latitudinal bands (e.g. Sunday et al., 2012), we compiled coldest month temperature measurements of each species fundamental thermal minima at their potential poleward limits along a transect (an isotherm; ~10-15 points) across the range of longitudes occupied within a species' native range. This method accounted for any continental and regional variation in global temperatures, which could be overlooked by collapsing temperatures into latitudinal bands. We also inferred cold-edge limits within a single hemisphere, which was the hemisphere with the species' largest range and most poleward-extending boundary. The realized cold-edge latitudinal limit of each species within only its native range was determined by consulting multiple published and online distribution records to confirm native range extents (e.g. Global Biodiversity Information Facility, http://www.  et al. (2011, 2012), to determine how biological traits and habitat differences influence thermal niche filling patterns while controlling for method-based variables. We evaluated two morphological traits that have been associated with dispersal and mobility in previous fishbased studies (Bower & Winemiller, 2019;Comte & Olden, 2018;Troia & Gido, 2015). First was body shape (BS), which was quantified as the ratio of total length to maximum body depth. High BS values indicate streamlining and are associated with species occupying high velocity environments. The second was swim factor (SF), which was quantified as the ratio of minimum caudal peduncle depth to maximum caudal fin depth. Low values of SF are associated with deep caudal fins capable of generating a high amount of thrust. Digital measurements were taken from online, vouchered photographs of specimens from the Global Biodiversity Information Facility, which were required to have fins properly spread and only photographs of adult specimens were considered. Anguilliform-swimming species (n = 5) were excluded from analyses due to their unique morphologies. In addition to the macro-habitat categories of marine or freshwater, we also broadly categorized fishes into four micro-habitat categories (pelagic, bentho-pelagic, demersal, reef-associated) to test for additional habitat differences in cold-edge niche filling. For example, we may expect species occupying deep-water habitats with colder temperatures (demersal traits) to exhibit different filling patterns than species occupying near-surface habitats (pelagic or reef-associated traits) because of temperature differences across depth levels (Brito-Morales et al., 2020;Pinsky et al., 2013). Finally, we included whether a species exhibited migratory behaviour as a predictor of underfilling or overfilling. All categorical traits were derived from FishB ase.org.
We also included an interaction term between macro-habitat and the species' native latitudinal median (°) to distinguish between temperate and tropical species. This interaction was included because high latitude species have been documented over-filling their coldedge limit more than low latitude species (Svenning & Skov, 2004), and this pattern varies for terrestrial versus marine organisms (Sunday et al., 2012). Method-based covariates included whether studies were performed on juvenile or adult life stages because juveniles tolerate a narrower thermal range than adults (Beitinger et al., 2000), acclimation history (acclimated or not), and the specific thermal minima metric used (LT min or CT min ). Variation in these co-variates could bias overfilling or underfilling patterns (i.e. in one species, an estimate from the juvenile life stage using CT min could produce more underfilling than an adult estimate using LT min ). We included a hierarchically nested random intercept term for taxonomy (Order (O), Family (F), Genus (G); like Sunday et al., 2011Sunday et al., , 2012 to account for non-independence among species because cold-limits and mobility traits are conserved across taxonomic families Pellissier et al., 2013). Latitudinal range size (°) was also considered as a potential co-variate to account for the relationship between range size and thermal tolerance (Gutiérrez-Pesquera et al., 2016), but was excluded due to high collinearity with latitudinal median (see Table S1; Figure S1). Our full model in R syntax was:

| Objective 2: Realized thermal niche between elevation and latitudinal gradients
To compare elevation and latitudinal limits, we used a second, independent dataset for stream fish species occurrence across the Missouri River basin of the Rocky Mountains -Great Plains region, U.S.A. Species elevation limits within their native range were determined for 28 species using a comprehensive dataset of 3189 stream surveys conducted across five states (see  for descriptions of the dataset). The dataset covered a large climatic gradient from the Great Plains and into the Rocky Mountains (elevation range: 310-3246 m above sea level; latitudinal range: 38.5-49.0°). We only included species that had >100 occurrences in the dataset (3.2%; similar to Siefert et al., 2015). For each species, we calculated the 0%, 2.5%, 5%, and 7.5% elevation quantiles from the stream survey dataset to estimate realized upper elevation limits for each species. While the 0% quantile may represent the "truest" high elevation limit, we did not use the 0% quantile to avoid potential outlier occurrences (e.g. misidentifications, introduced populations).
We present results based on analyses using the 2.5% quantile because all quantile estimates were highly correlated (all r > 0.95) and because the 2.5% quantile has been the most frequently applied in similar studies (e.g. Randin et al., 2013;Siefert et al., 2015).
We estimated minimum air temperatures at high elevation limits for each species using three climate sources. The first source was the WorldClim database used in objective 1 (period from 1950-2000), and air temperatures during the coldest month were determined based on the associated latitude and longitude of the location representing the 2.5% quantile elevation limit. The second was a linear regression model constructed from 133 air temperature stations across our study region of the Great Plains -Rocky Mountains, which were averaged for the period 1981-2010. The linear model was based on air temperatures during the coldest month, included predictor variables of elevation and latitude, and had high performance for predicting coldest month temperatures based on those variables (R 2 = 0.81, p < 0.001). The third source was the daily surface weather and climatological summaries database (Daymet; daymet.ornl.gov). Daymet estimates daily maximum and minimum air temperature measurements for gridded cells at a spatial resolution of 1 km 2 . Maximum and minimum daily temperature measurements based on the latitude and longitude of the 2.5% quantile elevation limit of each species were averaged for the coldest month, and the coldest months were averaged for the period 1980-2010. overfilling of elevation limits). Species whose elevation filling value did not cross zero when subtracting or adding the SD of the three cold air temperature sources were considered to be "substantially" overfilling or underfilling their realized elevation limits. Thus, species were classified into three groups: (1) those with matching elevation and latitudinal thermal limits (mean filling values overlap 0 ± SD), (2) those "substantially" overfilling their elevation limits (values >0 ± SD) and (3) those "substantially" underfilling their elevation limits (values <0 ± SD). Because latitudinal cold limits were colder than elevation limits, all latitudinal cold-edge limit values that exceeded the coldest possible elevation temperature (−11.0°C) were standardized to that minimum coldest elevation temperature (Table S2).
To assess patterns of elevation filling, we used a similar mixedeffects model as objective 1 (lmer function; R v. 3.4.2) with the following R syntax: using the same morphological traits as objective 1 (shape factor [BS] and swim factor [SF]). We also included a life history trait variable (LH), which was quantified from a principal components analysis (PCA) of four commonly used fish life history traits (total adult length, age at maturity, longevity, fecundity). We retained the first PCA axis as an important trait variable, which summarized 70.1% of the variation among traits and was the only axis with an eigenvalue >1.0. The axis distinguished large-bodied, late maturing, high fecundity species from smallbodied, short-lived, low fecundity species. Life history traits were considered important for this analysis because a complementary study found that these traits were strongly correlated with thermal gradients of elevation and latitude . All trait values were derived from that study.

| Objective 1: Latitudinal thermal niche filling
For fundamental thermal minima data collected on 126 fish species, 95 met the criteria for analyses. While there was similar representation between marine (n = 50) and freshwater (n = 45) species, thermal minima data was more abundant for tropical species (median latitudes <23°; n = 55) than temperate species (median latitudes >35°; n = 23) and subtropical species (median latitudes between 23° and 35°; n = 17). There was also a bias towards cold-edge limits being quantified in the northern hemisphere (n = 84) and a bias to CT min limits over LT min limits (n = 75 vs. 20, respectively). Finally, there was large variation in how the 95 species filled out their fundamental thermal niche across latitudes (mean = −1.1°; SD = 8.8°; range = −19.4 to 31.3°; Figure 1).
After accounting for non-significant method-based co-variates of life stage, acclimation history, and thermal minima metric (all coefficients [β] with p ≥ 0.203; Table 1), we did not observe support for our first hypothesis because freshwater species were not more likely than marine species to underfill the cold-edge of their latitudinal distribution (macro-habitat; β = −0.72, p = 0.873; Figure 1).
Similarly for hypothesis two, neither body shape nor swim factor were significant predictors (both p ≥ 0.248) and there were no differences among micro-habitat categories for explaining cold-edge

| Objective 2: Elevation and latitudinal niche filling comparisons
There was a positive correlation between a species' cold-edge elevation limit and cold-edge latitudinal limit (r = 0.47, p = 0.010), with high elevation species also having high latitude limits. However, elevation and latitudinal distribution limits were not significant predictors of elevation filling patterns (both p ≥ 0.100; Table 2). In contrast to our third hypothesis, neither body shape (β = 0.04, p = 0.944), swim factor (β = 13.28, p = 0.081), nor the life history trait variable (β = −0.59, p = 0.170) were significant predictors of whether species overfilled their cold-edge elevation limit relative to their latitudinal limit. Despite observing no significant relationships, fixed effects from the elevation filling model explained 31.1% of the variation in filling patterns (R 2 = 0.31).
Longnose sucker (Catostomus catostomus) and central stoneroller (Campostoma anomalum) were the only species with matching elevation and latitudinal limits (Figure 4). Only longnose sucker occupied the highest elevation limit predicted from its latitudinal limit, even though 16 species (57.1%) theoretically could have based on their latitudinal cold-edge limits (green species clump on left-side; Figure 3).  Guisan et al., 2014). Fewer assessments have tested the climate equilibrium hypothesis across thermal gradients in species' native ranges, especially for elevation gradients, likely due to the need for temperature and distribution data at fine spatial resolutions (Kambach et al., 2019;VanCompernolle et al., 2019).

F I G U R E 2 (a) Latitudinal
However, such studies reached similar conclusions on a lack of support for the climate equilibrium hypothesis. Sunday et al. (2012) compared fundamental and realized thermal niches across latitudes for a wide range of ectothermic groups, similar to our analysis in objective 1, and found thermal niche filling patterns varied widely among groups. Whereas terrestrial ectotherms were more likely to overfill their cold-edge distribution limit, marine ectotherms were more likely to underfill or match their cold-edge limit. Randin et al. (2013) demonstrated that 61.1% of 18 deciduous tree species did not have matching elevation and latitudinal cold-edge limits, similar to our analysis in objective 2. They also documented that the majority of tree species underfilled their latitudinal limit relative to their elevation limit, contrary to our results of fishes being more likely to underfill elevation limits relative to latitudinal limits. Similarly, Freeman (2016) and Slatyer & Schoville (2016) demonstrated that the lower thermal tolerance limits of birds and ground beetles, respectively, were not correlated with minimum temperatures experienced at high elevation limits.
An intuitive hypothesis is that the underfilling or overfilling of thermal niches in these studies is related to dispersal and mobility because larger range sizes are correlated with greater movement abilities (Calosi et al., 2010;Li et al., 2016;Luiz et al., 2013). While this hypothesis has been largely unexplored in the context of the climate equilibrium assumption, Siefert et al. (2015) demonstrated that tree species with unassisted seed dispersal and short heights Comte & Olden, 2017a; Rahel, 2007).
Swim factor and body shape were also not significant predictors of thermal niche filling along elevation or latitudinal gradients, which is surprising because morphometrics indicative of strong, streamlined body designs are associated with long distance dispersal behaviours and are favourable for persistence in high elevation environments (Comte & Olden, 2018;Lamouroux et al., 2002). Rather, migration behaviour was the only biological trait associated with thermal niche filling, with migratory species overfilling their latitudinal cold-edge boundary more than non-migratory species. Migratory species are more likely to overfill their cold-edge boundaries because (1) they are big-bodied species with strong swimming abilities and (2) migrations allow them to occupy seasonal habitats that are outside their non-migratory, thermal range (Comte & Olden, 2018). Hence, because migratory species have such high vagility, they are likely to track climatic changes better than smallbodied species with low vagility (Radinger et al., 2017;Sunday et al., 2015;VanCompernolle et al., 2019).
Instead of dispersal-related differences, the most prominent pattern we observed was related to latitudinal median, with temperate species more likely to overfill their cold-edge latitudinal limit and tropical species more likely to underfill. We propose several  Table S2 hypotheses to explain this pattern of high latitude overfilling and low latitude underfilling. First, the extensive levels of overfilling in our study may be due to inaccurate estimations of laboratory derived fundamental thermal minima's, especially if laboratory studies do not use low acclimation temperatures that would consistently overestimate lower thermal tolerance (Beitinger et al., 2000). The second hypothesis is that species at high latitudes have behavioural adaptions to persist under extreme cold (Sunday et al., 2012), such as entering dormant states (e.g. hibernation) to reduce metabolic demands that could allow organisms to survive in places colder than expected based on their fundamental thermal minima's (Ultsch, 1989). A third hypothesis is that the actual minimum temperature experienced by an organism is warmer than the temperature inferred from our macro-climate estimates, such as differences between surface water temperature and water temperatures across depth levels (Brito-Morales et al., 2020;Li et al., 2019). Such temperature differences are attributable to micro-climatic conditions that produce warmer minimum temperatures that more accurately reflect growth and survival compared with macro-climatic temperatures (e.g. the presence of ice cover for thermal insulation or warming from groundwater inputs; Cunjack, 1996;French et al., 2017). Hence, macro-climatic estimates of temperature for realized cold-edge limits may have weaker predictive power for accurately inferring thermal niche limits compared with fine-scale, micro-climatic conditions (Gutiérrez-Pesquera et al., 2016).
Indeed, fish species with the highest thermal tolerances also have the lowest thermal plasticity (Comte & Olden, 2017b). Hence, thermal niche filling in fishes is likely more reflective of evolutionary responses to climate variation than dispersal processes. In addition to tropical fish species having the highest sensitivity to climate warming (Comte & Olden, 2017a;Pinsky et al., 2019), low latitude species with small range sizes are the most likely to experience range contractions while not exhibiting range expansions from climate change (Botts et al., 2013;Jarić et al., 2019;Sunday et al., 2015).
However, these hypotheses regarding thermal niche filling could also be confounded by range size relationships, such as Rapoport's Rule, whereby high latitude species have larger range sizes than low latitude species (Li et al., 2016). Relationships between latitude and range size could be attributable to ecological and evolutionary processes unrelated to thermal tolerance, such as ecological strategies pertaining to habitat and trophic niche breadth (Hecnar, 1999;Slatyer et al., 2013) or the importance of biotic interactions for determining range limits (e.g. competition and predation; Louthan et al., 2015). We observed a positive relationship between latitudinal range size and filling patterns in our study confirming this confounding hypothesis ( Figure S1; see also Gutiérrez-Pesquera et al., 2016), but this relationship did not substantially alter the conclusions of our objective 1 analysis (Table S1). However, a recent study on marine fishes have revealed evidence for an inverse Rapoport's Rule (i.e. larger range sizes for tropical species than temperate species; Pie et al., 2021) that has challenged this long-standing tenet of macroecology. Overall, our conclusions emphasize the need for a better understanding of whether thermal niche filling has a causative or correlative association with geographic range size.
Despite elevation and latitude sharing similar thermal changes from cold to warm, they are not identical gradients. Each has different barriers to dispersal that could explain the overfilling and underfilling patterns observed in objective 2 .
Most prominently for latitude, relationships between deglaciation and post-glacial dispersal likely limited the poleward extent of many freshwater species. This has been coined the "post-glacial migration lag" hypothesis, whereby species with poor dispersal or low latitude distributions lag in filling out their cold-edge, high latitude limits (Siefert et al., 2015;Svenning & Skov, 2004). While not significant (p = 0.100), there was weak support that low latitude species were more likely to overfill their elevation limit relative to their latitudinal limit. Low latitude species in the Rocky Mountains-Great Plains region were late-dispersal species when the glaciers receded (red species on the far right of Figure 3) and have had more time to disperse to high elevations than high latitudes. In contrast, early-dispersal species have had more time to expand poleward (Hoagstrom & Berry, 2006). Indeed, glacial history and historical connectivity of water bodies often best explains the range patterns of freshwater fish (Carvajal-Quintero et al., 2019).
Because elevation corresponds to major habitat changes in running water ecosystems, other abiotic variables can limit the upstream elevation limit of freshwater fish besides temperature (Lipsey et al., 2005;Vannote et al., 1980). The most common habitat limitation constraining upstream distributions as prominently as temperature is stream size (Buisson et al., 2008;Rahel & Nibbelink, 1999).
Most of the 57.1% of species underfilling their high elevation limits (left-side; Figure 3) occupy large-stream or river habitats farther downstream (Rahel & Hubert, 1991). A second habitat limitation is that increasing channel slope results in faster hydraulic gradients, which requires species to have high burst swimming capabilities and adaptations for surviving energetically demanding flow conditions (Bower & Winemiller, 2019;Lamouroux et al., 2002). We observed a striking dichotomy for some species to support this distributional constraint from channel slope. Despite several species sharing similar cold-edge temperature limits, they differed profoundly in the maximum channel slope habitat they could occupy ( Figure S2).
Hence, steep channel gradients are likely preventing species from filling their elevation thermal niche. These changes in habitat size, channel slope, and geomorphology have also been suggested to limit climate-induced range shifts (Gibson-Reinemer et al., 2017).
Understanding the relationship between thermal tolerance limits and range limits is crucial for understanding the ability of species to track climate change Sunday et al., 2012), and a comparison of realized and fundamental thermal niche data serves as a proxy for estimating the sensitivity of species to warming (Freeman, 2016). In that regard, our idiosyncratic patterns of thermal niche filling support the literature documenting idiosyncratic range expansions (Freeman et al., 2018;Gibson-Reinemer & Rahel, 2015). These results support the observation that thermal constraints are acting on the elevation limits of these species because they are underfilling their latitudinal limits (i.e. they are low latitude species limited by "post-glacial migration lag"). Hence, the constraints of cold temperatures at elevation limits will be relaxed with climate warming and elicit expansions for these species. Quantitative studies testing whether thermal niche filling patterns actually correspond to climateinduced range shifts are scarce. In one such study, Freeman (2016) demonstrated that the magnitude of climate-induced range expansions for bird species was unrelated to climatic mismatches between their fundamental, thermal minima and their cold-edge elevation limit. Determining the role of temperature versus additional abiotic and biotic factors in setting elevation and latitudinal limits should improve our ability to explain idiosyncratic range shifts from climate change.

| CON CLUS IONS
The climate equilibrium hypothesis is a critical assumption of species distribution models attempting to predict range shifts under climate change. However, as concluded in similar studies on plants (Early & Sax, 2014;Kollas et al., 2014;Randin et al., 2013) and other ectothermic animals (Sunday et al., 2012), our results indicate that half of fish species do not have matching fundamental and realized thermal niches across thermal gradients. Hence, many species are likely limited by additional abiotic or biotic factors at their cold-edge limits (Slatyer & Schoville, 2016), and a lack of knowledge about these additional factors will hinder the accuracy of models for predicting climate-induced range shifts (Buisson et al., 2008). In particular, our results describing the underfilling of 78% of fish species along the elevation gradient of the Rocky Mountain-Great Plains suggest that distribution models focused on thermal constraints may overestimate the potential for expansions to high elevations in many fish species. Finally, we recommend future research comparing laboratory-derived, fundamental thermal minima data of freshwater species with their high elevation distribution limits (e.g. Freeman, 2016;Slatyer & Schoville, 2016). The ability to make this comparison in our study was limited by the lack of thermal minima data for many fish species, especially when compared with the abundance of data on maximum thermal tolerance (Beitinger et al., 2000).
Nonetheless, such studies will provide greater insight into the role of temperature versus other factors in setting cold-edge, high elevation limits.

ACK N OWLED G EM ENTS
Collection permits were not required for data collection in this publication. Funding was provided by numerous scholarships through the Department of Zoology and Physiology and a grant from the Biodiversity Institute at the University of Wyoming. We thank Amy Krist, Annika Walters, Daniel Laughlin, Michael Dillon, Shannon Albeke and two anonymous reviewers for providing comments on earlier drafts that greatly improved the manuscript.

CO N FLI C T S O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and the corresponding R code for the analyses of objective