Competition-induced starvation drives large-scale population cycles in Antarctic krill.

Antarctic krill (Euphausia superba) - one of the most abundant animal species on Earth - exhibits a 5-6 year population cycle, with oscillations in biomass exceeding one order of magnitude. Previous studies have postulated that the krill cycle is induced by periodic climatological factors, but these postulated drivers neither show consistent agreement, nor are they supported by quantitative models. Here, using data analysis complemented with modeling of krill ontogeny and population dynamics, we identify intraspecific competition for food as the main driver of the krill cycle, while external climatological factors possibly modulate its phase and synchronization over large scales. Our model indicates that the cycle amplitude increases with reduction of krill loss rates. Thus, a decline of apex predators is likely to increase the oscillation amplitude, potentially destabilizing the marine food web with drastic consequences for the entire Antarctic ecosystem.

Introduction

Abundances of Antarctic krill can fluctuate inter-annually over an order of magnitude on both local1–4 and regional scales5,6, constituting one of the world largest population cycles in absolute biomass. This cycle has important consequences for the entire Antarctic ecosystem, because krill occupies a key position in Southern Ocean food webs7,8, acting as a direct link between primary producers and apex predators.

A clear example of such oscillations is represented in the time series of summer krill abundance and length distributions collected over a period of 18 years4,9 on the Palmer long-term ecological research (Palmer-LTER) grid in the Western Antarctic Peninsula (WAP) region (see supplementary materials, Methods). During this time, krill population exhibited nearly four complete population cycles in length distribution, recruitment index, biomass and abundance with a period of 5-6 years (Fig. 1, left panel). An earlier monitoring program on Adélie penguin diet10, starting already in 1987, shows evidence of five such cycles. Typically, cycles start with the spawning of new strong cohorts within two successive years11 (Fig. 1a), followed in the next years by peaks in the recruitment index (Fig. 1c) and an abrupt increase in krill biomass and abundance4 (Fig. 1e). Although decreasing in abundance, these cohorts dominate the population during the following 3-4 years, until new strong cohorts can prevail and the cycle repeats.

Fig. 1

Anatomy of the krill population cycle, PAL-LTER vs. model.

a, b, Krill length density distributions in color coding, during summer. Vertical lines indicate spawning of strong krill cohorts4. c, d, Recruitment index; maxima occur one year after the successful spawning (grey shading shows ±1 standard deviation across the PAL-LTER grid lines). e, f, Total abundance (black) and biomass (red) of juvenile and adult krill (shading shows ±1 standard deviation). g, h, Negative effect of total krill biomass (adults + juveniles) on the abundance of juveniles in the following year, calculated in (g) for different lines of the PAL-LTER grid (grey) and for the whole grid (black). The insets show a box plot of the same data divided into two groups (Small, Large) with the total krill biomass either smaller or larger than 3 mg DW/m3 in (g) and 20 mg DW/m3 in (h). The simulations (right panel) assume a scenario without inter-annual variations in algal productivity. Field data are presented in Supplementary Table 5, see Methods for data analysis, and model description.

The mechanisms underlying this large amplitude population cycle remain a topic of strong debate. A schematic representation of the krill life cycle (Supplementary Fig. 1b) shows that krill reproduction and maturation depend on the food level, which in turn is driven either directly by climate variability or indirectly (in a feedback loop) by consumption. The origin of the krill cycle has usually been attributed to the direct link caused by periodic climatological influences affecting krill fitness, such as the El Niño-Southern Oscillation index4, Southern Annular Mode9, anomalies in primary production12, sea ice duration6,11 or sea surface temperature3.

Results

Climatological influences, however, do not explain the idiosyncratic pattern of the krill cycle: the occurrence of two successive strong year classes each followed by successful recruitment one year later, and the co-occurrence of the spawning of any new cohort with the extinction of the existing strong cohort. Nor can the climatic factors explain the observed negative effect of the total krill biomass, or of the biomass of reproductive females, on reproduction (Fig. 1g and Supplementary Fig. 1b). As shown in Fig.1g inset, an increase in the total krill biomass above 3 mg DW/m3 results in a small value (mean = 0.006 ind/m3) of juveniles from a narrow range (std = 0.009 ind/m3) in the following summer. In contrast, when total krill biomass is small (<3 mg DW/m3) juvenile abundance in the following summer is highly variable (mean = 0.1 ind/m3, std = 0.098 ind/m3), which is a significant difference in the mean value (t-test, p = 0.0038) and variance (F-test, p < 0.001). This suggests that recruitment is constrained by intraspecific competition when krill biomasses exceed a critical level, but depends on the environmental conditions when the krill biomass is small – an effect that stresses the influence of a feedback link of krill biomass on food level.

To evaluate to what extent intracohort and intercohort interactions contribute to the krill cycle, and how this biotic self-regulation works alongside bottom-up climate regulators, we propose and analyze a population dynamics model of Antarctic krill (Table 1, Supplementary Fig. 1a). We use a bioenergetic model to capture the effects of seasonality on reproduction and ontogenetic development of krill during its entire life cycle. The model relates the energetic demands of growth and fertility to the difference between ingestion and maintenance rates. Krill dominantly feed on pelagic phytoplankton and ice algae13, and our model assumes that all krill stages compete for a food resource (phytoplankton) during spring, summer and autumn, but that during winter (when primary productivity is limited by light) adult and one-year juveniles overwinter in deep waters, while larvae are associated with sea ice and feed on the sea ice community14 (Table 1). We also assume that the starvation risk is significant for larval krill due to their low body lipid levels15,16, while adults and juveniles can accumulate large lipid reserves to survive longer periods of starvation17. Environmental variability is implemented as stochastic inter-annual changes of phytoplankton carrying capacity. A detailed model description and parameters are provided in Supplementary Information and Supplementary Tables 1-4.

Table 1

Main stages of krill development (see Methods for the model description)

Stage	Embryos	Larvae	Juveniles	Adults
Age	0-30 days	30-365 days	>1 year	>2 years
Consumption	no	ice algae or phytoplankton	phytoplankton	phytoplankton
Starvation mortality	no	yes	yes, if size < 10mm	no
Reproduction				If size > 35 mm

Parametrized for WAP conditions our model provides results that are in excellent agreement with observations (Fig. 1, right panel). The model predicts population cycles with periods of about 5-6 years that occur even in the absence of inter-annual variability in plankton productivity. These cycles are related to cyclic increases in starvation mortality of larval krill (Fig. 2a) and are induced by the interplay between the seasonality of the Antarctic environment with the ontogenetic structure of krill populations (Fig. 2c, d). The seasonal succession of phytoplankton and ice algae – the two main food resources – leads to seasonal variations in the grazing impact of krill on phytoplankton, with two local maxima (up to 70% of primary production18) during spring and autumn when the ingestion rates of krill are high but primary production is small. This periodic depletion of resources gives rise to cyclic peaks in starvation mortality among larvae and juveniles. While during spring larvae can still profit from sea ice, during autumn there is typically a time gap between the formation of ice algae after the breakdown of phytoplankton, making autumn a potential bottleneck for krill survival. Thus, although krill affect phytoplankton concentrations during the whole year, the top-down control during autumn is the crucial mechanism of the density dependent regulation of reproduction.

Fig. 2

Mechanism of the multi-annual cycles in krill biomass.

a, b, Temporal changes over one period of the cycle, showing the simulated time course of (a), starvation mortality among larvae (orange) and juveniles (green) and (b), total krill biomass (red) and its food resources: phytoplankton (green) and ice algae (orange). c, d, Schematic representation of the annual dynamics for large (c) and small (d) krill biomass, showing phytoplankton (green), ice algae (orange), and krill grazing impact (percentage of primary production grazed by krill) on phytoplankton (red dashed line). The grey bars indicate starvation periods. When the total krill biomass is large, autumn becomes the main bottleneck for population survival.

Consequently, autumn phytoplankton concentrations, and hence the duration and intensity of the starvation periods, are strongly sensitive to krill biomass density. An abundant krill population accelerates the decline of the phytoplankton concentration in autumn and depletes phytoplankton to extremely low concentrations early, leading to a long starvation period (Fig. 2c). In contrast, a low krill population has a smaller impact on phytoplankton density leading to reduced starvation mortality (Fig. 2d), so that autumn phytoplankton concentrations are sufficient for larvae to survive and to reach an advanced stage before winter. After the first recruitment at the onset of a cycle, the population consists mainly of juveniles, which have limited impact on food levels, and a second larvae cohort can survive. However, with further increase of krill biomass, phytoplankton concentrations are reduced to a critical level, which elevates larval starvation mortality and results in recruitment limitation in the following year. A new cohort of larvae thus only survives when the biomass of the dominating cohorts has fallen below a critical level – giving rise to oscillations in biomass with periods corresponding to the lifetime of the dominating krill cohort (Fig. 2b).

This mechanism naturally explains the above described patterns of the krill cycle - two successive years of successful recruitment followed by 3-4 years of reduced recruitment (Fig. 1d); the appearance of a new strong cohort related to the extinction of an old one (Fig. 1b); and the negative effect of high post-juvenile krill biomass on the juvenile abundance one year later (Fig. 1h). The mechanism generating population cycles by food limitation is robust to different forms of grazing, including a scenario when the habitats of adults and juveniles decouple so that phytoplankton is grazed only by larvae (Supplementary Fig. 3). In this case a large cohort of adults affects the phytoplankton concentrations indirectly: They produce numerous larvae (up to 30% of the biomass), which have a strong grazing impact on the phytoplankton at the onset of autumn causing their starvation later.

It is often argued that the peaks in recruitment index correlate with elevated phytoplankton levels during the preceding summer9,12. We observe the same effect in the krill model (Supplementary Fig. 4b), however, here correlation does not imply causation: In the model, variations in phytoplankton concentrations are modulated by changes in population grazing rates during the population cycle, and not vice versa. A more detailed analysis shows that embryos and larvae reach high abundances in the years with low summer phytoplankton concentrations, because these periods match the maxima of the adult population. These cohorts, however, become extinct because survival during autumn, rather than fertility during summer, determines recruitment success (Supplementary Fig. 5).

Our model framework allows addressing the interplay between exogenous and endogenous factors in shaping the krill population cycle (Supplementary Fig. 1b). Model simulations incorporating among-year variations in algal productivity show that the characteristic population cycles are retained under different external disturbance modes, although the cycle dynamics can be strongly affected (Supplementary Fig. 4, left panel). The correlation between summer phytoplankton concentrations and krill abundance in the following year at first slightly decreases with increasing inter-annual variability, but starts to increase monotonically for larger disturbance levels (Supplementary Fig. 4, middle panel). The lack of food can lead to low recruitment even when the total krill biomass is small (Supplementary Fig. 4, right panel). However, irrespective of the summer conditions, the maximal abundance of juveniles always decreases with the total krill biomass in the previous year – clearly indicating the limitation of krill survival by resource competition.

Even though inter-annual environmental variability may not be the main driver of local cycles in krill abundance, it could well be a decisive factor for synchronizing them across a larger scale. If the oscillations were only induced by the population dynamics, we would expect them to be out of phase in different parts of the region. In contrast, field data indicate that the oscillations are synchronized across the entire PAL-LTER sampling grid, over more than 500 km from North to South4. One possible explanation is that populations in different regions are synchronized by random climatological factors (Moran effect)19,20. This hypothesis is confirmed by numerical simulation of krill populations subjected to inter-annual variability in algal productivity: the correlation between two separate populations, starting at different initial conditions, increases with the level of environmental disturbances, leading ultimately to a complete synchronization of two uncoupled populations (Fig. 3).

Fig. 3

The effect of inter-annual variations in algal productivity on krill population dynamics and synchrony.

a-c, The simulated dynamics of two krill populations, P1 (red) and P2 (black), subjected to identical environmental perturbations (levels of Kph are shown by green dotted lines). All plots compare simulations, starting from different initial conditions, after a transient of 80 simulation years. a, In the absence of external perturbations, ΔK = 0, the populations remain out of phase. b, c, With larger environmental variability, ΔK = 120 mg C/m3 (b) and ΔK = 210 mg C/m3 (c), the population cycles are increasingly synchronized, even though they started with different initial conditions (Moran effect). d, e, Response to different levels of external disturbances, indicated by the value of ΔK. d, The power spectrum of a single population (color coding) shows a maximum at a period of around 6 years for the whole range of noise levels. e, The average correlation coefficient Corr(P1, P2) between the two krill populations (±1 standard deviation shown in blue shading) increases with the level of external disturbances and reaches complete synchronization at ΔK > 200 mg C/m3. See Methods for model parameters.

Discussion

Our results demonstrate that the 5 to 6 year cycle between high and low krill abundance is a self-sustaining mechanism, driven by intraspecific food competition, whereas the climatology is a superimposed factor that can modulate its phase and synchronization over large scales. The suggested model was parametrized for WAP region, but can be used as powerful tool to test possible hypotheses on the krill life history in other regions of the Southern Ocean. According to our model, the stability of the krill cycle depends crucially on the delicate balance between krill reproduction and consumption, which is strongly influenced by ongoing changes in the Antarctic environment. In particular, any factor that increases future krill loss rates is likely to dampen the krill cycle (Supplementary Fig. 6). Thus, an increased predation pressure, caused for instance by recovering whale populations, paradoxically might stabilize krill populations. This effect can be further enhanced by the krill consumption by Baleen whales, which releases Fe locked up in the krill body. This released Fe could be a stimulant for phytoplankton growth21, thus potentially reducing the autumn bottleneck and further suppressing the krill cycle. On the other hand, future krill loss rates might well be reduced by a variety of factors, such as the anticipated decrease in other apex krill predators (e.g., seals, penguins, albatross)10,22,23, changes in harvesting of krill, or the separation between krill and predator habitats caused by climate change4. Such reductions in krill loss rates might, in turn, amplify the amplitude of the oscillations in the biomasses of krill and their predators, with the potential to destabilize the Antarctic marine food web24,25 and ecosystem. Possible consequences include increased likelihood of algal blooms, reproductive failure of krill consumers, and opening up a niche for krill competitors such as salps6.

Methods

Data

The Palmer Long Term Ecological Research (PAL LTER) program studies the Western Antarctic Peninsula (WAP) region4,12. A grid of PAL LTER sampling stations extends from Anvers Island (64.77°S, 64.05°W) in the north to approximately 700 km south near Charcot Island (69.45°S, 75.15°W), and from coastal waters to more than 200 km offshore. Grid lines perpendicular to the Peninsula are spaced 100 km apart with coordinates in km from Marguerite Bay (000) to Anvers Inland (600). The stations along each line are spaced 20km apart with coordinates also in km from the continent (000) to offshore (200 and further), see Supplementary Figs. 10-12 for the maps and a detailed description of the sampling methods in Steinberg et al.12 Only lines 200, 300, 400, 500, and 600 and standard stations with numbers -80, -60, … 260 (as annotated in the data file) were included in our study, as for these locations the data were available for the whole period of observations. Thus, the presented data are averaged over a grid, which extends for 400 km from South to North along the Atlantic Peninsula and for more than 200 km from West to East. We analyze data obtained from a series of 22 cruises: one in spring (November 1991) and 21 in summer (1 January–10 February) for krill size distributions (from 1991 to 2012) and for krill abundances (from 1993 to 2013). The sampling frequency across the grid changed with time: in 1991 only 9 stations were sampled, in the period from 1993 to 2008 the sampling intensity was the highest and 36 to 55 stations were sampled each year, and in the last five years from 2009 to 2013 the sampling intensity was strongly decreased and only 13 to 17 stations were sampled per year. Following Ross and Quetin4 we removed the outliers with extremely high krill densities (one measurement in 1993 and in 2002 and two measurements in 1998).

The average abundance and size distribution

The average abundance and variance was estimated as the mean across all grid stations based on the delta distribution26 as this approach provides a more efficient estimator for the mean abundance when sampling data contain a large proportion of zeros, as typically happens for swarming animals such as krill. The implementation of this method in deltadist function of fishmethods package in R was used, but ported into MATLAB. Krill size distribution was calculated with 1mm bin size and we used the delta distribution mean to find the average size distribution across the whole grid.

The recruitment index (Fig. 1c,d) was defined, for the field data, as the fraction of individuals in the population smaller than 30 mm and, for the model, as the fraction of 1 year old individuals among all juvenile and adult krill. The total dry biomass was calculated as the product of krill abundance and krill average weight < w > = ∫ w(L)f(L, t)dL, where w(L) = c L is the relationship between length L and dry weight (see below) and f(L, t) is the normalized krill length distribution at time t.

Environmental interannual disturbances (Fig. 3, Supplementary Fig. 4) were simulated as random inter-annual changes in phytoplankton carrying capacity K. The values K were uniformly distributed in the range [K], where K = 120 mg C/m3 − 0.24 ΔK, K = 120 mg C/m3 + ΔK and ΔK = 0 … 480 mg C/m3 sets the level of inter-annual variability in phytoplankton production. Figures 3a-c show results for K = 120 mg C/m3 and K between 91 … 240 mg C/m3 and 70 … 330 mg C/m3, respectively (corresponding to ΔK = 0,120 and 210 mg C/m3).

In all simulations the populations started with different initial conditions. Thereby, the initial phase difference between population P and P was set to three years (so that the two populations started approximately in antiphase). Figure 3 shows simulation runs after a transient of 80 years; that is only the last 20 years of the total simulation of 100 years are shown.

To calculate the correlation coefficients (Fig. 3e), we first modelled each population independently for a randomly chosen period of time from 0 to 12 years. Then the final state was used as a new initial condition and both populations were subjected to the same environmental signal, K(t), during 100 years. The correlation coefficient was found for the last 20 years of the simulation time and averaged over 60 pairs of populations with random initial phase shifts. Supplementary Figures 4 (top three rows) show results for K = 120 mg C/m3 and K between 77 … 350 mg C/m3 and 30 … 600 mg C/m3, respectively (ΔK = 0,230 and 480 mg C/m3). In Supplementary Figure 4 bottom row, we assume that phytoplankton carrying capacity, K, is driven by the Southern Annular Oscillations index9, namely, K = K,0 exp(−0.3 SAM).

Ontogenetic model of krill population dynamics (Supplementary Fig. 1). We model continuous growth and development of krill cohorts during their life cycle. We assume that growth, maturity and fecundity of individuals in each cohort depend on the food level and on their size and stage, where the food level is a balance between primary production, losses and consumption. Although krill are omnivorous, they feed primarily on phytoplankton and ice algae27. We assume that algal productivity is seasonally driven with a maximum of phytoplankton during summer and a maximum of ice algae during winter (Supplementary Fig. 7).

Krill development includes four functionally different stages: embryos, larvae, juveniles and adults (Table 1). During the embryonal stage, krill do not eat and develop using yolk reserves. This stage lasts approximately 30 days15. The second stage (larvae) lasts approximately one year. In most simulations, we assume that larvae compete with other krill stages for phytoplankton during summer (see, however, Supplementary Fig. 3), but only larvae can feed on ice algae during winter, when adults are inactive14. Larvae have high metabolic rates and cannot survive a long period of starvation15,16. After the first winter krill reach the juvenile stage. Juveniles feed on phytoplankton during summer and increase in size. The starvation mortality decreases with size and approaches zero for 20mm long krill. These individuals can compensate for maintenance during starvation by losing their weight17. Background mortality (including predation and fishing) decreases with size but rises upon reaching the individual age of around 6 years28. Krill metabolic rates depend on temperature, we model this dependence implicitly assuming that the metabolic rate follows a periodic function, and during winter juveniles and adults reduce their metabolic activity up to 20% of the summer level29 (Supplementary Fig. 7a). After the second winter, krill enters the adult stage. Adult krill can reproduce if resources are high enough and they have reached the reproductive body size. To model krill development we use a net-production type energy budget approach, which links the growth and fecundity rates to the difference between ingestion and maintenance rates.

Let w denote the dry weight of individuals, a their age and P the resource availability. Then the evolution of the age-weight distribution c(t, a, w) of krill can be described in terms of a partial differential equation30. where g(P, w) is the daily growth rate in weight and d(P, a, w) is the death rate, which includes predation and starvation mortality.

Since the number of new individuals (of weight w and age 0) equals the amount of eggs produced per unit of time by females above the reproductive weight w and age a we obtain the left boundary condition where μ(P, a, w) is the female fecundity, S is the egg hatchability and the factor 1/2 in front of the integral appears because females constitute approximately half of the entire population.

The growth of individuals is coupled with the resource dynamics. Denote P the density of phytoplankton in the water column and P the density of ice algae. Then the resource dynamics follow the equations

Here the periodic functions g (t) and g_ice (t) modulate the seasonal growth of phytoplankton and ice algae; K and K are the corresponding carrying capacity for each resource; I(P) and I(P) the ingestion rates; ρ a concentration factor relating the krill densities, when feeding on ice algae in a small volume below the ice, to their density within the water column; l the loss rate of the resource via sinking and grazing by other zooplankton species; δ and δ the input from external sources, such as adjacent patches. To simulate environmental disturbances we assume that different summers are characterized by random values of K.

The model was parameterized using the latest available data (Supplementary Tables 1-4) and simulated with a numerical method for the integration of physiologically structured models31.

Relationship between krill length L and dry mass

During the embryonal and the following non-feeding stages krill individuals increase their length, but the weight stays approximately unchanged, w. At later stages typically w~L, with an exponent σ that increases with krill size (Supplementary Fig. 8a)32–35. To match the relationship throughout the entire range, we assume a quadratic relation between the logarithm of weight and logarithm of length, ln w = a + b ln L + c ln2 L, which can also be written as

The inverse relation expresses the krill length as a function of krill weight

We use the last equation to compare the modeled krill traits, expressed as a function of krill weight, with field measurements typically expressed as a function of krill length.

Growth rate and fecundity are defined in terms of a bioenergetics approach. We assume that assimilated food after covering maintenance requirements is invested into either body growth or reproduction. Namely, a fraction K(w) of the difference between assimilated carbon and maintenance costs is invested into growth, and the rest (1 − K(w)) into maturation or reproduction processes, where K(w) is a decreasing function of krill size (Supplementary Fig. 8b). During starvation krill does not invest into maturation, and the maintenance is covered at the expense of body reserves.

To formulate the bioenergetics model we express the net production rate A(P, w) as the difference between the assimilation rate and maintenance costs where I(P, w) is the weight specific ingestion rate, ε the conversion coefficient from the assimilated ingested carbon into the dry body weight equivalents, and T(w) the weight specific maintenance rate.

Then the daily growth rate in weight equals Thus, krill grow in size when the assimilation net production is positive and shrink in size for small resource levels17.

The number of eggs produced by a female of the reproductive size and age equals where ε is the relative weight eggs in ovary tissue to the total dry weight of ovary tissue. Because A(P, w) depends on the food level, reproduction is possible only when food level is sufficient to cover the maintenance rate. Consequently, the combination of grazing impact and seasonal variation in plankton production preclude reproduction in winter.

For the function K(w) we use the logistic function36 where the parameter L defines length at which 50% of krill stock attains maturity and k characterizes the steepness of the decrease of K(L) with increasing L. For simplicity, we assume the same splitting function K(w) for male and female krill.

Ingestion rates

Laboratory experiments show that the ingestion rate likely follows a sigmoidal Holling type 3 functional response37,38, but field data indicate that the individual growth rate instead follows Holling type 2 functional response33,39. Using different functional responses for growth and ingestion is problematic because it implies e.g. that at small resource concentration krill consumption rate is small, but the growth rate can be relatively high. To resolve this we exploit the observation that the average grazing rate in a heterogeneous environment does not necessarily match the local grazing rate, measured in short-term experiments40, but should match on a global scale. Thus, to obtain the correct dependence of the growth rate on the resource availability we assume that the weight specific ingestion rate also follows a Holling type 2 function of the food concentration (which can be interpreted as an effective ingestion rate in a heterogeneous environment) where I is the maximal weight specific ingestion rate, P is the resource concentration, H is the half-saturation constant. For simplicity, we assume that I and H do not depend on krill weight and age41 (see Supplementary Table S2 for parameters).

For krill larvae, which can feed on ice algae or phytoplankton, we assume an active switching strategy42,43 where β = P/(P + P) are the dynamical preferences for resource i (with i = Ph or i = Ice).

Our model neglects the dynamic feedback on phytoplankton and ice algae from other grazers, such Thysanoessa or salps, which also play an important role in WAP region44. However, krill is likely to be the main grazer in its habitat, because krill typically dominate in the costal and shelf region, while salps dominate offshore, and krill and salp abundances negatively correlate12.

Background maintenance rate

According to equation (3) the change of weight of a starving adult equals where T is the weight-specific background maintenance rate for adults. Therefore, we use the relation w(t) = w(0)e and determine T to fit the results of starvation experiments17 (Supplementary Fig. 9d). Starving larvae have an approximately three times higher oxygen consumption rate in comparison with starving adults16,17. Therefore, we assume for larvae T = 3T and a gradual decrease of the maintenance rate from T to T with increasing krill size following the logistic function where L is the characteristic size of larvae, at which this transition should occur and k is the steepness of this transition.

Daily growth rate

Experimental data, typically, provide the daily increase in body length, dL/dt, and for our model we need to determine the increase in body weight, dw/dt. To convert the units we use the rule for the derivative of a composite function Substituting equations (1) and (2) into this expression and calculating all derivatives we obtain the maximal daily growth rate in length

Assimilation efficiency, ε, was chosen in such a way that the maximal daily growth rate (5) follows a hull above the experimental data points (the red line in Supplementary Fig. 9a). The other lines in this figure show the model predictions for the growth rate for different levels of food availability.

Fecundity as a function of body size can now be determined according to equation (4). Supplementary Figure S9b show a good match between the experimental data45 and model predictions. Note, Tarling at al.45 calculated the number of mature oocytes, which can be more than 15000 per female and can be released over 2 or 3 partial spawns. To estimate the daily fertility rate we divided the number oocytes by 45 days (one half of summer duration).

Krill mortality is a sum of background mortality d(w), which depends on krill weight, and starvation mortality, d(P, w), which depends on krill weight and on available resources,

We assume that the background mortality decreases with stage, and krill has a finite life span of 5.9 years46. We use a fractional number for the life span to avoid a resonance between annual cycle in the model and life span of krill.

Pakhomov (1995)46 found that background mortality rate per year ranges from 0.52, during the maturation period, to 1.1, during the first year. We assume that mortality changes stepwise from 1 year-1 for juveniles to 0.5 year-1 for adults. From the experimental data the background mortality can be found as , where θ is the percentage of the population survived after time t. In the experiments by Ross et al.15 after 60 days approximately 20% of the population survived, which corresponds to m = 7.3 year-1. However, these experiments included only the first feeding stages (Calyptopis 1 to 3) and the authors mentioned that at the end of the period the population stabilized, so that the mortality rate decreased. The highest field estimation of mortality 2.31 year-1, was found by Brinton and Townsend47 for one year old krill. To keep the background mortality between these two extremes we used the value of 5 year-1 for larvae.

The data on starvation mortality are even scarcer; in particular, no data is available on the dependence of starvation mortality on krill size and food availability. Thus, we parametrize the model to match the available data and use some additional assumptions to extend the dependences in the range, where no data is available. Furthermore, as our model assumes an exponential mortality rate, we can only specify a proportion (<100%) of the population which dies during a period of time, while in starvation experiments typically the entire population dies after some period. In particular, it was found that the first feeding stages (Calyptopys 1) die during 10-14 days when starving15,16. To match this result we assume that during 14 days approximately 80% of the population dies, which corresponds to the maximal starvation mortality rate, d = 40 year-1.

The maximal starvation mortality decreases with krill size, because adults can survive long periods of starvation. We assume that it decreases in the same way as the krill maintenance rate

Finally, if assimilation is positive, but not sufficient to cover the individual maintenance costs, we assume that the mortality rate follows a linear dependence on the relative difference between lack of resources for maintenance with respect to the maintenance rate (Supplementary Fig. 9c)

Metabolic rate of adults changes with the season29. We assume that the relative metabolic rate changes between f and 1, and is synchronized with the seasonality. To model this periodic dependence we use a modified von Mises distribution function where t is the summer midpoint where f reaches a maximum, T is the duration of the active period, and time is measured in years. A decrease of metabolic rate causes a proportional decrease of the maximal ingestion rate, which during winter reaches only 20% of the summer levels38. In the model the period of high metabolic rates matches approximately with the period of maximal phytoplankton concentrations.

Phytoplankton and ice algae

Behrenfeld48 noted that net growth rate of phytoplankton differs substantially between the ocean and lab experiments (0.06 day-1 in the ocean vs. 1 day-1 in lab experiments). Furthermore, the phytoplankton growth rate changes abruptly from positive during summer to negative during winter. Therefore, we model the growth rate by a stepwise function with a maximum in summer. To define it mathematically we denote T the duration of the period of phytoplankton growth and t the midpoint of this period (both measured as a fraction of a year). Then the growth is positive if the differences between the current moment of time and the closest midpoint is less than T/2. Mathematically it can be, for instance, formulated as Thus, during the period of positive growth (spring and summer) the net growth rate equals g − l and during autumn and winter the net growth rate equals − l (Supplementary Fig. 7b, black line). We choose g, t, T, l, K to match the experimental data obtained at the coastal station of the Palmer LTER (Supplementary Fig. 7b).

For ice algae we assume that the growth rate is a periodic function of time which reaches a maximum of g at the midpoint of the ice algae production, t, and T characterizes the duration of ice algae growth (see Supplementary Fig. 7b for an exemplary model outcome).

Note that our model shows a strong effect of krill on pelagic phytoplankton but not on ice algae. This is caused by the model parametrization, which assumes a relatively smaller net growth rate of the pelagic phytoplankton in comparison with ice algae.