Climate Alters the Migration Phenology of Coastal Marine Species Climate Alters the Migration Phenology of Coastal Marine Species

: Significant shifts in the phenology of life-cycle events have been observed in diverse taxa throughout the global oceans. While the migration phenology of marine fish and invertebrates is ex - pected to be sensitive to climate change, the complex nature of these patterns has made measurement difficult and studies rare. With continuous weekly ob - servations spanning 7 decades in Narragansett Bay, Rhode Island (USA), the University of Rhode Island Graduate School of Oceanography trawl survey provides an unprecedented opportunity to investigate the influence of climate on the migrations of marine species in the northwest Atlantic. Analyses of the survey observations of 12 species indicated that residence periods have changed by as much as 118


INTRODUCTION
The timing of life-cycle events can profoundly affect the fitness and reproduction of organisms.Often associated with environmental cues, the phenologies of such events frequently provide direct and easily measured evidence of the ecological effects of global climate change (Walther et al. 2002, Root et al. 2003).While many past efforts have focused on birds, butterflies, and plants due to their accessibility (Walther et al. 2002, Robinson et al. 2009), phenological shifts ABSTRACT: Significant shifts in the phenology of life-cycle events have been observed in diverse taxa throughout the global oceans.While the migration phenology of marine fish and invertebrates is expected to be sensitive to climate change, the complex nature of these patterns has made measurement difficult and studies rare.With continuous weekly observations spanning 7 decades in Narragansett Bay, Rhode Island (USA), the University of Rhode Island Graduate School of Oceanography trawl survey provides an unprecedented opportunity to investigate the influence of climate on the migrations of marine species in the northwest Atlantic.Analyses of the survey observations of 12 species indicated that residence periods have changed by as much as 118 d, with shifts in the timing of both ingress to and egress from the coastal zone.The residence periods of warmwater species expanded while those of cold-water species contracted.Dirichlet regressions fit to the annual presence−absence patterns of each species identified interannual temperature variations, fluctuations in ocean circulation, and long-term warming all as having a significant effect on migration phenology.Additionally, temperature gradients within Narragansett Bay were shown by generalized additive models to cause detectable shifts in local spatial distributions during coastal residency.These novel findings mirror results found in the spatial domain and therefore suggest that the studied species are adapting their spatiotemporal distributions to track their thermal niche in a changing climate.If so, characterizing the spatial and temporal aspects of climate responses across species will be critical to understanding ongoing changes in marine ecosystems and successfully managing the fisheries they support.
have been observed across a diversity of taxa (Robinson et al. 2009, Poloczanska et al. 2013, Cohen et al. 2018, Staudinger et al. 2019).Temporal changes in the timing of life-history events tied to climate processes acting at local to global scales impact breeding success, patterns of population abundance, migration routes and durations, and trophic interactions (Robinson et al. 2009, Socolar et al. 2017, Cohen et al. 2018).These responses can vary among and within species by feeding strategy, life history, geography, and sex or age (Walther et al. 2002, Mac Mynowski & Root 2007, Kluen et al. 2017).Migratory species appear to be particularly vulnerable to effects of climate change, as their reliance on seasonal availability of resources in habitats that can be separated by great distances may make adaptation difficult (Robinson et al. 2009).When species are able to alter the timing of life-cycle events, perhaps in concert with changes in spatial distribution, they may be able to track their preferred niche in a changing climate (Poloczanska et al. 2013, Socolar et al. 2017).Because the extent of niches among and between trophic levels may shift in space and time at different rates (Robinson et al. 2009, Cohen et al. 2018), changing phenologies can upset competitive interactions and predator−prey dynamics (Soberón & Nakamura 2009, Poloczanska et al. 2013) and complicate efforts to anticipate future changes.
Despite this breadth of knowledge, our understanding of phenological changes in marine fish and mobile invertebrates is limited.The complex life histories of many marine species make them vulnerable to climate change, but investigations of their responses are difficult and relatively rare (Sydeman et al. 2015).As ectotherms, fish and marine invertebrates can be impacted by climate change directly, via temperature changes (Poloczanska et al. 2013), or indirectly through food web interactions (Bal et al. 2017, Bell et al. 2017).Similar to results found for other taxa, phenological changes of life-cycle events in fish affect growth and body condition (Bal et al. 2017, Pershing et al. 2018), trophic dynamics (Bell et al. 2017, Pershing et al. 2018), survival and recruitment (Bell et al. 2017), migration route and duration (Sims et al. 2004, Pauly & Keskin 2017), genetic selection (Manhard et al. 2017), and vulnerability to harvest (Peer & Miller 2014).Many of these features are tied to changes in migration phenology, which can be particularly difficult to track in fish.The majority of the examples in the literature focus on salmonids (e.g.Kovach et al. 2013, Otero et al. 2014, Bell et al. 2017, Manhard et al. 2017), likely due to their accessibility, as they use rivers for spawning and nursery habitat.The rare studies of migration phenology in marine fish and invertebrates, meanwhile, have been limited to investigations of a single species and time series of less than 30 yr (e.g.Sims et al. 2001, 2004, Peer & Miller 2014, Pauly & Keskin 2017) or point-to-point changes among a group of species in only portions of the year over longer time series (e.g.van Walraven et al. 2017).While providing some insight, neither study design can fully provide the detailed understanding of how climate change influences diverse marine species to shift phenological patterns that is needed to fully assess the implications for ecological function and fisheries production.
Developing an understanding of shifts in patterns of migration phenology may be particularly important on the Northeast US (NEUS) continental shelf, a region experiencing extremely rapid rates of temperature change (Pershing et al. 2015, Chen et al. 2020).Supporting very productive and valuable commercial fisheries (NMFS 2018), the NEUS shelf has undergone significant ecological transformation in recent decades.The seasonal distributions of many species have changed in extent, depth, and/or latitudinal range in response to long-term warming, climate fluctuations, and interannual temperature variations (Nye et al. 2009, 2011, Henderson et al. 2017).In addition to distribution shifts, climate forcing has also affected species productivity (Xu et al. 2018, Tableau et al. 2019).As a result, the species assemblage of sub-regions of the shelf have changed significantly, each now resembling the historic biological community of the sub-region to its immediate south (Lucey & Nye 2010).Furthermore, model projections suggest that the thermal habitats for the species of the NEUS shelf will undergo additional dramatic poleward shifts by the end of the 21 st century (Kleisner et al. 2017).
At the latitudinal center of the NEUS shelf system, Narragansett Bay, Rhode Island, has similarly exhibited significant biological change in the past 60 yr.This large southern New England estuary provides a unique opportunity to investigate the effects of climate on marine phenological change.Since 1959, the University of Rhode Island Graduate School of Oceanography (URI GSO) has conducted a continuous weekly trawl survey of 2 distinct habitats within Narragansett Bay (Fig. 1), providing the sampling frequency necessary to measure changes in migration timing and their relationship to environmental processes that is often lacking in marine systems (Staudinger et al. 2019).The data collected by this survey indicate pronounced changes in the fish and invertebrate community, with a distinct assemblage shift detected in the early 1980s linked to warming temperatures and climate fluctuations (Collie et al. 2008).On the boundary between the colder Gulf of Maine and the warmer Mid-Atlantic Bight (Fig. 1; Thomas et al. 2017, Chen et al. 2020), this ecosystem has historically supported species from each region in different seasons.However, many of the coldwater resident species have been largely replaced by warm-water seasonal transients (Collie et al. 2008).Because many of the species captured in the URI GSO survey migrate seasonally to the mid-to outercontinental shelf (Cross et al. 1999, Packer et al. 1999, Steimle et al. 1999a,b,c, Jacobson 2005), examining variations in observed coastal occupancy can provide insight into integrated responses to conditions encountered along migratory routes spanning hundreds of kilometers.
In order to gain new insight into the relationships between the migrations of marine species and climate, we had 2 objectives for the dominant migratory fish and invertebrates captured in trawl surveys of Narragansett Bay: (1) to characterize changes in migration phenology and identify links to climatic processes acting at several spatial scales, and (2) to evaluate local distributional responses to temperature during the coastal residence period.In doing so, the results can be used to better understand how marine fish and invertebrate movement patterns are changing in a warming ocean and the implications for future shifts in distribution and productivity.

URI GSO survey
The weekly URI GSO trawl survey is conducted at 2 stations 15.3 km apart in Narragansett Bay: a midbay site with a mean depth of approximately 7 m and an outer-bay site with a mean depth of approximately 23 m (Fig. 1; Collie et al. 2008).At each site, a bottom otter trawl net with an effective opening of 6.5 m and headrope height of 1 m is towed at 3.7 km h −1 (2 nautical miles h −1 ) for 30 min.Because the net is made of 7.6 cm mesh changing to 5.1 cm in the codend, collections consist primarily of the adult life stages of fish and invertebrates.Auxiliary information during the tows, including date and surface and bottom water temperature, is recorded.Upon retrieval of the net, all catch is classified by species or the nearest discernible taxonomic level.Each species is then enumerated and weighed in aggregate.Four research vessels have been used to conduct the trawl survey since 1959, but the trawl nets, deployment, tow speed, and tow duration have been standardized throughout the time series.
Bottom temperature was inconsistently measured by the URI GSO trawl survey between 1964 and 1977, and occasional missing values are present throughout the remainder of the time series.Because the available temperature data suggest that Narragansett Bay typically has small, predictable vertical thermal gradients, the surface water temperature and day of the year were strong predictors of bottom temperature.To address these missing data, generalized additive models (GAMs) were fit to the measured bottom temperatures at each of the URI GSO trawl survey stations separately, with surface temperature, day of year, and year as independent variables.Fit with the 'mgcv' package (Wood 2013) in R version 1.3.959(R Core Team 2020), the GAMs explained over 97% of the residual deviance.These GAMs were then used to impute missing bottom temperature values throughout the URI GSO trawl survey time series in order to perform subsequent analyses on the temperature distribution of species observations at the depth of capture.

Species selection and classification
Candidate species for analysis were selected as those for which at least 50 individuals were caught over at least 10 sampling events per year, in at least 20 years between 1959 and 2016, in order to ensure suffi-cient detection frequency for measuring phenological changes.Of these candidates, any species that was observed for an average annual duration of ≥ 335 d was considered a year-round resident and excluded.Additionally, any taxa that had not been identified to the species level or were known to be non-migratory (e.g.gastropods) were also excluded.The species that met these criteria were then split into 2 groups based upon presence patterns during the seasonal temperature cycle: summer species were those that had been primarily observed at water temperatures above 15°C, and winter species were those that had been primarily observed below 10°C.In order to realign the calendar to match the continuity of the presence−absence patterns of each species group for analysis, Days 1−365 were set to 1 March to 28 February (29 February in leap years) for summer species and 16 August to 15 August for winter species.Winter species that were present in the study area in the fall of 2016 would not migrate offshore until spring 2017, and thus analysis of this group considered years through 2015.

Change in residence time
Annual residence time for each species was defined as the difference in days between the first and last observation in the URI GSO trawl survey during each adjusted calendar year.The first and last observations were first screened for outliers.Outliers could occur, for example, if rare individuals of a species overwinter in the coastal zone while the vast majority of the population migrates elsewhere, as is the case in juvenile summer flounder Paralichthys dentatus (Packer et al. 1999).Specifically, observations of a single individual of a species separated from the nearest observation by at least 1 standard deviation of the year-day abundance distribution, in a year in which at least 30 individuals were captured, was considered an outlier and removed.Changes in residence time between 1959 and 2016 were assessed by performing t-tests on the residence times between the first 10 and last 10 yr of the time series in which a given species was observed.We tested the hypothesis that summer species expanded their residence time between 1959 and 2016, while winter species contracted their residence time.Similarly, the bottom temperatures measured at the first and last observation and the minimum and maximum bottom temperatures at which each species was observed were compared with the same method and selected years to test the hypothesis that species migrate in response to temperature cues and exhibit a stable thermal niche over time.
The number of days before the first observation, the residence time, and the number of days after the last observation of each species in each year were then converted to proportions of the year (referred to herein as presence−absence components) and analyzed jointly by Dirichlet regression (Maier 2014) with the package 'DirichletReg' in R (Maier 2015) to evaluate the role of climate in shifting patterns of residency.Dirichlet regression is an analytical approach for modeling compositional data, which may exhibit heteroskedasticity or skewness, without need for transformation (Aitchison 1986, Maier 2014).The suite of potential covariates to be tested in the regression models included measures of abundance, temperature, and precipitation.Past work has suggested that fish populations may expand their spatial distribution when abundance increases (MacCall 1990, Bell et al. 2015).To control for the possibility that such a phenomenon may occur in time, in addition to the direct, saturating relationship between local species abundance and detectability in the trawl survey, the annual mean log catch-per-tow of each species was included in the central presence component of all Dirichlet regression fits regardless of statistical significance.
The long-term warming trend observed in Narragansett Bay was represented as the annual mean bottom temperature recorded in the URI GSO trawl survey.In order to evaluate the hypothesis that temperature conditions over the spatial extent of migratory routes influence the phenology of coastal ingress and egress, bottom temperature trends across the shelf were included.Specifically, the April (1968April ( − 2016) ) and October (1963October ( −2016) )  The temperature gradients between these regions were also considered, calculated as the difference in the seasonal mean temperatures between adjacent zones.The seasonal temperature gradients across the mouth of Narragansett Bay were calculated as the difference between the mean April and October bottom temperatures recorded by the URI GSO trawl survey and the seasonal mean bottom temperatures in the Rhode Island nearshore zone.In order to characterize spatiotemporal patterns in these region-specific temperatures, temporal trends were investi-gated using linear regressions, and differences between the 2 halves of the time series (1959−1987 and 1988−2016) were evaluated among months and regions using 2-way ANCOVAs, where year was included as a covariate, with Tukey's HSD tests.
Past work on the NEUS shelf has shown that climate fluctuations, including the North Atlantic Oscillation (NAO) and Atlantic Multidecadal Oscillation (AMO), and changes in ocean circulation patterns influenced spatial distributions of fish and invertebrates (Nye et al. 2009, 2011, Lucey & Nye 2010).Similarly, these oscillations have been linked to shifts in abundance of species at various trophic levels within Narragansett Bay (Collie et al. 2008, Borkman & Smayda 2009) and changing circulation patterns in the southern New England coastal zone (Luo et al. 2013).In order to test for such relationships with migration phenology, the winter (December−March mean) NAO index was obtained from the National Weather Service Climate Prediction Center (www.cpc.ncep.noaa.gov) as in Xu et al. (2015).The AMO index was cal culated as the annual mean of the unsmoothed monthly time series available from the NOAA Earth System Research Laboratory (www.esrl.noaa.gov/).Finally, the Regional Slope Water Temperature (RSWT) index, representing the coupled impacts of the Labrador Current and Gulf Stream on the temperature of shelf slope waters along the NEUS shelf, was calculated as described by Pershing et al. (2001).The AMO and RSWT index were considered at lags of up to 2 yr.
To investigate the effects of local-scale climatic variations, seasonal temperature anomalies and precipitation accumulations were calculated for Narragansett Bay.Here, mean bottom water temperatures for winter (January− March), spring (April− June), summer (July− September), and fall (October− December) were taken from the URI GSO trawl survey and de-trended using loess smoothers with a span and degree of 1. Seasonal precipitation values, calculated as the total accumulation at the Theodore Francis Green Memorial State Airport in Warwick, Rhode Island, were obtained from the NOAA National Centers for Environmental Information (www.ncdc.noaa.gov).For the winter season, rain and snow accumulations were tested separately.Trends in the seasonal bottom temperatures and precipitation accumulations were evaluated using linear regressions and generalized linear models with gammadistributed errors (gamma regressions), re spectively.
Forward model selection of the Dirichlet regression fits for each species was performed using a change in the Bayesian information criterion (ΔBIC) of 2 points for variable inclusion to maximize parsimony.Here, covariates were tested separately in each of 3 possible components of the annual presence−absence pattern.In order to minimize the incidence of multicollinearity, a single covariate was never loaded onto adjacent components.Further, covariates were only loaded onto presence−absence components that had temporal overlap or occurred later in the year.For example, the fall temperature anomalies were not tested in the component defining ingress for a species that enters Narragansett Bay in spring.For all temperature covariates, or climate fluctuations associated with temperature, we tested the hypothesis that summer species will exhibit positive correlations with residence time and winter species will exhibit negative correlations.Meanwhile, past work by Hanson (1982) suggested that precipitation events result in increased input of bioavailable nutrients into Narragansett Bay.Therefore, we tested the hypothesis that increased precipitation may result in augmented primary productivity and thus increase fish residence through bottom-up control.

Distribution patterns during coastal residence
Finally, we examined the distributional response of each species to temperature gradients within Narragansett Bay during the residence period.Evidence of shifting distribution patterns was evaluated by fitting GAMs with beta-distributed errors using the 'mgcv' package (Wood 2013).Here, the proportion of the total catch of a given species on each survey day was calculated for the mid-bay site.These proportions were then transformed to strictly fall between 0 and 1 using the transformation suggested by Smithson & Verkuilen (2006).Interannual and seasonal trends were controlled for with splines fit to the year and mean bottom temperature between the stations for each observation.The bottom temperature difference between the survey stations at the time of each observation was then used to test the hypothesis that summer (winter) species shift their local spatial distribution toward (away from) warm temperatures during the residence period.

Change in residence time
Out of the 154 species or species groups documented in the URI GSO trawl survey between 1959 and 2016 (a total of 5546 trawls across both stations), 26 surpassed the abundance thresholds to be considered for evaluation.Of these, 12 seasonal migrant species (11 fish and 1 squid) met the criteria for analysis of movement phenology (Table 1).Eight species were classified as summer species and the remaining 4 as winter species.Ten of the species were present throughout the 58 yr time series.Ocean pout (see  1. Species selected for analysis, their seasonal classifications, and the number of years in which residence time could be calculated.Significant p-values for the t-tests on changes in the residence time (days) between the first and last 10 yr in which each species was observed are designated as: *p < 0.05, **p < 0.01, ***p < 0.001 classified as outliers were only detected and removed in 4 species (longfin squid, northern searobin, scup, and summer flounder), and no single species had more than 4 outliers.Longhorn sculpin, ocean pout, red hake, scup, butterfish, summer flounder, striped searobin, and longfin squid exhibited statistically significant changes in residence time between the first and last 10 yr in the time series during which each was observed (t-tests, p < 0.05), with detected changes ranging from 31.8 to 118.0 d (Table 1, Fig. 2).Of those with detectable changes, all summer species expanded their residence period, and all winter species contracted theirs, as hypothesized (Figs. 2 & 3).Notably, red hake qualified as a year-round resident species prior to the mid-1980s, after which its annual residence period underwent a sharp decline into a strictly seasonal pattern of presence in Narragansett Bay.The contrast between the significant changes observed in the days of first and last observation each year, which define the residence time, and the results found for the interior quantiles of the abundance distribution by Henderson (2013) suggest that the largest shifts in temporal distributions have occurred at the extremes while the core remained relatively more stable (Figs. 2  & 4).
The 4 winter species (longhorn sculpin, ocean pout, Atlantic herring, and red hake) showed statistically significant increases in the minimum temperature of observation while having no detectable changes in the maximum temperature.Two of the 4 winter species had significant differences in the mean temperatures of first and last observation each year, with longhorn sculpin arriving at warmer temperatures and ocean pout departing Narragansett Bay at warmer temperatures (Table 2).The temperature at last observation of these species also shifted colder by more than 3°C over time.The 4 summer species that had exhibited the largest increases in residence time, i.e. butterfish, striped searobin, summer flounder, and longfin squid, were observed at significantly colder minimum temperatures (t-tests, p < 0.05) in the most recent decade (Table 3).All summer species, except tautog, also showed signi ficant increases in the maxi-Fig.2. Average ingress and egress timing (in year-days) of the studied species during the first 10 yr (blues, top) and last 10 yr (reds, bottom) in which each was documented in the URI GSO trawl survey.Summer and winter species groupings are designated on the righthand axis.The residence period (lightest), 80% abundance interval (moderate), and 50% abundance interval (darkest) are indicated by shades and transparency levels of each color, and the median observation is depicted by a colored diamond.Statistically significant (p < 0.05) changes in ingress or egress timing are denoted by an asterisk (*).Mean striped searobin ingress during the first period was calculated using 1963−1966, 1969−1973, and 1977.Mean ocean pout egress during the second period was calculated using 1989−1998.In all other cases, the first period corresponded to 1959−1968 and the second period corresponded to 2007−2016.Red hake qualified as a year-round resident species during the first period and thus is depicted by an unbroken interval for the entire year mum temperature of ob servation.Tautog, northern searobin, fourspot flounder, scup, butterfish, and summer flounder had statistically significant differences between the mean temperatures at first and last observation, where final observations were warmer on average for all but butterfish (Table 2).Four summer species, namely butterfish, summer flounder, striped sea robin, and longfin squid, also showed significant decreases in the temperature at first and/or last observation over time.

Changes in the physical environment of southern New England
The annual average bottom water temperature in Narragansett Bay increased by 1.50°C between 1959 and 2016 (linear regression, p < 0.001).The rate of warming varied by season, such that the time series mean increase was greatest in winter (2.17°C) and least in summer (1.06°C) (linear regressions, p < 0.05).The minimum winter temperature warmed by 1.53°C, and the maximum summer temperature warmed by 1.40°C (linear regressions, p < 0.01).The winter bottom temperature was also significantly positively correlated to the NAO (linear regression, p < 0.001).The effect of the NAO was therefore removed prior to calculating the winter temperature anomalies to eliminate collinearity with any environmental process acting across larger spatial scales.
The seasonal precipitation accumulation time series exhibited no significant trends (gamma regression, p > 0.05) or correlations with other covariates (Pearson's product moment correlation, p > 0.05).
Patterns of warming differed between seasons and among the studied regions on the continental shelf and in Narragansett Bay (Fig. 5).Bottom temperatures in October (0.77°C) warmed more than in April (0.58°C) between the first (1959−1987) and second half (1988−2016) of the time series (Tukey's HSD, p < 0.05), although region-specific warming trends were only significant in the fall (linear re gressions, p < 0.05).During April, the warmest bottom temperatures were found at the mid-and outer shelf, whereas the pattern was reversed in October.The observed temperature ranges across the shelf regions spanned approximately 9°C in both seasons.Compared to the 0.024°C yr −1 increase in October bottom temperature in Narragansett Bay between 1963 and 2016 (linear regression, p = 0.006), the shelf regions were found to be warming at an ap proximately 1.4−2.3times greater rate.Given that the URI GSO trawl survey sites in Narragansett Bay are much shallower than the designated shelf re gions, this difference suggests that more than at mospheric forcing is affecting October bottom temperatures offshore, in agreement with  2011), the RSWT index correlates strongly to the October warming trends in the Rhode Island nearshore, mid-shelf, and outer shelf (Pearson's product moment correlation, ρ = 0.77− 0.89, all p < 0.001).In April, meanwhile, correlations between the RSWT index and the Rhode Island nearshore, mid-shelf, and outer shelf are also positive but substantially weaker (ρ = 0.36−0.58,all p < 0.05).

Relating residence time to changes in the physical environment
The Dirichlet regression fits indicated that the residence patterns in Narragansett Bay of 10 of the 12 species, excluding Atlantic herring and tautog, were significantly influenced by climate variables (Table 4; see Table S1 in the Supplement at www. int-res.com/ articles/ suppl/ m660 p001 _ supp.pdf).All seasonal temperature anomalies in Narragansett Bay and precipitation accumulations between spring and fall were included in re gression fits for at least 1 species.Beyond these localscale processes, the NAO, AMO, and RSWT index, in addition to the seasonal bottom temperatures in the designated regions spanning migration routes across the continental shelf, were also selected in regression fits.In nearly all cases, the tested covariates representing water temperature, either directly or indirectly as in the case of the NAO or RSWT index, had the hypothesized direction of effect on species phenology.The only 2 exceptions occurred in the regression fits for ocean pout and striped searobin.However, both cases had significantly correlated covariates in adjacent presence−absence components (ocean pout: October outer-shelf temperature and annual mean Narragansett Bay bottom temperature, Pearson's product-moment correlation = 0.547, p < 0.001; striped searobin: annual mean Narragansett Bay bottom temperature and spring temperature anomaly, Pearson's product-moment correlation = 0.430, p < 0.001), raising the possibility that collinearity altered the estimated regression coefficients.Seasonal precipitation accumulations were selected in the re gression fits of longhorn sculpin, butterfish, striped searobin, fourspot flounder, and summer flounder, with the effects estimated to be positive in the first 3 species and negative in the latter 2 (Table 4).Finally, the fitted coefficients on mean log catch-per-tow were statistically significant (Diri chlet regression, p < 0.05) for longhorn sculpin, ocean pout, Atlantic herring, tautog, butterfish, summer flounder, and longfin squid.

Distribution patterns during coastal residence
Analysis of the bottom temperature data from the URI GSO trawl survey revealed that a positive gradient (warmer mid-Bay) was present on average from mid-March to mid-October (Fig. 6), after which the gradient inverted until the following spring.Three of the 4 winter species, namely longhorn sculpin, ocean pout, and red hake, as well as fourspot flounder were observed at the mid-bay station in less than 10% of observations, indicating they had little habitat flexibility within Narragansett Bay (Table 5).Of the remaining 8 species, the bottom temperature gradient between the trawl stations exhibited a significant relationship with distribution (GAM, p < 0.05) in all but northern searobin (Table S2, Fig. S1).Within those 7 species, the direction of the effects suggests that all of the summer species shift their distribution toward the warm end of the gradient (Fig. S1b,d−h) and the winter species (Atlantic herring) shifted its distribution toward the cold end of the gradient as hypothesized (Fig. S1a).The selected model fit for longfin squid included an interaction term between the temperature gradient and the mean bottom temperature of the estuary.Specifically, the fit indicated that this species favored the warm end of the temperature gradient until the estuary reached an average bottom temperature of 16°C, the approximate upper bound of its preferred thermal range (Jacobson 2005).Beyond this mean bottom temperature, squid shifted its distribution toward the cold end of the gradient (Fig. S1h).The one species that had been ob served at the mid-bay site in at least 10% of capture events and did not have a significant relationship with the temperature gradient, i.e. northern sea robin, underwent a significant habitat shift during the URI GSO trawl survey time series.Prior to 1990, 68.9% of northern searobin were caught at the mid-bay station.After 1990, this percentage dropped to 23.3%.A GAM fit for northern searobin catch data prior to 1990 (Fig. S1c) suggested a strong positive relationship with the thermal gradient as in the other summer species.

DISCUSSION
Many of the species included in this analysis have been identified as having high sensitivity to climate change (Hare et al. 2016).Consequently, the majority of the studied species, longhorn sculpin, ocean pout, red hake, scup, butterfish, summer flounder, striped searobin, and longfin squid, exhibited significant, and often nonlinear, changes in their residence time in Narragansett Bay between 1959 and 2016.As expected and consistent with previous research (Day et al. 2018), winter species decreased their time spent in Narragansett Bay, while summer species extended their residence time.After controlling for the effect of abundance on detectability, the observed phenologi-  2019) across the defined shelf regions (see Fig. 1B) during the periods of 1959−1987 (blue) and 1988−2016 (red).The median of each spatiotemporal group is designated by a horizontal black line.Lowercase letters indicate statistically significant differences among locations within each season identified by a 2-way ANCOVA (Tukey's HSD, p < 0.05).The average bottom temperatures during the latter period were 0.58°C warmer in April and 0.77°C warmer in October (Tukey's HSD, p < 0.05), but no significant differences in warming among locations were detected between halves of the time series cal patterns in a group of species representing 10 taxonomic families and a wide range of life histories were best explained by the effect of climatic processes along migratory routes on their respective thermal niches, as hypothesized.Seasonal migration into and out of Narragansett Bay appears to be influenced by temperature conditions from the inner to outer NEUS continental shelf.Supporting this result, fish and squid species exhibited distributional responses to kilometer-scale temperature gradients within Narragansett Bay during the coastal residence period.Bottom water temperatures on the southern New England shelf have warmed significantly since the 1960s due to a combination of climate change and ocean circulation (Nye et al. 2011, Thomas et al. 2017, Chen et al. 2020), while also varying seasonally and spatially.Although the warming rates of each region have differed, Table 4. Selected covariates in the Dirichlet regression fits for each species, where blank cells denote that no significant covariates were identified.The symbols (+) and (−) indicate effects that increased and decreased the residence period, respectively.The annual mean log catch per tow (CPT) was included in all regression fits as a control variable, and its associated pvalues are designated as: () p > 0.05, *p < 0.05, **p < 0.01, ***p < 0.001.BT: bottom temperature; TA: temperature anomaly; RNS: Rhode Island nearshore region; MS: mid-shelf region; OS: outer-shelf region; NBBT: annual mean bottom temperature of Narragansett Bay; RSWT: Regional Slope Water Temperature index; AMO: Atlantic Multidecadal Oscillation index; NAO: North Atlantic Oscillation Winter index; Apr: April; Oct: October Fig. 6.Bottom temperatures recorded at the mid-bay (orange) and outer bay (purple) (see Fig. 1C) by day of the year in the URI GSO trawl survey.The average seasonal cycle is shown by a loess smoother with a span of 0.5, and its 95% confidence interval (gray).The bottom temperature gradient is positive (warmer in the mid-bay) between mid-March and mid-October before inverting for the late fall and winter months the seasonal cross-shelf patterns in bottom temperature remained the same (Fig. 5).During the spring (fall) migration of the studied species, the Rhode Island nearshore is the coldest (warmest) region and Narragansett Bay is similar in temperature to the mid-shelf.Therefore, ingressing summer species will experience the least suitable (coldest) temperatures in the Rhode Island nearshore as a result of winter mixing of the water column (Li et al. 2015) during their cross-shelf migration.As fall egress begins, however, warm summer surface water mixed to the bottom by fall storms (Li et al. 2015) will be slow to cool and thus this region will gradually become more suitable than Narragansett Bay.The same suitability pattern would hold for winter species ingress and egress, but with the seasons reversed.Therefore, this region may act as a hurdle to migrating species during coastal ingress but provide suitable temperatures during the initiation of egress.Notably, however, the relative magnitudes of seasonal warming trends in Narragansett Bay were the opposite of what has been documented on the shelf (Thomas et al. 2017).This suggests that differential warming rates across space throughout the year may create heterogeneity in the evolution of temperature conditions of seasonal habitats.
As has been suggested for a number of other fish and squid species (e.g.Sims et al. 2001, 2004, Peer & Miller 2014, Kovach et al. 2015, van Walraven et al. 2017), our results indicate that changes in temperature significantly affected observed migration phenology.Significant changes in the timing of first and last observation of each species were all in the hy pothesized directions, whereby summer species exhibited earlier ingress and later egress and winter species exhibited the reverse, matching global patterns of predictability of phenological re sponses across taxa (Root et al. 2003, Poloczanska et al. 2013).The average magnitude of these changes (ingress: 0.58 d yr −1 ; egress: 0.77 d yr −1 ) was also in line with previously observed values in fish (van Walraven et al. 2017) and other organisms (Root et al. 2003, Sydeman et al. 2015), albeit with high variance.
Atlantic herring, tautog, northern searobin, and fourspot flounder did not exhibit a significant change in their residence times between the first and last 10 yr of the URI GSO trawl survey time series, although the residence times of Atlantic herring and fourspot flounder appeared to display larger shifts between these intervals.While tautog and Atlantic herring met the abundance thresholds to be included in this analysis, neither is well-sampled by bottom trawl gear (Jech & Sullivan 2014, Vidal et al. 2018).The majority of the annual catch of fourspot flounder, which exhibited a significant shift toward earlier ingress (Table 3), and tautog occurs in early summer (Fig. 2), after which observations are less frequent.Tautog are known to spawn during this time period (Vidal et al. 2018) and, based upon collections of young-of-the-year individuals in late fall in the URI GSO survey and in other surveys of the southern New England region (Sullivan et al. 2000), it may be the spawning season for fourspot flounder as well.Therefore, it is possible that variable catchability or a link between detectability and spawning behavior affected the phenology observations of these 3 species.If so, it could explain why no environmental covariates were selected in the regression fits of Atlantic herring and tautog.Northern searobin were also primarily caught during the spring and early summer before becoming less abundant in Narragansett Bay, as observed by McBride & Able (1994) it is possible that they are unrelated to temperature and therefore their phenology did not shift over time.
It may be that all 4 of these species did not shift their migration phenology due to characteristics of their thermal preferences, but further investigation is required in order to account for the influence of abundance, catchability, and behavior on the results found here.
Accompanying shifts in residence time were alterations in the realized thermal niche of species in Narragansett Bay.In all significant cases, the minimum observed temperature decreased in summer species and increased in winter species, counter to our expectation of stability in realized thermal niches.Interestingly, the increases in the minimum temperatures at which the winter species were observed were approximately 0.5−1.0°Cor more above the measured increase in minimum annual temperature and winter mean temperature in Narragansett Bay.Given that observed warming in the mean and minimum winter temperatures was similar between the mid-and outer-bay, this suggests that these species contracted the extent of their habitat use away from the mid-bay, where the coldest bottom temperatures are observed, as their residence in the estuary declined.Changes in the maximum observation temperature, meanwhile, were only observed in summer species and were similar to the increases recorded in summer mean and annual maximum temperatures in the estuary.This pattern agrees with the findings of Sandblom et al. (2016) and Stuart-Smith et al. (2017) that organisms in shallow marine waters have flexible thermal minima (summer species), but inflexible maxima (winter species).Furthermore, the majority of species left Narragansett Bay at temperatures equal to or warmer than those at which they were first observed in a given year (Table 2).This differential could be due to migration being influenced by conditions along the route traveled as observed in birds (Haest et al. 2018), or that temperatures in the Rhode Island nearshore were as or more suitable than Narragansett Bay at the time of egress.
Supporting the hypothesized influence of conditions along the migration route, temperatures across the shelf were frequently significantly correlated to the observed ingress timing in Narragansett Bay (Table 4).Similarly, the temperature of the Rhode Island nearshore was significantly correlated to the timing of the last observations of longfin squid and summer flounder, suggesting the thermal gradient between the estuary and the shelf influenced migration as reported by Sims et al. (2004).It is possible that such a gradient influences the migration timing of other species as well, but the temporal mismatch between the available temperature data on the shelf (mid-October) and egress timing (mid-October to mid-December, Fig. 2) may have made such relationships difficult to detect.This mismatch, and the use of seasonal average temperatures, could also explain why the estimated gradient covariates between the shelf regions were not significant in any cases.Additionally, lagged values of the RSWT index were correlated to migration phenology in 2 species.In the case of summer flounder, the lagged correlation with ingress timing agrees with a recently proposed link between the effects of slope water dynamics on shelf temperature and juvenile survival (O' Leary et al. 2019) and research suggesting young individuals are first to migrate into the coastal zone (Langan et al. 2019).While the cause of such patterns cannot be definitively identified from the available data, they do suggest that temperature may also act indirectly on migration through effects on the life cycle.
Temperature appears not to be the only process affecting migration.Seasonal precipitation positively influenced the length of residence time (through delayed egress) in 3 species, consistent with the hypothesis that migration timing may be influenced by precipitation-driven changes in prey availability.Yet, fall precipitation negatively influenced the egress timing of fourspot and summer flounder.Additional work is required to fully understand these links between precipitation and phenology.The crossshelf temperature patterns suggest that ingress would not occur as observed if temperature was the only driver.For example, recent observations made at the Ocean Observatories Initiative Coastal Pioneer Array (https://oceanobservatories.org/array/coastalpioneer/), located at the southern New England shelfbreak, indicate that bottom temperatures at the overwintering habitat of summer species do not begin to substantially increase until at least mid-April.The initiation of outer-shelf spring warming near to when these species are first documented in the coastal zone suggests that the onset of warming may not cue migration.Further, the low thermal suitability of the Rhode Island nearshore in spring for summer species and fall for winter species would prevent each group from crossing this region at the recorded migration times if movement was solely influenced by temperature.It is possible that an additional cue, such as photoperiod (Otero et al. 2014), may trigger migration, but temperature significantly modifies the actual timing and/or rate of movement.Although photoperiod cannot explain the egress phenology of the studied species exiting Nar-ragansett Bay, the potential role of the Rhode Island nearshore as a temporary thermal refuge may complicate our ability to characterize the initiation of egress from the coast zone.Other local physical and ecological characteristics of all seasonally occupied habitats are likely to play a role as well and merit further investigation.Other species, such as dolly varden Salvelinus malma, have been shown to synchronize migrations with those of their prey (Sergeant et al. 2015, Bell et al. 2017).While such links were outside the scope of this analysis, migration timing is likely correlated with trophic interactions that may also have differential effects across life stages.
During the residence period, the studied species altered their spatial distributions in response to thermal gradients between the mid-and outer-bay sites (15.3 km apart) in a similar manner to observations of smelt made by Power & Attrill (2007).For most species for which at least 10% of observations occurred at each of the URI GSO trawl stations, the catch between stations shifted with the difference in bottom temperature in the expected direction based on their seasonal classifications.However, the frequent observations of each species at both survey stations despite the presence of a temperature gradient indicates that temperature does not play an absolute role in setting distribution patterns.Other dynamics, such as habitat preferences or competition, appear to modify the response of individual organisms within their niche as suggested by Soberón & Nakamura (2009).When compared to migration phenology and the evolution of the estuarine thermal gradient throughout the year, the documented responses indicate that ingressing species (summer and winter) favor the mid-bay relative to their mean habitat distribution.Depending upon species-specific temperature preferences, each then shifts toward the outer bay as the residence period progresses.The species would then be able to 'feel' and respond to the temperature of the Rhode Island nearshore as a result of estuary−shelf exchange (Pfeiffer-Herbert et al. 2015) at the end of the residence period, supporting the influence on migration of shifting thermal gradients between Narragansett Bay and waters just offshore.
The findings of this first-of-its-kind study have many parallels to investigations conducted in the spatial domain.Observations of marine taxa suggest that their spatial distributions have changed in response to climate change, climate fluctuations, ocean circulation, and interannual variability (Nye et al. 2009, 2011, Lucey & Nye 2010, Henderson et al. 2017, Day et al. 2018), generally resulting in range expansions in warm-water species and con-tractions in cold-water species over time (Nye et al. 2009, Lucey & Nye 2010).The present work shows that warm-water species are expanding their temporal range in a coastal temperate habitat, while cold-water species are contracting their residence periods, in response to the same suite of climate processes.Phenological change appeared to be greatest at the edges of the temporal distribution, suggesting that the effects of environmental processes, and potentially other variables such as competition, intensify at the extremes, as noted for spatial distribution by past work (Boudreau et al. 2015).This similarity is important because research of climate responses in marine fish and invertebrates has primarily focused only on the spatial aspect of their ecology, without explicitly decoupling spatial and temporal signals in the data.On the NEUS shelf, for example, monitoring of shifts in spatial distribution and management of marine fisheries have re lied upon the seasonal (Spring/Fall) Bottom Trawl Surveys conducted by the NEFSC.The observed phenological patterns in Narragansett Bay suggest both that the NEFSC Bottom Trawl Surveys are conducted during the cross-shelf migrations of many species and that those migrations have shifted in time relative to the surveys.It is possible that the observed spatial distributions of the survey catch will have shifted significantly, with portions of the distribution residing in the under-sampled nearshore zone (ASMFC 2002), due only to changing migration timing.
In this manner, changes in migration phenology of the magnitude observed here may have profound effects on marine fisheries.In addition to altered species availability impacting the results of the fisheriesindependent surveys used to develop management measures and assess spatial distributions, shifts in seasonal patterns of habitat occupancy may influence population productivity (Bal et al. 2017, Manhard et al. 2017, Staudinger et al. 2019).Because it may be difficult to distinguish if population trends are due to changes in abundance or availability, this could create significant uncertainty in the management process.This is particularly important for species that are regulated by spatiotemporal management measures, in the form of spatial allocation of harvest or spatiotemporal harvest restrictions.Summer flounder, for example, are managed using seasonal closures and quota allocations to individual or groups of states along the NEUS shelf (Terceiro 2018).If changes in migration phenology alter the ability to assess the spatial abundance patterns used to evaluate quota distribution or impact the effective-ness of seasonal closures in restricting harvest, as has been observed for other species (Peer & Miller 2014), complications could arise in a very valuable fishery with a history of controversy between the fishing industry and fisheries management (Terceiro 2018).Furthermore, shifts in species availability to the fishery can also have significant economic consequences in cases where there is a seasonality in market demand or seafood processing capacity (Mills et al. 2013, Staudinger et al. 2019).In order to prevent such consequences from unaccounted changes in migration phenology, it is necessary that further effort is devoted to characterizing these patterns and incorporating them into assessments of fish stocks.The current paradigm of treating spatial and temporal responses to environmental forcing as separable in fisheries-independent data must change in order to successfully anticipate the effects of continued climate change on the composition and productivity of marine populations and the fisheries that they support.
The results of this analysis suggest that marine species respond to temperature gradients across their seasonal habitats and migratory routes, resulting from complex interactions of climate processes acting at several spatial and temporal scales.The unique nature of the URI GSO trawl survey both enabled the novel analyses presented here and limits the ability to compare the results to other ecosystems.However, observations of synchronous life history phenologies across subpopulations of sole (Fincham et al. 2013) and salmonids (Kovach et al. 2013(Kovach et al. , 2015) ) signal that the patterns observed in this estuary may be similar across the coastal zone of the NEUS shelf.The identified environmental correlations, meanwhile, suggest that the spatial extent of these similarities will be limited to the region in which the drivers of spatiotemporal temperature gradients on the shelf are similar.Given the significant potential impacts on ecosystem function and fisheries production, further efforts to understand spatiotemporal responses to climate by marine species will prove critical in providing insights to scientists, fishery managers, and industry stakeholders confronting marine climate change around the globe.

Fig. 1 .
Fig. 1. (A) Northeast US (NEUS) Shelf, with southern New England indicated by a black box.The 50, 90, and 150 m depth contours are shown as gray lines.(B) Defined southern New England shelf regions with respect to Narragansett Bay (NB), Rhode Island (RI).The 50, 90 and 150 m depth contours are shown as gray lines.(C) Locations of the northerly mid-bay (orange arrow) and southerly outer-bay (purple arrow) stations sampled by the University of Rhode Island Graduate School of Oceanography (URI GSO) weekly trawl survey in Narragansett Bay (1959−2016).The length and orientation of the arrows represent the distance and direction trawled interpolated grids developed by Friedland et al. (2019) from measurements recorded by the Spring and Fall Bottom Trawl Surveys conducted by the Northeast Fisheries Science Center (NEFSC) were used to calculate the seasonal mean temperatures in 3 regions: (1) the outer shelf (69.75−73.0°W,≤ 41°N, 90−150 m), (2) the midshelf (69.75−73.0°W,≤ 41°N, 50−90 m), and (3) the Rhode Island nearshore (70.5−72°W, > 41°N) (Fig. 1).

Fig. 3 .
Fig. 3. Annual observed residence times (points) of the 12 studied species in the URI GSO trawl survey.Summer (red) and winter (blue) species groupings are designated by the point colors.The central trends in species residence times are represented by a loess smoother (black line) with a span of 0.5 and its 95% confidence interval (gray)

Fig. 4 .
Fig. 4. Abundance distribution by (A,B) day of the year and (C,D) temperature for longfin squid (A,C) and ocean pout (B,D).Longfin squid expanded its temporal and thermal niche over the URI GSO trawl time series, while ocean pout contracted both and was eventually extirpated from the ecosystem (diamonds denote single observations).Day 1 for longfin squid represents 1 March, and Day 1 for ocean pout represents 16 August.There is a gap in the 1964 temperature distribution of ocean pout due to missing winter measurements

Fig. 5 .
Fig. 5. (A) April and (B) October bottom temperatures from the URI GSO fish trawl survey and the interpolated grids of Friedland et al. (2019) across the defined shelf regions (see Fig. 1B) during the periods of 1959−1987 (blue) and 1988−2016 (red).The median of each spatiotemporal group is designated by a horizontal black line.Lowercase letters indicate statistically significant differences among locations within each season identified by a 2-way ANCOVA (Tukey's HSD, p < 0.05).The average bottom temperatures during the latter period were 0.58°C warmer in April and 0.77°C warmer in October (Tukey's HSD, p < 0.05), but no significant differences in warming among locations were detected between halves of the time series Table1for species binomials) were highly abundant during winter in the first several decades of the URI GSO trawl survey before becoming very rare by the late 1990s.The 1998 season was the last in which more than a single individual was observed, making this the terminal year used in analysis.Striped searobin, on the other hand, first appeared in 1963 and was only observed sporadically until the late 1970s, after which it maintained a consistent summer presence in Narragansett Bay.Observations

Table 3
. Changes in the day of the first and last annual observation and annual minimum and maximum recorded bottom temperature (°C) of each species between the first and last 10 yr was observed in the University of Rhode Island Graduate School of Oceanography (URI GSO) trawl survey.For the changes in day of observation, positive (negative) values indicate a shift later (earlier) in the year.P-values for the t-tests on changes in day and temperature

Table 2 .
Mean temperature (°C) of the first and last observations of each studied species.Comparison of the temperatures of the first and last obser vation in each year was made with a 1-sample t-test on the measured dif ferences.Change in observated temperatures during the URI GSO trawl survey time series was assessed by comparing the first and last 10 yr of ob servations for each species with a 2-sample t-test.P-values for the t-tests on changes in day and temperature are designated as: ns: p > 0.05, *p < 0.05, **p < 0.01, ***p < 0.001

Table 5 .
, and shifting their distribution toward the outer-bay station of the URI GSO trawl survey.While it is unclear what drives these patterns of habitat use, Results of generalized additive model fits on daily distribution of observations of each species between the 2 stations of the URI GSO trawl survey.Positive effects of the temperature gradient, corresponding to a greater proportion of observations coming from the mid-bay station, are represented by '+,' while negative relationships are designated by '−.'