The model introduced here is “Coltrane” (Copepod Life-history Traits and Adaptation to New Environments) version 1.0. Matlab source code is available at https://github.com/neilbanas/coltrane. An overview of the model structure is shown in Figure 2.

Like many individual-based models, Coltrane represents the time-evolution of one cohort of a clonal population, all bearing the same traits and spawned on the same date t₀, with a set of ODEs. The state variables describing a cohort are relative developmental stage D, where D = 0 represents a newly spawned egg and D = 1 an adult; survivorship N, the fraction of initially spawned individuals that remain after some amount of cumulative predation mortality; structural biomass per individual S, and “potential” or “free scope” φ, which represents all net energy gain not committed to structure, or equivalently, the combination of internal energy reserves and eggs already produced. Combining reserves and eggs into one pool in this way lets us cleanly separate results that depend only on the fundamental energy budget (gain from ingestion, loss to metabolism, and energy required to build somatic structure) from results that depend on particular assumptions about egg production (costs, cues, and strategies). An alternate form of the model explicitly divides φ into internal reserves R and egg production rate E: the simpler model without this distinction will be called the “potential” or φ model and the fuller version the “egg/reserve” or ER model.

The φ and ER models take different approaches to generating population-level results from this cohort model, as explained in detail below (Section 2.4). In both cases, the logic changes from the simple forward time-integration at the cohort level: one runs the cohort model for all possible spawning dates t₀, retroactively determines which spawning dates would prove optimal or sustainable, and considers the cohort time series from those t₀ values, appropriately weighted, to constitute the model solution (Section 2.4). The biological logic here is similar to the backwards-in-time dynamical optimization method frequently used in studies of optimal annual routines (Houston et al., 1993; Varpe et al., 2007), although our solving method is quite different and less exact. This is a compromise with the eventual goal of coupling Coltrane to oceanographic models as a spatially explicit IBM.

Communities are generated in Coltrane 1.0 simply by running families of cases of the population-level model that vary one or more traits. Treating coexisting populations as uncoupled vastly simplifies the interpretation of the landscape of viable strategies in a given environment, or the fundamental niche of a particular trait combination, our primary modes of analysis. At the same time, it tightly restricts our choices regarding the formulation of predation mortality. In reality, coupling through shared predators can rival bottom-up effects as a determinant of community structure (Holt et al., 1994; Chesson, 2000; Record et al., 2013, 2014), and we expect that many potential applications of this model would require that this be better represented. In the present study, we have taken the minimalist, incrementalist approach of imposing the simplest possible form of predation mortality—a linear function, with scalings that closely mirror the growth and development functions (Section 2.2)—and restricting the terms of analysis. In particular, we will describe model output in terms of trait correlations, optimality, and viability, but not in terms of absolute copepod biomass or abundance. Likewise, while some plankton models resolve the process of adaptation explicitly (Clark et al., 2013), we address it only in the indirect sense of mapping the viable and optimal regions of the strategy landscape. This approach is less mechanistic but also helpfully agnostic about whether adaptation in the copepods arises through individual plasticity, species composition shifts, or natural selection per se.

An environment in Coltrane 1.0 is defined by annual cycles of three variables, total concentration of phytoplankton/microzooplankton prey P, surface temperature T₀, and deep temperature T_d. At present, these annual cycles are assumed to be perfectly repeatable, so that a “viable” strategy can be defined as a set of traits that lead to annual egg production above the replacement rate, given P, T₀, and T_d as functions of yearday t. The level of predation mortality (Section 2.2.4) might also be viewed as an environmental characteristic.

The evolution equations above Equations (1), (6), (7), (16) specify the development of one cohort in the φ model. If this model is elaborated with an explicit scheme for calculating total egg production over time E(t), then it is possible to define R(t), individual storage/reserve biomass, and interpret R as a state variable and φ as a derived quantity. The relationship between the two is

Thus, φ tracks the reserves that an animal would have remaining if it had not previously started egg production. This is a useful metric for optimizing reproductive timing, as we will show (Section 2.4).

Any explicit expression for E(t) allows Equation (25) to replace Equation (7). In one model experiment below (Section 3.4), we use the following scheme: E(t) is the sum of income egg production E_inc and capital egg production E_cap, which are 0 until maturity is reached (D = 1) and an additional timing threshhold has been passed (t > t_egg). Past those threshholds, they are calculated as

where E_max is a maximum egg production rate which we assume to be equal to food-saturated assimilation:

Thus, the trait t_egg determines whether egg production begins immediately upon maturation (if t_egg is prior to the date on which D reaches 1) or after some additional delay. Instead of t_egg, expressed in terms of calendar day, one could introduce the same timing freedom through a trait linked to light, an ontogenetic clock that continues past D = 1, or a more subtle physiological scheme. However, since we run a complete spectrum of trait values in each environmental case, it is not important to the results how the delay is formulated, provided we only compare model output, rather than actual trait values, across cases.

A population-level simulation (Figure 2) consists of integrating either the φ model (Equations 1, 6, 7, 16) or ER model (Equations 1, 6, 16, 25) for a full annual cycle of spawning dates t₀, and then identifying optimal and viable values of t₀ in terms of the egg fitness F, future egg production per egg (Varpe et al., 2007). Calculating a time series of F in the φ model requires an estimate of individual egg biomass W_e in order to convert φ(t) from carbon units into a number of eggs, and a similar issue arises in the ER population model. Thus, a digression on the determination of W_e is required.

Community-level predictions in Coltrane take the form of bounds on combinations of traits that lead to viable populations in a given environment (Figure 2). There are many copepod traits represented in the model that one might consider to be axes of diversity or degrees of freedom in life strategy: u₀, I₀, θ, D_s, K_s, W_e/W_a, and even m₀ to the extent that predation pressure is a function of behavior (Visser et al., 2008). Record et al. (2013) allowed five traits to vary among competitors in their copepod community model. We have taken a minimalist approach, where in the φ model we allow only one degree of freedom: variation in u₀ from 0.005 to 0.01 d⁻¹. Banas and Campbell (2016) showed from a review of lab studies that u₀ variations appear to be the primary mode of variation in adult size among large calanoids (W_a > 50 μgC) including Calanus and Neocalanus spp., with slower development leading to larger adult sizes. That study also suggests that variation in I₀ is responsible for copepod size diversity on a broader size or taxonomic scale (e.g., between Calanus and small cyclopoids like Oithona). However, variation in I₀ (energy gain from foraging) probably only makes sense as part of a tradeoff with predation risk or egg survivorship (Kiørboe and Sabatini, 1995) and we have left the formulation of that tradeoff for future work. We therefore expect Coltrane 1.0 to generate analogs for large and small Calanus spp. (~100–1000 μgC adult size) but not analogs for Oithona spp. or even small calanoids like Pseudocalanus or Acartia.

Choices regarding reproductive strategy require another degree of freedom. In the φ model, this does not require additional parameters, because the difference between, e.g., capital spawning in winter and income spawning in spring is simply a matter of the time t at which F is evaluated in postprocessing: each model run effectively includes all timing possibilities (Equation 33). In the ER model, however, diversity in reproductive timing must be made explicit. Under the simple scheme for egg production specified above (Section 2.3), this takes the form of running a family of cases varying t_egg for each t₀ and u₀.

This study comprises three complementary experiments (Table 2). The first of these is an idealized global testbed which addresses broad biogeographic patterns. The second is a testbed representing the Eastern Bering Sea shelf, which addresses time-variability in one population in one environment. The last is a testbed representing Disko Bay, West Greenland, which addresses trait relationships along the size spectrum in detail. The first two are evaluated entirely in terms of the φ model, while in the Disko Bay case we use the ER model to allow more specific comparisons with observations.

The global testbed consists of a family of idealized environments in which surface temperature T₀ is held constant, and prey availability is a Gaussian window of width δt′ centered on yearday 365/2:

We compare environmental cases in terms of T₀ and an effective season length

which rescales the δt′ cases in terms of the equivalent number of days of saturating prey per year. We assume that deep, overwintering temperature T_d = 0.4 T₀. The ratio 0.4 matches results of a regression between mean temperature at 0 and 1000 m in the Atlantic between 20 and 90°N, or 0 and 500 m in the Northeast Pacific over the same latitudes (World Ocean Atlas 2013: http://www.nodc.noaa.gov/OC5/woa13/). Over the same data compilation, the mean seasonal range in temperature is approximately 5°C at the surface (and approximately zero at 500–1000 m); an alternate formulation of the testbed that models T₀ as an annual sinusoid with this range gives results that are somewhat noisier but heuristically very similar to those shown in Section 3.2 below.

The Bering Sea testbed considers interannual variation in temperature, ice cover, and the effect of ice cover on in-ice and pelagic phytoplankton production (Stabeno et al., 2012b; Sigler et al., 2014; Banas et al., 2016). Variation between warm, low-ice years and cold, high-ice years has previously been linked to the relative abundance of large zooplankton including C. glacialis/marshallae (Eisner et al., 2014), and we test Coltrane predictions against 8 years of C. glacialis/marshallae observations from the BASIS program. Seasonal cycles of T₀, T_d, and P are parameterized using empirical relationships between ice and phytoplankton from Sigler et al. (2014) and a 42-year physical hindcast using BESTMAS (Bering Ecosystem Study Ice-ocean Modeling and Assimilation System: Zhang et al., 2010; Banas et al., 2016). Details are given in the Appendix in Supplementary Material.

The Disko Bay testbed represents one seasonal cycle of temperature and phytoplankton and microzooplankton prey, based on the 1996–1997 time series described by Madsen et al. (2001). We use this particular dataset not primarily as a guide to the current or future state of Disko Bay but rather as a specific circumstance in which the life-history patterns of three coexisting Calanus spp. (C. finmarchicus, C. glacialis, C. hyperboreus) were documented (Madsen et al., 2001). Details are given in Section 3.4 and the Appendix in Supplementary Material.

One case from the global experiment with u₀ = 0.007 d⁻¹, T₀ = 1°C, and δt = 135 is shown in detail in Figure 3 to illustrate the analysis method described in Section 2.4.2. In this case, out of cohorts spawned over the full first year, only those spawned in spring reached adulthood without starving (Figure 3B, blue–green lines; non-viable cohorts not shown). The fitness function F (Equation 42) declines during winter diapause and rises during the following summer when prey are available. There is no equivalent peak during the third summer, indicating that by this time cumulative predation mortality is so high that there is no net advantage to continuing to forage before spawning.

The maximum value of F for most cohorts (^*, Figure 3C) comes at ~1.5 year into the simulation, at the peak in prey availability following maturation. This point in the annual cycle, however, does not fall within the window of spawning dates at which maturation is possible (compare year 2 in Figure 3C with year 1 in Figure 3B), and thus is an example of “internal life history mismatch” (Varpe et al., 2007), the common situation in which the spawning timing that maximizes egg production by the parent is not optimal for the offspring. The long-term egg fitness corresponding to stable 1-year and 2-year cycles is marked for each cohort (Figure 3C, red, orange circles). Some but not all of the cohorts that reach maturity are able to achieve F > 1, egg production above the replacement rate, in these cyclical solutions (solid circles). The best 1-year and 2-year strategies achieve similar maximum fitness values (red vs. orange solid dots), although they require slightly different seasonal timing.

Note that although F can be described as the egg-fitness function, the lines in Figure 3C—time series of F for particular spawning dates t₀—are not the same as the seasonal curve of egg fitness that results from a backwards-in-time dynamic optimization (e.g., Figure 6F in Varpe et al., 2007). Rather, each curve of F(t₀; t) in our approach gives a series of possible values for egg fitness at t₀ depending on what future strategy is taken. The forwards and backwards calculations converge (at least qualitatively) once the internal life-history mismatch is resolved and a stable long-term cycle is found (red and orange circles, Figure 3C). As expected (Varpe et al., 2007), these stable values of egg fitness peak, for each generation length, somewhat prior to the bloom maximum (Figure 3C).

In the global experiment, populations like that shown in Figure 3 were run for a spectrum of u₀ values, across combinations of T₀ and δt from −2 to 16°C and 0 to 310 d (the latter corresponding to δt′ from 0 to 150 d). Across these cases, at a given u₀, the model predicts a log-linear relationship between adult size and temperature, which is not much perturbed by variation in prey availability (Figure 4). The slope of this relationship is equivalent to a Q₁₀ of 1.8–2.0, consistent with that predicted by Equation (32):

Field observations of size in relation to temperature in C. finmarchicus and C. helgolandicus across the North Atlantic show a similar relationship (Q₁₀ = 1.65, Wilson et al., 2015, with prosome length converted to carbon weight based on Runge et al., 2006). Somewhat surprisingly, even wide variation in prey conditions (clusters of gray dots, Figure 4) has only minor effects on this slope.

The intercept of the size-temperature relationship depends on u₀ (Figure 4), with u₀ = 0.005–0.01 d⁻¹ corresponding to the range of adult size from C. finmarchicus to C. hyperboreus at the cold end of the temperature spectrum (Disko Bay, ~0°C: Swalethorp et al., 2011). It is not always fair, however, to associate a particular u₀ value with a particular species over the full range of temperatures included. As Banas and Campbell (2016) discuss further, the temperature response of an individual species is often dome-shaped, a window of habitat tolerance (Møller et al., 2012; Alcaraz et al., 2014), whereas Coltrane 1.0 uses the monotonic, power-law response observable at the community level (Forster et al., 2011). C. finmarchicus, for example, is fit well by u₀ = 0.007 d⁻¹ at higher temperatures (4–12°C), whereas near 0°C in Disko Bay, it has been observed to be considerably smaller than extrapolation along the u₀ = 0.007 d⁻¹ power law would predict. Past studies have also found C. finmarchicus growth and ingestion to be suppressed at low temperatures, i.e., to show a very high Q₁₀ compared with the community-level value (Campbell et al., 2001; Møller et al., 2012).

With this caveat on the interpretation of u₀, we can observe a sensible gradation in life strategy along the u₀ axis (Figure 5). From u₀ = 0.01 d⁻¹ (C. finmarchicus-like at 0°C) to u₀ = 0.005 d⁻¹ (C. hyperboreus-like), the environmental window in which multi-year life cycles are viable (F_1/2 ≥ 1) expands dramatically. This window overlaps significantly with the window of viability for 1-year life cycles (F₁ ≥ 1; Figure 5, black vs. gray contours). In all u₀ cases, there is a non-monotonic pattern in maximum fitness as a function of either temperature or prey (Figure 5, color contours), as environments align and misalign with integer numbers of generations per year or years per generation.

The overall gradient from high to moderate F with increasing temperature (Figure 5) is largely an artifact of displaying F normalized to generation as opposed to per calendar year. In general, these results should not be taken as a quantitative prediction of annual production rates: the linear mortality closure that simplifies the analysis also omits the role of density dependence in stabilizing growth rates. Accordingly, in what follows, we consider only whether F in a given circumstance exceeds replacement rate, not whether it exceeds it modestly or dramatically.

The number of generations per year in the timing strategy that optimizes F for each (T₀, δt) habitat combination is shown in Figure 6 for u₀ = 0.007 d⁻¹. This u₀ value corresponds in adult size to Arctic C. glacialis and temperate C. marshallae populations in the Pacific (Figure 4), species which coexist and are nearly indistinguishable in the Bering Sea. In the lowest-prey conditions, no timing strategy is found to be viable. As prey and temperature increase, the model predicts bands proceeding monotonically from multiple years per generation to multiple generations per year. Validating these model predictions requires parameterizing places (in terms of T₀ and δt) in addition to parameterizing their inhabitants, and thus the meaning of either success of failure is ambiguous. Still, we can observe the following. Ice Station Sheba in the high Pacific Arctic (Figure 1) falls in the non-viable regime (Figure 6), consistent with the conclusion of Ashjian et al. (2003) that Calanus spp. are unable to complete their life cycle there. Disko Bay falls on the boundary of 1- and 2-year generation lengths, consistent with observations of C. glacialis there (Madsen et al., 2001). At Newport, Oregon, near the southern end of the range of C. marshallae, the model predicts multiple generations per year, consistent with observations by Peterson (1979).

These idealized experiments (Figures 5, 6) suggest that very short productive seasons place a hard limit on the viability of Calanus spp., regardless of size, temperature, generation length, or match/mismatch considerations (although these factors affect where exactly the limit falls). A decade of observations in the Eastern Bering Sea provide a unique opportunity to resolve this viability limit with greater precision. This analysis takes advantage of the natural variability on the Southeastern Bering Sea shelf described by the “oscillating control hypothesis” of Hunt et al. (2002, 2011): in warm, low-ice years, the spring bloom in this region is late (~ yearday 150; Sigler et al., 2014) and the abundance of large crustacean zooplankton including C. glacialis/marshallae is very low, while in colder years with greater ice cover, the pelagic spring bloom is earlier, ice algae are present in late winter, and large crustacean zooplankton are much more abundant. The task of replicating these observations serves to test the Coltrane parameterization, and situating them within a complete spectrum of temperature/ice cover cases also allows the model to provide some insight into mechanisms.

Mean surface temperature T0¯ was used to index annual cycles of surface and bottom temperature on the Eastern Bering Sea middle shelf (Appendix in Supplementary Material; insets in Figure 7). Date of ice retreat t_ice was likewise used to index phytoplankton availability over each calendar year (Appendix in Supplementary Material; insets in Figure 7). Coltrane was run for each (T0¯, t_ice) combination with u₀ = 0.007 d⁻¹, thus consistent with Figure 6 except for the more refined treatment of environmental forcing, and an adjustment to K_s to match results of Bering Sea feeding experiments (Campbell et al., in press). The maximum egg fitness F for a one-generation-per-year strategy is shown as a function of T0¯ and t_ice in the main panel of Figure 7. Coltrane predicts that one generation per year is the optimal life cycle length everywhere in this parameter space except for the cold/ice-free and warm/high-ice-cover extremes (white contours), combinations which do not occur anywhere in a model hindcast of middle-shelf conditions back to 1971 (Figure 7, red and blue dots).

Late summer measurements of C. glacialis/marshallae abundance (individuals m⁻²), averaged over the middle/outer shelf south of 60°N, are shown in Figure 7 for 2003–2010 (n = 364 over the 8 years; Eisner et al., 2014, 2015). Both these observations and the predicted maximum F from Coltrane show a dramatic contrast between the warm years of 2003–05 (t_ice = 0) and the cold years of 2007–2010 (t_ice = 100–130), with the transitional year 2006 harder to interpret. Eisner et al. (2014) found that there was less contrast between cold year/warm year abundance patterns on the northern middle/outer shelf, consistent with the model prediction that all hindcast years on the northern shelf fall within the “viable” habitat range for C. glacialis/marshallae (Figure 7, blue dots).

The viability threshhold that the Southeastern Bering Sea appears to straddle is qualitatively similar to that in the more idealized global experiment (Figures 5, 6), primarily aligned with the phenological index (horizontal axis) rather than the temperature index (vertical axis). The threshhold in the Bering Sea experiment (t_ice ≈ 90–100) falls somewhat beyond the dividing line imposed in the experiment setup between early, ice-retreat-associated blooms and late, open-water blooms (t_ice = 75: see Appendix in Supplementary Material, Sigler et al., 2014). This gap (whose width depends on the mortality level m₀: not shown) indicates that some period of ice algae availability is required by C. glacialis/marshallae in this system, in addition to a favorable pelagic bloom timing.

The experiments above test the ability of Coltrane 1.0 to reproduce first-order patterns in latitude and time but do not provide sensitive tests of the model biology. A model case study in Disko Bay, where populations of three Calanus spp. coexist and have been described in detail (Madsen et al., 2001; Swalethorp et al., 2011), allows a closer examination of the relationships among traits within the family of viable life strategies predicted by Coltrane.

The model forcing (Figure 8) describes a single annual cycle, starting with the 1996 spring bloom. This represents a cold, high-ice state of the system, compared with more recent years in which the spring bloom is earlier (e.g., 2008, Figure 8, Swalethorp et al., 2011) and the deep layer is warmed by Atlantic water intrusions (Hansen et al., 2012). This particular year was chosen because measurements of prey availability and Calanus response by Madsen et al. (2001) were particularly complete and coordinated. A simple attempt to correct the prey field for quality and Calanus preference was made by keeping only the >11 μm size fraction of phytoplankton and adding total microzooplankton, in μg C. The measured phytoplankton C:chl ratio was used to convert the sum to an equivalent chlorophyll concentration, and this time series was then slightly idealized for clarity (Figure 8, Appendix in Supplementary Material).

Sensible results were only possible after tuning the predation mortality scale coefficient m₀. It is likely that our simple mortality scheme introduces some form of bias, compared with the reality in this system of predation by successive waves of visual and non-visual predators, which will be considered in a separate study. Still, a sensitivity experiment using the φ model shows that varying m₀ has, as intended, a simple, uniform effect on fitness/population growth (Figure 9) that leaves other trait relationships along the size spectrum unaffected. The φ model predicts that copepods similar to C. finmarchicus in size have much greater fitness at a generation length of 1 year than at 2 years or more; that C. hyperboreus would be unable to complete its life cycle in 1 year, but is well-suited to a 2-year cycle; and that C. glacialis falls in the size range where 1- and 2-year life cycles have comparable fitness value. These results are consistent with observations (Madsen et al., 2001) and more general surveys of life strategies in the three species (Falk-Petersen et al., 2009; Daase et al., 2013).

These results naturally raise the question of whether even lower values of u₀—further reductions in development rate—would produce even larger copepods with even longer life cycles in this environment. C. hyperboreus has been reported to have a life cycle of up to 5 years in other systems (Falk-Petersen et al., 2009) and so the question is more than theoretical. In this version of Coltrane, the lower limit on development rate (and thus the upper limit on adult size) are set by the assumption that modulation of this rate is spread uniformly across the developmental period, rather than concentrated in late copepodid stages, as might be more realistic (Campbell et al., 2001). The largest viable adults in the Disko experiment are those that barely reach D = D_s, the start of reserve accumulation, before a first diapause is required.

For greater specificity, we switched from the φ to the ER model version, running a spectrum of t_egg cases (earliest possible egg-production date: see Section 2.5) along with a spectrum of u₀ (development rate) cases. The predicted “community,” then is the set of all combinations of u₀ and t_egg that lead to a viable level of lifetime egg production. An analog for each of the three Calanus spp. is constructed by averaging model results over the set of viable (u₀, t_egg) cases that predict an adult size within 30% of the average measured adult size for that species. The ER model imposes additional constraints on the model organisms—e.g., they are no longer allowed an infinite egg production rate—and to compensate we reduced m₀ from 0.08 d⁻¹ to 0.06 d⁻¹.

The relationship between generation length and adult size across all (u₀, t_egg) combinations is shown in Figure 10. Results are consistent with the φ model (Figure 9) only a 1-year life cycle is viable for C. finmarchicus in this environment, only a 2-year or longer cycle is viable for C. hyperboreus, and C. glacialis again lies near the boundary where the two strategies are comparable. Note that the model allows for a continuum of intermediate cases in the C. finmarchicus–C. glacialis size range, consistent with the observation of hybridization between these species (Parent et al., 2015).

The ER model also predicts a time series of egg production associated with each trait combination, which we can compare with observations for each species. The model predicts that C. finmarchicus analogs spawn in close association with the spring bloom, that C. hyperboreus spawns well before the spring bloom, and that C. glacialis is intermediate (Figures 11, 12A). These patterns are all in accordance with Disko Bay observations (Madsen et al., 2001; Swalethorp et al., 2011), although the absolute range is muted: Madsen et al. (2001) report C. hyperboreus spawning as early as February. As one would expect from these timing patterns, the model predicts a significant trend between size and the capital fraction of total egg production E_cap/(E_inc + E_cap) (Figure 12C). Again, the pattern is qualitatively correct but muted: Coltrane predicts 80% income breeding at the size of C. finmarchicus (a pure income breeder in reality) and 80% capital breeding at the size of C. hyperboreus (a pure capital breeder in reality). More notable than the error is how much of the income/capital spectrum can apparently be reproduced as a consequence of optimizing reproductive timing alone Varpe et al. (2009), without imposing the physiological difference between the two strategies as an independent trait (Ejsmond et al., 2015).

The model predicts (Figure 12B) that the largest model organisms, with the longest generation lengths, enter their first diapause near the boundary between copepodite stages C4 and C5 (D ≈ 0.75), whereas smaller organisms enter first diapause well into stage C5. Madsen et al. (2001) found that both C. glacialis and C. hyperboreus diapause as C4, C5, and adults in Disko Bay, suggesting that the model is biased toward fast maturation. The discrepancy could also be related to intraspecific variation in the real populations or non-equiproportional development in the late stages, i.e., a variable conversion scale between actual developmental stage and D.

Finally, the ER version of Coltrane allows an estimate of the fraction of individual carbon in the form of storage lipids R/(R + S) (Figure 12D). Averaging each model population from the first diapause-capable stage D = D_dia through adulthood, weighted by survivorship N, yields an overall range that compares well with the species-mean wax ester fractions measured by Swalethorp et al. (2011): ~30% for C. finmarchicus to ~60% for C. hyperboreus. In the middle of the size spectrum, reserve fraction is highly variable across viable 2-year strategies, a warning that the success of this final model prediction may be partly fortuitous. Still, taken as a whole, this experiment has yielded a striking result: that a small set of energetic and timing contraints is able to correctly predict, a priori, that Disko Bay should be able to support a spectrum of calanoid copepods from income breeders with an adult size ~100 μg C, a 1-year life cycle, and a wax ester fraction ~30% to capital breeders with an adult size ~1000 μg C, a two-or-more-year life cycle, and a wax ester fraction ~60%.

In the results above, whether prey availability is treated simply (Figures 5, 6) or with site-specific detail (Figure 7), it appears that the viability of the calanoid community near its high-latitude limit is more sensitive to prey abundance and phenology than to temperature. This result is heuristically similar to the conclusions of Ji et al. (2012) and Feng et al. (2016), although the exact physiological mechanisms differ. Alcaraz et al. (2014) suggested based on lab experiments that C. glacialis reaches an bioenergetic limit near 6°C, and Holding et al. (2013) and others have hypothesized that thermal limits will produce ecosystem-level tipping points in the warming Arctic. Our results, in contrast, suggest that thermal tipping points, even if present at the population level, do not generalize to the community level in copepods. Rather, the model predicts complete continuity between the life strategy of Arctic C. glacialis and temperate congeners like C. marshallae (Figure 6). It also suggests that even on the population level in the Bering Sea, warm/cold-year variation in prey availability is a sufficient explanation of variability in the abundance of C. glacialis/marshallae (Figure 7), without the invocation of a thermal threshhold.

Both the global and Bering experiments suggest, furthermore, that increasing water temperature per se is not necessarily a stressor on copepod communities, even high-latitude communities. In both cases, the low-prey viability threshhold actually relaxes (i.e., is tilted toward lower prey values) as temperature increases, indicating that in these testbeds, the positive effect of temperature on growth and maturation rate actually outweighs the effect of temperature on metabolic losses and overwinter survival (This result may be reliant on the model assumption that the stage of first diapause is highly plastic). In cases where deep, overwintering temperatures increase faster than surface temperatures (Hansen et al., 2012) this balance may not hold, and in the real ocean changes in temperature are highly confounded with changes in phytoplankton production and phenology. Still, it is notable that the model predicts that warming temperatures will have a non-monotonic effect on copepod populations (∂F/∂T₀ ≷ 0, Figures 5, 6) even when metabolic thermal threshholds sensu Alcaraz et al. (2014) and changes in prey availability are not considered. These results are a caution against overly simple climate-impacts projections based on temperature alone.

The biology in Coltrane could be refined in many ways, but two issues stand out as being both mechanistically uncertain and sensitive controls on model behavior. These correspond to the two parameters that it was necessary to tune among model experiments (Table 2) the obstacles to formulation of a fully portable scheme that could produce accurate results across the full range of environments considered here with a single parameterization.

The first of these is the perennial problem of the mortality closure. We modeled predation mortality as size-dependent according to the same power law used for ingestion and metabolism, a choice which is mathematically convenient and makes the effect of top-down controls, if not minor, then at least simple and easy to detect (Figure 9). This size scaling is consistent with the review by Hirst and Kiørboe (2002) but that study also shows that the variation in copepod mortality not explained by allometry spans orders of magnitude (cf. Ohman et al., 2004). Indeed, in some cases one might posit exactly the opposite pattern, in which mortality due to visual predators like larval fish increases with prey body size (Fiksen et al., 1998; Varpe et al., 2015). This latter pattern is one hypothesis for why in reality C. hyperboreus is confined to high latitudes, whereas the model predicts no southern (warm, high-prey) habitat limit to C. hyperboreus analogs based on bottom-up considerations (Figure 5). Merging Coltrane 1.0 with a light- and size-based predation scheme similar to Varpe et al. (2015) or Ohman and Romagnan (2015) would allow one to better test the balance of bottom-up and top-down controls on calanoid biogeography.

Second, our experience constructing the Bering Sea and Disko Bay cases suggests that the greatest uncertainty in the model bioenergetics is actually not the physiology itself—empirical reviews like Saiz and Calbet (2007), Maps et al. (2014), Kiørboe and Hirst (2014), and Banas and Campbell (2016) have constrained the key rates moderately well—but rather the problem of translating a prey field into a rate of ingestion. Within each of our model testbeds, the prey time series P remains subject to uncertainty in relative grazing rates on ice algae, large and small pelagic phytoplankton, and microzooplankton, despite a wealth of local observations and a history of work on this problem in Calanus specifically (Olson et al., 2006; Campbell et al., 2009, in press). The precision of each testbed, and even moreso the ambition of generalizing across them, is also limited by uncertainty in the shape of the functional response (Frost, 1980; Gentleman et al., 2003), here represented by a half-saturation coefficient, which has not been found to be consistent across site-specific studies (Campbell et al., in press; Møller et al., 2016) or well-constrained by general reviews (Hansen et al., 1997). This ambiguity is perhaps not surprising when one considers that ingestion as a function of chlorophyll or prey carbon is not a simple biomechanical property, but in fact a plastic behavioral choice. Accordingly, it might well be responsive not only to mean or maximum prey concentration but also to the prey distribution over the water column, the tradeoff between energy gain and predation risk (Visser and Fiksen, 2013), prey composition and nutritional value, and the context of the annual routine. These issues are fundamental to concretely modeling the effect of microplankton dynamics on mesozooplankton grazers. Addressing them systematically in models will require novel integration between what could be called oceanographic and marine-biological perspectives on large zooplankton.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer ØF and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.