Spatiotemporal patterns of emergence phenology reveal complex species‐specific responses to temperature in aquatic insects

Climate change is broadly affecting phenology, but species‐specific phenological response to temperature is not well understood. In streams, insect emergence has important ecosystem‐level consequences because emergent adults link aquatic and terrestrial food webs. We quantified emergence timing and duration (within‐population synchronicity) of insects among streams along a spatiotemporal gradient of mean water temperature in a montane basin to assess the sensitivity of these phenological traits to heat accumulation from mid‐winter through spring emergence periods.


| INTRODUC TI ON
The consistent warming trend currently affecting natural systems has stimulated strong research interest in understanding thermal controls on life history traits, including the timing (phenology) and duration of life cycle events within and among populations (Parmesan & Yohe, 2003;Schwartz et al., 2012). At the population level, phenological concerns include the possibility of shifts in life history events to suboptimal timing for resource availability (Thackeray et al., 2016;Van Dyck et al., 2015), in addition to the more extreme possibility that increasing temperatures push metabolic rates and/or by-products beyond tolerable ranges (e.g. Chou et al., 2018;Galen & Stanton, 1991;Shah et al., 2017). Consequences may also occur at community and ecosystem levels, including development of temporal mismatches between trophic groups in food webs (Edwards & Richardson, 2004;Renner & Zohner, 2018;van Asch & Visser, 2007) or between plants and pollinators (Forrest & Thomson, 2011;Hegland et al., 2009).
Climate-associated phenological shifts are strong possibilities for ectotherms, plants, and other organisms whose biological rates are directly influenced by ambient temperature (Cohen et al., 2018;Van Dyck et al., 2015). For example, simple measures of accumulated "heat units," typically measured as degree-days, from the onset of growing seasons in the temperate zone are often broadly predictive of the timing of insect life cycle events such as pupation and emergence (Bartomeus et al., 2011), salmon fry emergence (Kaylor et al., 2021), and budbreak and flowering in plants (Ward et al., 2018). However, many ectotherms and plants also have periods of diapause or quiescence, during which the biological clock associated with accumulation of heat units is temporarily paused. During such periods, growth and development stop or slow substantially regardless of temperature, making thermal effects on life cycle phenology more difficult to predict. Animals for which we currently have the best physiological understanding of the interplay of thermal regime, diapause and phenological traits are mostly terrestrial insects of economic importance, typically pests or pollinators. Common environmental cues for both initiating and ending periods of diapause in insects include temperature thresholds (high or low) and photoperiod (Bean et al., 2007;Bentz et al., 2014;Forrest, 2016). Hence, the combination of degree-day accumulation rate during periods of growth and activity, and the timing and length of diapause, together are major variables useful for predicting phenological response to changes in temperature. Cues and responses can vary substantially among species, however, such that an understanding of phenological response to temperature in one species often does not translate to another (Li et al., 2011;Marshall et al., 2020).
In aquatic insects, the timing and duration of life cycle events such as adult emergence have well-studied ecological upshots including trophic subsidies from aquatic to terrestrial systems (Anderson et al., 2019;Nakano & Murakami, 2001;Schindler & Smits, 2017). Although emergence phenology has widely appreciated consequences at the ecosystem level, the underlying physiology associated with phenological response to temperature is not well understood in aquatic insects (Woods et al., 2021). General patterns of faster development rates and earlier emergence timing with warmer conditions have been documented either in natural systems within a limited spatial extent (e.g. Baranov et al., 2020;DeWalt & Stewart, 1995;Finn & Poff, 2008;Peckarsky et al., 2000) or in controlled experiments in simplified systems (e.g. Harper & Peckarsky, 2006;Li et al., 2011;Sweeney et al., 2018). As such, general predictions associated with warming in aquatic ecosystems typically assume a simple and direct link between warmer conditions and earlier timing of life cycle events (Durance & Ormerod, 2007, Finn et al., 2014, Larsen et al., 2016McCulloch & Waters, 2018).
However, as noted for the better studied terrestrial insects, aquatic insect phenological response to thermal regime is likely not so simple. For example, Newbold et al. (1994) used a model to show that both high-and low-temperature diapause thresholds were necessary to explain the maintenance of univoltine life cycles and nearly synchronous emergence in several mayflies along a lengthy thermal gradient in streams spanning 16° latitude. Although this model was able to explain phenological outcomes with temperature patterns as the sole input, thermal regime probably interacts with photoperiod or other environmental cues to affect phenological traits in some stream insects (Cabanita & Atkinson, 2006;Mendez & Resh, 2008).
Even in some mayfly species modelled by Newbold et al. (1994), a high-temperature diapause could not be induced under laboratorycontrolled thermal regimes (Funk et al., 2019). Funk et al. (2019) suggested that photoperiod likely interacts with water temperature to cue high-temperature diapause in natural streams.
Thermal regime might also influence the duration of emergence periods of some aquatic insects, with longer emergence periods associated with populations experiencing warmer conditions during growth and development (Richardson & Clifford, 1986;Li et al., 2011;Cheney et al., 2019;but see DeWalt & Stewart, 1995). A longer emergence period at the population level indicates increasing life cycle asynchrony and may be a response to unstable environments (Calabrese & Fagan, 2004;Li et al., 2011;Lytle, 2001). Alternatively, life cycle stages. We hypothesize a trade-off between complex phenological response that synchronizes emergence among heterogeneous sites and other traits such as adult longevity and dispersal capacity.

K E Y W O R D S
aquatic insects, headwater streams, life history traits, phenology, temperature in early-emerging insects, longer emergence periods may be an adaptive response that extends adult activity periods in years when a longer growing season is likely (Buckley et al., 2015).
We still have limited understanding of the direct impact of thermal regime on phenological traits in aquatic insects, but streams are ideal systems for detailed investigation. Stream networks are natural mosaics of habitat conditions, including temperature (Ebersole et al., 2003;Johnson, 2004;Uno & Power, 2015;Ward & Stanford, 1982). Mountain stream networks in particular contain a mosaic of local, reach-scale thermal environments thanks to elevational gradients, shading, substrate, and heterogeneity of hydrological sources to headwaters (Hotaling et al., 2017;Rupp et al., 2020;Ward, 1994).
Furthermore, the topographical complexity of mountain landscapes coupled with interannual heterogeneity of weather results in thermal conditions in one year that are rarely predictive of thermal conditions in the next (Bales et al., 2006;Reiter et al., 2020). As such, mountain stream networks should be useful settings for documenting phenological response along thermal gradients, both spatially and temporally, while controlling for photoperiod and other local/ regional-scale environmental cues.
Here, we documented emergence timing and duration of four spring-emerging aquatic insect species from six headwater streams within a fifth-order stream network of the Cascade Range, USA, for 6 consecutive years. High spatial and temporal variability in thermal conditions allowed us to address the overarching question: Is emergence phenology of common stream insects predictable based strictly on differences in heat accumulation from mid-winter, across both space and time? We focused efforts on spring-emerging insects because environmental conditions in the study region are predicted to be highly responsive to changing climate during this season (Cook et al., 2012;Safeeq et al., 2013;Wu et al., 2012). Additionally, springemerging insects are likely to be directly responsive to thermal influences during the transition between winter (when diapause or quiescence is most likely to occur) and the peak growing season. This study was conducted in coordination with other research into spring phenological events across trophic levels in the same basin, including plant budbreak (Ward et al., 2018), terrestrial insect activity (Schmidt, 2019), and bird arrival and dispersal .
Because we had minimal understanding of the potential physiological drivers of emergence phenology for our study taxa, for each species we tested two basic and plausible hypotheses addressing the influence of water temperature on (1) timing of emergence and (2) duration of the emergence period ( Figure 1). H1: Adult emergence timing responds directly to water temperature, such that warmer streams in the network and warmer years yield predictably earlier emergence. H2: Duration of emergence period is longer in warmer streams in the network and in warmer years. For both hypotheses, we quantified "warmth" as heat accumulation in degree-days (DD) from a mid-winter period of reduced ecosystem-scale metabolic rates through the early summer growing season.
F I G U R E 1 Hypotheses and our definitions of two traits: peak emergence timing (t) and duration (d) of emergence period. Blue curves and variables indicate cooler (c) streams or years; orange curves and variables indicate warmer (w) streams or years. Panels (a and b) represent Hypothesis 1 with distributions presented in two ways: (a) absolute, and (b) cumulative emergence through time. Vertical broken lines illustrate our definition of "peak" emergence timing as the mean ordinal date from fitted Gaussian curves. Panels (c and d) represent Hypothesis 2, again with distributions presented as absolute (c) and cumulative (d) emergence. Here, both curves in each panel have the same peak emergence timing, but emergence duration as the number of days between the 5th and 95th percentiles of the Gaussian models is longer in the warmer than in the cooler stream or year. Vertical dotted lines represent 5th and 95th percentiles for the cool stream/year; vertical double-lines represent 5th and 95th percentiles for the warm stream/year 2 | ME THODS

| Study sites
This study took place in six headwater streams of the Lookout Creek basin in the H.J Andrews Experimental Forest (Oregon, USA; Figure 2). Dominant overstorey vegetation is Douglas fir (Pseudotsuga menziesii), with riparian woody vegetation of Douglas fir, western red cedar (Thuja plicata), red alder (Alnus rubra) and big leaf maple (Acer macrophyllum) (Frady et al., 2007). Each of the six focal streams was first or second order at the 1:24,000 scale, and streams spanned a range of elevations (457-1010 m) and two major categories of surrounding forest age ( Table 1). Annual precipitation in the Lookout Creek basin averages 2285 mm, with the majority occurring as winter rainfall (October through April). Higher-elevation locations in the basin experience winter snowfall that typically melts during the season, and lower elevations experience winter snow only sporadically. Thermal regimes vary among streams with differences in basin area, elevation and hydrology. The six study streams spanned a wide range of environmental conditions in the basin (Table 1), including one with spring-fed base flow generating stable, cold temperature and flow year-round (Stream E, "Cold Creek," Figure 3) and one with a small basin area and bedrock substrate, often drying down to disconnected pools by late summer (Stream F, "South Creek"). We collected data from 30-m reaches at each of the study streams and will henceforth refer to these locations as "streams," which are coded alphabetically (A-F) in order of increasing elevation ( Table 1).

| Focal species
We selected the four most abundant, spring-emerging species from the orders Ephemeroptera, Plecoptera, and Trichoptera ("EPT") across the six study streams. Prior studies in the region have shown spring-summer adult emergence phenology for each.
We assumed all four species to be univoltine because each had single narrow or broad annual emergence periods in other studies from western Oregon (Anderson et al., 1984;Harper et al., 1995;Kerst & Anderson, 1974;Progar & Moldenke, 2009). However, just one of the four species (the mayfly Neoleptophlebia temporalis) has been confirmed to be univoltine (Lehmkuhl & Anderson, 1971).
Moselia infuscata has also been considered univoltine (Banks et al., 2007;Muchow & Richardson, 2000), but no published life cycle information is available to confirm. Diapause has not been evaluated for any of our focal species but is common in Plecoptera in egg or nymph stages (Bogan, 2017;DeWalt & Kondratieff, 2019), and some Trichoptera diapause in egg, larval or adult stages (Holzenthal et al., 2015). Ephemeroptera are not known to diapause as nymphs or adults, but egg diapause has been documented in a few species (Sartori & Brittain, 2015). Henceforth, we refer to the four focal species by their generic names.

| Data collection
Temperature dataloggers (HOBO Pro V2 U22-001, accuracy 0.2°C, Onset Computer Corporation) were programmed to record instantaneous temperature at half-hour intervals, placed in protective flowthrough housings, and secured in each of the six study streams. All data were downloaded twice yearly, archived, and compiled at the end of the study into a single file for each stream. Following the statistical procedure detailed by Frey, Hadley, Johnson, et al. (2016), we quality checked all data from 1 January 2009 through 31 July 2014 to flag any outliers and gaps. Unflagged half-hourly data were then averaged to an hourly time step, and linear regression relationships among all six streams were quantified. We also calculated regressions with hourly water temperature data from eight other nearby streamflow gauges in the Lookout Creek basin (https:// doi.org/10.6073/pasta/ 9437d 16030 44f5b 92189 110dd 8343763; Gregory & Johnson, 2019). Data flagged as anomalous or missing from our six study streams for the 6-year duration of the study (<2% of all values) were then replaced or filled according to the resulting best-fit regressions (see Appendix S2 for an example). We then converted the hourly temperatures to daily means to calculate degreedays (as described below). All final hourly temperatures are publicly available (https://doi.org/10.6073/pasta/ 7dc5e 1f888 28108 1de71 81808 2bd0c78; Johnson, 2017).
We collected emergent insects in each stream/year throughout the spring and early summer season with four emergence traps deployed per stream (Frady et al., 2007). Each trap covered 0.3 m 2 of stream surface and consisted of a pyramidal PVC frame draped with no-see-um netting. At the top of each trap frame, we suspended an insect collection cup containing propylene glycol with a few drops of unscented biodegradable soap. Traps remained in place from a mean start date of 18 April to a mean stop date of 1 July (Appendix S1) in each of the six years (2009-2014), with collections typically retrieved every 3-4 days for the duration of the study period. Annual start date was delayed at the two highest-elevation streams in some years when late-season snow prevented access, particularly in the coldest year of study (2011 ,   Table 1). On each retrieval date, we used a 250μm sieve to filter insects from the preservative and aspirated any additional individuals remaining in the trap mesh. Specimens were then preserved in 95% ethanol and transported to the laboratory for identification and enumeration.

| Data analysis
We used mean daily water temperature for each stream throughout the study period to calculate cumulative degree-days (DD) from ordinal date 1 (1 January) to 189 (8 July in nonleap years). Cumulative DD to ordinal date 189 represented how warm each stream/year was, in terms of heat-unit accumulation, from mid-winter throughout the observed spring emergence period. We also calculated the total DD accumulated to peak emergence date (described in the following paragraph) for each species/stream/year combination.
Many insects have emergence patterns that can be described by a normal (Gaussian) distribution of abundance of emergent adults through time (e.g. Harper & Peckarsky, 2006;Li et al., 2011). For each of the four focal species at each stream/year, we used a maximum likelihood routine in MATLAB (MathWorks, Inc.) to fit a parametric model to our empirical observations of numbers of emergent adults through time. Code for the model is available from Warren et al. (2012). From each fitted model, we used the mean ordinal date to represent the "peak" emergence timing for each species/stream/ year, and we extracted the 5th and 95th percentiles of each distribution to represent endpoints for describing a standardized duration, in days, of emergence period for each species/stream/year ( Figure 1).
Data were discarded for any individual species/stream/year combination that did not fit the model, typically due to low sample size (too few emergent adults collected). If peak emergence timing modelled for a particular species/stream/year fell later than the final date of empirical field sampling, we also removed that data point from further analysis. This latter issue was primarily faced during the coolest F I G U R E 2 Map of the stream network in the Lookout Creek basin with six focal headwater stream reaches (Streams A-F) marked. Broken line for Stream F indicates surface water is usually reduced to disconnected pools by late summer. See Table 1 for additional information about study streams sample year (2011) for the stonefly Moselia (which exhibits a relatively long emergence period, Table 2).
We used analysis of variance (ANOVA) with each of the four focal species to summarize and address whether overall patterns of three phenological traits varied consistently among the six streams in the Lookout Creek basin. The first set of ANOVAs used the modelled ordinal date of peak emergence as the response variable (timing).
The second set of analyses addressed whether duration of the emergence period varied, and the third asked if DD accumulated from ordinal date 1 through the peak emergence timing varied among streams, for each species. In each of these ANOVAs, replicates were sample years within streams. We assessed significance at α = 0.05.
We then used Analysis of Covariance (ANCOVA), with ordinal date of peak emergence for each species/stream/year as the phenological response variable to ask if there were consistent negative relationships between peak emergence timing and temperature. In a first set of ANCOVA, we used stream as the categorical independent variable and DD accumulated though ordinal date 189 as the covariate, for each of the four focal species. If slopes of the relationships were significant and negative, we interpreted this as evidence in support of H1 (earlier emergence in warmer conditions). Under the same hypothesis, we also expected a similar species-specific response among streams; therefore, we did not expect ANCOVA intercepts to be different. If slopes of the relationships between DD accumulation and ordinal date of emergence were not statistically different from zero (tested with 95% confidence intervals around slope values), we could not reject the null hypothesis and could infer that temperature alone was not strongly predictive of emergence timing.
Following an unexpected observation of significant slopes within streams but significantly different intercepts among streams for two of the four focal species (Alloperla and Neoleptophlebia), we ran a second set of ANCOVA in which we again used ordinal date of peak emergence for each species/stream/year as the response variable and DD accumulated through ordinal date 189 as the covariate, but with year instead of stream as the categorical independent variable. Under H1, we expected slopes of the relationships between DD accumulation and ordinal date of emergence to be significantly negative both within sample years (among streams of different temperature) and within streams (as above, among years with different temperature). Between the two ANCOVAs, if just one of the two sets of slopes was significant, this result would indicate that temperature influenced emergence timing but in a more complex manner than stated in H1. (19%) and fourth most abundant was Dolophilodes (7%) (Appendix S1). Emergent Dolophilodes were rarely trapped at the two highestelevation streams but were collected in sufficient numbers at the remaining four streams to be included in the statistical analyses.  Figure 3). Temperature patterns also varied among the 6 years of the study, with 2011 the coolest year and 2014 the warmest ( Table 1).

| RE SULTS
Gaussian models fitted 87% (125) of the 144 possible species/ stream/year emergence datasets. These models allowed us to extract ordinal dates of peak emergence and duration of emergence ( Of the four species, Moselia had the longest emergence duration (mean 50.6 days between the 5th and 95th percentiles of the distributions; Table 2, Figure 4b). Dolophilodes (25.0 days) and Alloperla (  Moselia was missing 2011 data for streams C-F because modelled peak emergence occurred later than the empirical collection period in the coolest year, and those data were discarded between temperature and peak emergence timing. The caddisfly Dolophilodes exhibited an overall significant negative relationship between temperature and emergence timing temporally (year-toyear) and spatially among streams in the basin; that is, Dolophilodes emerged predictably earlier both in warmer years and within years in warmer streams (Figure 5a). The second observed pattern was for the stonefly Moselia, which showed a weak but significant overall negative relationship between temperature and emergence timing but no clear effects either within streams or within years due to substantial scatter of the data around the stream-or year-specific lines of best fit (Figure 5b, Table 3). Pattern 2 might be applicable to any species with lengthy emergence periods that have naturally high variation around the peak timing.
The third and final species-specific pattern was a strong response to interannual (temporal) variation in temperature but limited spatial response among streams within years. Pattern 3 was observed for both the stonefly Alloperla and the mayfly Neoleptophlebia ( Figure 5c,d). The strong effect of stream but negligible significance of the interaction terms in ANCOVA (Table 3) Table 4), but simple linear regression applied to each single year of data revealed only two of six slopes to be significantly nonzero for Alloperla and just one of six slopes as nonzero for Neoleptophlebia (Table 5). Furthermore, the magnitudes of slopes were substantially lower within years than within streams (the first ANCOVA).

| DISCUSS ION
Contrary to an expectation of a simple relationship of earlier emergence timing with warmer stream conditions (hypothesis H1), we observed three distinct phenological patterns among the four aquatic insect species in our study. Of the three patterns conceptualized in Figure 6, Pattern 1 (Figure 6a) depicts emergence timing responding to water temperature differences across both space and time and supports H1. We observed this pattern in just one of our four focal species, the caddisfly Dolophilodes dorca. Overall, our observations suggest that temporal (interannual) variation in temperature, more so than spatial variation, influences emergence timing. However, emergence duration appears to be more tightly associated with spatial than temporal variation in temperature, such that populations in warmer streams have longer emergence periods (supporting H2; Figure 6). Perhaps surprisingly, the answer to the general question of whether emergence phenology of spring-emerging aquatic insects is predictable based solely on differences in heat accumulation is a resounding "no."

| Timing of emergence
All four species showed some evidence of earlier emergence timing with warmer conditions (Figure 6a-c), but two of the four exhibited the intriguing Pattern 3, which implicates a stronger phenological response to temporal than spatial differences in water temperature (Figure 6c). Describing a spatially near-synchronized response to interannual temperature differences, Pattern 3 reflects minimal phenological response to among-stream thermal differences despite relatively strong interannual variation. To our knowledge, Statistical significance denoted by ***=.001, **=.01, *=.05. Interaction terms are associated with the null hypothesis that slopes were not different among streams. Results plotted in Figure 5.
Abbreviation: NS, not significant. a Model 1: Includes interaction between stream and DDrate. b Model 2: Does not include interaction between stream and DDrate. this pattern has never been documented nor hypothesized for any aquatic insect. It is particularly noteworthy because the amongstream range of mean water temperature (3.4°C between coldest and warmest streams) was more than double the interannual range (1.6°C between coldest and warmest years of study; Table 1).
Among-stream differences in heat accumulation required to reach peak emergence for the Pattern 3 stonefly Alloperla fraterna and mayfly Neoleptophlebia temporalis approached 500 degree-days between the coolest and warmest streams. These differences allowed the near-synchronous emergence timing observed across heterogeneous streams in the Lookout Creek basin within years. In some insects, photoperiod cues drive spatially coordinated response (e.g. Bean et al., 2007), but our Pattern 3 species differed in peak emergence timing by nearly a month between the coolest and warmest years of observation. These interannual differences rule out photoperiod cues and point instead to a strong role of diapause or quiescence periods induced by common environmental variables experienced similarly in streams throughout the Lookout Creek basin.
Diapause and quiescence are different types of dormancy in the life cycle, during which development is paused and metabolic activity drastically slowed, leading to similar phenological responses (Diniz et al., 2017). Diapause is under genetically determined hormonal control, whereas quiescence is a plastic response to environmental cues and can therefore be "switched on and off" much more quickly than diapause. Most databases for aquatic insect traits (e.g. Poff et al., 2006;Tachet et al., 2010) include diapause but not quiescence, possibly because the capacity for quiescence is ubiquitous in insects and typically initiated in response to temperature extremes (e.g. Newbold et al., 1994). For example, diapause at any life stage is rare in mayflies, but most Ephemeroptera probably have minimum and/or maximum temperature thresholds for nymphal growth and development (Newbold et al., 1994;Tokeshi, 1985;Ward & Stanford, 1982). Importantly, combinations of high-and/or low-temperature thresholds to dormancy can allow escape from extreme environmental conditions and synchronize population life cycles across heterogeneous environments. In insects, such synchronization has mainly been recognized along latitudinal gradients (Bentz et al., 2014;Newbold et al., 1994). Our results now suggest that similar drivers might synchronize life cycle events across much smaller spatial extents.
Conditions affecting dormancy response during preceding winters might play a role in synchronizing some spring-emerging  (Forrest, 2016;Marshall et al., 2020). Winter temperatures can influence both the timing and depth of dormancy periods in complex and sometimes counterintuitive ways. For example, many terrestrial insects require low temperatures, known as the "chill factor," to induce a deep diapause, after which spring development can then proceed more rapidly than if winter temperatures remained warm enough to forego diapause (Forrest & Thomson, 2011;Stålhandske et al., 2017). Furthermore, temperature thresholds to induce coldtemperature dormancy can vary among instars within populations (Bentz et al., 2014;Logan & Bentz, 1999), a mechanism that helps synchronize population phenology. More research on the effects of both summer (Sweeney et al., 1992) and winter thermal regimes on complex phenological responses in aquatic insects is strongly needed.

| Duration of the emergence period
In support of H2, populations in warmer streams tended to have longer emergence periods. In contrast to emergence timing, emergence duration was more responsive to spatial variation than to temporal variation. Our measures of emergence duration can be translated directly as the degree of synchrony among individuals completing the life cycle transition to the adult stage within a population. Mechanistic drivers of intraspecific differences in emergence duration/synchrony are not well understood but are hypothesized to be influenced by environmental stability at the scale of the local habitat (Corkum et al., 2006;Li et al., 2011). In our study, the observation of longer emergence duration (decreased local synchrony) at lower-elevation streams suggests that water temperature is a cue.
Other variables thought to affect emergence period duration, including predation pressure (Peckarsky et al., 2000), flow disturbance (Lytle, 2001), and quality of conditions for adult activity (Schultheis et al., 2008) are unlikely to vary consistently with elevation within the limited spatial extent of Lookout Creek or among its trophically similar headwaters (Zatkos et al., 2021). Although the mechanism is unclear for how water temperature influences emergence period duration, the pattern has been observed elsewhere (e.g. Cheney et al., 2019;Li et al., 2011;Richardson & Clifford, 1986) and merits further investigation. As with thermal drivers of emergence timing, responses are likely complex and variable among species. For example, the "chill factor" required in some terrestrial insects to induce deep diapause has also been linked to greater emergence synchrony in spring-emerging species (Forrest, 2016). It is also possible in some species that extended duration of emergence period following warmer winters confers a selective advantage in anticipation of longer breeding seasons, such as in spring-emerging montane grasshoppers (Buckley et al., 2015).

| Climate change implications
Among-year differences in peak emergence timing of nearly 1 month associated with differences in mean water temperature of <2.0°C in our Pattern 3 species is a cause for concern associated with climate change. Both earlier emergence timing and increasing duration of emergence periods have been documented in a European stream that warmed during a 42-year study period by 1.9°C (Baranov et al., 2020), and similar patterns are expected in other regions. Lookout Creek has experienced a recent trend of decreasing late-summer streamflow (Kaylor et al., 2019), but despite overall warming of air temperature, multi-decadal trends in water temperature appear to be highly variable among streams of the topographically and hydrologically complex basin (Arismendi et al., 2012). We argue here that thermal heterogeneity within and among streams in a naturally complex basin, interactions of other species-specific biological traits with phenological traits, and the potential for phenological traits to adapt to changing conditions might each contribute to ecosystemscale resilience in the face of ongoing climate change.
A value of our study in comparison to multi-year studies at single locations (e.g. Baranov et al., 2020;Finn & Poff, 2008;Harper & Peckarsky, 2006) is that we incorporated the natural mosaicism of thermal conditions among streams in a topographically complex basin. Concurrent with our study, terrestrial investigations in the Lookout Creek basin also demonstrated spatially heterogeneous microclimates associated with air temperature (Ward et al., 2018), and habitat use by birds  and terrestrial insect flight patterns (Schmidt, 2019) were linked to microhabitat-scale thermal regimes, which shifted spatially from year to year. Like these terrestrial taxa, most emergent aquatic insects fly and can therefore respond behaviourally to environmental heterogeneity across streams in a basin by dispersing among tributaries. Similarly, thermal heterogeneity typically exists at the local scale within streams (Johnson, 2004;Kaylor et al., 2021), to which immature stages of aquatic insects can also respond behaviourally (Uno & Power, 2015), particularly when a hyporheic zone is present and accessible (Dorff & Finn, 2020). For example, Moselia infuscata likely occupies the cooler and more thermally stable hyporheic zone for at least part of its nymphal development period (Gill et al., 2015), which could help explain the highly variable patterns of emergence timing and dura- But caddisfly adults are relatively long lived (on the order of weeks) and typically are strong flyers. The combination of asynchronous emergence among streams with relatively long-lived and TA B L E 5 Slope values from simple linear regression between heat accumulation to ordinal date 189 (in degree-days) and ordinal date of peak emergence for the stonefly Alloperla and the mayfly Neoleptophlebia (the two "Pattern 3" species). Stream is the categorical variable (with sample years as replicates) for the first two columns; year is the categorical variable (with streams as replicates) in the final two columns. The relationship between temperature and timing of emergence was overall stronger within streams (among years) than within years (among streams) for these two species F I G U R E 6 Conceptual diagrams representing three major patterns of emergence timing and duration observed in this study. (a) Pattern 1: relatively consistent spatial and temporal response of emergence timing to temperature, and longer duration of emergence period in warmer years and at warmer streams (Dolophilodes); (b) Pattern 2: some evidence for response to heat accumulation rate but overall lengthy emergence periods make patterns difficult to distinguish among streams and years (Moselia); (c) Pattern 3: temporal (among-year) response to temperature but spatial near-synchrony of emergence timing across cooler and warmer streams within years, with longer emergence period in warmer years (Alloperla and Neoleptophlebia) Overall, the phenological responses we quantified in this study were highly variable among species and among populations within species. There was no evidence that phylogenetically more closely related species shared similar phenological traits, an observation that suggests a high degree of evolutionary lability ("flexibility") to adapt relatively quickly to changing conditions (Poff et al., 2006).
Although greater lability implies that these traits are difficult to predict and must be measured on a species-by-species basis, it also presumably confers greater potential for adaptation to relatively rapid environmental change. Evidence for rapid evolution of phenological traits exists for some terrestrial insects, including pine beetles that have intraspecific genetic differences associated with dormancy cues along latitudinal gradients (e.g. Bentz et al., 2014). Overall, both high levels of plasticity in some traits (e.g. duration of emergence) and the potential for relatively rapid evolution of others might allow flexibility of phenological response into the future and contribute to ecosystem-scale resilience with climate change.

| CON CLUS IONS
Aquatic ecology has a rich heritage of mechanistically linking thermal regimes to distributions and habitat associations of invertebrates (e.g. Vannote & Sweeney, 1980;Ward & Stanford, 1982), but the physiological mechanisms associated with phenological response to specific thermal cues remain little understood (Woods et al., 2021). Our detailed spatial and temporal observations among thermally different headwater streams shed new light on the complexity of phenological responses of aquatic insects to water temperature. Among just four species, we observed three distinct patterns of emergence timing associated with temperature. The intriguing Pattern 3 suggests complex physiological response to thermal cues generating near-synchronous emergence from thermally distinct streams despite strong differences in peak timing among years. Variability in peak emergence timing by up to a month between years with <2℃ mean temperature difference can help explain interannual variation in taxa present in benthic communities sampled at the same time each year (e.g. Frady et al., 2007). Another practical implication is that because many insects likely respond to winter thermal conditions, it is important to rear individuals throughout the life span for studies measuring phenological response to experimental treatments (e.g. Funk et al., 2019;Li et al., 2011).
The sensitive and variable responses documented among species underscore the importance of mechanistic explanations linking phenological traits to water temperature and improved predictive ability associated with climate change in aquatic ecosystems. For aquatic ectotherms in general, we recommend closer collaboration between physiologists and ecologists (Chmura et al., 2019) and more attention to the terrestrial insect literature for promising research angles, particularly associated with dormancy cues and the potential role of winter chill factors. It will also be important to improve our understanding of other biological traits that could interact with phenological traits (Forrest, 2016;Nelson et al., 2020;Thomas & Bacher, 2018). We hypothesized a trade-off between complex response mechanisms associated with phenological synchronization and traits such as adult longevity and/or dispersal ability. Studies designed to test this hypothesis might identify predictable "syndromes" of strongly associated traits that can improve communityand ecosystem-level prediction of phenological response to climate change in aquatic ecosystems.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The stream temperature data that support the findings of this study are openly available from the EDI data portal at https://doi. org/10.6073/pasta/ 7dc5e 1f888 28108 1de71 81808 2bd0c78. All insect emergence data are provided in supplementary material and in