Abstract
Warming temperatures caused by climate change have the potential to impact spawning phenology of temperate marine fish as some species have temperature-dependent gonadal development. Inter-annual variation in the timing of Atlantic cod (Gadus morhua) spawning in the northern North Sea, central North Sea and Irish Sea was estimated by calculating an annual peak roe month (PRM) from records of roe landings spanning the last three decades. A trend towards earlier PRM was found in all three regions, with estimates of shifts in PRM ranging from 0.9 to 2.4 weeks per decade. Temperatures experienced by cod during early vitellogenesis correlated negatively with PRM, suggesting that rising sea temperatures have contributed to a shift in spawning phenology. A concurrent reduction in the mean size of spawning females excluded the possibility that earlier spawning was due to a shift in size structure towards larger individuals, as large cod spawn earlier than smaller-sized individuals in the North Sea. Further research into the effects of climate change on the phenology of different trophic levels within the North Sea ecosystem should be undertaken to determine whether climate change-induced shifts in spawning phenology will result in a temporal mismatch between cod larvae and their planktonic prey.
Introduction
Sea temperatures are warming at unprecedented rates as a result of global climate change (IPCC, 2013). These changing environmental conditions have the potential to impact marine animal and plant populations, including fish, through changes in distribution, phenology, and physiology (Hughes, 2000). Although climate change-induced shifts in phenology are now well documented for terrestrial plant and animal species (e.g. Roy and Sparks, 2000; Roberts et al., 2015), research focusing on the impact of climate change on phenology of marine fish is more limited. Phenology can impact the fitness of marine fish by influencing whether life-history events such as spawning occur under optimal conditions. For example, the match/mismatch hypothesis (Cushing, 1975) relates survival rates of fish larvae to the match in timing between a “critical period” of larval first-feeding (Hjort, 1914) and the peak of spring plankton blooms in temperate marine environments (Cushing, 1990). Cod larvae depend on copepod prey for survival and changes in the plankton community during the main period of larval cod presence has been linked to variation in recruitment, highlighting this as an important mechanism of bottom-up control (Beaugrand et al., 2003). Changes in spawning phenology therefore have the potential to impact fish recruitment, if the phenology of prey and predator become decoupled (e.g. Platt et al., 2003).
The physiological mechanism through which temperature influences spawning time of Atlantic cod (Gadus morhua, hereafter referred to as cod) has been well described through experimental studies. Vitellogenesis is initiated at the autumnal equinox but the subsequent rate of oocyte development is influenced by temperature (Kjesbu et al., 2010). Warmer temperatures result in a faster rate of oocyte development and therefore an earlier spawning time (Kjesbu et al., 2010). Changes in temperature during early vitellogenesis have a greater effect on spawning time than changes in temperature during late vitellogenesis (Kjesbu, 1994). In addition to the effects of temperature, a physiological link between body size and spawning time can also contribute to the determination of cod spawning phenology (Kjesbu et al., 2010). Larger cod generally spawn earlier than smaller individuals in warmer environments such as the North and Irish Seas, possibly due to their higher thermal sensitivity (Kjesbu et al., 2010).
Despite the potential for the reproductive phenology of commercially exploited temperate marine fish stocks to be influenced by climate change, with potential implications for fitness and recruitment, the majority of research on the response of fish populations to rising temperatures has so far focused on changes in spatial distribution and productivity (e.g. Clark et al., 2003; Perry et al., 2005; Brander, 2007; Dulvy et al., 2008; Engelhard et al., 2011). One exception is the recent study by Fincham et al. (2013) which examined 40 years of maturity-stage data for sole (Solea solea), collected through an intensive monitoring scheme of fisheries landings from waters around the United Kingdom. A shift in spawning phenology of ca. 1.5 weeks per decade across four populations of sole was detected, and this shift correlated with rising winter sea surface temperature (SST) (Fincham et al., 2013). Previous research on the impacts of temperature on spawning time of cod has also revealed a relationship between spawning time and temperature in wild cod off the coast of Newfoundland (Hutchings and Myers, 1994) and in the Baltic Sea (Wieland et al., 2000). Subsequent research investigating the link between spawning time and sea temperature for cod in the North Sea did not find evidence to support this hypothesis, but these results may have been constrained by small sample sizes and a restricted time series (Morgan et al., 2013).
Investigations into phenological shifts in marine fish populations are uncommon due partly to the limited availability of long-term datasets with a high degree of intra-annual temporal resolution. The resolution of scientific sampling programmes is limited both temporally and spatially (e.g. Carscadden et al., 1996; Greve et al., 2005) although exceptions do exist (e.g. Hutchings and Myers, 1994; Fincham et al., 2013). There is potential for systematic, industry-based data sources to provide long-term datasets for commercial species. For example, catch records from the fishing industry have been used to successfully explore variation in the timing of annual migrations of tuna (Thunnus thynnus) (Dufour et al., 2010) and mackerel (Scomber scombrus) (Jansen and Gislason, 2011). The collection of industry data relevant to determining spawning time is hindered by the fact that, unlike sole, in Scotland and Ireland gadoids such as cod are gutted at sea, precluding the direct assessment of maturity status when landed. However, for species for which there is a roe market there is potential for fisheries statistics on timing and quantity of roe landings to provide long-term data relating to spawning time. As cod roe is of highest commercial value when cod ovaries are ripe or vitellogenic (Katsiadaki et al., 1999) peaks in roe landings can be used to infer timing of spawning. Pederson (1984) demonstrated how the relationship between cod and cod roe landings could be used to estimate the timing of spawning of Barents Sea cod. However, the use of fisheries data can introduce bias into analysis which must be accounted for when interpreting results. In particular, variation in fishing effort (e.g. Dufour et al., 2010) and fluctuations in demand for products such as roe (e.g. Pederson, 1984) can potentially influence the results. A detailed understanding of the fishery is therefore required when using fisheries-dependent data to investigate the ecology of an exploited species.
In the North Sea, sea temperatures have been warming at a rate of 0.2–0.6 °C per decade over the past thirty years (Baudron et al., 2014). North Sea temperatures have been increasing at rates faster than global averages, with summer temperatures increasing at particularly rapid rates (Belkin, 2009). The North Sea is therefore well-suited to investigating how warming temperatures are affecting reproduction dynamics of wild fish stocks. Although synchronous changes in individual growth rate across a varied range of fish species in response to climate change have already been demonstrated in the North Sea, North Sea cod was a notable exception (Baudron et al., 2014). North Sea cod from different regions can tolerate a wide range of temperatures (Righton et al., 2010) and generally display site fidelity to historic spawning sites (e.g. Neat et al. 2014). They will not necessarily migrate to cooler waters that are available even when experiencing temperatures that are apparently super-optimal (Neat and Righton, 2007). This evidence suggests that adult cod in the North Sea may not yet have made large-scale distributional shifts to cooler waters. Given that there is little evidence for putative impacts of warming temperatures on growth and distribution of North Sea cod the possibility of a climate change-induced shift in cod spawning phenology deserves investigation.
The objectives of this study were firstly to investigate whether spawning phenology of cod in the North Sea has shifted over time, and secondly whether shifts in spawning phenology are related to warming trends in sea temperature. Based on the experimental work by Kjesbu et al. (2010), cod spawning phenology was hypothesized to have shifted to earlier in the year in response to warming temperatures caused by climate change due to the physiological link between rate of oocyte maturation and water temperature. To test this, roe landings from the North Sea to Scotland recorded over a 30-year period were used to calculate peak roe month (PRM) which was used as a proxy for the commencement of spawning. PRM was calculated from peaks in roe landings and from peaks in a roe index (RI), an aggregate measure of roe production that is conceptually analogous to a gonadal somatic index (GSI) estimated at the individual-level. Comparable data were also available for the Irish Sea, a region where cod stocks already experience ambient temperatures at the upper end of the thermal tolerance of the species (Brander, 1995) and which are therefore expected be negatively impacted by future warming (Drinkwater, 2005). These data were included in analysis, so that the potential shifts in phenology of cod from different regions could be compared. A survey of demersal fishermen was carried out to test the validity of the assumption that roe landings were composed mainly of cod. The overall rates of change in PRM over the 30-year period were calculated through regression analysis and compared across the upper North Sea, central North Sea and Irish Sea. Trends in average autumn (corresponding to early cod vitellogenesis) and average winter (corresponding to late cod vitellogenesis) SST in the three regions were investigated, and correlations between PRM and seasonal temperature trends were explored. Finally, inter-annual trends in size of spawning females over the study period were estimated, to investigate whether a concomitant change in size of cod could potentially confound the interpretation of the rate of change in spawning time.
Methods
Data sources
Landings data
Records of the landed weight of roe (LWR) to Scottish ports from the North Sea for the years 1985–2014 were obtained from the Marine Scotland Fisheries Statistics of the Scottish Government. Roe landings data were available as the total monthly landings to each port. This restricted the temporal resolution of subsequent analysis to month. The International Council for Exploration of the Sea (ICES) fishing area where the roe originated was the only available information relating to the location of catches. The roe landed to Scotland from the North Sea originated from either ICES fishing area IVa, in the northern North Sea, or from ICES fishing area IVb, the central North Sea (hereafter referred to as IVa and IVb) (Figure 1). Roe data from IVa and IVb were analysed separately, in an attempt to distinguish between the separate populations of cod which spawn in the northern and central North Sea (Neat et al., 2014). Data on live weight of cod landed monthly to Scottish ports from IVa and IVb were obtained from the same source and used to calculate the RI described below.
The location of ICES fishing grounds IVa and IVb in the North Sea, and VIIa in the Irish Sea, from which records of landings of roe were available.
Records of daily LWR per vessel from ICES fishing area VIIa (the Irish Sea, hereafter referred to as VIIa) (Figure 1) were obtained from the Centre of Environment, Fisheries and Aquaculture Science (Cefas). Records were available for the years 1982–2014. Although there is evidence suggesting that the cod spawning in the eastern and western Irish Sea comprise two separate populations (Armstrong et al., 2001), it was not possible to distinguish between separate cod stocks due to the spatial resolution of the data. Although landings data were available per day for VIIa, the index of spawning phenology was calculated to resolution of month for consistency with the North Sea data.
The annual weight of roe landed from all regions has decreased over time (Figure 2). This decrease in roe landings may be attributable to a combination of declining market for roe, a reduction in cod quotas over the study period (ICES, 2015a), and particularly in the case of the Irish Sea due to the decline in cod abundance (ICES, 2014). A sharp reduction in LWR from VIIa is noticeable in 2000, when the Irish Sea cod fishery was closed (ICES, 2005) (Figure 2). In the North Sea, monthly LWR from IVa was consistently higher than from IVb (Figure 2).
The total annual LWR, landed during months January to April from regions IVa (a) and IVb (b) and during February to May from region VIIa (c).
Information on the species of roe landed was not recorded in the landings dataset, so an assumption was made that all of the roe landed was cod. To reduce the potential contamination from other species, only roe data from the months during which cod are likely to be maturing and spawning were used. Consequently, only roe landings from December to April were extracted from the North Sea dataset, and landings from January to May were extracted from the Irish Sea dataset, to encompass the spawning seasons of North Sea cod and Irish Sea cod, respectively (Brander, 1994). The data extraction excluded herring (Clupea harengus) roe, which spawn between August and October (Rogers and Stocks, 2001). Roe from haddock (Melanogrammus aeglefinus), ling (Molva molva), saithe (Pollachius virens), whiting (Merlangius merlangus), and hake (Merluccius merluccius) could potentially be mixed in with cod roe given that the spawning time of these species overlaps with cod (Rogers and Stocks, 2001). Of these species, only haddock roe could be landed in large enough quantities to contaminate the dataset significantly.
Questionnaires
To examine the validity of the assumption that cod dominate the landed roe in the months used for the analysis (see above), a questionnaire-based survey of North Sea demersal fishermen who land to Scottish ports was conducted. The survey was undertaken at a skipper’s exposition in Aberdeen (27 and 28 May 2016), an event attended by ca. 89 demersal fishermen, many of whom are based in Scotland and fish regularly in the North Sea. Although a fully randomized design could not be implemented, this event nevertheless provided a convenient opportunity to question many active and retired fishermen from across the country. Of the attendees approached, 32 fishermen who had landed roe from the North Sea to Scottish ports between the years 1985 and 2014 agreed to complete a questionnaire. The questionnaire consisted of a mix of closed- and open-ended questions (Supplementary Material S1). Participants were asked a series of simple questions to assess the relevance of their experience to this study. They were then asked which species of roe they landed, and what proportion of roes landed were cod. Answers gathered from the survey were used to qualitatively inform interpretation of results that are based on assumptions regarding the species composition of the roe landings data.
Temperature data
As a demersal species, sea bottom temperature would be most appropriate for investigating the impacts of temperature on cod phenology. However, a complete sea bottom temperature data set which covered the full spatial and temporal scale of this study was not available. The North Sea, which can be described as stratified throughout autumn in some areas, becomes well-mixed in winter (Connor et al., 2006). During autumn, the stratification of the Irish Sea breaks down completely, resulting in a well-mixed water column (Connor et al., 2006). SST therefore provides an approximation of temperatures at deeper depths in these regions during autumn and winter, assuming that the temporal trend in SST is the same as the temporal trend in sea bottom temperature. SST has been widely used as an indicator of ocean climate in studies relating aspects of cod ecology to temperature (e.g. Planque and Frédou, 1999; Armstrong et al., 2004; Engelhard et al., 2014). Therefore, SST was considered appropriate here as a relevant proxy for the ambient temperature experienced by cod.
Median monthly SSTs were extracted from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST1) (http://www.metoffice.gov.uk/hadobs/hadisst/data/download.html). The HadISST1 dataset is maintained by the UK Met Office and is based on opportunistic sampling from merchant ships and research vessels, but with extensive processing applied to minimize the potential impact of sampling biases (Rayner et al., 2003). This dataset provides complete global coverage of monthly SST estimates on a 1° latitude-longitude grid, dating back to 1871 (Rayner et al., 2003). For each year of this study, mean autumn (September, October, November) and mean winter (December, January, February) SST estimates were calculated for each fishing area, to match the spatial resolution of the roe landings data (see Figure 1). Temperatures during these seasons were assumed to represent the ambient thermal environment experienced by cod during early (autumn) and late (winter) vitellogenesis.
Length data
Data on the length and maturity stage of North Sea cod were obtained from the North Sea international bottom trawl survey (NS-IBTS) co-ordinated by ICES. These data were downloaded from the DATRAS database (https://datras.ices.dk) as a Sex Maturity Age Length Key dataset. Data from Quarter 1 surveys (conducted during January-February) were used to estimate the length of spawners. Data for female fish classed as maturing, mature, spawning, or spent, and which were sampled from ICES sampling area 1 (ICES, 2015b) were extracted from the DATRAS dataset. This restrictive definition of spawners excluded male and immature fish, and fish which were caught outside of fishing area IVa. Lengths of spawners in fishing area IVb were obtained using the same dataset. As IVb crosses the boundaries of several sampling areas, and encompasses no single area in its entirety, the lengths of spawners from sampling areas 2, 3, 4, 6, and 7 (ICES, 2015b) were used to infer average length of spawners in IVb. There were no comparable length data available for Irish Sea cod which covered the extent of the study period.
Indices of spawning phenology
Peak roe month
The calculation of PRM assumes that annual peaks in LWR are related consistently to spawning phenology over time. Total monthly roe landings were explored graphically to confirm that they formed a roughly bell-shaped distribution each year. It is possible that temporal variation in fishing effort may influence the timing of peaks in roe landings. To address this potential source of bias, a RI was calculated (described below) and used as an alternative weighting factor in estimation of PRM.
Roe index
The use of RI as a weighting factor for calculating PRM assumes that peaks in RI are related consistently to phenology of spawning. Cycles in the analogous GSI can be used to indicate spawning cycles in fish populations over time (Jons and Miranda, 1997). RI increases as adult cod direct energy allocation to gonad growth prior to spawning, peaks around the commencement of spawning, and decreases after the initiation of spawning as roes become unsuitable for human consumption (Pederson, 1984). The RI was therefore explored graphically to check that it followed a roughly bell-shaped distribution each year.
Statistical analysis
Trends in PRM
Before temporal trends in the timing of PRM were analysed, the presence of serial correlation in PRM values was tested for using the autocorrelation function (ACF) in R Core Team (2015). The value of the ACF at different time lags provides an indication of whether there is any autocorrelation in the data. No significant autocorrelation (p > 0.05) was detected in the time series of PRM calculated for each region, using either RI or LWR as a weighting factor. Temporal trends in PRM for each region were then investigated through linear regression. The slope fitted to the regression of PRM against year provided estimates of rate of change in months per year. These rates were used to calculate rate of change in spawning time in weeks per decade, so that results could be compared with results from similar research on sole (Fincham et al., 2013).
Regional differences in PRM, and in temporal trends in PRM, were analysed through analysis of covariance (ANCOVA), with year, region and the interaction between year and region included as explanatory variables. Analyses were carried out with PRM weighted by RI for IVa and IVb, and by LWR for IVa, IVb, and VIIa.
To test whether PRM weighted by RI and PRM weighted by LWR provided similar estimates of temporal trends in spawning time, an ANCOVA was fitted to compare the difference in PRM, and trends in PRM over the years, when calculated using the two different weighting factors. The model was fitted with PRM as response variable, year, weighting factor, and the interaction between year and weighting factor as explanatory variables.
Trends in temperature
Prior to analysis of temporal trends in temperature, the datasets of seasonal SST averages were tested for presence of autocorrelation using the ACF function in R. As no significant autocorrelation was detected (p > 0.05), analysis of temporal trends was conducted using simple linear regression. Linear regressions were fitted for each region with average seasonal SST as the response variable and year as the explanatory variable. This analysis was conducted separately for autumn and winter. The slopes estimated from the linear regressions were used to calculate the rate of change in temperature over time. To test whether directional trends in SST over time differed significantly between regions, an ANCOVA was carried out with seasonal SST as the response variable and year, region and the interaction between year and region as fixed effects. This analysis was repeated for each season.
Relationship between spawning phenology and temperature
To test for correlations between PRM and autumn and winter SST, regressions were carried out with PRM as the response variable and autumn or winter SST as the explanatory variables. As the seasonal SST averages encompass a range of temperature estimates across each fishing area over 3 months there is considerable error associated with each average which simple linear regression would not account for. The exact location and thermal experience of sampled cod is unknowable, so including the error associated with estimates of seasonal SST averages in the model accounts for some of the uncertainty associated with this variable. Deming regression, a form of type II regression analysis, was therefore carried out. In a Deming regression model, the regression line is fitted by minimizing the sum of squares for the perpendicular distance between the data points and the regression line, therefore producing a single line regardless of which of the two variables is used as the response variable (Cornbleet and Gochman, 1979). The advantage of Deming regression analysis is that error in the explanatory as well as the response variable can be included in the model (Cornbleet and Gochman, 1979). The estimates of s.e. calculated for each estimate of average seasonal SST and PRM were incorporated into the model, using the “deming” package in R. The significance of the slope estimate was assessed through inspection of the 95% CIs (Nakagawa and Cuthill, 2007). These analyses were repeated for each region.
Trends in spawner length
Trends in average body length of female spawners were analysed through linear regression, with fish length as the response variable and year as explanatory variable. This analysis was repeated for IVa and IVb.
All statistical analyses were carried out using R v3.2.5 (R Core Team, 2015).
Results
Trends in PRM
There was a negative relationship between PRM and year in all three regions when using either LWR (IVa: F1,28 = 45.96, p < 0.001; IVb: F1,27 = 60.09, p < 0.001; VIIa: F1,31 = 18.92, p < 0.001) or RI (IVa: F1,28 = 11.56, p = 0.002; IVb: F1,27 = 17.65, p < 0.001) as the weighting factor (Figure 3).
The relationship between PRM and year in regions IVa, IVb, and VIIa, when weighted by either LWR or RI. Regression lines are included to indicate a significant (p < 0.05) trend in PRM over time, based on linear regression analysis. Standard error bars are included to illustrate the s.e. associated with each estimate of PRM, but were not included in regression analyses.
Neither estimates of average PRM (IVa: F1, 56 = 1.243, p = 0.27; IVb: F1,54 = 0.001, p = 0.97), nor trends in PRM over the years (IVa: F1,54 = 0.589, p = 0.45, IVb: F1,54 = 1.52, p = 0.22) differed significantly when weighted by either weighting factor, suggesting these weighting factors are interchangeable for calculating PRM using this dataset. However, the rate of change of PRM calculated for both IVa and IVb was slightly higher when weighted by LWR than when weighted by RI (Table 1).
Average PRM and rate of change in PRM over the study period, for each region and calculated using each weighting factor (LWR and RI).
| Region | Weighting factor | Average PRM | Slope ± s.e. | Rate of change in PRM (weeks per decade) |
|---|---|---|---|---|
| IVa | LWR | 2.34 | −0.024 ± 0.0035 | 1.03 |
| IVa | RI | 2.41 | −0.019 ± 0.0056 | 0.81 |
| IVb | LWR | 2.11 | −0.057 ± 0.0073 | 2.44 |
| IVb | RI | 2.11 | −0.041 ± 0.0099 | 1.76 |
| VIIa | LWR | 3.21 | −0.022 ± 0.0051 | 0.94 |
| Region | Weighting factor | Average PRM | Slope ± s.e. | Rate of change in PRM (weeks per decade) |
|---|---|---|---|---|
| IVa | LWR | 2.34 | −0.024 ± 0.0035 | 1.03 |
| IVa | RI | 2.41 | −0.019 ± 0.0056 | 0.81 |
| IVb | LWR | 2.11 | −0.057 ± 0.0073 | 2.44 |
| IVb | RI | 2.11 | −0.041 ± 0.0099 | 1.76 |
| VIIa | LWR | 3.21 | −0.022 ± 0.0051 | 0.94 |
Presented is the slope and associated error estimated by the regression models, interpreted as an estimate of change in months per year. The rate of change in PRM in weeks per decade calculated from this slope is also presented.
Average PRM and rate of change in PRM over the study period, for each region and calculated using each weighting factor (LWR and RI).
| Region | Weighting factor | Average PRM | Slope ± s.e. | Rate of change in PRM (weeks per decade) |
|---|---|---|---|---|
| IVa | LWR | 2.34 | −0.024 ± 0.0035 | 1.03 |
| IVa | RI | 2.41 | −0.019 ± 0.0056 | 0.81 |
| IVb | LWR | 2.11 | −0.057 ± 0.0073 | 2.44 |
| IVb | RI | 2.11 | −0.041 ± 0.0099 | 1.76 |
| VIIa | LWR | 3.21 | −0.022 ± 0.0051 | 0.94 |
| Region | Weighting factor | Average PRM | Slope ± s.e. | Rate of change in PRM (weeks per decade) |
|---|---|---|---|---|
| IVa | LWR | 2.34 | −0.024 ± 0.0035 | 1.03 |
| IVa | RI | 2.41 | −0.019 ± 0.0056 | 0.81 |
| IVb | LWR | 2.11 | −0.057 ± 0.0073 | 2.44 |
| IVb | RI | 2.11 | −0.041 ± 0.0099 | 1.76 |
| VIIa | LWR | 3.21 | −0.022 ± 0.0051 | 0.94 |
Presented is the slope and associated error estimated by the regression models, interpreted as an estimate of change in months per year. The rate of change in PRM in weeks per decade calculated from this slope is also presented.
ANCOVA revealed regional differences in the average PRM when weighted by LWR (F2,86 = 150.44, p < 0.001) or RI (F1,55 = 9.61, p = 0.003), with average PRM of cod earliest in IVb and latest in VIIa (Table 1). The trends in PRM over time also varied between regions when weighted by LWR (F2, 86 = 11.77, p < 0.001) or RI (F1,55 = 3.98, p = 0.05), with shifts in PRM fastest in IVb when weighted by either factor (Table 1).
Roe composition
The results from the questionnaires supported the assumption that the majority of roes landed by Scottish fishermen were composed of cod. 94% of participants reported landing cod roe and 54% of the participants reported that >75% of the roe they landed was cod. Participants reported landing cod, haddock, hake, whiting, saithe, and ling roe. After cod, haddock was the most commonly landed species, as 75% of participants reported landing haddock roe. 21% of participants reported landing more haddock roe than cod roe. This suggests that haddock is the species most likely to contaminate the roe landings dataset. Of the respondents that reported landing more haddock than cod roe, 29% reported that they rarely landed any roe. This leaves only 12% of participants who claimed to land a considerable amount of haddock roe. Anecdotal evidence from discussions with fishermen suggested that some fishermen would not land haddock roe as it is smaller in size compared to cod roe (owing to smaller individual body size of individual haddock caught) and more difficult to handle than cod roe. Furthermore, in Scottish fish markets demand is highest for cod roe, creating an incentive for landing predominantly cod gonads (D. Anderson, Aberdeen Fish Producers Organization, pers. comm, January 2017).
Temperature trends
There was a significant, positive relationship between average autumn SST and year in all three regions (Figure 4). Average autumn SST increased at a rate of 0.35 °C per decade in IVa (F1,29 = 21.91. p < 0.001), 0.48 °C per decade in IVb (F1,29 = 20.86, p = 0.001) and 0.36 °C per decade in VIIa (F1,31 = 26.8, p < 0.001). ANCOVA revealed a significant effect of year (F1,88 = 55.10, p < 0.001) and region (F2,88 = 246.67, p < 0.001) on the average autumn SST. Across the study period, average SST was 11.28 °C in IVa, 12.86 °C in IVb and 13.64 °C in VIIa. There was no interaction between region and year (F2,88 = 0.68, p = 0.51) in relation to autumn SST, suggesting that the rates of autumn SST warming did not differ significantly between regions.
The trends in SST over time, for autumn and winter averages in regions IVa, IVb, and VIIa. Standard error bars are included to illustrate the s.e. for each of the calculated seasonal means, but were not included in regression analyses. Regression lines are included where there was a significant trend (p < 0.05) in the seasonal mean SST over time.
There was no relationship between winter SST and year in IVa (F1,30 = 2.094, p = 0.16) or VIIa (F1,32 = 2.75, p = 0.11) (Figure 4). There was a positive relationship between winter SST and year in IVb, with winter SST increasing at a rate of 0.24 °C per decade (Figure 4). Similarly to autumn SST, ANCOVA revealed a significant effect of year (F1, 95 = 6.062, p = 0.016) and region (F2, 95 = 163.871, p < 0.001) on winter SST. Across the study period, average winter SST was 7.94 °C in IVa, 7.44 °C in IVb and 9.42 °C in VIIa. There was no interaction between region and year in relation to winter SST (F2,95 = 1,134, p = 0.33), again suggesting no significant difference between the long-term temperature trends in each region.
Relationship between PRM and temperature
Owing to the lack of a significant difference between PRM weighted by LWR and PRM weighted by RI, analysis was conducted using PRM weighted by LWR to represent trends in spawning phenology. This allowed results to be compared across all three regions.
PRM correlated negatively with autumn SST in all regions (Figure 5, Table 2). No relationships between PRM and average winter SST were revealed for any region (Figure 5, Table 2). For regions and seasons where a significant relationship between PRM and average seasonal SST was observed, the slope estimates from regression analyses were used to calculate predicted shifts in spawning phenology in response to increases in SST (Table 2). The rate of change in spawning time in relation to autumn temperature rise ranged from 1.16 weeks per 1 °C change in autumn SST in IVa, to 3.92 weeks per 1 °C rise in autumn SST in VIIa.
The relationship between LWR weighted PRM and average seasonal SST. Dashed lines indicate significant relationships as calculated through Deming regression (95% CIs do not contain zero).
The relationship between PRM weighted by LWR and seasonal SST, analysed using Deming regression.
| Season | Region | Weighting factor | Deming slope estimate ± s.e. | 95% CI for Deming slope estimate (lower, upper) | Rate of change (week per °C) estimated from Deming slope |
|---|---|---|---|---|---|
| Autumn | IVa | LWR | −0.29 ± 0.059 | −0.41, −0.18 | 1.16 |
| IVb | LWR | −0.41 ± 0.14 | −0.67, −0.14 | 1.64 | |
| VIIa | LWR | −0.98 ± 0.30 | −1.57, −0.38 | 3.92 | |
| Winter | IVa | LWR | −0.041 ± 0.11 | −0.26, 0.18 | |
| IVb | LWR | −0.35 ± 0.36 | −1.05, 0.36 | ||
| VIIa | LWR | −1.55 ± 0.79 | −3.10, 0.01 |
| Season | Region | Weighting factor | Deming slope estimate ± s.e. | 95% CI for Deming slope estimate (lower, upper) | Rate of change (week per °C) estimated from Deming slope |
|---|---|---|---|---|---|
| Autumn | IVa | LWR | −0.29 ± 0.059 | −0.41, −0.18 | 1.16 |
| IVb | LWR | −0.41 ± 0.14 | −0.67, −0.14 | 1.64 | |
| VIIa | LWR | −0.98 ± 0.30 | −1.57, −0.38 | 3.92 | |
| Winter | IVa | LWR | −0.041 ± 0.11 | −0.26, 0.18 | |
| IVb | LWR | −0.35 ± 0.36 | −1.05, 0.36 | ||
| VIIa | LWR | −1.55 ± 0.79 | −3.10, 0.01 |
The slope estimate can be interpreted as rate of change in month per °C change in temperature. Statistically significant slope estimates (i.e. 95% CIs do not contain zero) are emboldened. For statistically significant relationships the regression slope for the model has been used to estimate a rate of change in weeks per °C.
The relationship between PRM weighted by LWR and seasonal SST, analysed using Deming regression.
| Season | Region | Weighting factor | Deming slope estimate ± s.e. | 95% CI for Deming slope estimate (lower, upper) | Rate of change (week per °C) estimated from Deming slope |
|---|---|---|---|---|---|
| Autumn | IVa | LWR | −0.29 ± 0.059 | −0.41, −0.18 | 1.16 |
| IVb | LWR | −0.41 ± 0.14 | −0.67, −0.14 | 1.64 | |
| VIIa | LWR | −0.98 ± 0.30 | −1.57, −0.38 | 3.92 | |
| Winter | IVa | LWR | −0.041 ± 0.11 | −0.26, 0.18 | |
| IVb | LWR | −0.35 ± 0.36 | −1.05, 0.36 | ||
| VIIa | LWR | −1.55 ± 0.79 | −3.10, 0.01 |
| Season | Region | Weighting factor | Deming slope estimate ± s.e. | 95% CI for Deming slope estimate (lower, upper) | Rate of change (week per °C) estimated from Deming slope |
|---|---|---|---|---|---|
| Autumn | IVa | LWR | −0.29 ± 0.059 | −0.41, −0.18 | 1.16 |
| IVb | LWR | −0.41 ± 0.14 | −0.67, −0.14 | 1.64 | |
| VIIa | LWR | −0.98 ± 0.30 | −1.57, −0.38 | 3.92 | |
| Winter | IVa | LWR | −0.041 ± 0.11 | −0.26, 0.18 | |
| IVb | LWR | −0.35 ± 0.36 | −1.05, 0.36 | ||
| VIIa | LWR | −1.55 ± 0.79 | −3.10, 0.01 |
The slope estimate can be interpreted as rate of change in month per °C change in temperature. Statistically significant slope estimates (i.e. 95% CIs do not contain zero) are emboldened. For statistically significant relationships the regression slope for the model has been used to estimate a rate of change in weeks per °C.
Trends in spawner length
There was a negative trend in length of female spawner over time in fishing area IVa, with average length decreasing at a rate of 20 mm per decade (F1, 2252 = 17.69, p < 0.001), and in fishing area IVb, with average length decreasing at a rate of 79 mm per decade (F1, 2472 = 175.2, p < 0.001).
Discussion
The results suggest that over three decades cod spawning has been occurring progressively earlier in two adjacent regions of the North Sea and in the Irish Sea. During approximately the same period a shift to earlier spawning was observed for sole in both the North and Irish Sea (Fincham et al., 2013). The study on sole estimated the long-term trend towards earlier spawning at 1.5 weeks per decade (Fincham et al., 2013). In this study, cod were estimated to have shifted their spawning time at a rate of roughly 0.9 weeks per decade in the Irish Sea, between 1.8 and 2.4 weeks per decade in the central North Sea, and between 0.8 and 1 week per decade in the northern North Sea. The average rate of change in spawning time of sole therefore falls within the estimated rates of change for cod. Differences between the two species likely result from the different environmental conditions they experience given that sole is a shallow-water, coastal species (Fincham et al., 2013), whereas cod is a deeper-water, demersal species. The rate of change presented for sole is furthermore an average across several spawning stocks, and therefore may include variation in rates between different local populations. Overall, the shift towards earlier spawning time observed in these two species occupying a warming regional sea can be interpreted as a common, cross-species response to climate change.
Increasing autumn SSTs were observed for all regions in this study between 1985 and 2015. This is consistent with previous research which has demonstrated rapid increases in SST in the North Sea in the last few decades (MacKenzie and Schiedek, 2007). In this study, significant decadal trends were not found in winter SST in all regions. Previous research has demonstrated that of all the seasons, summer temperatures in the North Sea have increased most rapidly and consistently (MacKenzie and Schiedek, 2007). The lack of a clear temporal trend in winter temperatures in all regions might be attributed to the increase in extreme winter weather events which have occurred due to rapid Arctic warming (Cohen et al., 2014). The occurrence of intermittent colder winters, illustrated in Figure 4, prevents a clear warming trend in winter SST from being observed over the 30-year period of this study.
The trend towards earlier spawning time in all three regions, which was negatively correlated with a warming trend in autumn SST, supports the hypothesis that warming temperatures in the North Sea and Irish Sea are associated with a shift towards earlier spawning of cod. The negative relationship between seasonal SST and PRM is consistent with previous experiments which have shown that cod held at warmer temperatures spawn earlier (Kjesbu et al., 2010). Therefore, the most parsimonious explanation of the observed correlation between spawning time and temperature is that it reflects the causal, relationship between temperature and oocyte development rate and subsequently spawning date (Kjesbu et al., 2010). Stronger relationships were found between autumn SST and PRM than between winter SST and PRM. Model predictions have suggested that changes in temperature will have the greatest effect on spawning time the sooner after vitellogenesis these changes occur. These results therefore agree with these predictions, as vitellogenesis occurs around the autumnal equinox (Kjesbu, 1994).
Alternative mechanisms which could produce the temporal trends in PRM as calculated from the roe landings dataset must be considered when interpreting the results. The opportunistic roe landings dataset was collated from records from the fishing industry, therefore variation in fishing effort could bias the data. For example, the management measures which have been implemented in the North Sea and Irish Sea throughout the time period of this study to aid the recovery of cod stocks may have affected the distribution of fishing effort both spatially and temporally. Measures such as real-time closures in Scottish waters since 2008 (ICES, 2015c), closures of spawning areas in the North Sea in 2001 (Anon, 2001), closure of spawning areas in IVb in 2010 (Kraak et al., 2013) and the continuous closure of spawning grounds in the Irish Sea since 2000 (ICES, 2005) may each have acted to direct fishing pressure away from spawning cod in the study areas, potentially affecting roe landings. The potential bias associated with temporal variation in fishing effort was addressed through the calculation of a RI. The RI represents roe landed as a proportion of total cod landed, therefore adjusting for fishing effort. Estimates of PRM and trends in PRM calculated using RI or LWR did not differ significantly, suggesting variation in fishing effort did not confound analysis. However, the estimated rates of change in PRM were lower when weighted by RI than by LWR. It is possible that this is due to an overestimation of the rate of change in PRM when weighted by LWR, which would suggest that there was some non-biological factor influencing the timing of roe landings. However, it is not obvious how the management measures mentioned previously could produce a gradual but systematic shift in fishing effort to earlier in the year throughout the northern and central North Sea.
Timing of peaks in roe landings could also be affected by market demand for roe. When estimating average spawning times from roe landings, Pederson (1984) reported that during years when there was a lack of demand for roe an earlier average spawning time was calculated, as only the highest quality roe was retained and partially spent roe later in the season was discarded. As there is a negative trend in reported roe landings from all regions in this study, it could be hypothesized that this is due to a decreasing trend in demand for roe. However, decreasing roe landings may simply be driven by reductions in catches of cod in response to decreasing quotas (ICES, 2015a), rather than by a decline in demand for roe.
The limited spatial resolution of the roe landings dataset prevented the delineation of location of spawning cod more precisely than to fishing area. Each of the fishing areas for which data was available likely encompass more than one cod spawning ground. Several spawning aggregations of cod can be found in IVa, including concentrations west of Orkney, around Shetland and within the Moray Firth (Wright et al., 2006). Anecdotal evidence from surveyed fishermen suggested that cod spawning time differs between spawning areas within the North Sea. Variation in timing of spawning within the central North Sea, related to distance from the coast, has been reported (Brander, 1994). Cod spawning in the western Irish Sea appears to occur earlier than in the eastern Irish Sea (Brander, 1994). There is a possibility that spawning cod respond to rising temperatures through spatial shifts in spawning habitats, as has been recorded for Northeast Arctic cod (Sundby and Nakken, 2008). However, previous studies using data storage tags to track the movement of spawning cod in the North Sea have demonstrated that individuals appear to show site fidelity to spawning sites (Neat et al., 2014). Without spatially resolved data it is not possible to assess whether a shift in fishing effort to areas where cod spawn earlier has occurred during this study period. However, this shift would have to have occurred in all three regions studied here to result in the consistent, earlier trends in PRM recorded.
It is well known that fishing can truncate the age and size structure of exploited fish stocks (Brander, 2007). Warming sea temperatures can also result in cod maturing at smaller sizes and younger ages (Drinkwater, 2005). Consistent with these observations, the average length of spawners in the North Sea has decreased over the study period. Given that larger cod in the North Sea begin spawning earlier than smaller individuals (Kjesbu et al., 2010) the observed shift towards earlier spawning cannot be attributed to shift in size structure of the stocks. The trend towards reduced size of spawners would be expected to result in an opposite trend in spawning phenology than was observed. This suggests that spawner length cannot be the only factor influencing the spawning phenology of cod in the North Sea. If the trends in spawning time recorded here were indeed driven by rising temperatures, the magnitude of the rate of change may actually have been inhibited by concurrent demographic changes.
Another factor which has been proposed as influencing spawning time of cod is dietary lipid composition. Specifically, a delay in spawning time of eastern Baltic cod has been linked to a reduction in the availability in cod prey of arachidonic acid (ARA), a lipid associated with ovary development. This change in dietary lipid composition was hypothesized to result from changes in the plankton composition of the Baltic Sea (Røjbek et al., 2012). The phytoplankton and zooplankton communities in the North and Irish Seas have undergone dramatic changes in phenology and biogeography in recent decades (Beaugrand et al., 2002; Edwards and Richardson, 2004; Alheit et al., 2005). In 2003 North Sea cod ovaries were found to contain high levels of ARA (Tomkiewicz et al., 2010). As cod around the United Kingdom are generalist feeders, with a diet dependent on prey abundance and which can exhibit size-dependent temporal and spatial prey switching behaviour (Floeter and Temming, 2003; Pinnegar et al., 2003; Trenkel et al., 2005), it would be a complex task to assess whether changes in the plankton communities have resulted in changes in the lipid and fatty acid composition of cod that in turn impacted spawning time.
Of the various factors discussed above which could result in the trend towards earlier PRM, rising temperature is currently the only variable for which there is consistent evidence throughout the three regions. Furthermore, the trend is consistent with the putative causal relationship between temperature, oocyte development rate and spawning time (Kjesbu et al., 2010). The finding that advancing spawning phenology appears to be correlated to warming autumn sea temperatures in three different regions has implications for understanding the impacts of current and future climate change on cod populations. The rates of change calculated here suggest that an increase in SST of 1 °C may result in a shift in cod spawning time of up to 3.92 weeks. Sea temperatures experienced by cod in the North and Irish Seas are expected to continue to increase throughout the 21st century, potentially warming by 2 °C by 2100 (Drinkwater, 2005). In addition to temperature, spawning phenology of cod is believed to be influenced by a combination of abiotic and biotic factors, namely the age and size of mature fish (Morgan et al., 2013) and photoperiodic cues (Kjesbu et al., 2010). Age and size at maturity may be influenced by rising temperatures (Drinkwater, 2005) whereas photoperiodic cues could remain fixed. The responses of spawning phenology to future climate change in this species, known to be behaviourally flexible and with a wide thermal tolerance (Righton et al., 2010) are therefore difficult to predict. However, the trend towards earlier spawning concurrent with rising temperatures suggest that consideration of reproductive phenology is essential when assessing the population-level response of cod to climate change.
Consequences of climate change-induced shifts in spawning phenology of cod include the potential for a mismatch in timing with vital prey. Cod larvae rely particularly on copepods as a major food source (Beaugrand et al., 2003), hence the timing of peaks in larval abundance must overlap with peaks in copepod abundance to ensure survival of larvae. Research on the impact of climate change on plankton phenology in the North Sea has demonstrated that the timing of spring peaks in copepod abundance has no relationship with SST (Edwards and Richardson, 2004). Between 1958 and 2002 the average timing of the spring bloom did not change, likely because diatom phenology is controlled by photoperiod rather than temperature (Edwards and Richardson, 2004). It has already been demonstrated that climate change-induced changes in the composition and succession within the plankton ecosystem in the North Sea has impacted cod recruitment (Beaugrand et al., 2003). As spawning time is necessarily linked to timing of peaks in larvae, the evidence presented in this study demonstrates how a temporal mismatch between the arrival of cod larvae in the plankton and the spring plankton blooms could occur and potentially further impact recruitment. It may be possible that such a mechanism has already been impacting cod recruitment, contributing to the slow recovery of North Sea cod (ICES, 2015c) and the lack of recovery observed for Irish Sea cod stocks (Kelly et al., 2006).
Conclusion
To date, no change in the growth rate of cod in the North Sea has occurred (Baudron et al., 2014), and the evidence suggesting that there have been shifts in spatial distribution of North Sea cod is equivocal (Righton et al., 2010). Therefore, climate change-induced shifts in the reproductive phenology of cod may currently provide the strongest evidence that rising temperatures have had a detectable impact on the life history of cod in the North Sea. Owing to the coarse temporal and spatial resolution and opportunistic nature of the datasets used in this study, the rates of change in spawning phenology calculated here should be considered as preliminary estimates of overall trends. Further work to validate these preliminary findings could involve tagging studies to quantify the temperatures experienced by cod undergoing vitellogenesis and so validate SST as a representative indicator of thermal experiences of cod. In addition, egg surveys could be used to confirm that estimates of spawning time calculated through cumulative egg curves correlate with those calculated through roe landings from the study regions (e.g. Pederson, 1984). However, the agreement in observed trends in temperature and spawning time across three different regions, which is contrary to the concurrent change in size composition towards smaller spawners, provides consistent support for attributing the change in phenology to warming temperatures. As a shift in spawning time has now been suggested for two species in the North and Irish Sea, with divergent life-strategies and exploitation histories, shifts in reproductive phenology may represent a synchronous response to climate change across marine fish species in these regions.
Supplementary data
Supplementary material is available at the ICESJMS online version of the article.
Acknowledgements
Thanks to Mike Armstrong at Cefas for providing landings data from the Irish Sea and to Gail Burns at Marine Scotland Fisheries Statistics for providing the landings data from the North Sea required for this project. Thanks also to Peter Witthames for advice and suggestions for methodology and Rodrigo Wiff for his advice on statistical analysis.
Funding
The Scottish Fisherman’s Trust provided financial support for this research.
References
