Temperature-dependent adaptation allows fish to meet their food across their species’ range

See allHide authors and affiliations

Science Advances  25 Jul 2018:
Vol. 4, no. 7, eaar4349
DOI: 10.1126/sciadv.aar4349


In seasonal environments, timing is everything: Ecosystem dynamics are controlled by how well predators can match their prey in space and time. This match of predator and prey is thought to be particularly critical for the vulnerable larval life stages of many fish, where limited parental investment means that population survival can depend on how well larvae match the timing of their food. We develop and apply novel metrics of thermal time to estimate the timing of unobserved stages of fish larvae and their prey across the north Atlantic. The result shows that previously identified life-history strategies are adaptive in that they allow parents to “predict” a beneficial environment for their offspring and meet larval fish food timing that varies by 99 days across a species’ range.


Many environments are seasonal, requiring species to match growing season conditions in time and space to ensure population persistence (1). This match is particularly important in vulnerable early life-history stages of fish, where population abundance is thought to be linked to survivorship and growth during critical larval periods where individuals are “first feeders,” switching to exogenous food sources after using up yolk reserves (2). It is this match (or mismatch) of the first-feeding stages with their prey in time and space that affects the vulnerability of first feeders to starvation and predation (3), consequently shaping future population abundance [that is, the match-mismatch hypothesis (MMH) illustrated in fig. S1].

Larval predator (first-feeder) timing is expected to be adapted to the timing of their prey, but adult predators must start the process of reproduction long before the chemical or visual cues of direct food encounter (4). Predators must respond to proximate cues to time reproduction and development to match first-feeder timing with variable prey timing (3). In the case of many fish species in seasonal environments, this process requires the timing of spawning (annual reproduction) and development from egg (where they rely on endogenous yolk resources) to first-feeding larvae (where they are dependent on exogenous food resources) to match the availability of larval zooplankton stages (for example, copepod nauplii) in the environment. Larval fish growth and survival will change relative to their ability to match with their food in time and space, and thus, life-history (spawning) timing and development rates could be adapted to coincide first-feeding predators with median prey timing conditions (3).

Estimating the degree of match-mismatch between animals and their prey is challenged by a lack of observations of the relevant stages (for example, in the case of fish, the first-feeding larvae, and their copepod naupliar food) and, until now, has been confined to studies of lower trophic levels [for example, shrimp and their phytoplankton prey (4)]. Because of sampling difficulties, observations of the timing of relevant stages are often not made on a temporal and/or spatial resolution sufficient to resolve growing season adaptation across large spatial ranges. However, recent advancements on the mechanisms controlling the development of aquatic ectotherms (for example, the roles of light and temperature) allow for the creation of novel, physiologically relevant timing metrics that one can use to estimate the timing of unobserved stages. In particular, the recent successes using integrated temperature metrics (for example, degree days) to represent a “theory of relativity” for ectotherms (time scaled by temperature) has led to our development of the “thermal development fraction” (TDF; dimensionless), a flexible metric that allows one to scale development progress within a stage with observed temperature (and beyond, for example, food). In this way, we are able to estimate the timing of unobserved stages in dynamic environments, even allowing for nonlinearities in the temperature-development relationship found at extreme temperatures [for example, living beyond tolerance limits (5)]. The timing of both sides of the match-mismatch equation (predators and prey) are shown to respond to temperature (6, 7), with larval fish potentially “using” temperature-dependent development as a predictor of median growing season timing (8).

Here, we use novel metrics of thermal time to estimate timing match-mismatch for Atlantic cod (Gadus morhua) and their prey across their species’ range (~40° to 80°N; Fig. 1). In doing so, we test a century-old hypothesis and demonstrate that temperature-independent life-history variability previously identified for this species (8) is adaptive in that it allows first-feeding offspring to “match” high concentrations of prey that will support high growth and survival. There is a 1:1, coincident relationship in expected mean annual timing between fish larvae and their prey field across the species’ range due to adaptive variations in reproductive timing across populations. We use these mechanism-driven results to discuss implications for population abundance in the face of changing climate and, in particular, changing life-history cues that may lead to decoupling of predators and their prey.

Fig. 1 Chart showing the North Atlantic Ocean with locations of predator (Atlantic cod, G. morhua, acronyms) and prey (copepods, filled circles and shaded areas, corresponding colors) observations used in this study.

Populations are Georges Bank (GEO), Gulf of Maine (GOM), spring-spawning western Scotian Shelf (WSS1), northern Gulf of St. Lawrence (NSL), south Newfoundland (SNL), Grand Banks (GB), Flemish Cap (FC), southern Labrador and eastern Newfoundland (LAB), west Greenland offshore (WGO), west Greenland Inshore (WGI), Iceland (ICE), Faroe Plateau (FP), northeast Arctic (NEA), western Baltic Sea (WBS), North Sea (NS), Irish Sea (IRS), and Celtic Sea (CEL). Fish population positions were approximated from (62). Copepod sampling locations as per (10, 21, 22).


Variability in spawning time shows an 86-day change across populations of Atlantic cod, a pattern shown previously to be independent of environmental controls on gonadal maturation (Fig. 2) (8). On the basis of thermal controls of early cod development (via the TDF method, see Materials and Methods), we show that this pattern translates to annual arrivals of vulnerable, first-feeding stages ranging 104 days across populations (Fig. 2). We also find shifts in larval copepod timing, with a 99-day change in the prey-field timing for the young cod across the north Atlantic (Fig. 2). As a result, we estimate the larval cod to match their food in time across their species’ range (~40° of latitude), with first-feeding cod larvae estimated to match in timing the occurrence of their larval copepod food (Fig. 3; SMA regression, P < 0.001; R2 = 0.68) with a slope and intercept no different than 1 (P = 0.19) and 0, respectively (P = 0.66; similar fit for individual sources given in fig. S2). These patterns hold when we test our assumption of using Calanus finmarchicus to represent the predator’s copepod diet (see Materials and Methods). Results are similar when the populations are restricted to only the northern populations that predominantly eat C. finmarchicus as their dominant prey source (fig. S3) and when information on warmer water Pseudocalanus sp. temperature-dependent development is included for the southern populations (fig. S4).

Fig. 2 Predator-prey timing across a species’ range.

Timing of spawning (red triangle), predators (fish larvae, red), and prey (larval copepods, blue) for Atlantic cod populations across the species’ range (population acronyms on dependent axis as described in Fig. 1). Predator timing and prey timing are estimated via TDF metrics (see Materials and Methods) using population-specific mean temperature phenology estimates (symbols, fish larvae, red circle; larval copepods, blue square) and estimates of mean temperature phenology ± one standard deviation (lines; fish larvae, red; larval copepods, blue).

Fig. 3 Predator-prey match across a species’ range.

Estimates of larval copepod (N3 stage) and fish larvae (first feeders) timing for populations of Atlantic cod (G. morhua) across the north Atlantic. Estimates of N3 stage copepods and first-feeding fish larvae timing are made via the TDF method (see Materials and Methods) using observations of copepodite timing (table S1) (10, 21, 22) and spawning time (42), respectively. We made estimates using SDs based on population-specific mean daily temperature estimates (data points with error bars ± one standard deviation from mean temperature phenology; data labels refer to populations in Fig. 1). Variability in larval fish timing due to uncertainty around the relationship between temperature and time to yolk absorption (fig. S8) is shown by open symbols. Also given is the standard major axis (SMA) line fits (solid line, with 95% confidence intervals around the slope as dashed lines; P < 0.001; R2 = 0.68; slope not different from 1, P = 0.19; intercept not different from 0, P = 0.66). The 1:1 line (that is, slope of 1, intercept of 0) is given in the solid gray line.

The importance of estimating unobserved stages via the TDF metric that allows for temperature-dependent, time-varying stage durations (SDs) is demonstrated by comparing the results in Fig. 3 to those found when we either (i) compare observed stages or (ii) assume constant SDs. First, observed stages for predator and prey include eggs (via spawning time estimates) and juvenile copepods (copepodite observations). While these stages are observed directly, they are not the stages directly involved in the predator-prey interaction and, thus, not the players involved in possible food-driven adaptive evolution of life-history timing. The limitations of relying on indirect evidence from observed developmental stages can be found when we compare the timing of these observed stages (fig. S5), where predator and prey proxies are correlated in time (SMA regression, P < 0.001, R2 = 0.70; slope not different from 1, P = 0.21) but (unsurprisingly) not coincident (that is, elevation differs significantly from 0; P < 0.001) and estimates of predator-prey temporal match are not possible. Second, using constant SDs based on population-specific average temperature and ignoring temperature variability during the development of each stage (versus including it via the TDF metric) result in estimates of predator (larval fish) timing that are again correlated with those of their prey (larval nauplii; SMA regression, P < 0.001, R2 = 0.72) but not coincident (slope different than 1, P = 0.022), with first-feeding larvae predicted to arrive up to ~41 days before prey appears in some populations (fig. S6). It is only when we estimate the stages directly involved in the predator-prey interaction based on SDs that are environmentally dependent and allowed to vary over the development period via the TDF metric that accurate estimates of stage timing are obtained and the strong (1:1 and coincident) agreement of predator and prey timing are found (Fig. 3).


Here, we show that timing of first-feeding cod larvae and that of their prey are not only correlated but coincident (1:1) for populations spanning the north Atlantic. It is this coincidence (slope of 1 and intercept of 0) that describes predators and prey that are overlapping in time. This finding provides an adaptive explanation to previous observations of varying temperature-independent life-history constants (reproductive timing) across the species’ range (8). In our previous work, we found that cod populations in cold regions of the species’ range had much lower thermal requirements for gonadal development than in warmer regions and therefore demonstrated counter-gradient adaptation of a key life-history trait with decreasing thermal needs to develop gonads for spawning with increasing latitude (8). This adaptation is reflected in parallel latitudinal clines in thermal time to spawning found on both sides of the Atlantic and mirrors other evidence of adaptation including the evolutionary history of the species (8). Here, we show that the patterns in thermal requirements of cod reproduction allow for temporal overlap of first-feeding larvae with the local prey field, thereby likely increasing growth and survival potential of their vulnerable first-feeding larvae. In the absence of this local adaptation, many cod populations would produce larvae at times when the concentrations of their preferred prey would likely affect survival either directly or indirectly through growth.

These results confirm predictions of modeled larval fish production indicating that first-feeder timing can be selected for by the timing of prey abundance (9) and for Atlantic cod in particular where spawning time variations were previously observed (1012) and unexplained but hypothesized to allow for a match with prey timing (1113). Here, we provide findings consistent with this hypothesis across the species’ range by presenting and using a method to estimate the timing of the relevant predator-prey players from observed timing of other stages. The use of the TDF metric allows us to scale time with temperature in a manner that includes nonlinearity of the temperature–development rate relationship at temperature extremes, similar to the motivation behind metrics being used to estimate and explain the timing of, for example, flowering in terrestrial systems (14), as well as fish growth and development (15, 16). Our timing estimates agree with the limited field observations available: Our prey timing estimates are consistent with observations of C. finmarchicus nauplii (various stages) timing for Faroe Shelf waters [near FP, April to June 2004 (17)], and our predator timing estimates overlap with available observations of larval cod timing in the field for the Newfoundland [SNL, June to September 1998 (18)] and WGI [May to July 2010 (19)] populations. The power of being able to directly consider the relevant stages involved in the predator-prey interaction through the TDF metric (versus being limited to observed players) can be found by comparing results to indirect measures of predator and prey timing. In the latter, we observe departures from coincidence (this study, see Results and fig. S5) or lack of correlations at all (12).

While prey-timing adaptation is hypothesized as the driver for life-history timing in many species [for example, tuna (20)], cod represents a prime candidate for tests of the MMH (spring-spawning populations in environments with distinct growing seasons). Still, the relationship observed in Fig. 3 includes some interesting among-population variability, particularly concerning those populations furthest from the coincident line (for example, WBS, NS, and NEA). Where MMH is at play, large residuals from the coincident line may indicate data limitations of the current study, either through limited observations (for example, copepod timing estimates only available at the edge of the WBS area; Fig. 1) or departures from the species-wide, temperature-dependent development rates (fig. S8). The latter highlights the need for population-specific estimates of temperature effects on development (also discussed in “Challenges when estimating timing” and “Life in a warmer ocean” sections). More broadly, quantifying match-mismatch dynamics depends on identifying the resource measure that will limit survival either directly or indirectly through growth. Here, we use estimates of both onset (10) and maximum (21, 22) copepod timing, while the link between food and survival may be different in other populations [for example, autumn spawners shown in previous work to be outliers to thermal life-history constant patterns (8)] or species [for example, herring (23)]. For example, the degree of match with food in time may affect survival indirectly, through growth, with a need to maximize size before overwintering [for example, herring (23)]. Beyond MMH, factors could include abiotic conditions that optimize survival of particular stages [for example, egg hatch success and temperature (24)] or the avoidance of high abundances of predators at certain times of year. In all cases, our method could be used to help disentangle the various candidate hypotheses and determine what, if any, role the MMH may play in determining population abundance for a given population or species.

Challenges when estimating timing

Because of limitations in life-history characterizations for fish and zooplankton populations across the species’ range, a number of assumptions were necessary to estimate predator timing and prey timing. First, we assume that the relationships between temperature and larval development for cod are similar across populations. Temperature is the dominant factor controlling egg development time (24), with negligible [for example, anchoveta, Engraulis ringens (25)] to measurable [for example, pilchard, Sardinops sagax (26)] variation in temperature-dependent development among populations in a number of fish species. However, observed temperature-driven variability in development times may as easily be explained by differences in experimental protocol, with as much variation within and between studies on the same population as across populations (27). In addition, development time will also depend on maternal effects such as egg quality (27) and size (6) but will be dominated by temperature effects, particularly when comparisons are made over large latitudinal (and temperature) gradients (8, 15).

For copepods, we use a temperature-dependent development rate for C. finmarchicus sampled from waters between Georges Bank and Cape Cod in the northwest Atlantic (7). This parameterization includes three assumptions: (i) Temperature affects development in similar ways across the species’ range, (ii) temperature effects overwhelm possible food effects on development, and (iii) C. finmarchicus is an appropriate representative of the larval cod prey field across the species’ range. First, while spatial variability in temperature-dependent copepod development may lead to errors in the estimates of larval timing, there is evidence for similar temperature-dependent development rates across space [for example, Pseudocalanus spp. (28)]. Second, our assumption that food is not limiting copepod development is arguable here where we consider long-term mean development timing across large spatial areas where variability in long-term mean temperature among locations will be much higher than that of food. Effects of food availability on naupliar timing could prove important for future studies needing to estimate naupliar timing variation within a population, among years for areas where food availability is likely to be below satiating levels (7). In these cases, we can easily extend the TDF metric to include food effects via the dynamic SD estimates in Eq. 1 (see Materials and Methods). Finally, our assumption that C. finmarchicus represents the prey field of the first-feeding cod may not hold for all populations, particularly in warmer waters (for example, WSS, GOM, GEO, NS, IRS, and WBS populations in Fig. 1) (10). We show here that our results are robust to these assumptions with similar predictions of predator-prey match across space when (i) only northern C. finmarchicus feeding populations are included (see Results and fig. S3) or (ii) developmental rates for Pseudocalanus spp. are included for southern populations (see Results and fig. S4). Future studies that characterize spatial variability in temperature-dependent development rates for both fish and their prey will further refine these estimates of predator-prey overlap, allowing for intrapopulation estimates of match among years and more robust predictions of the spatially varying response of individual populations under changing climate conditions (for example, helping identify populations that are most vulnerable). Alternatively, direct observation (monitoring) programs of the relevant stages could be implemented.

Implications for abundance

Here, we provide a meta-analysis of the adaptive significance of spawning time variations for cod populations. Investigations of the relative match-mismatch between cod and their prey may also be examined within a population, among years, to get an estimate of the implications of match-mismatch timing for population abundance using similar approaches as here (where data are available). Such an endeavor would quantify the fish’s ability to “predict” changes in prey timing that occurs among years (29). The new approaches to investigate life-history processes and links to environmental conditions based on thermal time could provide changes in the relative timing of critical stages and their food for each year-class to begin to quantify variability in interannual match-mismatch and relate this mismatch to year-class abundance at older ages for applications of population and ecosystem dynamics, as well as fisheries management. This method will have particular promise if we are able to further characterize the abundance and magnitude of the prey field to make quantitative estimates of the implications of match-mismatch over variable conditions where processes may become nonlinear (1). Such an analysis would be a direct test of a long-held hypothesis concerning fish abundance variability—MMH stating that survival of juvenile fish is correlated with the degree of match between vulnerable first-feeding stages and their prey (1, 2). This analysis can also help (i) explore variability in timing (for example, spawning time and duration) within a population along with population structure (age and size) and (ii) disentangle potential abundance drivers, quantifying the influence of predator-prey match versus other factors affecting survival [for example, growth, egg or overwintering survival, etc. (23, 24, 30)] and how these vary among populations and species (23).

Life in a warmer ocean

Advances in life-history timing are predicted for ectotherms in warming waters and are already found for, for example, fish spawning phenology (8, 31) and their prey (32). The survival effects of these life-history timing shifts will be minimized if prey timing and predator timing respond to environmental change in similar directions and magnitudes. In these cases, shifts in predator and prey life-history timing may act to maintain the predictive strength of the adaptive cues currently used by fish to match their prey, resulting in little change in the match between fish and their food and small consequences to abundance. In contrast, the implications for population abundance will be highest when temperature changes and/or responses differ across trophic levels and trophic levels risk becoming decoupled (32, 33). The potentials for predator-prey mismatch may occur because of temperature-dependent development timing if/when (i) temperature dependencies of development rates vary in shape across trophic levels, (ii) temperature changes vary between habitats used by the different trophic levels, and (iii) dominant forcing factors (for example, light versus temperature) vary across trophic levels. First, how animals at different trophic levels are predicted to respond to temperature changes will depend on the shape of their respective development-temperature relationships. Predator and prey SDs explored here both show declines in SD with increasing temperature. When predator and prey responses are compared, relative decreases in SD with increased temperature are similar at higher temperatures (>4°C) but higher for prey (copepods) versus predators (fish) at lower temperatures (<2°C). This may mean that the potential for decoupling of predator and prey timing is highest for the cold water populations. Characterization of the species- and spatial-specific variability in temperature-dependent rates (including tolerance limits) will be necessary to predict whether predators and prey timing will respond to temperature changes in similar ways or at all. Second, temperature changes may vary in direction and magnitude across the habitats used by predators and prey with the potential for future predator-prey match depending on, for example, the relative changes in sea surface versus bottom temperature (4). Third, the risk of decoupling may be particularly important for trophic levels lowest on the food chain where one major controlling factor (temperature) is shifting, while the other (light) is held relatively constant. Beyond temperature and light, changes in circulation, stratification, and mixing may also influence predator-prey timing, match, and resulting larval survival.

Changes in climate may mean that the signals predators “use” to foreshadow prey timing lose predictive power, resulting in shifts in predators relative to their prey that are no longer adaptive and an increased mismatch between predators and their prey under future conditions (32). These effects may be mitigated in a number of ways. First, generalist feeders may switch their prey preferences to available food (34). While major forcing factors are expected to act in similar directions across species of the same trophic level, variation in prey characteristics (for example, tolerance limits and life-history strategies) may encourage one prey type over another (32), with a predator’s success depending on its ability to switch along with dominant prey types. Second, fish may adapt their temporal distribution (via spawning times or development rates) to match temporal shifts in their prey. Potential for adaptive variability in reproductive timing is seen in space [(8) and this study], but how this relates to the potential for future adaptation over time [that is, through space-for-time substitutions (35)] will depend on the relative paces of adaptation and climate change. Third, fish may change their spatial distribution to maintain temporal match with their larval-stage prey by, for example, moving poleward and/or deeper when waters warm. Migration between cooler versus warmer waters may delay/speed up development to coincide with prey timing. Spatial shifts in community players have been observed and correlated with warming waters (29), with examples of spatial and (in some cases) temporal shifts in larval fish prey correlated with temperature in the northeast Pacific (32).

In conclusion, we have implemented new approaches to address a century-old hypothesis in marine science: Are patterns in reproductive timing adapted to the prey-field timing of vulnerable larval stages? By developing metrics that accommodate thermally dependent physiology to estimate timing of relevant players (predator and prey), we find a 1:1 overlap between widely distributed fish populations and their prey. Coupled with our previous finding that cod reproduction throughout its north Atlantic range appears locally adapted to thermal regimes (8), this finding not only resolves a long-studied issue in ecology but also demonstrates how a species can evolve its life history, so reproduction can track the environment to time larvae with the conditions (for example, the production of suitably sized zooplankton) that will increase offspring survival and ultimately overall fitness. This work demonstrates how both environmental (for example, temperature) and genetic cues (34) shape predator-prey phenology across space and time and provides tools to explain and predict abundance and distribution timing for wide-ranging species now and in the future.


Assessing predator-prey match requires timing estimates of both first-feeding fish larvae and larval (naupliar) zooplankton. These stages are often not observed directly (previously limiting mechanistic tests of the MMH), necessitating backward/forward calculation to estimate timing of the relevant stages (fig. S7). Here, we estimate the timing of unobserved stages using species-specific TDF metrics.

Thermal development fraction

For ectotherms, environmental temperature acts to scale time [cf. (36)] and thus observed growth and development rates [including those of zooplankton (21) and fish (6)]. This has traditionally complicated the estimate of stage timing in dynamic field environments. Here, we present the TDF metric, which allows us to scale time with temperature (15) to back-calculate (or forward-calculate) the timing of unobserved stages in dynamic environments when relevant temperature measures are available (box S1).

With this method, an animal’s progress through a stage is tracked as TDF (dimensionless), a relative measure, between 0 (new to stage) and 1 (achieved SD and predicted to develop to the next stage)TDF(t+Δt)= TDF(t)+ΔtSD(t)(1)where Δt is the time step (for example, day), SD(t) is the SD based on current environmental (for example, temperature) conditions [for example, SDyolk(t) or SDstage(t) in Eqs. 2 and 3; days], and TDF is the TDF at the beginning (t) and end (t + Δt) of the time step (dimensionless). The TDF metric can be used to estimate the timing of unobserved stages by following time forward (or backward) from the observed stage and incrementing development rate estimates until TDF = 1 and the animal is predicted to occupy the next (or previous) stage. For example, at each time step (Δt, for example, 1 day), an increment of development (ΔtSD(t)) is added to the TDF sum based on the time-step duration (for example, 1 day) and the temperature-specific SD estimated for that day [SD(t)]. The day on which TDF ≥ 1 is the estimated date of the unobserved next stage (easily estimated for multiple stages, for example, with a “while” loop when programming). The TDF method is illustrated in box S2, where it is used to estimate timing of larval fish for three populations of cod.

TDF models were used to estimate the occurrence of unobserved predator and prey stages (first-feeding cod larvae and larval copepods) from measures of observed stages (cod spawning and juvenile copepods; details below). The TDF model extends previous work, exploring thermal constants of development (8, 37) to include information on the nonlinear responses of ectotherm development rates to temperature [for example, beyond tolerance limits (5)], building on the molt cycle fraction (38) and incremental growth-at-temperature models (39), which can accurately represent development rates in fluctuating temperatures (39). The TDF metric presented here could easily be extended to allow for the effect of other factors (for example, food) on development rate through the dynamic SD estimates.

Estimating larval cod timing

The average timing of larval (first-feeding) Atlantic cod (G. morhua) was estimated for 17 populations across the species’ range (Fig. 1). The first-feeder stage was used as the relevant predator stage rather than hatching, as fish may not hatch at the same development stage (40), and it is the timing of yolk absorption rather than hatching that indicates the need to switch to exogenous food sources (41). Because of the limited availability of in situ estimates of first-feeding cod, we forward-calculate first-feeder timing from population-specific spawning times (8, 42).

Using the TDF metric requires estimates of SDs over time. As above, SD from egg to yolk absorption is temperature-dependent across ectothermic fish species (fig. S8) (43), with temperature as the only factor to consistently and significantly affect the development rates of Atlantic cod eggs [for example, versus egg size (6) and versus salinity (44); note also that eggs are a nonfeeding stage]. Thus, we used the TDF metric (Eq. 1) to forward-calculate the timing of first-feeding cod from information on temperature-dependent SDs and population-specific spawning time and temperature estimates. Cod egg development occurs in temperatures from <0°C to 12°C, with most populations developing in waters between 2° and 7°C (13, 27, 42). We consolidated published estimates of development time from egg to first feeder (yolk absorption) to fit a predictive model of temperature-dependent yolk absorption timing for Atlantic cod (fig. S8) asSDyolk(t)=eayolkT(t)+byolk(2)where SDyolk(t) is the SD (days) from spawning to yolk absorption estimated for the current time step t, that is, at a temperature T(t) (degree Celsius), and ayolk = −0.125°C−1 (standard error = 0.0086°C−1) and byolk = 3.72 (standard error = 0.055). The model was fit using linear regression with logn-transformed SDyolk (P < 0.001; R2 = 0.87). Statistically similar coefficient estimates were obtained using nonlinear least-squares model fitting. The data underlying the fit in Eq. 2 include development rates from cod in the Gulf of Maine (45), Newfoundland (6, 46), northeast Arctic (47), and western Scotian Shelf (48) populations—populations spanning 30° (43°N to 73°N) of the cod’s 40° latitudinal range. Because of a lack of population-specific data on yolk absorption at variable temperatures, we are limited to the use of this species-wide relationship. This may not capture the full picture if cod egg development demonstrates thermal adaptation across populations (see Discussion). However, no difference in temperature-dependent development of early cod stages was previously found for east versus west or north versus south populations, at least across temperatures of 0° and 12°C [for example, hatch timing (27)]. Moreover, temperature was found to explain between 90 and 98% of variation in egg development time within an individual fish species (24). The potential for limited variability in thermal-dependent development time across populations may also be reflected in the observation of temperature-dependent growth rates shown to be relatively fixed across the species’ range (15).

We combined information on population-specific spawning time (8, 42) and temperatures to estimate timing of first feeders via Eqs. 1 and 2. Temperatures were extracted for the period 1913–2011 from the International Council for the Exploration of the Sea and Fisheries and Oceans Canada databases for population-specific areas given in table S1. In addition, inshore Greenland temperature estimates were obtained for Nuuk fjord (8). Temperatures were restricted to those corresponding to the locations of the cod populations (table S1) and depths between 5 and 100m to correspond to waters in which the cod eggs, larvae, and copepods are developing (4952). Sample years for temperature data spanned 1913 through 2011, with most of sampling occurring in the last ~30 years of this range. All sample years available were used in the analysis, as we are concerned here with long-term average predator timing and prey timing. In all cases, results were near identical when temperatures used were restricted to only the last 30 years available (1982–2011). A population-specific temperature time series (average monthly temperature ± one standard deviation) was estimated for each population, and monthly temperature estimates were linearly interpolated to approximate daily temperature estimates (fig. S9). These population-specific temperature time series were used to give the daily estimates of SD needed in Eq. 1 to forward-calculate first-feeder timing from spawning time for each population in a dynamic environment (see box S2 for an example of this method). First-feeder timing was estimated using population-specific mean temperature estimates and estimates of mean temperature phenology ± one standard deviation to capture variability in temperature observations across years.

Estimating larval zooplankton timing

Cod feed primarily on calanoid copepods with C. finmarchicus dominating cod diets for populations in the north and Paracalanus and Pseudocalanus spp. dominating diets for populations in the south (10, 5355). Calanoid copepods develop through six larval or naupliar stages (N1 to N6) before a further six juvenile or copepodite stages (C1 to C6). First-feeding cod larvae begin feeding on naupliar stages of the copepods (beginning at stage N3) and include later stages as they develop (29, 54, 56). Hence, calanoid copepod larvae (primarily Calanus and Pseudocalanus spp.) at stage N3 were chosen to represent the prey field of first-feeding cod.

Estimates of prey-field timing across the species’ range were limited to observations of older copepod stages with data on population-specific observed timing of copepodite stages (C1 to C3) obtained from Anderson [timing of maximum C. finmarchicus stage C2 percent composition (21)], Melle et al. [timing of maximum C. finmarchicus stages C1 to C3 (22)], and Heath and Lough [first appearance of copepod biomass above the annual average observed postspawning including Calanus spp. C1 to C4, C. finmarchicus C5 to C6, and Paracalanus/Pseudocalanus spp. stages (10)]. Heath and Lough (10) reported copepod timing for both an early (1961–1978) and late (1991–2002) period. We used estimates for the early period (1961–1978) as these patterns were consistent with those reported by Myers et al. (12) for the same data but discuss temporal shifts in prey timing in Discussion.

Development rates of zooplankton were shown in the laboratory to be highly sensitive to temperature, with longer SDs estimated at lower temperatures (7) and zooplankton timing related to environmental temperature in the field (32). Similar to above, we used a TDF metric to estimate the unobserved larval naupliar stages from observed timing of the C1 to C3 stages via environmental temperatures. Daily temperature estimates were used to estimate SDs for C. finmarchicus, a well-studied prey species for cod from Campbell et al. (7) asSDstage(t)=astage(T(t)+9.11)2.05(3)where SDstage(t) is the SD (days) from the previous stage estimated for the current time step t, that is at a temperature T(t) (degrees Celsius), and astage is a stage-specific parameter estimate (7). We assume that this response of development to temperature is similar throughout the C. finmarchicus species’ range, that (as above) C. finmarchicus represents a main prey species for the cod larvae and/or temperature-dependent development of other copepod species respond similarly with temperature. We test the effects of these assumptions by comparing results to those obtained when we (i) restrict our analysis to cod populations known to have C. finmarchicus–dominant diets (fig. S3) (10, 57) and (ii) include temperature-dependent development estimates for Pseudocalanus spp. based on temperature-dependent SD estimates in Lee et al. (58) as well (fig. S4). Here, we also assume that food dependencies on SDs are negligible. Food is known to influence copepod SDs for feeding stages (N3 and older), but population specific estimates of copepod food timing are not available over this temporal and spatial resolution. Given that we are interested here in macroecological patterns across a species’ range, we argue that the large difference in temperature among stocks (~10°C range) will dominate developmental variation (also shaping the copepod food timing). However, the methods we present here can easily be extended to include information on temperature-independent (for example, food) effects on SDs in the future (via the TDF metric) should estimates of copepod prey fields over space and time become available.

We combined Eq. 3 with population-specific daily temperature estimates (as above) to back-calculate timing from stages C2 to N3 using iterative application of the TDF method. The TDF metric is set at 1 for the initial stage (for example, C2), and the zooplankton are virtually moved backwards through each day by subtracting a thermal increment equal to the ratio between the time step (Δt, 1 day) and the SD estimated for the temperature experienced during that time step [SD(t), as per Eq. 1]. This was continued through time while TDF > 0, and the timing of the previous stage (for example, C1) was the day when TDF = 0. The procedure was then repeated for other stages (for example, N6 through N4) to estimate larval zooplankton timing. As above, larval copepod timing was estimated using population-specific mean temperature estimates and estimates of mean temperature phenology ± one standard deviation to capture variability in temperature observations across years.

Comparing timing of predator and prey

The relationship between larval zooplankton and first-feeder fish timing was made using SMA line fitting (59). SMA line fitting is preferred to standard linear regression for estimating correlations and slope parameter values when both variables may contain measurement error (59). To remove the effect of correlation between the estimated slope and intercept parameters and allow appropriate inference, timing data were centered on the overall mean (that is, a single value representing the mean of range in both predator timing and prey timing was subtracted) before subsequent analysis. A correlation between zooplankton and first-feeder timing was established using a Pearson correlation coefficient before fitting the line, as per Warton et al. (59). Where significant, the resulting slope parameter was compared to a value of 1 by comparing residual and fitted axis scores using 1 as the slope (59). The fitted slope value is not significantly different from 1 when residual and fitted axis scores are uncorrelated (59). In addition, the intercept parameter was compared to 0 via a t statistic comparing the difference between the estimated intercept and 0 to the standard error of intercept (59). All analyses were performed in R with base and smatr packages (5961).


Supplementary material for this article is available at

Table S1. Observed and estimated predator timing and prey timing used in this study.

Box S1. A description of the TDF model.

Box S2. Examples of using the TDF model to estimate predator timing.

Fig. S1. Illustrating MMH.

Fig. S2. Predator-prey match across a species’ range with source-specific estimates.

Fig. S3. Predator-prey match across a species’ range for northern populations.

Fig. S4. Predator-prey match across a species’ range using alternate prey species.

Fig. S5. Predator-prey match across a species’ range using observed predator-prey stages.

Fig. S6. Predator-prey match across a species’ range using timing estimates based on constant SD predictions.

Fig. S7. Illustration of method to estimate timing of unobserved stages via the TDF metric.

Fig. S8. Temperature-dependent SDs for Atlantic cod.

Fig. S9. Intra-annual variation in population-specific temperatures estimated from temperature observations between 5- and 100-m depth.

References (63, 64)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank S. Tsoukali for the discussions regarding egg development rates for marine fish. We thank R. Hedeholm for providing the Greenland temperature data. The study has been supported by the European Union 7th Framework Programme (FP7 2007-2013) under grant agreement number 308299 (NACLIM) and the Horizon 2020 research and innovation programme under grant agreement number 727852 (Blue-Action). This work was completed during an AIAS-COFUND Fellowship to A.B.N. at the Aarhus Institute of Advanced Studies, which receives funding from the Aarhus University Research Foundation (Aarhus Universitets Forskningsfond) and the European Union’s Seventh Framework Programme, Marie Curie Actions (grant agreement 609033). We thank the Center for Macroecology, Evolution and Climate, where some initial concepts and planning of this work was conducted by A.B.N. and B.R.M., and which has received funding from the Danmarks Grundforskningsfond. Author contributions: A.B.N. conceived the study with B.R.M and M.R.P. A.B.N. developed and implemented the analysis. A.B.N. wrote the manuscript with contributions from B.R.M. and M.R.P. Sequence of authors reflects declining contribution. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article