An ecothermal paradox: bull trout populations diverge in response to thermal landscapes across a broad latitudinal gradient

Maintaining natural thermal regimes in montane stream networks is critical for many species, but as climate warms, thermal regimes will undoubtedly change. Mitigating impacts of changing thermal regimes on freshwater biodiversity requires knowledge of which elements of the thermal regime are limiting factors for aquatic biota. We used full-year stream temperature records sampled across a broad latitudinal gradient to describe the diversity of the thermal landscapes that bull trout (Salvelinus confluentus) occupy and identify potential divergences from thermal regimes where this species has been studied previously. Populations of bull trout occupied stenothermic, cold thermal niches in streams that exhibited low to moderate thermal sensitivity throughout the species’ range. However, winter thermal regimes in the central and northernmost streams were colder and more stable than in the southernmost streams, reflecting differences in sensitivity to air temperature variation and contributions of perennial groundwater to baseflow. In the southernmost streams, bull trout distributions appeared to be regulated by warm summer temperatures, whereas in northern streams, unsuitably cold temperatures may be more limiting. Our results also suggest that local differences in the extent of complete freezing during winter among northern streams may further limit the distributions of suitable habitats. Contrasts in limiting factors at bull trout range extents would suggest differential responses to climate warming wherein northern populations extend their range while southern populations contract, and an overall change in species status that is less dire than previously anticipated.

temperature is an important dimension of an organism's ecological niche, the geographic distributions of taxa across the globe are determined, at least in part, by the range of temperatures species can occupy-i.e., realized thermal niche (Angilletta 2009;Sexton et al. 2009;Peterson 2011). Therefore, temperatures that species experience at their north-south distributional limits provide a high-level constraint often associated with physiological limitations (Pörtner and Farrell 2008;Somero 2010;Bates and Morley 2020). When organisms experience temperatures outside their thermal niche, they must disperse to track thermally suitable habitats, adjust behaviorally or phenotypically (acclimate), adapt (over the course of generations), or risk local extirpation (Parmesan 2006;Sunday et al. 2012;Sunday et al. 2014). As climate warms, populations of species occupying warm-edge niche boundaries in the Northern Hemisphere are likely to experience range contractions as individuals are forced to move to higher elevations or poleward to track suitable thermal habitats (Parmesan 2006;Deutsch et al. 2008;Sunday et al. 2014). Less studied, however, are the range extensions and associated mechanisms that may occur in populations of the same species that occupy cold-edge niche boundaries near a northern range terminus (Chamaillé-Jammes et al., 2006;Clarke and Zani 2012;Campana et al. 2020). The discrepancy in the effects of climate warming at range extents could lead to outcomes wherein no net loss occurs for a species or the loss that occurs is smaller than otherwise estimated by studies focusing on range subsets where losses are most likely to occur.
The above considerations are particularly relevant for lotic ectotherms that are strongly controlled by temperatures and inhabit linear networks where dispersal opportunities are limited (Fagan 2002;Somero 2010;Sunday et al. 2014). Streams and rivers have a demonstrated sensitivity to air temperature increases and have been warming in association with regional climate trends (Isaak et al. 2017a;Michel et al. 2020;Su et al. 2021) that are often most pronounced at higher latitudes (Prowse et al. 2006a, b;Heino, 2020). This has precipitated numerous regional climate risk assessments for societally important cold-water fishes such as salmon, trout, and char (Isaak et al. 2015;Lynch et al. 2016;Wenger et al. 2011), but most assessments focus on southern subsets of species' ranges in mid-latitude areas where the preponderance of researchers also live. The vulnerability of these populations has been confirmed by observations showing range contractions and phenological adjustments in populations that are attempting to avoid unsuitably warm conditions (Lynch et al. 2016;LeMoine et al. 2020). The importance of temperature in driving these changes highlights the critical nature of maintaining or restoring natural thermal regimes for preserving aquatic biota (Olden and Naiman 2009). Doing so, however, requires a better understanding and description of thermal regimes, especially from high latitude areas where temperature records have traditionally been lacking. The advent of inexpensive, miniature temperature sensors and robust data collection protocols over the last decade (Stamp et al. 2014) have reduced this constraint, and broad comparisons with a focus on the ecologically relevant aspects of stream thermal regimes are becoming possible.
Bull trout (Salvelinus confluentus) is a char species that has experienced significant population declines throughout much of its range in western North America. Populations in the USA and Canada have protected status under the Endangered Species Act and Species at Risk Act, respectively (USFWS 1999;COSEWIC 2012), and the species is classified as being highly sensitive to climate change (Rieman et al. 2007). Bull trout typically occupies montane watersheds and is patchily distributed in headwater streams that are most conducive to adult spawning and subsequent survival of natal and juvenile life stages (Rieman and McIntyre 1995;Isaak et al. 2017b;Mochnacz et al. 2021). The occurrence of bull trout populations has been strongly linked with the distribution of the coldest water temperatures within landscapes (Benjamin et al. 2016;Isaak et al. 2017b;Kovach et al. 2017). However, this view is derived almost entirely from research conducted at the southern extent of the species' range where the primary focus has been on limitations associated with maximum summer stream temperatures (Isaak et al. 2015;Kovach et al. 2017). Maximum temperatures, however, represent only one component of annual thermal regimes, and it is well documented that there is substantial spatial-temporal heterogeneity in thermal regimes across stream networks (Steel et al. 2016;Isaak et al. 2020). One element that has not been examined for bull trout is whether streams near the northern geographic range extent may be cold-limiting. For example, studies have shown two ways in which cold stream temperatures act to limit fish survival. First, some streams do not provide enough thermal units throughout the growing season for juveniles to attain sufficient size and lipid stores to survive throughout the winter (Finstad et al. 2004b;Coleman and Fausch 2007a;Berg et al. 2009). Secondly, some streams may freeze completely during intense winters or experience high variability in habitat availability during this period and therefore provide little viable year-round habitat to support populations (Cunjak 1988;Cunjak et al. 1998). Here, we explore these issues using full-year temperature data records compiled from three montane areas distributed across the range-wide distribution of bull trout in North America. We use these datasets to describe the annual thermal regimes that these populations experience and gain a broader understanding of which elements might act as factors limiting distributions (Fig. 1). Our objectives were to (1) compare annual temperature records and juvenile distribution data from representative streams across the bull trout range and (2) compare biologically relevant aspects of thermal regimes across the range.

Study areas
The study area encompasses three montane stream networks situated in the southern (45° N, 124° W), central (55° N, 124° W), and northern (61° N, 124° ~ 2000 km Sources: Esri, GEBCO, NOAA, National Geographic, DeLorme, NAVTEQ, Geonames.org. This map was created using ArcGIS® software by Esri. ArcGIS® and ArcMap TM are the intellectual property of Esri and are used herein under licence. Copyright © Esri. All rights reserved. Fig. 1 Locations of the three sampling areas from Idaho (ID, red), Alberta (AB, orange), and Northwest Territories (NT, blue) where stream temperature records (n = 15/site) were used for this study. The black line depicts the approximate bull trout distribution W) regions of the bull trout range, spanning approximately 2000 km (Fig. 1). The three areas have topographically complex terrain but experience different climatic conditions. The southern area (Idaho, USA) has cold, wet winters, and dry, hot summers, whereas the central (Alberta, Canada) and northern (Northwest Territories, Canada) areas experience colder and longer winters, and shorter, hot, dry summers (Holland and Coen 1983;Halliwell and Catto 2003;Isaak et al. 2018; see differences in air temperatures, Table 1). The dominant vegetation in all three areas is mixed coniferous forests in higher elevation areas and shrubs, grasses, and willows in lower elevation areas, but the overall density of vegetation is far less in the northern area than the two others. The geology in the south consists of resistant granites and volcanics (Isaak et al. 2018), whereas the central and northern locations are composed of limestone, dolomite, shale mantled by till, and sandy fluvioglacial drift (Holland and Coen 1983;Halliwell and Catto 2003). Unpaved roads/ trails are present in all three areas but are least extensive in the north. All three stream networks have self-sustaining, healthy bull trout populations in multiple streams (M. Taylor unpublished data; Isaak et al. 2015;Mochnacz et al. 2021). The southern area has the most diverse fish assemblage with 12 species (Isaak et al. 2017a), followed by the central and northern areas which each have four species (Schindler 2000;Babaluk et al. 2015).

Stream temperature datasets
Hourly temperature data were collected from 15 sites across bull trout streams of similar size and gradient in each of the three areas using Tidbit temperature sensors (Onset Computer Corporation, Pocasset, Massachusetts, USA; Table 1). These sensors have measurement accuracies of 0.21 °C and resolutions of 0.02 °C. Most site records had temperature recordings on at least 60% of the days (but average completeness of most records was ≥ 80%) during a 3-year period. Records from northern and southern areas covered the period of 1 August 2013 to 31 August 2016, whereas the central area data ran from 1 July 2016 to 26 July 2019. Although it would have been desirable to have monitoring sites that ran concurrently over a longer time span to inform thermal regime research (Jones and Schmidt 2018), 2 or 3 years of data have been shown sufficient for representing many key aspects of thermal regimes . In areas where site monitoring records were spatially dense at the stream scale, only records from streams that were ≥ 2.5 km apart (with one exception in the northern area; Fig. S1) were selected to achieve spatial balance and minimize the probability of spatial dependency that often occurs in temperature records that are close to one another (Isaak et al. 2010). Because some study streams freeze during winter, we only used records with monthly mean daily water temperatures that were warmer than − 1.2 °C. This threshold was based on both experimental and field data, which show that supercooling temperatures in rivers range between 0.07 and − 1.0 °C (Devik 1949;Nafziger et al. 2013) and, because the accuracy of our temperature sensors was ± 0.21 °C, we assumed that any values < − 1.2 °C were not from flowing waters.
In some instances, records were missing daily stream temperature data because loggers were lost due to high flows or removed from sites earlier than anticipated for logistical reasons. Missing daily stream temperature values were imputed using the missMDA package (Josse and Husson 2016) in R (R Development Core Team 2018). This technique uses correlations among sites in time-series records to accurately estimate missing values by first applying standard principal components analysis to the incomplete dataset where missing values have been replaced with column means. Data are then reconstructed from the principal components and the initial analysis step repeated but with missing values replaced using estimates from the reconstructed data (Josse and Husson 2016). Similar to what others have shown (Isaak et al. 2018;Johnson et al. 2021), we found that the quality of imputed data based on temporal covariation from nearby stream sites was good, and all of the correlations between daily observed records and predictions from the imputation were high (r ≥ 0.99). After imputation, all stream site records consisted of 1127 mean daily temperature records across 3 years. To provide measures of the climatic variability that streams experienced, mean daily air temperature data were downloaded from local monitoring stations within 100 km of sites in each area (Idaho: Cooperative Observer Network, https:// www. ncdc. noaa. gov/ data-access, last access, 01 July 2020; Canada: Environment Canada, https:// clima te. weath er. gc. ca/ histo rical_ data/ search_ histo ric_ data_e. html, last access, 01 July 2020) (Environment Canada 2020; NOAA 2020). Mohseni et al. (1998) has shown that air temperature stations up to 250 km from stream monitoring sites can be useful in this regard.

Bull trout thermal metrics
Metrics were calculated to describe the thermal regimes of streams from each region based on magnitude, variability, and timing (Table 2). Although dozens of metrics are available, many are strongly correlated and redundant, so we focused on a smaller subset that was most relevant to key elements of bull trout biology (Fig. S2) (Chu et al. 2010;Arismendi et al. 2013;Isaak et al. 2018). For example, August mean temperature and winter mean temperature were used as magnitude metrics. The former has been used to define the thermal niche that bull trout occupy in southern latitudes (Dunham et al. 2003;Isaak et al. 2015Isaak et al. , 2017b, but August mean temperature has not been reported from mid-and high-latitude streams across the range. Juveniles and adults in the southern and central regions of the range prefer to occupy streams with mean summer stream temperatures < 11 °C and, although adults can survive in warmer water (> 16 °C; Isaak et al. 2015;Parkinson et al. 2016), growth is typically poor at these temperatures (Selong et al. 2001). Accounts of winter stream thermal regimes that bull trout experience are rare, but the winter season is considered by many as a survival bottleneck for freshwater salmonids as it influences development and growth of eggs, timing of hatching and juvenile emergence (i.e., free swimming fish), and survival of adults and juveniles (Cunjak 1988;Shuter et al. 2012). During winter, early-life stages of salmonids often experience high mortality due to exhaustion of energy reserves (Finstad et al. 2004b).
Temperature is also an important determinant of the phenology of bull trout life history events, such as hatching and emergence timing, as well as growth rates that may reflect either local phenotypic adjustments, or adaptations to unique thermal regimes that act to maximize individual fitness (Sparks et al. 2017;Austin et al. 2019;Campbell et al. 2019). Moreover, the timing of phenological events varies across the species' range because both climate and the length of development and growing seasons differ across latitude-i.e., winter begins earlier and is longer in the north than it is in the south (Reist et al. 2006b;Shuter et al. 2012). Once eggs are deposited in the gravel, the rate and magnitude of thermal units that accumulate throughout the incubation and initial growth period determine when hatching occurs and juveniles emerge-i.e., where the yolk sac is absorbed and juveniles can swim on their own (Neuheimer and Taggart 2007; Fuiman and Werner 2009). Accumulated thermal units (ATU) are often used to quantify how temperature influences developmental rates, and is quantified by calculating cumulative thermal units over time, where 1 °C for 24 h = 1 thermal unit (Neuheimer and Taggart 2007). Laboratory and field studies on bull trout show that individuals require approximately 800 ATU from egg deposition to 50% juvenile emergence (Gould 1987;Bowerman et al. 2014), and Bebak et al. (2000) show that survival of juvenile Arctic char (Salvelinus alpinus) is highest when fish experience an additional 800 ATU after emergence. Together, these data suggest that chars require 1600 ATU from egg deposition to the onset of winter. Based on this assumption, we used 1600 ATU as a guideline for estimating the minimum number of thermal units that bull trout require to survive the winter.
ATU are also important for understanding differences in fish growth across latitudinal gradients and during important life stages, such as the summerearly fall (Coleman and Fausch 2007b; Neuheimer Table 2 Metrics used to compare different elements of thermal regimes across montane river networks spanning the distributional range of bull trout a Cumulative thermal unit calculations started and ended on different dates to coincide with spawning timing windows (October 01, September 15, September 01) and estimated emergence dates (March 07, June 25, May 29) associated with southern, central, and northern watersheds, respectively b The start of the growing period coincided with the estimated emergence date (March 22, May 20, April 14) and the end of the growing season was defined by the date when water temperatures were ≤ 2°C on two consecutive days (November 01, October 15). These dates correspond to southern, central, and northern watersheds, respectively. Note that the end date was the same for the central and northern watersheds c These are general start-end dates for the growing period. For calculations, the start of the growing period coincided with the estimated emergence date (March 07, June 25, May 29) and the end of the growing season was defined by the date when water temperatures were ≤ 2°C on seven consecutive days (November 01, October 15). These dates correspond to southern, central, and northern watersheds, respectively. Note that the end of the growning season was the same for the central and northern watersheds (October 15) d Calculations for analyses during this period started and ended on different dates to coincide with spawning timing windows and estimated emergence dates specified above where E is an effective value (range of 0-1) describing the relative daily contribution to development; log e a = 5.59 and b = 0.126 are model coefficients, based on thermal relationships and emergence timing for bull trout (Austin et al. 2019); and T is the daily mean water temperature on each day of incubation. Because fish eggs accumulate E over the course of the incubation period, the model predicts 50% emergence when the sum of E = 1. The initial model developed by Beacham and Murray (1990) requires an estimate of mean water temperature during incubation but, when this is unknown, the Sparks et al. (2019) model allows one to predict hatch timing by using daily mean water temperature and each day's respective contribution towards development. The growing period was an estimate of emergence date for each respective area through to the onset of winter, defined as the period when stream temperatures were ≤ 1.5 °C for 7 consecutive days (1 November in the south and (1) 15 October in the central and northern areas). We did not choose these timing windows to predict exactly when timing of key events happens, but rather to compare differences in the estimated hatch dates and magnitude of thermal units accumulated across developmental periods and among regions.
Finally, we calculated thermal sensitivity, which is a measure quantifying how streams respond to air temperature variation and an important metric for understanding the thermal stability of streams (Snyder et al. 2015;Bolduc and Lamoureux 2018). Streams with low to moderate thermal sensitivity (0.10-0.45) are defined as thermally resilient because they have more stable stream temperatures throughout the year due to minimal effects imposed by changes in air temperature. Conversely, streams with thermal sensitivities greater than 0.55 are defined as thermally reactive, where the thermal response to air temperature is more intense (Kelleher et al. 2012;Mayer 2012;Piccolroaz et al. 2016). Additionally, perennial groundwater is an underlying mechanism driving low thermal sensitivity in streams, and areas associated with perennial groundwater constitutes high-quality spawning and rearing habitat for bull trout (Baxter and McPhail 1999;Baxter and Hauer 2000). To identify variations in water temperature relative to changes in air temperature, we calculated the thermal sensitivity of all sites across years in each area. Thermal sensitivity was expressed as the slope of the linear regression relationship between weekly water temperature (T w ) and weekly air temperature (T w ) records across a given year. Weekly time steps were used because they typically provide more precise thermal sensitivity relationships than daily time steps, and full-year records were used to capture the broadest range of temporal variability in this relationship (Kelleher et al. 2012). Negative air temperatures were not included in this calculation because linear relationships below this threshold poorly predict stream temperatures and less accurately represent the influence of groundwater buffering on thermal sensitivity (Morrill et al. 2005;Kelleher et al. 2012;Mayer 2012).

Data analyses
We used a combination of generalized linear (binomial) and linear mixed models to test for differences in thermal regime metrics among regions because these models can handle unbalanced data and account for differences associated with random effects among datasets (Zuur et al. 2009). Year was set as a fixed effect in all models to account for potential differences between years because data were not collected across the same periods in all regions. Site was specified as a random intercept to account for repeated measures across years for models examining winter stream temperature, thermal sensitivity, emergence date, and accumulated thermal units. Variation in winter mean temperature, thermal sensitivity index, estimated egg hatch date, and ATU were examined using linear mixed models with location (i.e., region-as factor), elevation, their interaction (location × elevation), and year (factor) included as fixed effects. Seasonal differences in ATU were compared among regions across incubation, growing, and annual periods. We modeled mean August stream temperature using a combination of linear and generalized linear models (see details below). Model selection was performed using backwards stepwise regression with marginal F tests. No data transformations were performed for linear mixed models but, for the generalized linear model, continuous variables were standardized to a mean of 0 and a standard deviation of 1. The assumptions of models were tested following the methods of Zuur and Leno (2016). Tukey pairwise post hoc multiple comparisons tests were used for among-region comparisons when location was found to be an influential fixed effect. For linear mixed models, marginal (R 2 m) and conditional (R 2 c) coefficients of determinations were used to quantify the proportion of variance explained by fixed factors and fixed and random factors, respectively. All analyses and figures were completed in R (R Development Core Team 2018) and significance was assessed at the 0.05 level. Analyses were performed using the following R packages: Tukey tests with lsmeans (Lenth 2016), LMM with nlme (Pinheiro et al. 2016), and R 2 m and R 2 c with MuMIn (Barton 2016).
Because juvenile bull trout distributional data were available across all three areas (ID, n = 180; AB, n = 185; NT, n = 415) and was collected using similar methods-i.e., two spatial/temporal replicates across sites allocated using a stratified random design (Isaak et al. 2017b;Mochnacz et al. 2021)we combined these data with mean August stream temperatures to portray the available and occupied thermal niche in each respective area. To provide a consistent means of comparison across the fish survey sites which lacked co-located temperature sensors, spatially explicit temperature models were built for streams in each of the three study areas that predicted mean August temperatures (Isaak et al. 2017a;Mochnacz 2021). By modeling August stream temperatures across these three areas, we were able to precisely predict temperatures (≤ 1.0 RMSPE) across a broader spatial scale than we could have using the limited number of full-year temperature records available in each area (n = 15). These data were analyzed as a two-step process. First, full-year temperature records (n = 15) were used to test for regional differences in mean August stream temperature (dependent variable), at sites occupied by bull trout, using a linear model with location (as factor), elevation, and mean August air temperature included as fixed effects, and site as a random effect. Second, the distributional datasets were input into a generalized linear model (GLM) to define the realized summer thermal niche that juvenile bull trout occupy in each respective region (south, central, north) as well as a global model to represent a range-wide thermal niche. For each GLM, juvenile presence-absence data (dependent variable) was regressed against mean August stream temperature. This model was fit with both linear and quadratic terms for mean August stream temperature because both relationships have been shown to explain variation in occupancy of stream-dwelling salmonids elsewhere (Isaak et al. 2017a). Thermal response curves were plotted as the probability of occupancy versus mean August stream temperatures across the range of values from the dataset. Modeled thermal response curves for each region were plotted together to visualize the degree of similarity in the realized summer thermal niche among populations. Because prevalence differed across each region and resulted in differences in peak probability of thermal response curves (i.e., height of response curves), probabilities were rescaled to the maximum value observed (0.80) for visual comparison. In the northern area, most streams freeze completely during the winter; therefore, full-year temperature records are sparse. Consequently, small sample sizes of other metrics associated with fullyear records (e.g., mean winter stream temperature, ATU) precluded integration with distributional data and modeling, as described above in step 2.
As latitude increases, winters become colder and longer, which translates into more extensive freezing of freshwater rivers in montane systems (Prowse et al. 2006a, b;Crites et al. 2020). However, at higher latitudes, perennial groundwater discharge in some streams attenuates the magnitude of freezing, creating areas that are ice free or do not completely freeze to the bottom (Utting et al. 2013;Crites et al. 2020). This phenomenon is reflected in the thermal sensitivity of streams, whereby as perennial groundwater contributions increase, thermal sensitivity declines (Kelleher et al. 2012;Bolduc and Lamoureux 2018;Hare et al. 2021). Given that perennial groundwater is the primary mechanism preventing shallow (< 1.5 m) streams from freezing at higher latitudes (Utting et al. 2013), it follows that streams which remain unfrozen in the winter should have lower thermal sensitivity. To investigate whether this relationship was present in our dataset, we randomly selected 80% of the available full-year temperature records from each region (ID, n = 314; AB, n = 84; NT, n = 68), and calculated thermal sensitivity. Sites were classified as frozen or unfrozen, based on mean daily winter temperatures being below or above the flowing water threshold of − 1.21 °C. Pairwise t-tests were used to determine if thermal sensitivity of frozen versus unfrozen sites differed within regions.

Results
A comparison of daily mean air temperatures to daily mean water temperatures showed that stream temperatures from all three regions exhibited a dampened response to air temperature fluctuations during the summer months (Fig. 2). Results of our linear mixed model support these trends and showed that the mean (± SD) thermal sensitivity was low to moderate for all three areas (ID: 0.43 ± 0.06; AB: 0.35 ± 0.09; NT: 0.35 ± 0.12) but differed between southern and central-northern regions (location: F 2,42 = 4.03, p = 0.02; Tukey test: ID-AB, Z = 0.08, p < 0.0001; Tukey test: ID-NT, Z = 0.09, p < 0.0001; Tukey test: AB-NT, Z = 0.005, p = 0.955; Fig. 3(A)) and varied across years, but the year effect size was relatively small (year: F 6,129 = 8.03, p < 0.0001, coefficient range 2014-2019 = − 0.004 to − 0.10, SE = 0.01-0.02). Elevation was removed from the model because it was not significant (F 2,39 = 1.10, p = 0.34). Location accounted for 20% of the variation in thermal sensitivity, whereas the random effect of site accounted for 65%.
Thermal response curves indicated that the realized thermal niches of these populations overlapped, but the northern population occupied a colder and narrower thermal niche than both the central and southern populations (Fig. 4A). The curves for the central and southern populations showed that occurrence probabilities peaked at 8.3 °C and 6.2 °C, Fig. 2 Time series of mean daily air (blue line) and stream (pink line) temperatures from Idaho (ID), Alberta (AB), and the Northwest Territories (NT). The water temperature values are from all sites in each respective area for a given year and the thick black line is the daily mean of all sites respectively, and both exhibited cold-and warmedge transition boundaries. The thermal curve for the southern population did not show a cold-edge boundary because sufficiently cold streams were not available to populations in this portion of the range, but this curve did display a warm-edge transition boundary at temperatures ≥ 11 °C (Fig. 4B).
Estimated number of days to emergence (mean ± SD) differed among regions with days to emergence in the south (157 ± 34.9) being fewer than the north (270 ± 7.0) and central (283 ± 5.9) areas (location: F 2,42 = 41.8, p = 0.0003; Table 3). Interestingly, the length of time to emergence was longest in the central region followed closely by the northern  (Fig. 5). ATU differed among regions across all developmental periods, but differences were smallest during the incubation period (Table 3). In the south, ATU during the annual and growing periods were more than double the central region and 1.8 to 2.1 times greater than the northern region. Conversely, differences in mean ATU between central and northern locations were much smaller (Table 3). Time series plots of ATU are illustrated in Fig. 5 for southern, central, and northern areas. These data show that in all years, ATU from all sites show an immediate increase during the initial part of the incubation period shortly after spawning, followed by a plateau during the winter months, an increase during the early spring, and then a more pronounced increase throughout the growing period. During the growing period, the northern and central areas showed a slower rate of increase in ATU than the southern location (Fig. 5).
Results of the linear mixed models for the incubation, growing, and annual development periods support the patterns shown in Fig. 5, as ATU differed between southern and central-northern locations during all periods (incubation: location-F 2,41 = 24.5, p < 0.001; growing: location-F 2,41 = 120.0, p < 0.0001; annual: location-F 2,41 = 120.0, p < 0.0001), but intercepts did not differ between the central and northern locations (incubation: Tukey-T 41 = − 0.12, p = 0.97; growing: Tukey, NT-AB-T 41 = 53, p = 0.44; annual: Tukey, NT-AB-T 41 = 99.8, p = 0.12). The effect of year was significant for all periods and elevation was significant only during the growing and annual periods (incubation: elevation-F 1,41 = 1.8, p = 0.19; year-F 4,86 = 15.1, p < 0.0001; growing: elevation-F 1,41 = 6.89, p = 0.01; year-F 3,72 = 45.3, p < 0.0001; annual: elevation-F 1,41 = 7.21, p = 0.01; year-F 3,72 = 46.0, p < 0.0001). Although year was a significant effect in the models, it had a relatively small effect on ATU across all three periods with model coefficients (± SE) ranging from − 0.003(0.004) to − 0.006(0.004). The interaction between location and elevation was not significant and was removed from all models. The random effect of repeated measures across sites accounted for 37%, 10%, and 11% of the variation in ATU during the incubation, growing, and annual periods, respectively, while the fixed effects accounted for 50%, 82%, and 83%, respectively. , and for all sites across the three areas (Global). The black arrow shows where peak occurrence occurs (8.4 °C) on the global response curve. Variation in juvenile occupancy was best described by a quadratic relationship with mean August stream temperature in northern and central areas and a negative linear relationship in the southern area (ID). The hatched blue line represents the theoretical curve extending beyond data used to build this model. Untransformed regression equations for each respective curve are shown. Models were built using standardized values for temperature terms and values were back transformed for the probability plot. Aug_Temp, is the mean August stream temperature and Aug_Temp 2 , is mean August stream temperature raised to the second power for use in the quadratic model

Discussion
We show that bull trout occupy a relatively cold thermal niche, relative to other cold-water fishes, across a broad latitudinal gradient (mean annual stream temperatures across the three study areas ranged from 1 to 5 °C). Further results of a range-wide species thermal response curve model predicted site occupancy peaked at 8.4 °C mean August stream temperatures. In the southern portion of the species' range, however, bull trout appeared to be regulated by warmer temperatures and occurrence probabilities declined to < 0.1 once summer temperatures exceeded 13 °C (Isaak et al. 2017b;Selong et al. 2001). Temperatures that warm were rarely observed in the central and northern study streams where bull trout appeared instead to be regulated by unsuitably cold temperatures and a lack of thermally suitable habitat in many areas throughout the year. This finding further reinforces that some populations in the southern distributional range extent occupy streams that are at the species' warm-edge range boundary where site extirpation is more likely as streams continue to warm (Eby et al. 2014;LeMoine et al. 2020). Regardless of the mechanism controlling population persistence (i.e., warm-vs. cold-limiting habitat), all the streams we examined exhibited high thermal stability, which suggests this is a key abiotic property of core habitat across the species' range (Luce et al. 2014;Isaak et al. 2016). However, in northern areas, this thermal property acts to stabilize stream temperatures to prevent freezing during the winter, whereas in the south, high thermal stability mediates air temperature warming effects in the summer because streams rarely freeze in the winter (Isaak et al. 2016).
These results highlight an important dichotomy in the prevailing thermal mechanisms governing the distributions of bull trout populations at southernnorthern range limits and may have important implications for the species as climate warms. In the south, populations inhabiting streams that warm beyond acceptable upper thermal limits will need to disperse upstream to colder, higher elevation streams, or adjust or adapt to occupy warmer thermal regimes. In the Bitteroot basin in Montana, both Eby et al. (2014) and LeMoine et al. (2020) report range contraction of bull trout through abandonment of previously occupied sites that have warmed, suggesting either dispersal or extirpation as possible mechanisms. In most instances, fish will probably move from warmer downstream habitat upstream to colder habitats at higher elevations; however, in some systems, these habitats may be at their carrying capacity or occupied by competitors (e.g., brook trout Salvelinus fontinalis;  Paul and Post, 2001). In the north, our results suggest that the distribution of populations is currently limited by streams being too cold in both the summer and winter (i.e., streams freeze to the bottom). However, as streams warm in these very cold systems, there is potential for expansion of thermally suitable habitat whereby streams that are currently too cold to support populations warm enough to provide new colonization opportunities. In addition, productivity may increase as the severity (e.g., warmer, less ice cover) and length of winters decrease (Prowse et al. 2006a, b), resulting in greater growth potential and better recruitment via higher survival throughout the winter. Campana et al. (2020) show that as climate warms in the north, lake trout (Salvelinus namaycush) populations will realize gains in productivity and experience potential range expansion by colonization of new habitats. Species geographic ranges are an artifact of physiological thermal safety margins, which are typically smaller at upper thermal limits (Sunday et al. 2014;Sandblom et al. 2016). Given this and our findings, we would expect that, as streams warm, bull trout will continue to experience range contraction in the south and potential expansion in the north. However, the overall net effect on the species' distribution will probably be driven by the proportion of cold-limiting areas versus warm-limiting areas found across the geographic range.

Potential for local adaptation
To our knowledge, the northern bull trout populations from this study occupy the coldest and narrowest summer thermal niche on record for the species, and represents a downward shift from thermal niches reported elsewhere across the range (Benjamin et al. 2016;Isaak et al. 2017b). However, this is not particularly surprising given this is the most northerly situated area for which thermal regimes have been described, and it is near the species' northernmost range boundary . The summer stream temperature of 6.2 °C, where occurrence peaked, aligns well with the thermal optimum (≤ 10 °C) for fishes in the Arctic thermal guild . The colder streams that the northern populations occupy is partially a reflection of a colder summer thermal niche that is available compared to others we examined. This implies that northern populations may be forced to occupy a colder summer thermal niche compared to those in central and southern areas. However, it is interesting that despite northern populations having access to relatively warm streams (7.0-9.0 °C), occurrence peaks at a much lower temperature of 6.2 °C, which is 2 °C below the peak temperature of the global response curve. This suggests that northern populations have either adapted to occupy a colder thermal niche than those found further south, possess the phenotypic plasticity to acclimate to a colder thermal regime, or occupancy of this cold thermal niche is driven by other factors (e.g., predation by older bull trout, competition with other species, access to resources). Local adaptation to different environmental conditions is common in salmonids (Eliason et al. 2011;Narum et al. 2013) and arises more often as geographic distance between population increases (Fraser et al. 2011). However, chars are also renowned for exhibiting exceptional phenotypic plasticity; therefore, it would not be surprising if acclimation was the prevailing mechanism (Klemetsen 2010;Chavarie et al. 2016).

Thermal regime contrasts
Thermal regimes across all three areas exhibit similar insensitivity to seasonal fluctuations in air temperature and highlights a common element of thermal stability in core habitat that these populations require. Similar stability in thermal regimes has been reported for bull trout in other regions (Isaak et al. 2016) but, to our knowledge, this is the first account of disparate populations experiencing similar thermal stability across a wide geographic area. Occupying habitat with stable thermal regimes is important for persistence of coldwater salmonids because streams with these properties are projected to mediate effects of climate warming (Lisi et al. 2015;Snyder et al. 2015;Mauger et al. 2017) and could promote positive growth during winter (French et al. 2017). The average thermal sensitivity in the southern area was slightly higher than that in the central and northern areas, but it was still well below the 0.55 threshold differentiating insensitive versus sensitive streams (Kelleher et al. 2012;Piccolroaz et al. 2016). Differences in thermal sensitivity could be related to the characteristics of these streams (i.e., elevation, sun angle, riparian canopy) and the prominent mechanism driving thermal sensitivity. Very stable water temperatures throughout the summer in the central and northern areas produced flatter air-water temperature hysteresis cycles (i.e., a cyclical lag between input and output) and much lower thermal sensitivities that are consistent with systems where perennial groundwater is the prevalent underlying mechanism (Kelleher et al. 2012). Conversely, the streams in the south exhibit slightly broader and higher air-water temperature hysteresis cycles that are reflective of systems where thermal sensitivity is driven more by a combination of precipitation, snow melt, riparian shading, and, to a lesser extent, groundwater (Kelleher et al. 2012;Lisi et al. 2015). Although it is not possible to define the exact mechanism(s) responsible for the patterns we report, a first-order understanding of the most probable mechanism(s) will be useful for predicting how climate warming may affect these systems differently. For example, Lisi et al. (2015) show that, as climate warms, the heterogeneity pattern of thermal sensitivity in high elevation snow-melt stream systems may be lost, resulting in more uniform warming across all streams.
Our analysis also demonstrates that in the north, only streams which do not freeze completely during the winter have low thermal sensitivities, suggesting groundwater as the primary mechanism maintaining warm enough environments for fish to survive throughout the winter. It follows that the distribution of suitable spawning and rearing habitat in this area appears to be driven by a balance between streams being warm, and presumably productive enough in the summer, but do not freeze completely during the winter due to perennial groundwater input. Similarly, the distribution of Dolly Varden (Salvelinus malma) spawning and overwintering habitat is also associated with prevalence of groundwater (Dunmall et al. 2016), and suggests the pattern we report reflects an important mechanism governing suitability of riverine char habitat north of the 60th parallel. Further, Crites et al. (2020) mapped the distribution of perennial groundwater sources in the western Arctic, and these align well with the broader distributions of riverine chars in this region, further reinforcing the importance of groundwater prevalence for sustaining riverine char populations at higher latitudes .
The finding that central and northern populations experience colder thermal regimes than the southern population also represents a divergence in patterns reported by others (Benjamin et al. 2016), and highlights a potential recruitment bottleneck associated with growth potential and emergence phenology. Our results show how populations exposed to colder thermal regimes will emerge later and have access to a shorter open water season, thereby reducing growth potential before the onset of winter. Variability in emergence timing across streams with different thermal regimes was also reported by Austin et al. (2019) for bull trout in Washington State and Sparks et al. (2019) for sockeye salmon (Oncorhynchus nerka) in Alaska. However, our results are contrary to those of Campbell et al. (2019), who reported high synchrony in hatch dates for coho salmon (O. kisutch) across streams with vastly different thermal regimes. This discrepancy could be related to an indirect interaction (i.e., not statistical) between latitude and winter thermal regime that are present in our results but not considered by Campbell et al. (2019), a missing co-variate, or interspecific differences in life history strategies that represent local adaptations to unique thermal regimes (Shuter et al. 2012). Determining hatch and emergence dates using otoliths would refine our understanding of hatch phenology for populations from these areas and validate emergence model estimates (Fuiman and Werner 2009;Campbell et al. 2019). In addition, the fact that interannual variation (i.e., year effect) had a relatively small influence on most metrics we examined suggests some level of stability in the thermal regimes we report and may also reflect a broader temporal signature in these systems. It is worth noting, for models of thermal sensitivity and winter stream temperatures, the random effects explained more variation than fixed effects, and may suggest that these two metrics vary more within regions (i.e., site effects) rather than among them (i.e., distributional effects). However, high sitelevel variability could also be present because few fixed effects were included in our models and important covariate controls associated with subsurface hydrology, and the amount of groundwater inputs, are not well represented by things estimable from digital elevation models, such as stream elevation.
The divergence in growth potential through exposure to different annual thermal regimes presents some interesting hypotheses that warrant consideration. Results suggest that limits on habitat carrying capacity currently force northern populations to occupy streams near cold-edge thermal boundaries. So how do these fish manage to survive in such cold and presumably unproductive environments? Northern populations may possess a compensatory mechanism to grow more efficiently in colder environments during shorter growing seasons (i.e., countergradient variation) (Conover and Schultz 1995). All of the streams in the northern area, where we describe fullyear thermal regimes, have abundant, self-sustaining bull trout populations ), yet several of these streams accumulate very few thermal units throughout the annual growing season (i.e., ≤ 800). However, we would not expect populations to occupy these streams if they could not survive under these growing conditions; therefore, they appear to have adjusted (i.e., genetically or via phenotypic plasticity/acclimation) to be successful in these colder thermal regimes. Such adaptations to cold winter thermal regimes have been observed in other salmonids (Finstad et al. 2004a). Furthermore, the acclimation time to readjust efficiency of physiological energy pathways is reduced as a species' fundamental thermal niche gets closer to ambient temperature (Pörtner and Farrell 2008). In this context, it is worthwhile noting that, in the northern area, the difference between the mean annual stream temperature (2.4 °C) and optimal temperature where summer occurrence probability peaks (6.2 °C) is approximately 3.8 °C.
Others report similar compensatory growth in Arctic char (Chavarie et al. 2010;Sinnatamby et al. 2015), and we think this observation warrants further investigation in bull trout.

Conclusion
Examination of full-year thermal regimes across a broad latitudinal gradient reinforces what others report from the southern range extent, which is that bull trout populations occupy a cold, narrow thermal niche relative to other cold-water salmonids (Benjamin et al. 2016;Isaak et al. 2017b). However, we show how divergent elements of annual thermal regimes appear to limit populations at southern and northern range extents (i.e., warm-versus coldlimiting streams). Consequently, as air temperatures continue to rise, managers will need to consider how these two areas differ if they wish to manage thermal regimes effectively to ensure long-term survival of core populations (Isaak et al. 2016). Even if managers carefully monitor and manage thermal regimes (e.g., mitigate flow disruption and maintain riparian zones), it seems inevitable that streams will continue to warm (Isaak et al. 2017a), which makes it important to understand the adaptive capacity that populations possess to adjust to potential changes in stream temperature. Given that local adaptations to unique thermal regimes occur in other salmonids (Eliason et al. 2011;Sparks et al. 2019), we view this as an important avenue of future research for conserving bull trout populations. A key element in this process is broadening understanding of thermal regimes that bull trout populations experience across a wide latitudinal gradient by further expanding stream temperature monitoring networks in both central and northern Canada. Doing so will determine if the patterns we report are a local area effect or are reflective of a broader classification of cold-water thermal regimes in montane areas across the geographic range.