Shift in the larval phenology of a marine ectotherm due to ocean warming with consequences for larval transport

Because environmental temperature has an important influence on developmental rate and physiology, marine ectotherms are vulnerable to phenology changes due to ocean warming. Identifying changes to phenology, the timing of biological events, and understanding their effect on recruitment and abundance is of critical importance to establish potential population effects. We examined the larval phenology of the commercially important Norway lobster (Nephrops norvegicus) and used a larval transport model to examine its effect on simulated transport patterns. Using a model to estimate annual larval release dates based on temperature‐dependent embryo incubation, an earlier shift of 17.2 d occurred between 1982–1995 and 2000–2010 in the Irish Sea, similar to an observed empirical shift in phenology of 19.1 d using historical zooplankton data sets. Despite this earlier phenology, temperature‐dependent pelagic larval durations were unchanged because the water column to which larvae were released earlier had also warmed. Larval transport simulations in the western Irish Sea indicated that the phenology shift had minimal effects on larval retention and advection distance overall, because major variations were observed only at very early or late stages of the larval season, that is, times when lower proportions of larvae were present. As the western Irish Sea grounds exports small but consistent quantities of larvae to nearby populations, especially off Scotland, it may act as an important source of larvae, especially when retention of native larvae is low. Overall, larval transport tools may indicate grounds that are periodically vulnerable to recruitment failures and offer potentially valuable information in fishery management.

Major alterations in phenology, or the seasonal timing of biological events, have been documented in marine, limnological, and terrestrial habitats worldwide (Walther et al. 2002;Root et al. 2003;Edwards and Richardson 2004). Global sea temperatures in the upper ocean have warmed by 0.11 C decade −1 between 1971 and 2010 (Rhein et al. 2013), while in the Northeast Atlantic, surface temperatures have undergone an increasing trend of 0.1-0.5 C decade −1 from 1983 to 2012 (Dye et al. 2013). Since the physiology of ectotherms is closely coupled to environmental temperature (Wear 1974), variations in sea temperature have been linked to changes in the timing of spawning, larval hatching, and migrations of marine species including flatfish, crustaceans, and squid in the North Atlantic (Sims et al. 2001;Richards 2012;Fincham et al. 2013). Identifying altered seasonalities and understanding how variations in phenology affect abundance and recruitment is of critical importance to establish potential negative population effects, particularly to commercially harvested species. For example, in the Wadden Sea, increased winter temperatures influenced the recruitment of the bivalve Macoma balthica by lowering reproductive output, while earlier spawning led to a timing mismatch with the phytoplankton bloom and increased predation risk (Philippart et al. 2003). Interestingly, the reverse may also be true: decreasing temperatures were implicated in delayed spawning in Baltic cod in the mid-1990s, coinciding with lower than expected recruitment (Wieland et al. 2000). Despite several examples of phenology shifts in the marine environment (Poloczanska et al. 2013), less is known on how phenology changes can affect recruitment, even in commercially harvested species.
Ocean warming can lead to a contraction of embryo incubation periods and an earlier release of meroplankton larvae. Larval phenology changes can potentially affect mortality and settlement patterns of larvae. A timing mismatch with optimal food levels may result in increased likelihood of starvation (Cushing 1990). Due to seasonal oceanographic features, early-vs. late-hatching larvae may be subjected to contrasting conditions that influence advection patterns. In springhatching larvae, temperature increases over the larval season so early hatchers encounter relatively cooler temperatures, therefore, we might expect these to develop slower and take a longer time to settle, potentially increasing dispersal away from suitable habitat and reducing local retention. In contrast, we expect retention to be highest when larvae are released late in the season, with higher temperatures leading to faster development and earlier settlement. So, does this mean that an earlier shift in larval phenology, for example, earlier release of larvae into cooler waters, could prolong development and lead to lower retention? To examine this, we also must consider the relative change in temperature at the sea bottom compared to the surface, which determine benthic and pelagic warming impacts, respectively. Bottom temperature affects the timing of larval release in certain benthic species due to its influence on embryonic development, while surface temperature affects pelagic duration of surface-dwelling planktonic larvae. Hence, the rate at which surface and bottom temperature change relative to one another affects the processes that influence larval retention or loss.
Along with biological factors (hatching time, larval duration, and behavior), transport in the pelagic larval stages is dependent on physical factors such as ocean circulation patterns which will affect recruitment to the adult population (Ospina-Alvarez et al. 2018;Blanco et al. 2019). Considerable interannual variability can arise in oceanographic transport patterns (Espinasse et al. 2017), therefore we need to understand how important this interannual variation is for larval retention within grounds/habitat, or exports of larvae to neighboring grounds. A spatial influence according to release location may also occur, if the grounds where larvae are produced is large, larvae being released from one area may have a higher likelihood of retention than those from another area. Unfortunately, although the pelagic phase is crucial, both for survival of larvae and enabling exchange of larvae between distant populations, larval transport is difficult to observe in situ due to the dispersal of vast numbers of miniscule larvae in very large water bodies. Larval transport models have been developed to address this difficulty and identify the processes that influence retention and connectivity (Metaxas and Saunders 2009).
Larval recruitment in certain crustacean populations is tied to the confines of suitable habitat, which make these ideal cases to examine how changes to larval phenology affect recruitment success via larval retention. Nephrops norvegicus is a small lobster and commercially important throughout its range in the Northeast Atlantic and Mediterranean, with European landings of 44,000 t valued at approximately 360 million EUR in 2016 (EUROSTAT, ec.europa.eu/eurostat/ web/fisheries/data/database). The spatial distribution of benthic juvenile and adult N. norvegicus is patchy due to their dependence on suitable muddy sediment to construct burrows. Adults spend much of their time within or close to their burrows, but the larval phase is pelagic: larvae hatch into the water column after first experiencing a period of temperaturedependent embryo incubation during brooding by adult females (Dunthorn 1967;Farmer 1974). Spawning, that is, extrusion of eggs to the female abdomen, takes place in late summer/early autumn and in the Irish Sea, embryo incubation is 7-9 months (Farmer 1974), whereas in warmer regions such as the Mediterranean it is approximately 6 months (Mori et al. 1998). Experiments by Farmer (1974) showed a 50% reduction in embryonic development duration as a result of a 10 C increase in temperature. In a second temperaturedependent process, the planktonic larvae pass through three developmental stages with pelagic larval duration (PLD) lasting between 1 and 2 months at temperatures of 8.5-13 C (Smith 1987;Thompson and Ayers 1989;Dickey-Collas et al. 2000) and larval transport is heavily influenced by local oceanographic conditions (O'Sullivan et al. 2015;Phelps et al. 2015). Larvae inhabit the upper 40 m of the water column and they perform a vertical migration over the diel cycle (McGeady et al. 2019). When larval development is complete, settlement on suitable muddy sediment is essential for survival and recruitment to the population. Subsequent development of the benthic juvenile is influenced by temperature, food availability, sediment type, and population density (Bell et al. 2006).
The Irish Sea is a semi-enclosed waterbody between Ireland and Great Britain, where ocean circulation is mainly driven by tides; however, relatively weak tidal currents in the western Irish Sea promote stratification during spring and summer (Simpson 1971). Mud patches in the western and eastern Irish Sea provide habitat to N. norvegicus populations, which are assessed by the International Council for the Exploration of the Seas (ICES) as separate functional units (Ungfors et al. 2013). The western Irish Sea population of N. norvegicus supports a highly productive fishery with yields of 5000-10,000 t annually and stock size has remained stable over the past 17 yr (Lundy et al. 2019). It has high adult densities, although individuals are smaller in size relative to other grounds due to density-dependent growth suppression Merder et al. 2020). Most larvae are present from April to June in the Irish Sea (Nichols et al. 1987). A seasonal near-surface gyre in the western Irish Sea was thought to aid recruitment to the population by acting as a retention mechanism (Hill et al. 1996). However, the effectiveness of the gyre in retaining larvae has since been shown to be significantly reduced because larvae perform a diel vertical migration, limiting daily exposure to gyral flows (Phelps et al. 2015;McGeady et al. 2019). Due to the potentially large numbers of larvae produced in the western Irish Sea each year, the grounds may act as an important, but variable (O'Sullivan et al. 2015), source of larval exports to other grounds, and exchange of larvae may be important to sustain the viability of individual patches/populations in the long term (Levins 1970).
In the present study, we tested the hypothesis that ocean warming has led to a contraction of the N. norvegicus incubation period and resulted in earlier larval phenology (release of larvae) in the Irish Sea. We tested this hypothesis using two complementary approaches. The first used experimental data to create a predictive model of larval release timing and was fed with empirical temperature data. The second examined the timing of the N. norvegicus larval season empirically using historical zooplankton data sets. We also tested the hypothesis that a shift in larval phenology affects the transport patterns of pelagic larvae. This was tested using a larval transport model to examine the effect of release date on larval retention, advection distance (including annual and spatial influences), exports, and PLD. A final aim was to identify grounds which are recipients of larvae from the western Irish Sea, as these grounds may depend on larval exchange to maintain their populations.

Temperature-dependent larval release date estimation
Two empirical temperature data sets were acquired from the British Oceanographic Data Centre (bodc.ac.uk/data/ bodc_database/, accessed 28 June 2019). The data sets were both part of the Isle of Man Long-term Environmental Time Series project, initiated by the University of Liverpool's Port Erin Marine Laboratory and taken over by the Isle of Man Government Laboratory in 2006. The first data set (Port Erin surface temperature) was used to examine long-term surface temperature trends in the Irish Sea, data were obtained from 1904 to 2010 at the Port Erin Breakwater, Isle of Man (54.08522, −4.76805; Fig. 1). From 1904 to 2006, a sample of water was collected from the sea surface twice daily, once in the morning and afternoon, and temperature was recorded using a thermometer. From 2006 to 2010, daily averages were calculated using logged data taken 2 h either side of high water. For the second data set (Cypris station bottom temperature), temperature data were available at five depths (0, 5, 10, 20, and 37 m) from 1954 to 2010 at the Cypris station (54.09167, −4.83333; Fig. 1) approximately 5 km west of Port Erin. Data were collected on a weekly to monthly basis depending on boat availability, season, and weather. Temperature data collected from 37 m (sea bottom) were used to examine long-term bottom temperature trends in the Irish Sea. The Cypris station is situated approximately 6 km from the north eastern edge of the western Irish Sea N. norvegicus grounds, making it a good approximation of bottom temperature on the grounds (Fig. 1).
To investigate temperature as a predictor of larval phenology, average temperature for the period 15 th August, that is, the approximate start of incubation (Farmer 1974), until 31 st March, that is, the approximate date that marks the end of incubation/start of peak larval hatching (Nichols et al. 1987) was used to estimate larval release timing. Incubation duration was estimated by applying temperature data from the Port Erin surface temperature time series to regression parameters ( Fig. 2) from previous laboratory studies describing the relationship between temperature and embryo incubation (Dunthorn 1967;Farmer 1974). Daily surface temperature data from the Port Erin Breakwater sampling station was used to calculate incubation duration for 1905-2010, due to the lack of adequate temporal coverage in the Cypris station bottom  temperature data set. However, inspection of Cypris station bottom temperature data indicated that the water column is mostly mixed from September to March meaning little difference (< 0.3 C) between surface and bottom temperatures during the incubation period (Fig. 3b,d). Between 2011 and 2018, surface temperatures from the Hybrid Coordinate Ocean Model (HYCOM) from a point source at the Port Erin Breakwater sampling site was used, due to a lack of empirical data and good agreement between HYCOM and Port Erin surface temperature was observed (Supporting Information Fig. S1). Therefore, the annual date of temperature-dependent larval release was defined as the estimated embryo incubation duration in days ( Fig. 2) from 15 th August, that is, warmer incubation temperatures led to an earlier larval release date estimate.

Larval phenology
To examine the larval phenology of N. norvegicus in the Irish Sea with empirical observations, zooplankton data sets were obtained from the CEFAS data hub (https://www.cefas. co.uk/cefas-data-hub/, accessed 17 April 2017). Data on N. norvegicus larvae (Stage I, II, and III) were available for 16 yr between 1982 and 2010 (Table 1). Due to individual research cruises having distinct objectives, temporal coverage in sampled years did not always overlap. The temporal range for the entire data set was day of year (DOY) 12-180 (i.e., January-June). Larval densities were extracted from 5177 samples collected on 60 research cruises (Fig. 1). Samples were considered part of the western or eastern Irish Sea grounds, respectively, if collected west or east of −4.5 longitude. Samples were collected from several vessels using a variety of plankton samplers, mesh size was 250-270 μm except in 1993, when an 800 μm mesh was used (Table 1). Plankton samplers were fitted with a CTD and flowmeter to record temperature, salinity, depth, and water volume filtered. Nets were deployed on a double oblique tow from the surface to within 2 m of the seabed. Hauling speeds were adjusted to ensure equal sampling of all depths. Upon recovery, nets were washed down and samples were preserved in a buffered 4% formaldehyde solution. Samples were transported to the laboratory where N. norvegicus larvae were sorted and staged. Larval densities were calculated as abundance m −2 , where the number of larvae in each sample was divided by the volume of water filtered and multiplied by the sampled depth.

Larval transport
To examine the influence of larval phenology on retention, advection distance, exports, and PLD of larvae, including annual and spatial influences, a biophysical larval transport model was used. The HYCOM model (Chassignet et al. 2007) was used to simulate oceanographic conditions during the N. norvegicus larval season (i.e., February-July when larvae were observed in empirical zooplankton data sets). HYCOM is a hybrid isopycnal coordinate model, that is, isopycnal in the open stratified ocean, but terrain-following in shallow coastal areas, and has fixed depths in the mixed layer and/or unstratified seas. The model is data assimilative and receives data from satellite observations and in situ temperature and salinity from expendable bathythermographs, Argo floats, and moored buoys. HYCOM model output was available at a 3-hourly temporal resolution and 1/12 th degree spatial resolution. HYCOM output was coupled with particle-tracking software (Ichthyop v3.3; Lett et al. 2008) to simulate the transport of N. norvegicus larvae released from the western Irish Sea grounds (Fig. 1). Simulations were conducted for each of the years that HYCOM output was available, that is, 25 yr from 1994 to 2018. For each simulated year, 5000 particles were released into the ocean domain every 5 d from DOY 46-166 (i.e., 25 release dates), resulting in a total of 3,125,000 released particles (125,000 per year). Particles were released from three spatial sections on the grounds, that is, western (−6.20 to −5.58 longitude), center (−5.58 to −5.26), and eastern (−5.26 to −4.76) zones, in order to examine spatial influences on larval transport. A diel vertical migration was incorporated to model larval behavior in the water column, where particles made an ascent to 10 m at dusk, a midnight sink to 20 m followed by a subsequent ascent and another descent to 20 m at dawn, following previous empirical observations ( McGeady et al. 2019). Simulations were forced using a Runge-Kutta 4 th order numerical scheme. Particle positions were calculated every 5 min and recorded every 3 h.
The PLD was estimated using parameters from laboratory studies showing the temperature-dependent developmental rate of N. norvegicus pelagic larval stages I-III (Smith 1987;Thompson and Ayers 1989;Dickey-Collas et al. 2000) combined in Dickey-Collas et al. (2000). The temperature experienced by each particle was used to determine its PLD, and when the PLD had elapsed, the particle was considered ready to settle. Settlement positions were used to calculate larval retention (%) on the western Irish Sea grounds, as well as exports to other nearby grounds. Particles settling outside of suitable habitat were considered "lost." In addition, the straight-line distance from initial release location to settlement position was calculated as advection distance (km).
To examine the influence of phenology on larval retention and advection distance, generalized linear models (GLMs) were used. Larval retention (%) and advection distance (km) for each release event (25 yr × 25 release dates × 3 release zones = 1875 release events) were used as response variables in the respective models. Categorical predictor variables were: year, release day (grouped into categories: , release zone on the grounds (west, center, and east) and their one-way interactions. Variance inflation factors (VIFs) were calculated to check for collinearity between predictors. It was not necessary to remove any of the predictors due to VIFs of < 10 (Montgomery et al. 2012). Initially a GLM with larval retention expressed as a proportion as the response variable and a binomial error distribution was used; however, due to overdispersion, a quasibinomial distribution was used. Advection distance as a response variable was log-transformed to attain residual normality and used in a GLM with a Gaussian error distribution. A simple linear regression was also used to examine the effect of release day on PLD.

Warming temperature trend
Surface temperature at Port Erin Breakwater displayed an increasing annual trend from 1904 to 2010. The maximum annual temperature (11.7 C) was observed in 2007 and 9 of the 10 warmest years were from 1998 to 2010 (temperature anomalies relative to the 1904-2010 average are shown in Fig. 3a). Annual temperature was above the 1904-2010 average for 14 consecutive years from 1997 to 2010 (Fig. 3a). For the time series 2000-2010, average annual surface temperature was warmer by 0.8 C compared to 1982-1995 (Fig. 3b). June (1.1 C), May (1.0 C), and September (1.0 C) had the largest increases and each month increased by a minimum of 0.6 C except for December (0.3 C; Fig. 3b). Bottom temperature at the Cypris station from 1954 to 2010 also showed a warming trend (Fig. 3c). Due to gaps in sampling dates, the data set was split into periods of 5 or 6 yr to ensure adequate coverage of each month. A sharp increase occurred in the late 1980s/early 1990s and reached a peak in the early 2000s (Fig. 3c). An annual bottom temperature increase of 0.9 C was observed between the earlier (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) and later (2000-2010) time series, with the largest increases of > 1.1 C in the months May, August, September, and October (Fig. 3d).

Temperature-dependent larval release date estimation
The estimated date of larval release, calculated using temperature-dependent incubation rates (Fig. 2) and the Port Erin surface temperature times series, showed an increasingly earlier trend over time between 1905 and 2018 (i.e., negative correlation p < 0.001, R 2 = 0.23) at a rate of 1.5 d decade −1 (Fig. 4a). This trend was stronger from 1982-2010 (p < 0.001, R 2 = 0.56), when larval release timing showed an earlier trend of 10.1 d decade −1 (Fig. 4a). Using the same temperaturedependent incubation relationship, a significant difference (Student's t-test, t = 5.8, p < 0.001) was observed between average larval release dates for 1982-1995 (DOY 119.8 AE 2.2, mean AE SE) and 2000-2010 (102.6 AE 1.8; Fig. 4b), resulting in an estimated difference of 17.2 d. The 17.2 d earlier larval release date was associated with a 0.9 C increase in temperature during the incubation period.

Larval phenology
In total, 85,654 N. norvegicus larvae were observed in 5177 zooplankton samples collected in the Irish Sea. In the western Irish Sea, the earliest date that N. norvegicus larvae were observed was 31 January 2006 (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) and later (2000-2010) time series. A well-distributed spatial coverage of sampling locations across the western and eastern Irish Sea grounds was available for each time period (Fig. 5). It was evident that high densities appeared earlier in the larval season on both grounds in the 2000-2010 time series (Fig. 5). This pattern of earlier peak densities in 2000-2010 was also apparent for Stage II (Supporting Information Fig. S2) and Stage III larvae (Supporting Information Fig. S3). Stage I larvae in the western Irish Sea for DOY 85-104 were more abundant (Mann-Whitney U-test, n 1 = 238, n 2 = 221, p < 0.05) in 2000-2010 (4.86 AE 0.70 abundance m −2 AE SE) compared to 1982-1995 (1.43 AE 0.35 abundance m −2 ) indicating an earlier release of larvae (Fig. 5). Meanwhile, larval abundance for DOY 125-144 was significantly higher (Mann-Whitney U-test, n 1 = 472, n 2 = 120, p < 0.001) in the 1982-1995 time-series (5.81 AE 0.44) than in the 2000-2010 time series (2.14 AE 0.46), indicating that a high proportion of larvae had already progressed to the next stage of development in 2000-2010 (Fig. 5). An earlier larval release was also apparent in the eastern Irish Sea as Stage I abundance was higher (Mann-Whitney U-test, n 1 = 122, n 2 = 348, p < 0.001) from 2000 to 2010 (0.47 AE 0.09) than from 1982 to 1995 (0.05 AE 0.02) for DOY 85-104 (Fig. 5). Larvae also progressed to the next stage earlier as lower densities were observed for DOY 124-145 (Mann-Whitney U-test, n 1 = 132, n 2 = 208, p < 0.001) in the later 2000-2010 (0.09 AE 0.04) compared to the earlier 1982-1995 (0.79 AE 0.28) time series. The SP of Stage I larval abundance from empirical observations was estimated as DOY 121.9 for 1982-1995 and DOY 102.8 for 2000-2010. Therefore, the SP of Stage I N. norvegicus larval abundance was 19.1 d earlier in the later time series. This empirical observation closely matched the difference in estimated temperature-dependent larval release date between the two time series (17.2 d-see above).

Larval transport
Simulations were completed for 25 yr in the western Irish Sea between 1994 and 2018 with virtual larvae released from three zones on 25 dates spaced 5 d apart from DOY 46 (15 th February) to DOY 166 (15 th June). Larval retention and advection distance was significantly influenced by release year, day, spatial zone, and their one-way interactions (p < 0.001 in each case; Table 2). Average annual larval retention was 24.2% and fell below 20% only in 2002 (14.2%), 2014 (17.7%), 2016 (15.8%), and 2018 (17.7%; Fig. 6a). Year explained 14.7% and 18.8% of deviance in retention and advection distance, respectively, indicating that interannual variability is an important factor (Table 2). A sustained period of high retention (> 27%) occurred from 2005 to 2009 and reached a maximum of 31.9% in 2006. However, a gradual decline was observed from 2009 to 2012, although high year-on-year variation was also apparent (2013-2018), showing a pattern of high retention followed by low retention until the end of the time series (Fig. 6a).
The day of release (p < 0.001) explained 12.8% deviance in both retention and advection distance (Table 2). Retention was lowest (17.4%) when larvae were released early in the year between DOY 46 and 66 (15 th February-7 th March) and increased substantially when released late in the year (32.0%) between DOY 146 and 166 (26 th May-15 th June; Fig. 6b). Although a trend was observed with gradually higher larval retention as the season progressed, little difference (< 3.5%) was observed between DOY 71 and 141 (12 th March-21 st May).
Spatial zone of release also affected retention and advection distance respectively by 5.4% and 8.9% (p < 0.001; Table 2), that is, retention decreased, and advection distance increased when larvae were released from the east of the grounds. Releases from the west had 28.5% retention, on average, but on the western Irish Sea grounds estimated using temperature-dependent incubation durations with surface temperature from Port Erin Breakwater (triangles) and HYCOM (circles). Lines indicate significant trend for 1905-2018 (red) and 1982-2010 (green) and shading specifies 95% confidence intervals. (b) Boxplot with difference in estimated temperaturedependent larval release date between the time series 1982-1995 and 2000-2010. Central lines in boxplots represent the median, box extremities indicate 1 st and 3 rd quartiles and whiskers specify range. decreased to 23.4% in the center zone and further reduced to 20.8% to the east (Fig. 6c). The interaction between year and day of release explained the most deviance in retention and advection distance (27.3% and 32.1%, respectively), indicating that interannual and intraseasonal variation are very important factors in larval transport (Table 2; Fig. 6d). In most years, retention improved with a later release (Fig. 6d), 15 of 25 yr had the highest retention in the latest release period, DOY 146-166. By contrast, lowest retention was observed in the early season (DOY 46-66) for 13 of 25 simulated years.
The interaction between year and zone also had a significant influence on retention (p < 0.001) and advection distance (p < 0.001; Table 2). Larvae released from the western zone had the highest retention for 21/25 yr and only had the lowest retention in 2001, 2010, and 2014; in these years, the eastern zone had highest retention. Release day and zone interaction explained only 1.4% deviance in retention and 0.4% in advection distance (p < 0.001; Table 2). The difference in retention between western and eastern releases was greatest for the earliest release dates 8.7%) in comparison to the latest releases 4.6%).

Exports from western Irish Sea
The western Irish Sea grounds exported larvae to several nearby grounds to the north, east, and south (Fig. 7a). The Minch grounds received the highest average exports (2.0%) followed by the Clyde and Jura (0.8%), eastern Irish Sea (0.3%), Celtic Sea (0.2%), and Stanton grounds (0.1%) (Fig. 7). Exports were subject to large interannual variability, with the Minch grounds receiving the highest exports in 1996, 2002, and 2018, yet receiving no exports in 2004 and 2007 (Fig. 7b). An earlier release date was associated with higher exports, maximum exports to the eastern Irish Sea, Minch and Celtic Sea grounds occurred between DOY 46 and 91 (Fig. 7c). Larvae released to the east or center of the grounds also had a higher likelihood of export to the Clyde and Jura, eastern Irish Sea, Minch and Stanton grounds, while those originating from the western zone were less likely to be exported (Fig. 7d).
With an overall average of 3.4% of larvae exported and 24.2% retained, the majority of larvae failed to settle on suitable habitat and were "lost" (72.4%). Due to high retention, the proportion of "lost" larvae was lowest late in the season, DOY 146-166 (66.9%) compared to early, DOY 46-66 (77.1%). This means that more larvae were lost early-on, when retention was lower and exports to other grounds were at their highest (Fig. 7c).

Discussion
The results of the present study showed a marked shift in the larval phenology of N. norvegicus in the Irish Sea and linked this to temperature. We used experimental data to create a predictive model of larval release timing using temperature data from Port Erin which agreed well with empirically observed shifts in larval phenology from historical datasets in the Irish Sea. Larval transport simulations indicated that larvae released very early in the season had lower retention compared to later releases; however, limited change in retention was apparent between larval phenology shift dates. Instead, Table 2. Analysis of deviance for GLM examining the associations between the response variables larval retention (%) and log advection distance (km) and the predictor variables release year, day, and zone from larval transport simulations on the western Irish Sea grounds. large interannual variability in hydrography-forced larval transport was apparent, which can result in strong changes to retention of larvae within grounds, and can alter exports from the high-density western Irish Sea population to neighboring grounds.

SS
Temperature was considered a major contender as a causal factor in the larval phenology shift due to its importance in embryonic development (Dunthorn 1967;Farmer 1974). Our predictive model of larval release date based on temperaturedependent incubation rates predicted a 17.2 d earlier shift between 1982-1995 (DOY 119.8 AE 2.2) and 2000-2010 (DOY 102.6 AE 1.8) due to a 0.9 C increase in temperature. This was very close to the empirical observation from historical larval data sets of a 19.1 d shift in the SP of Stage I larval abundance between 1982-1995(DOY 121.9) and 2000-2010, and provides support for the hypothesis that temperature is the causal mechanism. This phenology shift was observed in all three developmental stages and in both the western and eastern Irish Sea (Fig. 5, Supporting Information Figs. S2, S3). The empirical larval phenology shift coincided with a bottom temperature increase of 0.9 C at the Cypris station between both time series (1982-1995 and 2000-2010); with the largest temperature increases (> 1.1 C) occurring in August, September, and October, that is, at the early stages of embryo incubation (Fig. 3d). Eiríksson (2014) observed an increased ratio of posthatching : prehatching ovigerous N. norvegicus females in May from the 1970s-1990s to the 2000s, suggesting that the larval hatching season has also advanced in Icelandic waters. Earlier phenology due to temperature has similarly been recorded in other marine systems, for example, earlier timing of the biomass maximum of copepods in the North Pacific (Mackas et al. 1998) and the hatch timing of northern shrimp in the Northwest Atlantic (Richards 2012). However, far less is understood how ocean warming-driven phenology shifts affect recruitment dynamics.
Bottom temperature has increased decadally in the Irish Sea (Fig. 3c,d) and brought larval release dates earlier (Figs. 4, 5). This raises the question of whether temperature-dependent PLDs and transport patterns have also been altered? The surface temperature minimum occurs in February/March in the Irish Sea, after which time, temperature increases (Fig. 3b). If pelagic larvae are released earlier (due to a contraction of the incubation period as a result of increased bottom temperatures) when surface temperatures are cooler, then PLDs may be extended as development is slowed at lower temperatures, potentially resulting in increased dispersal leading to lower larval retention. Despite this expectation, the Port Erin surface temperature time series showed that temperatures on the empirical larval phenology shift dates, that is, DOY 121.9 (1982-1995) and DOY 102.8 (2000-2010 is relatively unchanged. Average surface temperature from Port Erin for the 1982-1995 time series was 8.49 C at DOY 122 and 7.75 C at DOY 103 (0.74 C apart), for the 2000-2010 time series it was 9.63 C at DOY 122 and 8.6 C at DOY 103 (1.03 C apart). Therefore, despite an earlier larval release, larvae would have experienced similar temperatures at the time of hatching, likely resulting in an unaltered larval duration. Regardless, in future, larvae may increasingly hatch when temperatures are at their lowest. The annual minima for surface and depthaveraged temperatures in the Irish Sea in February/March are projected to get later by 12 d from 1980 to 2100 (Olbert et al. 2012). Such a scenario could see increased overlap between a later-shifting annual temperature minima and an earlier-shifting larval season, with larvae increasingly hatching when surface temperatures are at their lowest, potentially prolonging development and increasing losses of larvae. The Irish Sea is forecasted to increase in surface and depthaveraged temperature by 1.89 C and 1.79 C, respectively, from the 1980s-2090s, with higher increases projected for autumn and winter (embryo incubation) than for spring (larval season; Olbert et al. 2012).
As HYCOM model availability was restricted to 1994-2018, it was not possible to simulate larval transport for the time series 1982-1995 and directly compare it to 2000-2010, when the earlier larval phenology occurred. Instead, we simulated the effect of larval timing across all available years (1994Fig. 6d) to examine a general effect of release date on larval transport. Differences in retention were greatest at opposite ends of the larval season and more subtle in the intervening periods. For example, an extremely late larval release (DOY 146-166) resulted in 14.6% higher retention, 51.6 km lower advection distance, and 30.4 d shorter PLD than an extremely early release (DOY 46-66; Fig. 6b). Extreme early hatchers, such as larvae observed on 31 January 2006 in empirical data, would have a lower likelihood of retention and higher dispersal potential than late-hatching larvae (Fig. 6b). Thus, larvae released on the fringes of the larval season, be it the first or last hatchers, would be most affected by a phenology shift. Nevertheless, the significance of this is reduced since the relative abundance of larvae in "fringe" periods is low. We can examine what the effect on the majority of larvae would be by applying the 19.1 d earlier shift, that is, from DOY 121.9 to DOY 102.8, which results in limited changes to retention and advection distance from larval transport simulations in the western Irish Sea. Thus, this study found that the empirical shift in phenology did not result in a pronounced reduction in larval retention.
More significant for the future may be the fact that larval retention was subject to a high degree of interannual variability, particularly toward the end of the time series (Fig. 6d). Although average retention (24.2%) in the western Irish Sea was generally high across the time series (Fig. 6a), this is not the case on all grounds (O'Sullivan et al. 2015). The enclosed nature of the Irish Sea and relatively weak currents, particularly near the coast, proved conducive to retention. In the western Irish Sea, vast quantities of larvae produced from a large population coupled with high retention may contribute to the stable and high density adult stock observed in underwater TV surveys (Lundy et al. 2019). In contrast, on isolated grounds with low densities and less potential for larval imports, interannual variability may have a significant effect on recruitment. For example, the relatively isolated Aran grounds off the west coast of Ireland has undergone stock fluctuations in recent years, with a decreasing density trend, and low modeled retention levels may render it vulnerable to periods of low recruitment (McGeady et al. 2019). Interannual variability has also been recognized as an important contributor to the transport patterns and ultimately recruitment of ichthyoplankton larvae (Ospina-Alvarez et al. 2015;Kvile et al. 2018).
Retention was highest for releases from the west of the western Irish Sea grounds and lowest to the east. Currents were weak in coastal areas to the west, resulting in short dispersal distances and enhanced retention. Releases to the east, particularly for earlier releases of larvae, were prone to advection through the North Channel by strong northward currents resulting in transport to grounds off the west of Scotland (Fig. 7). For this reason, the spatial distribution of spawning females, which determines the "zonation" of larval release densities on the grounds, is another factor affecting larval transport patterns. Adult densities are not uniform across the grounds and the 2019 underwater TV survey indicated that highest burrow densities are to the east; however, spatial distribution is subject to interannual variation (Lundy et al. 2019). Additionally, a spatial gradient of 1.3 C in bottom temperature was observed across the grounds during the incubation period: lowest temperatures occurred in deeper channels to the center, potentially resulting in prolonged incubation and later larval releases, while higher bottom temperatures near the coasts of Ireland and the Isle of Man could lead to an earlier release. Overall, if most spawning females are near the coasts where temperatures are higher, the majority of larvae would be released earlier. This temperature gradient is projected to become stronger in the future as coastal areas warm faster than deeper channels of the Irish Sea (Olbert et al. 2012). In theory, fishing pressure could indirectly affect the timing of larval release by altering the distribution of females on the grounds. For example, if fishing pressure was focused near the coast, resulting in diminished densities and proportionally higher densities to the center, where incubation temperatures are lower, a shift in the timing of larval release from earlier to later, and spatially from coastal to central zones, could occur. In order to promote a long larval release that protects against poor retentive conditions in the short term, a harvesting strategy that promotes a spread of female N. norvegicus across the temperature gradient of the grounds may be beneficial.
In terms of oceanography, the most favorable conditions for larval retention were low current velocities and high temperatures, leading to shorter PLDs. In scenarios with high current velocities, larvae could potentially be moved great distances in a short length of time. Therefore, despite shorter PLDs, associated with a later release, unfavorable oceanography led to larvae being transported away from suitable habitat. No consistent seasonal differences were observed in current speeds as a high degree of interannual variability was evident and instances of poor retention were recorded both early (2017) and late (2012) in the season (Fig. 6d). Circulation patterns in the Irish Sea are predicted to change in the future, stronger southward currents and a strengthening of the western Irish Sea gyre as a result of sharpened stratification (Olbert et al. 2012) may result in increased southward advection of larvae.
Empirical observations of larvae beyond the boundaries of the western Irish Sea grounds (Fig. 5) demonstrated the clear potential for larvae to be transported to other grounds. Simulations showed a significant number of larval exports to large grounds off the west coast of Scotland, particularly the Minch grounds, even though this was not the most adjacent area. The eastern Irish Sea and Celtic Sea grounds also received small but consistent exports (Fig. 7). Exports were highest early in the season, and larvae released from the east of the grounds were more likely to be exported to northern grounds due to their proximity to the North Channel, where strong northward currents often persisted. Due to its high-density biomass and large spatial extent, the number of reproducing females on the western Irish Sea grounds is vast. Therefore, larvae originating on the grounds being exported elsewhere could be very important to the annual larval settlement and recruitment of other grounds, particularly if retention is poor. The relatively high level of interground connectivity may be the reason that N. norvegicus in the Atlantic Ocean exhibit low levels of genetic differentiation, reflecting this species' capacity for high gene-flow (Maltagliati et al. 1998;Stamatis et al. 2004). Some genetic differences have been detected between Atlantic and Mediterranean populations, however (Gallagher et al. 2019).
Although a first estimation, larval retention may be a good proxy for estimating recruitment in N. norvegicus and other species. Other ways to measure recruitment are problematic: burrows of juveniles are not easily quantified in underwater TV surveys (juveniles are also known to share burrows with adults; Tuck et al. 1994), length-frequency distributions on research surveys or in fisheries catches tend to be rather unimodal and hence difficult to separate into year-classes (Farmer 1973), and there is no reliable aging method (Sheridan et al. 2016). Since successful recruitment depends on larvae settling on suitable habitat, hydrography plays a particularly important role due to this species' potential for larval transport beyond the boundary of mud habitat which is critical for survival (Bell et al. 2006). In the absence of other data, estimating retention and receipt of larvae from other grounds by modeling larval transport could provide an index of recruitment. Indeed, a strong recruitment relationship was demonstrated in other systems, for example, in the Mediterranean, between modeled anchovy larvae reaching a nursery area and densities of age-0 individuals from acoustic surveys (Ospina-Alvarez et al. 2015). Further improvements to the N. norvegicus larval retention index may include applying a release schedule that accurately reflects the timing of peak larval release, using the spatial distribution of adults across the grounds to determine the initial distribution of hatching larvae and estimating larval mortality before settlement. While potentially useful as an index, natural mortality due to parasites, predation or starvation, and so on (Farmer 1975;Cushing 1990) could also be significant before juveniles recruit to the fishery. However, larval retention models could potentially be used as an early warning signal of low recruitment due to spells of poor retention and/or low larval immigration. Early identification of poor recruitment could be accounted for in the management process of commercial species by reducing fishing mortality early-on, particularly in cases lacking alternative sources of recruitment data.
The observed phenology shift in the Irish Sea has relevance for N. norvegicus populations across its range which may have also been subjected to similar changes to phenology. Modeled larval retention may be useful as a proxy for recruitment in N. norvegicus and other species with a strong requirement for suitable habitat to survive. Such a retention/recruitment index could be highly beneficial in fisheries management and future work should examine whether modeled larval retention is linked to stock size estimates.