Vertically structured prokaryotic community can control the efficiency of the biological pump in the oceans

Abstract

Despite recent advances in molecular and metagenomic approaches, it is still unclear how spatiotemporal variations in microbial communities influence the biological pump, exporting organic carbon from the surface to the deep oceans. We address this important problem by constructing a simple model of a prokaryotic metacommunity in which two generalist ecotypes compete for two resources. The first and the second ecotypes have, respectively, a higher preference for particulate organic carbon (POC) and dissolved organic carbon. Sinking of POC would lead to a higher abundance of the first ecotype in the deep layer compared to the surface layer, but vertical mixing in turn generates a net upward flux of this ecotype to the surface layer. This upward movement accelerates the shifts in the community composition during the phytoplankton bloom, contributing to a higher efficiency in POC remineralization at the surface layer and reducing the carbon flux to the deep layer.

Appendices

Appendix 1

Local stability analysis of the surface-layer system with m = 0

We describe the stability analysis of our system by first examining the local equilibria at the surface layer. The Jacobian matrix of the variables corresponding to the surface layer is given by

$${\mathbf{J}}_{{\mathbf{SS}}} = \left[ {\begin{array}{*{20}c}a & 0 & b & c \\0 & d & e & f \\{ - g} & { - h} & { - i} & 0 \\{ - j} & { - k} & 0 & { - l} \\\end{array} } \right],$$

(11)

where $a = - {\left( {\kappa _{B} + \gamma } \right)} + \varepsilon {\left( {\omega _{1} \rho ^{*}_{{1{\text{P}}}} + \delta _{1} \rho ^{*}_{{1{\text{D}}}} } \right)}$, $b = \varepsilon \omega _1 B_{11}^*$, $c = \varepsilon \delta _1 B_{11}^ *$, $d = - {\left( {\kappa _{B} + \gamma } \right)} + \varepsilon {\left( {\omega _{2} \rho ^{*}_{{1{\text{P}}}} + \delta _{2} \rho ^{*}_{{1{\text{D}}}} } \right)}$, $e = \varepsilon \omega _2 B_{12}^*$, $f = \varepsilon \delta _2 B_{12}^*$, $g = \omega _1 \rho _{1{\text{P}}}^ *$, $h = \omega _2 \rho _{1{\text{P}}}^ *$, $i = \kappa + \beta + \omega _1 B_{11}^* + \omega _2 B_{12}^*$, $j = \delta _1 \rho _{1{\text{D}}}^*$, $k = \delta _2 \rho _{1{\text{D}}}^*$, $l = \beta + \delta _1 B_{11}^* + \delta _2 B_{12}^*$, with b, c, e, f, g, h, i, j, k ≥ 0.

At E _S0, b = c = e = f = 0, and

$$\begin{array}{*{20}c} {\rho ^{ * }_{{1{\text{P}}}} = \frac{1}{{\kappa + \beta }}s_{{\text{P}}} \equiv \rho ^{ * }_{{1{\text{P}}}} {\left( {E_{{{\text{S}}0}} } \right)},} \\ {\rho ^{ * }_{{1{\text{D}}}} = \frac{1}{\beta }s_{{\text{D}}} \equiv \rho ^{ * }_{{1{\text{D}}}} {\left( {E_{{{\text{S}}0}} } \right)}.} \\ \end{array} $$

The eigenvalues of the Jacobian matrix at E _S0 are a, d, −i, and −l, from which it follows that E _S0 is locally stable if and only if a < 0 and d < 0 or, equivalently, if and only if

$$- \left( {\kappa _B + \gamma } \right) + \varepsilon \left[ {\omega _1 \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}0} } \right) + \delta _1 \rho _{1{\text{D}}}^ * \left( {E_{{\text{S}}0} } \right)} \right] <0,$$

(12)

$$- \left( {\kappa _B + \gamma } \right) + \varepsilon \left[ {\omega _2 \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}0} } \right) + \delta _2 \rho _{1{\text{D}}}^ * \left( {E_{{\text{S}}0} } \right)} \right] <0.$$

(13)

The above relations correspond to the invasibility conditions for each bacterial ecotype at equilibrium: ecotype 1 (or 2) cannot invade the system when condition 12 (13) holds, that is to say, when $\left( {{1 \mathord{\left/ {\vphantom {1 {B_{1j} }}} \right. \kern-\nulldelimiterspace} {B_{1j} }}} \right)\left( {{{{\text{d}}B_{1j} } \mathord{\left/ {\vphantom {{{\text{d}}B_{1j} } {{\text{d}}t}}} \right. \kern-\nulldelimiterspace} {{\text{d}}t}}} \right) <0$ with j = 1 (or j = 2).

For the equilibrium when only one bacterial ecotype persists (E _Sj with j = 1 or j = 2), we have H _j  = H ₃ = H ₄ = 0 with $B_{1j}^* > 0$ (j = 1 or 2), resulting in a quadratic equation for $B_{1j}^ *$ whose solution gives

$$\begin{array}{*{20}l}{B_{1j}^ * \left( {E_{{\text{S}}j} } \right) = \frac{{M_{{\text{S}}j} + \sqrt {M_{{\text{S}}j}^2 + 4L_{{\text{S}}j} N_{{\text{S}}j} } }}{{2L_j }},} \hfill \\{\rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}j} } \right) = \frac{{s_P }}{{\omega _j B_{1j}^ * \left( {E_{{\text{S}}j} } \right) + \kappa + \beta }},} \hfill \\{\rho _{1{\text{D}}}^ * \left( {E_{{\text{S}}j} } \right) = \frac{{s_{\text{D}} }}{{\delta _j B_{1j}^ * \left( {E_{{\text{S}}j} } \right) + \beta }},} \hfill \\\end{array} $$

wherein

$$\begin{array}{*{20}l}{L_{{\text{S}}j} = \left( {\kappa _B + \gamma } \right)\omega _j \delta _j ,} \hfill \\{M_{{\text{S}}j} = \varepsilon \omega _j \delta _j \left( {s_{\text{P}} + s_{\text{D}} } \right) - \left( {\kappa _B + \gamma } \right)\left( {\omega _j b + \delta _j b + \delta _j \kappa } \right),} \hfill \\{N_{{\text{S}}j} = \varepsilon \omega _j \beta s_{\text{P}} + \varepsilon \delta _j \left( {\kappa + \beta } \right)s_{\text{D}} - \beta \left( {\kappa + \beta } \right)\left( {\kappa _B + \gamma } \right).} \hfill \\\end{array} $$

Focusing only on the stable equilibria, we can conclude that ecotype 1 (or 2) can invade (i.e., the condition 12 (or 13) does not hold at E _S0), if N _Sj  > 0.

At E _S1, e = 0 and f = 0. By applying the Routh–Hurwitz conditions to the characteristic equation of the Jacobian matrix (Murray 2002), we can show that E _S1 is locally stable if and only if d < 0 or, equivalently, if and only if

$$\omega _2 \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}1} } \right) + \delta _2 \rho _{1{\text{D}}}^ * \left( {E_{{\text{S1}}} } \right) <\varepsilon ^{ - 1} \left( {\kappa _B + \gamma } \right).$$

(14)

Similarly, the condition for local stability of E _S2 is written as

$$\omega _1 \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}2} } \right) + \delta _1 \rho _{1{\text{D}}}^ * \left( {E_{{\text{S2}}} } \right) <\varepsilon ^{ - 1} \left( {\kappa _B + \gamma } \right).$$

(15)

Eq. 14 (or 15) is equivalent to the invasibility conditions for ecotype 2 (or 1), e.g., ecotype 2 cannot invade the system when condition 14 holds because $\left( {{1 \mathord{\left/ {\vphantom {1 {B_{12} }}} \right. \kern-\nulldelimiterspace} {B_{12} }}} \right)\left( {{{{\text{d}}B_{12} } \mathord{\left/ {\vphantom {{{\text{d}}B_{12} } {{\text{d}}t}}} \right. \kern-\nulldelimiterspace} {{\text{d}}t}}} \right) <0$ at the equilibrium E _S1.

At E _SC, we obtain the equilibrium values of each variable as follows:

$$\begin{array}{*{20}c} {{\left( {\begin{array}{*{20}c} {{\rho ^{ * }_{{1{\text{P}}}} {\left( {E_{{{\text{SC}}}} } \right)}}} \\ {{\rho ^{ * }_{{1{\text{D}}}} {\left( {E_{{{\text{SC}}}} } \right)}}} \\ \end{array} } \right)} = \frac{{\varepsilon ^{{ - 1}} {\left( {\kappa _{B} + \gamma } \right)}}}{{\omega _{1} \delta _{2} - \omega _{2} \delta _{1} }}{\left( {\begin{array}{*{20}c} {{\delta _{2} - \delta _{1} }} \\ {{\omega _{1} - \omega _{2} }} \\ \end{array} } \right)}} \\ {{\left( {\begin{array}{*{20}c} {{B^{ * }_{{11}} {\left( {E_{{{\text{SC}}}} } \right)}}} \\ {{B^{ * }_{{12}} {\left( {E_{{{\text{SC}}}} } \right)}}} \\ \end{array} } \right)} = \frac{1}{{\omega _{1} \delta _{2} - \omega _{2} \delta _{1} }}{\left( {\begin{array}{*{20}c} {{\delta _{2} {\left( {{s_{{\text{P}}} } \mathord{\left/ {\vphantom {{s_{{\text{P}}} } {\rho ^{ * }_{{1{\text{P}}}} }}} \right. \kern-\nulldelimiterspace} {\rho ^{ * }_{{1{\text{P}}}} }{\left( {E_{{{\text{SC}}}} } \right)} - \kappa - \beta } \right)} - \omega _{2} {\left( {{s_{{\text{D}}} } \mathord{\left/ {\vphantom {{s_{{\text{D}}} } {\rho ^{ * }_{{1{\text{D}}}} }}} \right. \kern-\nulldelimiterspace} {\rho ^{ * }_{{1{\text{D}}}} }{\left( {E_{{{\text{SC}}}} } \right)} - \beta } \right)}}} \\ {{ - \delta _{1} {\left( {{s_{{\text{P}}} } \mathord{\left/ {\vphantom {{s_{{\text{P}}} } {\rho ^{ * }_{{1{\text{P}}}} }}} \right. \kern-\nulldelimiterspace} {\rho ^{ * }_{{1{\text{P}}}} }{\left( {E_{{{\text{SC}}}} } \right)} - \kappa - \beta } \right)} + \omega _{1} {\left( {{s_{{\text{D}}} } \mathord{\left/ {\vphantom {{s_{{\text{D}}} } {\rho ^{ * }_{{1{\text{D}}}} }}} \right. \kern-\nulldelimiterspace} {\rho ^{ * }_{{1{\text{D}}}} }{\left( {E_{{{\text{SC}}}} } \right)} - \beta } \right)}}} \\ \end{array} } \right)}.} \\ \end{array} $$

The choice of ecotype 1 as a better consumer of POC and ecotype 2 as a better consumer of DOC implies that ω ₁ > ω ₂ and δ ₂ > δ ₁ and that $\rho _{1{\text{P}}}^ * \left( {E_{{\text{SC}}} } \right)$ and $\rho _{1{\text{D}}}^ * \left( {E_{{\text{SC}}} } \right)$ are always positive. From the existence condition of E _SC, i.e., $B_{11}^ * \left( {E_{{\text{SC}}} } \right) >0$ and $B_{12}^ * \left( {E_{{\text{SC}}} } \right) >0$, it follows that

$$\begin{array}{*{20}l} {{B^{ * }_{{11}} {\left( {E_{{{\text{SC}}}} } \right)} > 0} \hfill} \\ {{ \Leftrightarrow s_{{\text{D}}} < \frac{{\delta _{2} }}{{\omega _{2} }}\frac{{\omega _{1} - \omega _{2} }}{{\delta _{2} - \delta _{1} }}s_{{\text{P}}} + \frac{{\varepsilon ^{{ - 1}} {\left( {\kappa _{B} + \gamma } \right)}{\left( {\omega _{1} - \omega _{2} } \right)}}}{{\omega _{2} {\left( {\omega _{1} \delta _{2} - \omega _{2} \delta _{1} } \right)}}}{\left\{ {\omega _{2} b - \delta _{2} {\left( {\kappa + \beta } \right)}} \right\}}.} \hfill} \\ \end{array} $$

(16)

$$\begin{array}{*{20}l} {B_{12}^ * \left( {E_{{\text{SC}}} } \right) > 0} \hfill \\ { \Leftrightarrow s_{\text{D}} > \frac{{\delta _1 }}{{\omega _1 }}\frac{{\omega _1 - \omega _2 }}{{\delta _2 - \delta _1 }}s_{\text{P}} + \frac{{\varepsilon ^{ - 1} \left( {\kappa _B + \gamma } \right)\left( {\omega _1 - \omega _2 } \right)}}{{\omega _1 \left( {\omega _1 \delta _2 - \omega _2 \delta _1 } \right)}}\left\{ {\omega _1 b - \delta _1 \left( {\kappa + \beta } \right)} \right\}.} \hfill \\ \end{array} $$

(17)

By applying Routh–Hurwitz conditions to the characteristic equation of the Jacobian matrix at E _SC, we can show that E _SC is locally stable if and only if

$$\Delta = \left| {\begin{array}{*{20}c}{\alpha _1 } & {\alpha _3 } & 0 \\{\alpha _0 } & {\alpha _2 } & {\alpha _4 } \\0 & {\alpha _1 } & {\alpha _3 } \\\end{array} } \right| >0.$$

(18)

where α ₀ = 1, $\alpha _1 = i + l$, $\begin{array}{*{20}l} {\alpha _2 = il + eh + kf + bg + cj} \hfill \\ \end{array}$, $\alpha _3 = ehl + kfi + bgl + cij$, $\alpha _4 = \left( {gk - jh} \right)\left( {bf - ce} \right)$, with α ₀, α ₁, α ₂, α ₃, α ₄ > 0 when the coexistence region E _SC is present.

With our choice of parameters (Table 1), numerical evaluation of this condition within our range of s _P and s _D values suggests that this equilibrium is locally stable if Eqs. 16 and 17 are satisfied.

Appendix 2

Local stability analysis of the deep-layer system with m = 0

The Jacobian matrix of the deep-layer system evaluated at equilibrium is given by

$${\mathbf{J}}_{{\text{DD}}} = \left[ {\begin{array}{*{20}c}a & 0 & b & c \\0 & d & e & f \\{ - g} & { - h} & { - i} & 0 \\0 & 0 & 0 & { - j} \\\end{array} } \right],$$

(19)

where $a = - \gamma + \varepsilon \left( {\omega _1 \rho _{2{\text{P}}}^ * + \delta _1 \rho _{2{\text{D}}}^ * } \right)$, $b = \varepsilon \omega _1 B_{21}^ *$, $c = \varepsilon \delta _1 B_{21}^ *$, $d = - \gamma + \varepsilon \left( {\omega _2 \rho _{2{\text{P}}}^ * + \delta _2 \rho _{2{\text{D}}}^ * } \right)$, $e = \varepsilon \omega _2 B_{22}^ *$, $f = \varepsilon \delta _2 B_{22}^ *$, $g = \omega _1 \rho _{2P}^ *$, $h = \omega _2 \rho _{2{\text{P}}}^ *$, $i = \beta + \omega _1 B_{21}^ * + \omega _2 B_{22}^ *$, $j = \beta + \delta _1 B_{21}^ * + \delta _2 B_{22}^ *$, with b, c, e, f, g, h, i ≥ 0.

Note that $\rho _{2\operatorname{D} }^ * = 0$ since H ₈ ≤ 0 when m = 0, irrespective of the equilibrium states at the surface layer. Among all possible equilibria, we just have to focus on a part of them. For example, when E _S0 is locally stable, we just have to study $E_{{\text{S0D0}}} \left( {B_{21}^ * = 0,{\text{ }}B_{22}^ * = 0} \right)$ and $E_{{\text{S0D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$ because $E_{{\text{S0D2}}} \left( {B_{21}^ * = 0,{\text{ }}B_{22}^ * >0} \right)$ and $E_{{\text{S0DC}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * >0} \right)$ cannot be locally stable. In such a situation, ecotype 1 is a better competitor for POC (ω ₁ > ω ₂) and excludes ecotype 2 when DOC is absent ${\left( {\rho ^{ * }_{{2\operatorname{D} }} = 0} \right)}$. Using default parameter values (“Appendix 2”), the existence and stability conditions for the other equilibria can be studied with a combination of analytical and numerical techniques. When E _S0 is locally stable, we consider only the two equilibria $E_{{\text{S0D0}}} \left( {B_{21}^ * = 0,{\text{ }}B_{22}^ * = 0} \right)$ and $E_{{\text{S0D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$ (see the main text). When E _S1 is locally stable, $B_{21}^ *$ is always positive due to the sinking of ecotype 1 $\left( {\kappa _B B_{11}^ * } \right)$ from the surface layer and dB ₂₂/dt is always negative at equilibrium since ω ₁ > ω ₂. In such a case, we consider only $E_{{\text{S1D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$. Similarly, when E _S2 is locally stable, we focus on the equilibrium $\left( {E_{{\text{S2D2}}} \left( {B_{21}^ * {\text{ }} = {\text{ }}0,{\text{ }}B_{22}^ * >0} \right)} \right)$ and $E_{{\text{S2DC}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * >0} \right)$ since $E_{{\text{S2D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$ is always unstable. When E _SC is locally stable, we consider only $E_{{\text{SCDC}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * >0} \right)$ since individuals of two ecotypes are supplied from the surface layer and the other equilibria are always unstable.

At E _S0D0,

$$\rho _{2{\text{P}}}^ * \left( {E_{{\text{S0D0}}} } \right) = \frac{\kappa }{{\beta \left( {\kappa + \beta } \right)}}s_{\text{P}} ,\quad \rho _{2{\text{D}}}^ * \left( {E_{{\text{S0D0}}} } \right) = 0,$$

and b = c = e = f = 0. The eigenvalues of the Jacobian matrix at E _S0D0 are simply its diagonal elements resulting in E _S0D0 being locally stable if and only if a < 0 and d < 0. These relations are equivalent to

$$\omega _1 \rho _{2{\text{P}}}^ * \left( {E_{{\text{S0D}}0} } \right) <\varepsilon ^{ - 1} \gamma ,$$

(20)

which implies that ecotype 1 and ecotype 2 (when ω ₁ > ω ₂) cannot survive at steady state.

For the equilibrium $E_{{\text{S0D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$, we have H ₅ = H ₇ = 0, from which we obtain

$$\begin{array}{*{20}l}{\rho _{2{\text{P}}}^ * \left( {E_{{\text{S0D1}}} } \right) = \varepsilon ^{ - 1} \omega _1 ^{ - 1} \gamma ,} \hfill \\{B_{21}^ * \left( {E_{{\text{S0D1}}} } \right) = \omega _1 ^{ - 1} \left( {\frac{{\kappa s_{\text{P}} }}{{\kappa + \beta }}\frac{1}{{\rho _{2{\text{P}}}^ * \left( {E_{{\text{S0D1}}} } \right)}} - \beta } \right),} \hfill \\\end{array} $$

where $B_{21}^ *$ is positive if the condition (20) does not hold.

At E _S0D1, a = e = f = 0, by Routh–Hurwitz conditions, it is found that E _S0D1 is locally stable if and only if d < 0 (when ω ₁ > ω ₂).

For the equilibrium $E_{{\text{S1D1}}} \left( {B_{21}^ * >0,{\text{ }}B_{22}^ * = 0} \right)$, we have

$$\begin{array}{*{20}l}{\rho _{2{\text{P}}}^ * \left( {E_{{\text{S}}j{\text{D}}j} } \right) = \frac{{\kappa \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}j} } \right)}}{{\omega _j B_{2j}^ * \left( {E_{{\text{S}}j{\text{D}}j} } \right) + \beta }},} \hfill \\{B_{2j}^ * \left( {E_{{\text{S}}j{\text{D}}j} } \right) = \frac{{M_{{\text{D}}j} + \sqrt {M_{{\text{D}}j}^2 + 4L_{{\text{D}}j} N_{{\text{D}}j} } }}{{2L_{{\text{D}}j} }},} \hfill \\\end{array} $$

where

$$\begin{array}{*{20}l}{L_{{\text{D}}j} = \gamma \omega _j ,} \hfill \\{M_{{\text{D}}j} = \kappa _B \omega _j B_{1j}^ * \left( {E_{{\text{S}}j} } \right) + \varepsilon \kappa \omega _j \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}j} } \right) - \beta \gamma ,} \hfill \\{N_{{\text{D}}j} = \kappa _B \beta B_{1j}^ * \left( {E_{{\text{S}}j} } \right),} \hfill \\\end{array} $$

and j = 1.

At E _S1D1, a < 0, e = 0, and f = 0. By Routh–Hurwitz conditions, it is found that E _S1D1 is locally stable if and only if d < 0 (when ω ₁ > ω ₂). A similar procedure applied to E _S2D2 (d < 0, b = c = 0) shows that E _S2D2 is locally stable if and only if a < 0, which is equivalent to the condition of ecotype 1 not being able to invade (dB ₂₁/dt < 0):

$$\omega _1 \rho _{2P}^ * \left( {E_{S2D2} } \right) <\varepsilon ^{ - 1} \gamma .$$

(21)

At E _S2DC, H ₅ = H ₆ = H ₇ = 0, and one obtains

$$\begin{array}{*{20}l}{\rho _{2{\text{P}}}^ * \left( {E_{{\text{S2DC}}} } \right) = \varepsilon ^{ - 1} \omega _1 ^{ - 1} \gamma ,} \hfill \\{B_{22}^ * \left( {E_{{\text{S2DC}}} } \right) = \frac{{\kappa _B B_{12}^ * \left( {E_{{\text{S}}2} } \right)}}{{\gamma \left( {1 - {{\omega _2 } \mathord{\left/{\vphantom {{\omega _2 } {\omega _1 }}} \right.\kern-\nulldelimiterspace} {\omega _1 }}} \right)}},} \hfill \\{B_{21}^ * \left( {E_{{\text{S2DC}}} } \right) = \omega _1^{ - 1} \left\{ {{{\kappa \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}2} } \right)} \mathord{\left/{\vphantom {{\kappa \rho _{1{\text{P}}}^ * \left( {E_{{\text{S}}2} } \right)} {\rho _{2{\text{P}}}^ * \left( {E_{{\text{S2DC}}} } \right)}}} \right.\kern-\nulldelimiterspace} {\rho _{2{\text{P}}}^ * \left( {E_{{\text{S2DC}}} } \right)}} - \beta - \omega _2 B_{22}^ * \left( {E_{{\text{S2DC}}} } \right)} \right\},} \hfill \\\end{array} $$

and a = 0, d < 0.

Numerical evaluation of the expression (18) with

$$\alpha _0 = 1,\quad \alpha _1 = - d + i + j,\quad \alpha _2 = - di - dj + ij + eh + bg,\quad \alpha _3 = - dij + ehj + bgj - bdg,\quad \alpha _4 = - bdgj,$$

suggests that E _S2DC is locally stable if ecotype 1 can invade ecotype 2 in the deep layer (Eq. 21 does not hold).

At E _SCDC,

$$B_{2j}^ * \left( {E_{{\text{SCDC}}} } \right) = \frac{{\kappa _B B_{1j}^ * \left( {E_{{\text{SC}}} } \right)}}{{\gamma - \varepsilon \omega _j \rho _{2{\text{P}}}^ * \left( {E_{{\text{SCDC}}} } \right)}}\;{\text{for}}\;j = 1\;{\text{and}}\;2,$$

where $\rho ^{ * }_{{2\operatorname{P} }} {\left( {E_{{{\text{SCDC}}}} } \right)}$ is the positive root of the following cubic equation:

$$\beta \varepsilon ^2 \omega _1 \omega _2 x^3 - \left[ {\varepsilon \kappa _B \omega _1 \omega _2 \left\{ {B_{11}^ * \left( {E_{{\text{SC}}} } \right) + B_{12}^ * \left( {E_{{\text{SC}}} } \right)} \right\} + \beta \gamma \varepsilon \left( {\omega _1 + \omega _2 } \right) + \varepsilon ^2 \kappa \omega _1 \omega _2 \rho _{1{\text{P}}}^ * \left( {E_{{\text{SC}}} } \right)} \right]x^2 + \left[ {\kappa _B \gamma \left\{ {\omega _1 B_{11}^ * \left( {E_{{\text{SC}}} } \right) + \omega _2 B_{12}^ * \left( {E_{{\text{SC}}} } \right)} \right\} + \beta \gamma ^2 + \gamma \varepsilon \kappa \left( {\omega _1 + \omega _2 } \right)\rho _{1{\text{P}}}^ * \left( {E_{{\text{SC}}} } \right)} \right]x - \kappa \gamma ^2 \rho _{1{\text{P}}}^ * \left( {E_{{\text{SC}}} } \right) = 0$$

Although the above equation may have three positive real roots, with our choice of parameters only one root satisfies $B_{2\operatorname{j} }^ * \left( {E_{{\text{SCDC}}} } \right) >0$.

Numerical evaluation of the expression (18) with

$$\alpha _0 = 1,\quad \alpha _1 = - a - d + i + j,\quad \alpha _2 = - di - dj + ij + eh + bg + ad - ai - aj,\quad \alpha _3 = - dij + ehj + bgj - bdg + adi + adj - aij - aeh,\quad \alpha _4 = adij - aehj - bdgj.$$

suggests that E _SCDC is always locally stable if E _SC is locally stable.

Appendix 3

Dependence of equilibrium values on the mixing rate

The Jacobian matrix of the full system (H ₁–H ₈; J _F) is as follows:

$$J_F = \left[ {\begin{array}{*{20}c}{{\mathbf{J}}_{{\text{SS}}} } & {\mathbf{0}} \\{{\mathbf{J}}_{{\text{DS}}} } & {{\mathbf{J}}_{{\text{DD}}} } \\\end{array} } \right],$$

(22)

where J _SS and J _DD are shown in “Appendix 1” and “Appendix 2,” and

$${\mathbf{J}}_{{\text{DS}}} = \left[ {\begin{array}{*{20}c}{\kappa _B } & 0 & 0 & 0 \\0 & {\kappa _B } & 0 & 0 \\0 & 0 & \kappa & 0 \\0 & 0 & 0 & 0 \\\end{array} } \right].$$

By simple calculations, the following inverse matrix of J _F ,

$${\mathbf{J}}_F^{ - 1} = \left[ {\begin{array}{*{20}c}{{\mathbf{J}}_{{\text{SS}}}^{ - 1} } & {\mathbf{0}} \\{ - {\mathbf{J}}_{{\text{DD}}}^{ - 1} {\mathbf{J}}_{{\text{DS}}} {\mathbf{J}}_{{\text{SS}}}^{ - 1} } & {{\mathbf{J}}_{{\text{DD}}}^{ - 1} } \\\end{array} } \right],$$

(23)

is obtained. Application of the implicit function theorem leads to Eq. 23. With small values of m, local stability is unchanged because the Jacobian matrix does not have eigenvalues with purely imaginary part when m = 0.

Appendix 4

How to use observed data and empirically derived model

We use time series of ¹⁴C primary production (mg C m⁻³ day⁻¹) at depth of 1 m in BATS station and those of primary production under light condition for 12 h (mg C m⁻³) at depth of 5 m in HOT station. Both datasets are available from the website of the US Joint Ocean Flux Study: http://www1.whoi.edu/research/bats_hot.html). For each time series (1989–2006), we conduct nonlinear regression (by the function “nls” in the software R) using the following expression for ten different n values (n = 1, 2,...10)

$${\text{PP}} = \operatorname{pp} _0 + \operatorname{amp} \cdot \sin ^{2n} \left( {\pi \cdot \left( {{\text{DY}} + \operatorname{phase} } \right)} \right),$$

where PP is primary production; DY is decimal year (0 ≤ DY ≤ 1).

For HOT and BATS datasets, respectively, the n value is chosen by Akaike Information Criterion scores by assuming Gaussian error to calculate the log-likelihood. The best model results in (see also the plot in Fig. 7):

$$\begin{array}{*{20}l}{{\text{PP}} = 7.58873 - 1.67250\sin ^2 \left( {\pi \cdot \left( {{\text{DY}} - 0.50762} \right)} \right),} \hfill \\{{\text{PP}} = 4.57232 + 6.81340\sin ^{12} \left( {\pi \cdot \left( {{\text{DY}} + 0.31582} \right)} \right).} \hfill \\\end{array} $$

Fig. 7

Seasonal fluctuations of primary production in HOT and BATS, and approximation by periodic functions

Based on this estimation, we especially use estimated factors n = 2 and n = 12 for representative of temporal changes in POC supply with moderate and strong seasonality, respectively, for Fig. 6. We set the two following functions as representative of seasonal changes of primary production (PP(t)) at the North Pacific Ocean and North Atlantic Ocean, respectively:

$$\begin{array}{*{20}l}{{\text{PP}}\left( t \right) = {\text{PP}}_{\min } + {\text{PP}}_{\max } \sin ^2 \left( {\frac{{\pi t}}{{365}}} \right),} \hfill \\{{\text{PP}}\left( t \right) = {\text{PP}}_{\min } + {\text{PP}}_{\max } \sin ^{12} \left( {\frac{{\pi t}}{{365}}} \right).} \hfill \\\end{array} $$

We also use a well-established model that allows one to estimate the production of POC and DOC from total primary production and temperature (Dunne et al. 2005). When we ignore the effect of temperature, the fraction of primary production by large phytoplankton (ϕ _L) is estimated to be (modified from Dunne et al 2005):

$$\phi _L \left( {{\text{PP}}} \right) = \frac{{ - 1 + \sqrt {1 + {{8{\text{PP}}} \mathord{\left/{\vphantom {{8{\text{PP}}} {{\text{PP}}_0 }}} \right.\kern-\nulldelimiterspace} {{\text{PP}}_0 }}} }}{{1 + \sqrt {1 + {{8{\text{PP}}} \mathord{\left/{\vphantom {{8{\text{PP}}} {{\text{PP}}_0 }}} \right.\kern-\nulldelimiterspace} {{\text{PP}}_0 }}} }},$$

where PP₀ (= 0.74 [mmol C m⁻³ day⁻¹]) corresponds to the total primary production at which production by small phytoplankton is equal to that by large phytoplankton (Dunne et al. 2005). The above equation represents how the contribution of large phytoplankton increases with the total primary production. It is possible to relate s _P and s _D to the total primary production via

$$\begin{array}{*{20}l}{s_{\text{P}} \left( {{\text{PP}}} \right) = \phi _{{\text{POC}}} \left( {{\text{PP}}} \right){\text{PP}},} \hfill \\{s_{\text{D}} \left( {{\text{PP}}} \right) = \left( {1 - \phi _{{\text{POC}}} \left( {{\text{PP}}} \right)} \right){\text{PP}}.} \hfill \\\end{array} $$

where $\phi _{{\text{POC}}} \left( {{\text{PP}}} \right) = \phi _{{\text{POCS0}}} \left( {1 - \phi _{\text{L}} \left( {{\text{PP}}} \right)} \right) + \phi _{{\text{POCL}}0} \phi _{\text{L}} \left( {{\text{PP}}} \right)$, linking in this way the fraction of particle production to ϕ _L (PP) and to the fraction of particle production by small and large phytoplankton (respectively, ϕ _POCS0 = 0.14 and ϕ _POCL0 = 0.74: Dunne et al. 2005). We set (PP_min, PP_max) = (1.0, 10.0) and thus (s _Pmin, s _Dmin, s _Pmax, s _Dmax) = (0.470, 0.530, 6.352, 3.648) where E _S2DC and E _SCDC are, respectively, locally stable with m = 0 (see dashed line in Fig. 2b).