Modeling polar marine ecosystem functions guided by bacterial physiological and taxonomic traits

Heterotrophic marine bacteria utilize organic carbon for growth and biomass synthesis. Thus, their physiological variability is key to the balance between the production and consumption of organic matter and ultimately particle export in the ocean. Here we investigate a potential link between bacterial traits and ecosystem functions in the rapidly warming West Antarctic Peninsula (WAP) region based on a bacteria-oriented ecosystem model. Using a data assimilation scheme, we utilize the observations of bacterial groups with different physiological traits to constrain the group-specific bacterial ecosystem functions in the model. We then examine the association of the modeled bacterial and other key ecosystem functions with eight recurrent modes representative of different bacterial taxonomic traits. Both taxonomic and physiological traits reflect the variability in bacterial carbon demand, net primary production, and particle sinking flux. Numerical experiments under perturbed climate conditions demonstrate a potential shift from low nucleic acid bacteria to high nucleic acid bacteria-dominated communities in the coastal WAP. Our study suggests that bacterial diversity via different taxonomic and physiological traits can guide the modeling of the polar marine ecosystem functions under climate change.

1 Introduction

Microbes regulate many key ecosystem functions in the marine food web. Unicellular primary producers fix organic carbon (i.e., an ecosystem function termed primary production), while heterotrophic marine bacteria and archaea (hereafter bacteria) utilize the fixed organic carbon for growth and biomass synthesis (i.e., an ecosystem function termed bacterial production, or BP; Azam et al., 1983). Thus, the variability in the abundance and activity of bacteria is central to understanding the balance between production and consumption of organic matter and ultimately particle export in the ocean. In flow cytometric analyses, bacteria cluster into two groups of cells with different nucleic acid content, including high nucleic acid (HNA) and low nucleic acid (LNA) cells (Bouvier et al., 2007; Gasol et al., 1999). These two groups are suggested to represent lineages (Schattenhofer et al., 2011; Vila-Costa et al., 2012) or physiological states (Bowman et al., 2017), in which HNA cells are generally larger in both cell and genome size compared to LNA cells (Bouvier et al., 2007; Calvo-Díaz and Morán, 2006). The significance of HNA versus LNA cells in determining distinct ecosystem states and functions has been investigated, but much is still unknown. In a recent study along the West Antarctic Peninsula (WAP), the high dimensionality of the bacterial community structure data was reduced via emergent self-organizing maps and subdivided into a small number of bacterial modes associated with specific taxonomic and functional traits (Bowman et al., 2017). Bowman et al. (2017) demonstrated that a combination of taxonomy, physiological structure (i.e., HNA and LNA cells), and abundance of bacterial communities explained up to 73 % of the variance in bulk BP. Their findings imply that physiological and taxonomic traits of bacteria may inform a predictive ecosystem model to further explore ecologically important questions including the following: would these bacterial traits reflect other important ecosystem functions such as the net primary production and particle sinking flux? If so, what would be the potential mechanisms and how will the relationship between bacterial traits and ecosystem functions be impacted by climate change?

The WAP is a rapidly warming marine ecosystem, with resulting changes in physical, ecological, and biogeochemical processes (Clarke et al., 2009; Cook et al., 2005; Ducklow et al., 2007; King, 1994; Meredith and King, 2005; Stammerjohn et al., 2008; Vaughan et al., 2003; Vaughan, 2006; Whitehouse et al., 2008). Routine monitoring through the Palmer Long-Term Ecological Research project (Palmer LTER; since 1991) has revealed climate-driven variations in seasonal phytoplankton accumulation (Saba et al., 2014; Schofield et al., 2017), bacterial dynamics (Bowman and Ducklow, 2015; Ducklow et al., 2012a; Kim and Ducklow, 2016; Luria et al., 2017, 2014), nutrient drawdown (Kim et al., 2016), and micro- and macrozooplankton dynamics (Garzio and Steinberg, 2013; Steinberg et al., 2015; Thibodeau et al., 2019). The wealth of Palmer LTER observations has enabled the construction of a numerical marine ecosystem model for the coastal WAP region (i.e., the WAP-1D-VAR v1.0 model; Kim et al., 2021) by adapting the regional test-bed models of other ocean basins (Friedrichs, 2001; Friedrichs et al., 2006, 2007; Luo et al., 2010, 2012). The WAP-1D-VAR v1.0 model is compared against roughly bi-weekly time-series observations over the growth season (October–March) near Palmer Station (64.77^∘ S, 64.05^∘ W; the mean depth ∼ 65 m, the bottom depth ∼ 75 m) that record seasonal variations in ecological processes modulated by variations in surface light, mixed layer depth, and surface sea-ice cover. The WAP-1D-VAR v1.0 model utilizes a data assimilation scheme to minimize the misfits between model results and observational data via a variational adjoint method (Lawson et al., 1995) that assimilates the available Palmer LTER data to objectively adjust the model parameters. Serving as a mechanistic model, assimilation of the Palmer LTER observations enables the model to constrain poorly measured bacterial processes (e.g., respiration, viral and grazing mortality, growth efficiency, carbon demand, and utilization of dissolved organic matter with varying lability) and to predict microbial system states in changing environments. Yet, incorporating molecular observations into a numerical ecosystem model is a challenge because of the difference in how levels of biological organization are treated in observations and models (Hellweger, 2020), as well as the high dimensionality of microbial molecular observations. One argument is that molecular-level changes in microbial dynamics may not directly translate into a clear picture of changes in community structure or resulting changes in bulk ecosystem functions.

In this study, we explore a potential link between bacterial traits and ecosystem functions in the warming coastal WAP using a bacteria-oriented ecosystem model modified from the WAP-1D-VAR v1.0 model (Kim et al., 2021). The bacterial traits examined in this study include physiological and taxonomic traits. For physiological traits, our model explicitly simulates the time-evolving dynamics of two ubiquitous bacterial groups with differing nucleic acid contents, the HNA group and the LNA group, by directly assimilating the group-specific carbon biomass observations estimated by flow cytometry. For taxonomic traits, taxonomic modes derived from bacterial 16S rRNA gene sequence data (Bowman et al., 2017) are compared to final model results at the corresponding time points with the assumption that bacterial taxonomy would provide information about bacterial ecosystem processes and structures. Our study indirectly incorporates bacterial molecular information into ecosystem-level dynamics, in contrast to genome-scale or gene-centric models predicting the time-evolving dynamics of microbial molecular processes.

Figure 1 Model structure. The model is forced by five different physical forcings, denoted as a horizontal row across the top of the schematic. As the ecosystem component, heterotrophic bacteria are divided into two groups of differing physiological states, including high nucleic acid (HNA) and low nucleic acid (LNA) bacterial compartments. The flows between the prognostic state variables with the name of the numbered flows in the legend only represent these two bacterial compartments.

2 Material and methods

2.1 Bacteria-oriented ecosystem model

The model is originally derived and modified from the one-dimensional (1-D) variational data assimilation planktonic ecosystem model for the coastal WAP region called the WAP-1D-VAR v1.0 model (Kim et al., 2021). The WAP-1D-VAR v1.0 model tracks C, N, and P of its state variables with flexible stoichiometry. For our study, we modify the original model's single bacterial compartment into HNA and LNA bacterial compartments and only discuss their C stocks and rates, as well as other state variables that remain the same as the WAP-1D-VAR v1.0 model. The state variables in this modified model include HNA and LNA bacteria, diatoms, cryptophytes, microzooplankton, krill, labile dissolved organic carbon (LDOC), semi-labile DOC (SDOC), ammonium (NH₄), nitrate (NO₃), phosphate (PO₄), and particulate (carbon) detritus (Fig. 1). Refractory DOC (RDOC) and higher trophic levels are implicit to serve as model closure terms (i.e., they are source or sink terms of other explicit state variables, and their time derivatives are not calculated in the model). The model is forced by mixed layer depth (MLD), photosynthetically active radiation (PAR) at the ocean surface, surface sea-ice concentration, water-column temperature, and eddy diffusivity (Fig. S1 in the Supplement) using a constant time step of 1 h and a second-order Runge–Kutta scheme (Texts S1 and S2 in the Supplement). The model allows both LDOC and SDOC as substrate sources for bacteria, and it is the nutrient quota of bacteria that allows the lability of SDOC to vary. In contrast to LDOC pool that is entirely available for bacteria, a parameter controlling the lability of SDOC (r_SDOC, Tables S2–S6 in the Supplement) regulates the size of the portion of SDOC that can be made available for bacterial uptake. Bacterial C growth is determined by their cellular nutrient quota, as well as available LDOC and SDOC concentrations (Kim et al., 2021).

Table 1 Initial values of bacterial model parameters. Different initial values are assigned to the model parameters of high nucleic acid (HNA) and low nucleic acid (LNA) bacterial groups to simulate their distinct physiological processes and trophic interactions.

The time derivative of C biomass for each bacterial group is determined as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{C}_{\text{HNA}}}{\mathrm{d}t}}=& G_{\text{HNA},\text{LDOC}}^{\mathrm{C}}+G_{\text{HNA},\text{SDOC}}^{\mathrm{C}}-R_{\text{HNA}}^{\mathrm{C}}\\ & -E_{\text{HNA},\text{RDOC}}^{\mathrm{C}}-E_{\text{HNA},\text{SDOC}}^{\mathrm{C}}\\ & -G{Z}_{\text{HNA}}^{\mathrm{C}}-M_{\text{HNA}}^{\mathrm{C}},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{C}_{\text{LNA}}}{\mathrm{d}t}}=& G_{\text{LNA},\text{LDOC}}^{\mathrm{C}}+G_{\text{LNA},\text{SDOC}}^{\mathrm{C}}-R_{\text{LNA}}^{\mathrm{C}}\\ & -E_{\text{LNA},\text{RDOC}}^{\mathrm{C}}-E_{\text{LNA},\text{SDOC}}^{\mathrm{C}}\\ & -G{Z}_{\text{LNA}}^{\mathrm{C}}-M_{\text{LNA}}^{\mathrm{C}},\end{array}\end{array}$$

where C is bacterial C biomass, $G_{\text{LDOC}}^{\mathrm{C}}$ is LDOC consumption, $G_{\text{SDOC}}^{\mathrm{C}}$ is SDOC consumption, R^C is respiration, $E_{\text{RDOC}}^{\mathrm{C}}$ is RDOC excretion, $E_{\text{SDOC}}^{\mathrm{C}}$ is SDOC excretion, GZ^C is the amount of C biomass grazed by microzooplankton, and M^C is mortality caused by viral attack (unit: mmol C m⁻³ for C and $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ for the rest of the terms). The sum of the first three terms on the right-hand side in Eqs. (1) and (2) is defined as BP for each bacterial group (i.e., ${\mathrm{BP}}_{\mathrm{HNA}}=G_{\text{HNA},\text{LDOC}}^{\mathrm{C}}+G_{\text{HNA},\text{SDOC}}^{\mathrm{C}}-R_{\text{HNA}}^{\mathrm{C}}$, ${\mathrm{BP}}_{\mathrm{LNA}}=G_{\text{LNA},\text{LDOC}}^{\mathrm{C}}+G_{\text{LNA},\text{SDOC}}^{\mathrm{C}}-R_{\text{LNA}}^{\mathrm{C}}$). Both C_HNA and C_HNA are constrained by the group-specific C biomass data estimated via flow cytometry, while only bulk BP (i.e., $\text{BP}={\mathrm{BP}}_{\mathrm{HNA}}+{\mathrm{BP}}_{\mathrm{LNA}}$), not the group-specific BP, is constrained by the observations due to the lack of the group-specific BP data at the study site. With the initial parameter values distinct to each bacterial group (Table 1), the model incorporates the observations of C_HNA, C_LNA, BP, and other state and rate variables to its built-in data assimilation scheme (Sect. 2.2) that optimizes the parameters and calculates the resulting C stock and flows of each bacterial group (Eqs. 1 and 2). In other words, the partitioning of BP into HNA and LNA groups is purely determined by parameter optimization using the information about other ecosystem observations.

Figure 2 Modeling framework. The modeling framework in this study consists of the mechanistic part as the processes associated with the bacteria-oriented ecosystem model (Fig. 1) and the data-driven part that represents how bacterial taxonomic modes are associated with the modeled bacterial and other key ecosystem functions. A variational adjoint method is employed for the parameter optimization and data assimilation processes (adapted from Glover et al., 2011). Gradient: the sensitivity of the total cost function with respect to the model parameter from optimization. Optimized model results are interpreted as a function of bacterial taxonomic and physiological traits.

2.2 Modeling framework

The modeling framework consists of the mechanistic and data-driven parts (Fig. 2). The mechanistic part represents prognostic, time-evolving microbial processes in the model (Fig. 1, Sect. 2.1) with its built-in data assimilation scheme via a variational adjoint method (Lawson et al., 1995). The data assimilation scheme minimizes the misfits between observations (i.e., assimilated data; Sect. 2.3) and model results by objectively optimizing a subset of model parameters (Friedrichs, 2001; Spitz et al., 2001; Ward et al., 2010; details in Text S3 in the Supplement). The data-driven part represents the pairing of the optimized model results from the mechanistic part and the bacterial taxonomic “modes” derived from 16S rRNA gene sequence abundance data (Bowman et al., 2017). We call this latter part data-driven because the high dimensionality of the bacterial community structure data is reduced to a single taxonomic mode using an unsupervised machine learning algorithm called Kohonen's self-organizing maps (Kohonen, 2001). Each taxonomic mode (mode hereafter) represents its specific taxonomic traits and is expressed as a single categorical variable without linear progression in between. For example, mode 1 is not necessarily closer to mode 3 than it is to mode 7. Modes are not necessarily correlated to the physiological traits of bacteria (i.e., HNA and LNA C biomass) despite being derived from the same samples. In other words, the taxonomic and physiological traits are independent of each other.

We select the nearshore Palmer LTER Station B (64.77^∘ S, 64.05^∘ W; ∼ 65 m) in the coastal WAP as the modeling site. The Station B datasets consist of roughly bi-weekly physical, chemical, and biological profiles collected via a profiling conductivity–temperature–density (CTD) rosette. Flow cytometric data for HNA and LNA C biomass and 16S rRNA gene amplicon data for taxonomic modes come from Arthur Harbour Station B at 10 m depth (situated 1 km from Station B) or Palmer Station seawater intake at 6 m depth (Bowman et al., 2017). Three upper-ocean depth levels, 0, 10, and 20 m (with the layer thickness of 2, 16, and 4 m, respectively), are modeled for four consecutive growth seasons, including November 2010–March 2011 (2010–11 hereafter), 2011–12, 2012–13, and 2013–14. However, only the results from 10 m are presented in detail because of the availability of the bacterial trait data at that depth. Despite the advantage of simulating the full water-column layers, it would be best to exclude the depth levels without bacterial trait observations, yet to include an adequate number of depth levels to simulate seasonal MLD and light impacts on bacterial dynamics. Thus, we choose to model three layers, that is, 0, 10, and 20 m, in a 1-D (vertical) water column.

The 1-D modeling of the coastal WAP region is justifiable given that the WAP shows relatively weak net advection compared with the Antarctic Circumpolar Current (ACC) or the subpolar gyres (Meredith et al., 2008, 2013). In addition, the CTD observations at Palmer Station do not show abrupt changes in physical and biogeochemical tracers as a result of lateral advection, with fairly homogeneous temperature and salinity distributions for the years and depths modeled in our study (Kim and Ducklow, 2016). There is a 6-month sampling gap in the austral autumn and winter months, so we optimize the model each year separately only for the austral spring to summer months. This results in each year possessing its own uniquely optimized parameter set that drives the minimized model–observation misfits for the given year. We also optimize the model for the climatological year, referred to as the climatological model, constructed by averaging 4-year observations (2010–11 to 2013–14; Text S4 in the Supplement). We do not average the whole Palmer LTER multi-decadal period (since 1991) because of the lack of HNA and LNA C biomass data except for those 4 years. Other modeling aspects (e.g., model initialization, spin-up, and bottom boundary conditions) are detailed in Text S4 and Kim et al. (2021).

2.3 Assimilated data

We assimilate the Palmer LTER observations from 0, 10, and 20 m that correspond to the compartments and flows in the model, including NO₃, PO₄, phytoplankton taxonomic-specific chlorophyll (Chl) for diatoms and cryptophytes (Schofield et al., 2017), microzooplankton C biomass (Garzio et al., 2013), bulk primary production (PP), bulk BP, HNA bacterial C biomass, LNA bacterial C biomass, SDOC, particulate organic carbon (POC), and particulate organic nitrogen (PON). Though not available in 2011–12, because of the importance in constraining the group-specific phytoplankton dynamics, the 4-year climatological value of the group-specific Chl is assimilated for 2011–12. NO₃ is not assimilated in 2010–11, while POC, PON, and SDOC are not assimilated in 2012–13 and 2013–14 because of the lack of observations in those years. Krill C biomass is not assimilated due to the strong patchiness of their distribution with many zero values that may hinder proper model optimization, while microzooplankton C biomass (2010–11) from a single year's measurement is assimilated for all 4 model years to at least provide constraints on the parameter values for phytoplankton grazing. The model–observation misfits for microzooplankton are not examined because of the discrepancy in the timing and location of those data assimilated compared to our study.

SDOC is calculated by subtracting the background (RDOC) concentration (40.0 mmol C m⁻³) from climatological total DOC concentration. POC (PON) is assimilated to represent the model detrital pool, but its measurements contain living biomass from bottle filter experiments. Climatological observations show that living phytoplankton and bacterial biomass account for 26 % of total POC and 29 % of total PON, so these fractions are used to exclude living biomass from the bulk particulate material pool. When converting Chl to phytoplankton C biomass, the maximum Chl / N ratio is used along with the reference (Redfield) $\mathrm{C}/\mathrm{N}$ ratio of 0.15. BP ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) is derived from ³H-leucine incorporation rate ($\mathrm{pmol}{\mathrm{l}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) data using the conversion factor of 1.5 kg C mol⁻¹ leucine incorporated (Ducklow, 2000). The bacterial group-specific C biomass (mmol C m⁻³) is estimated from bacterial abundance measured by flow cytometry (i.e., bulk bacterial C biomass multiplied by the fraction of each physiological group, fHNA or fLNA, with the conversion factor of 10 fg C cell⁻¹; Fukuda et al., 1998).

2.4 Cost function and portability index

The total cost function is calculated to represent the misfit between observations and model results as follows:

$$\begin{array}{}\text{(3)}& J={\sum }_{m=\mathrm{1}}^M{\displaystyle \frac{\mathrm{1}}{N_m}}{\sum }_{n=\mathrm{1}}^{N_m}{\left({\displaystyle \frac{{\widehat{a}}_{m,n}-a_{m,n}}{{\mathit{\sigma }}_m}}\right)}^{\mathrm{2}},\end{array}$$

where m and n represent assimilated data types and data points, respectively, M and N_m are the total number of assimilated data types and data points for data type m, respectively, σ_m is the target error for data type m, a_m,n is observations, and ${\widehat{a}}_{m,n}$ is model output. Hereafter, we present the total cost function as the total cost function normalized by M ($J^{\prime }=J/M)$ and normalized costs of individual data types ($J_m^{\prime }=J_m/M$) as the model–observation misfit equivalent to a reduced Chi-square estimate of the model goodness of fit (i.e., $J^{\prime }=\mathrm{1}$ as a good fit from optimization, J^′ ≫ 1 as a poor fit due to underestimation of the error variance or the fit not fully capturing the data, and J^′ ≪ 1 as an overfitting of the data, fitting the noise, or overestimation of the error variance). The base-10 logarithm of Chl and PP is used in Eq. (3) to account for high productivity of the WAP waters and the approximate log-normal distribution of those data types (Campbell, 1995; Glover et al., 2018). The target error σ_m is calculated for each data type m as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_m=\stackrel{\mathrm{‾}}{a_{m,n}}\cdot {\text{CV}}_m,\end{array}$$

where $\stackrel{\mathrm{‾}}{a_{m,n}}$ is the climatological mean of the observations, and CV_m is the adjusted coefficient of variation (CV) of the observations of each data type over 0, 10, and 20 m (due to observational error and seasonal and interannual variations). CV_m values for the 4 modeled years in our study are higher than those across every measured depth within the mixed layer for an extended year period in the original WAP-1D-VAR v1.0 model (2002–03 to 2011–12; Kim et al., 2021) and are therefore reduced to the levels in the mixed layer to avoid an overestimated target error of each data type (Text S5 in the Supplement). The rationale behind using the adjusted CV in the target error calculation is based on Luo et al. (2010), in which all properties should be completely mixed in the mixed layer, a perfect measurement without significant errors should generate similar values at every measured depth within the mixed layer, and the average CV of all depth profiles can be used as CV in the target error calculation. The standard deviation is used as target errors of the log-converted data types. The CV of the log-converted data type is estimated as the average of ± 1 standard deviation in log space converted back into normal space (Doney et al., 2003; Glover et al., 2018).

We compute the portability index (Friedrichs et al., 2007) to evaluate the broader applicability of the optimized model parameter set for each year in predicting the dynamics of the other year as follows:

$$\begin{array}{}\text{(5)}& \text{Portability index}=J_{\mathrm{c}}^{\prime }/J_{\mathrm{x}}^{\prime },\end{array}$$

where $J_{\mathrm{x}}^{\prime }$ is the normalized cross-validation total cost function when a model parameter set optimized for a given year is used to simulate another year, and $J_{\mathrm{c}}^{\prime }$ is the normalized total cost function of the climatological model. A portability index value close to 1 indicates a more portable model or a system that is not particularly sensitive to year-to-year variations in optimized parameters, while an index value ≪ 1 indicates a less portable model or a system that is sensitive to year-to-year variations in optimized parameters.

2.5 Uncertainty analysis

The uncertainties of the optimized parameters are estimated using a finite difference approximation of the complete Hessian matrix during the iterative data assimilation process (i.e., the second derivatives of the cost function with respect to the model parameters). When computed at the minimum of the cost function value, the square root of a diagonal element in the inversed Hessian matrix represents the logarithm of the relative uncertainty of the corresponding optimized parameter. The absolute uncertainty of the optimized parameter is calculated as $p_f\cdot e^{\pm {\mathit{\sigma }}_f}$, where p_f is the value of the optimized parameter, and σ_f is its relative uncertainty. We denote an optimized parameter with σ_f larger than 50 % as an “optimized” parameter, while an optimized parameter with σ_f smaller than 50 % is denoted as a “constrained” parameter (Text S3, Tables S2–S6). We then conduct Monte Carlo experiments to examine the impact of the uncertainties of the constrained parameters on the modeled fields. The Monte Carlo experiments consist of (1) creating an ensemble of parameter sets (N=1000) by randomly sampling values within the uncertainty ranges of the constrained parameters and (2) then performing a model simulation using each parameter set. All uncertainty estimates are calculated following standard error propagation rules and presented herein as ± 1 standard deviation.

Table 2 Data types, observed means, coefficient of variation, target errors, and costs before and after optimization. The observed mean ($\stackrel{\mathrm{‾}}{a}$), coefficient of variation (CV), and target error (σ) of each assimilated data type used for calculating the normalized cost function (unitless; Eq. 3) before ($J_{\mathrm{0}}^{\prime }$) and after optimization ($J_{\mathrm{f}}^{\prime }$). Data type unit: mmol N m⁻³ for nitrate (NO₃); mmol P m⁻³ phosphate (PO₄); mg m⁻³ for diatom chlorophyll (Chl_DA) and cryptophyte chlorophyll (Chl_CR); mmol C m⁻³ for HNA and LNA bacterial biomass, SDOC, and POC; mmol N m⁻³ for PON; and $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ for primary production (PP) and bacterial production (BP). NO₃ was not assimilated in 2010–11, while SDOC, POC, and PON were not assimilated in 2012–13 and 2013–14 (denoted as ”–” in the table).

3 Results

3.1 Model skill assessment

The iterative optimization procedure reduced by 24 %–93 % the misfits between observations and model results for each year and for the climatological year compared to those obtained using the initial parameter values (Table 2). The optimized parameter sets satisfied the pre-set convergence criteria, including the local minima achieved by the total costs, low gradients of the total costs with respect to each optimized parameter, and positive eigenvalues of the Hessian matrix (details in Kim et al., 2021). The total costs were reduced by optimizing only a subset of the parameters: five–seven constrained and three–six optimized parameters (Tables S2–S6). The optimized parameters in common across all years were α_DA (the initial slope of photosynthesis versus irradiance curve of diatoms, $\mathrm{mol}\mathrm{C}\left(\mathrm{g}\mathrm{Chl}\mathit{\text{a}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}\right(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$), μ_HNA (the maximum HNA bacterial growth rate, d⁻¹), μ_LNA (the maximum LNA bacterial growth rate, d⁻¹), and g_CR (the half-saturation density of cryptophytes in microzooplankton grazing, mmol C m⁻³). g_HNA (the half-saturation density of HNA bacteria in microzooplankton grazing, mmol C m⁻³), g_MZ (the half-saturation density of microzooplankton in krill grazing, mmol C m⁻³), and μ_KR (the maximum krill growth rate, d⁻¹) were the next most frequently optimized, at least for 4 years out of a total of 5 modeled years including the climatological year.

Figure 3 Model skill assessment. A Taylor diagram using a polar-coordinate system summarizing the model–observational correspondence for each model stock and flow for individual annual simulations for the 4 modeled years together (2010–11 to 2013–14; a) and for the climatological year (b). The angular coordinate denotes the Pearson correlation coefficient (r), the distance from the origin denotes the standard deviation normalized by the standard deviation of the observations, and the distance from point (1,0), marked as REF on the x axis, describes the centered (bias removed) root-mean-square difference (RMSD) between model results and observations.

Because of this study's focus on the modeled bacterial and other ecosystem functions as a function of bacterial traits (Sect. 3.3) rather than of year (Figs. S2–S5 in the Supplement), we combined the observations and model results from all 4 years together for model skill assessment. According to the Taylor diagrams, model skills were overall similar among the 4 study years (Fig. 3a) and the climatological year (Fig. 3b). Three core variables in this study, including HNA biomass, LNA biomass, and BP, had better model–observation agreements than other data types, with relatively high correlations, a low centered (bias removed) root-mean-square difference (RMSD), and the normalized standard deviation closer to 1. These variables also had better fits to the 4-year seasonal cycles of the observations than other data types (Fig. S7 in the Supplement). However, the model skill for HNA biomass slightly degraded in the climatological model (Fig. 3b), with the insignificant correlation (p=0.61 versus r=0.53 and p=0.003 in Fig. 3a), lower normalized standard deviation, and higher RMSD than the 4 years together (Fig. 3a). After optimization, the models tended to underestimate PP with relatively larger errors than for other data types (Fig. S7), while its temporal and spatial (depth) variability was captured well as shown by high correlations (Fig. 3). By contrast, there were slight positive model biases for POC and PON (Fig. S7), and their variability was not well captured as shown by their negative correlations (Fig. 3).

Table 3 Cross-validation cost and portability index. $J_{\mathrm{c}}^{\prime }$ is the normalized optimized cost from the climatological model (equivalent to $J_{\mathrm{f}}^{\prime }$ under the climatological model parameter set in Table 2), and $J_{\mathrm{x}}^{\prime }$ is the normalized cross-validation cost (Eq. 5), in which, for example, $J_{\mathrm{x},\text{2011–12}}^{\prime }$ under the row “2010–11 model parameter set” in Table 3 indicates the normalized cross-validation cost from simulating the 2010–11 model parameter set against 2011–12.

NA: not available.

Cross-validation cost analyses showed the increased model–observation misfits when a set of parameters optimized for a given year was applied to simulate another year's dynamics (Tables 2 and 3), suggesting that each year was best modeled using its own unique set of the optimized parameters. The magnitude of increase in the cost function varied by year pair, with the average portability index values indicating that the optimized model parameters for 2012–13 was most portable (0.76 ± 0.11), followed by those for 2013–14 (0.73 ± 0.17), 2010–11 (0.68 ± 0.08), and 2011–12 (0.61 ± 0.12; Table 3), though the differences were not always significantly different among the years.

Figure 4 Seasonal progression of modeled HNA and LNA bacterial carbon stocks and rates and key ecosystem functions across years. Seasonal patterns of HNA and LNA bacterial carbon stocks and flows, NPP and POC sinking flux at 10 m depth over the growth season (November–March) for each of the 4 simulation years (a), and coefficient of variation (the Monte-Carlo-derived standard deviation divided by each data point from Fig. 4a) from 1000 Monte Carlo experiments (b).

Figure 5 Annual mean carbon stocks and flows. Carbon stocks (mmol C m⁻³) and flows ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), particle sinking flux ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and other stocks (e.g., nutrients, mmol m⁻³) averaged over the growth season in each year are denoted as the numbers on the first row, while the numbers in the second row or in the parentheses are the standard deviation propagated from averaging over the growth season and the Monte-Carlo-derived uncertainties. Numbers around the arrows represent intercompartmental flows and do not necessarily balance to zero due to the build-up or loss in a compartment over the growth season. The magnitude of the N and P flows, as well as the flows smaller than 0.01 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$, are omitted. RDOC and higher levels are implicit.

3.2 Bacterial carbon stocks and flows

C stocks and flows for each bacterial group represented significant seasonal and interannual variability (Figs. 4a and S8 in the Supplement). Across years, HNA bacteria had significantly higher seasonal maximum values than their LNA counterparts when normalized by the group-specific biomass. These so-called cell-specific, seasonal maximum rates of the HNA group ranged from 0.10 ± 0.004 to 0.59 ± 0.24 d⁻¹, 0.03 ± 0.001 to 0.18 ± 0.12 d⁻¹, 0.07 ± 0.003 to 0.18 ± 0.08 d⁻¹, 0.05 ± 0.002 to 0.57 ± 0.26 d⁻¹, and 0.07 ± 0.03 to 0.36 ± 0.17 d⁻¹ for LDOC uptake, SDOC uptake, respiration, BP, and grazing rates, respectively (Fig. 4). For the LNA group, the maximum cell-specific rates ranged from 0.01 ± 0.002 to 0.12 ± 0.02 d⁻¹, 0.004 ± 0.002 to 0.03 ± 0.01 d⁻¹, 0.01 ± 0.001 to 0.02 ± 0.002 d⁻¹, 0.01 ± 0.003 to 0.13 ± 0.02 d⁻¹, and 0.02 ± 0.0004 to 0.17 ± 0.03 d⁻¹ for LDOC uptake, SDOC uptake, respiration, BP, and grazing rates, respectively (Fig. 4). For each year, C stocks and flows averaged over the growth season (Fig. 5) and those normalized by NPP (net primary production; normalized by NPP in 1 d for C stocks; Fig. S9 in the Supplement) summarized an annual snapshot of the group-specific bacterial dynamics. The annual mean LNA biomass was ∼ 17 times larger than that of HNA biomass in 2011–12 (Fig. 5b), in contrast to relatively similar average biomass values of both groups in other years (Fig. 5a, c, and d). Bacterial carbon demand (BCD; i.e., BCD = BP + bacterial respiration; blue arrows in Fig. 5) was mostly supported by LDOC (67 %–81 %) for both bacterial groups.

The rest of the modeled C stocks and flows fell into one of the following categories: (1) the variable for which a single year's values were assimilated (i.e., microzooplankton C biomass), (2) the variables for which observational values for the given year were assimilated (i.e., nutrients, POC or detritus, and SDOC), and (3) the variables that were not assimilated at all (i.e., krill C biomass, LDOC, NH₄, and particle sinking flux). Compared to bacterial variables, there was little interannual variability in the average microzooplankton C biomass (Fig. 5). Even in the years when NO₃, POC, and SDOC were not assimilated, their values were modeled similarly to those modeled in other assimilated years (Fig. 5). Modeled LDOC and NH₄ were also within the reasonable ranges of their typically small values (< 1 µM).

Figure 6 Bacterial physiological and taxonomic association with modeled ecosystem functions. The results from linear regression of the key modeled ecosystem functions on a categorical predictor of the observed mode (a–c) and on the observed fraction of HNA cells (d–f). Regression statistics: (a) number of observations (N)=43, error degrees of freedom (df)=35, root-mean-square error (RMSE) = 0.68, r² = 0.39, adjusted r² = 0.27, F statistic value = 3.22, p value = 0.01; (b) N=43, df = 35, RMSE = 2.88, r² = 0.48, adjusted r² = 0.37, F statistic value = 4.55, p value = 0.001; (c) N=43, df = 35, RMSE = 0.03, r² = 0.81, adjusted r² = 0.77, F statistic value = 20.7, p value < 0.001; (d) N=43, df = 41, RMSE = 0.65, r² = 0.36, adjusted r² = 0.34, F statistic value = 22.8, p value < 0.001; (e) N=43, df = 41, RMSE = 2.57, r² = 0.51, adjusted r² = 0.50, F statistic value = 43.0, p value < 0.001; (f) N=43, df = 41, RMSE = 0.04, r² = 0.57, adjusted r² = 0.56, F statistic value = 53.5, p value < 0.001.

3.3 Bacterial physiological and taxonomic association with ecosystem functions

Each mode was dominated by unique bacterial taxa, thereby representing taxonomic traits (Fig. S10 in the Supplement). Candidatus Pelagibacter was the most abundant in mode 6 (Fig. S10c), Dokdonia sp. MED134 in mode 7 (Fig. S10d), Candidatus Thioglobus singularis PS1 in mode 1 (Fig. S10e), Owenweeksia hongkongensis DSM 17368 in mode 2 (Fig. S10f), Rhodobacteraceae in mode 5 (Fig. S10g), and Planktomarina temperata RCA23 in mode 4 (Fig. S10h). To explore a potential link between the bacterial taxonomic traits and the key ecosystem functions, we first extracted the modeled NPP, POC sinking flux, and BCD from the ecosystem model (i.e., the “final optimized output” in Fig. 2) at the time the bacterial samples and depth (10 m) were placed into a single mode derived from the observations. We then performed a linear regression with the mode as a factor, in which the mode is a categorical predictor with eight modes rather than an ordinal or continuous variable (i.e., equivalent to a one-way ANOVA with eight different categories). 27 %, 37 %, and 77 % of the total variance in the modeled NPP, POC sinking flux, and BCD were explained by the bacterial taxonomic mode (Fig. 6a–c). In particular, modes 3, 5, and 7 were associated with 2–3 times higher NPP, POC sinking flux, and BCD compared to when mode 4 dominated (two-sample t test with unequal sample size, p=0.02 for NPP and p < 0.001 for POC sinking flux and BCD) or to when mode 6 dominated (p=0.03 for NPP, p=0.003 for POC sinking flux, and p < 0.001 for BCD).

The observed mode was also positively correlated to the observed fHNA (r² = 0.52, p < 0.001; not shown). Thus, we examined a potential link between the bacterial physiological traits and the key ecosystem functions as described above, using a linear regression with the observed fHNA as a predictor and the modeled ecosystem functions as dependent variables. The observed fHNA was positively correlated to the modeled NPP (r² = 0.34, p < 0.001; Fig. 6d) and to a stronger extent to the modeled POC sinking flux (r² = 0.50, p < 0.001; Fig. 6e) and to the modeled BCD (r² = 0.56, p < 0.001; Fig. 6f). The stepwise addition of one predictor variable to the other predictor variable (i.e., fHNA adding to mode or vice versa) did not improve the model performance (not shown). These results suggest a clear link between the modeled ecosystem functions and the bacterial taxonomic (modes) and physiological (fHNA) trait observations.

Figure 7 Climate change experiments. Seasonal progression of the modeled HNA and LNA bacterial carbon stocks and rates and key ecosystem functions under observed physical forcing and climate change experiments (a) and the percent change in the corresponding variable compared to observed fields in the second and third row of each panel, with the first row of each panel as zero to represent base states (b). For example, the percent anomaly of HNA biomass in (b) = (HNA biomass under +1 ^∘C and −10 % − HNA biomass under observed forcing) ⋅ 100$/$HNA biomass under observed forcing.

3.4 Climate change experiments

We explored the responses of the modeled bacterial dynamics and other ecosystem functions (Sects. 3.2 and 3.3) to changing climates along the WAP (Fig. 7). Due to the varying portability of the optimized parameter sets among the 4 study years, we used the optimized parameter set for the climatological year (Table S6) to simulate an overall WAP system response under perturbed ocean temperatures (i.e., +0.5 ^∘C and +1.0 ^∘C relative to observed temperatures) and sea-ice forcing fields (i.e., 5 % and 10 % loss of sea-ice concentrations relative to observed sea-ice concentrations). These experiments were conducted under each perturbed condition separately (i.e., warming alone in Fig. S11 in the Supplement versus melting alone in Fig. S12 in the Supplement), as well as simultaneously (i.e., climate change; Fig. 7). We only analyzed the results from the climate change experiments, given that despite different impacts of each forcing change (i.e., the impact of warming on rate processes versus the impact of melting on light and photosynthesis but not MLD in our model), climate change would cause simultaneous changes in sea ice and water temperature along the WAP.

The climate change experiments resulted in a combination of changes in overall bacterial C stocks and rates, as well as the key ecosystem functions, and shifts in their seasonal timing (Fig. 7a) compared to the base state (the first row as the base state in Fig. 7a and b, while the second and third rows as anomalies under perturbed conditions in Fig. 7b). HNA bacterial C stock and rates responded more strongly to the perturbed climate conditions compared to their LNA counterparts. Under combined warming and melting (+1.0 ^∘C and −10 %) conditions, there were the maximum increases in the HNA C stock and rates by 19 %–35 % (29 ± 89 % for biomass, 22 ± 67 % for LDOC uptake, 35 ± 111 % for SDOC uptake, 26 ± 79 % for respiration, 25 ± 78 % for BP, 29 ± 89 % for viral mortality, 19 ± 26 % for grazing, and 29 ± 89 % for RDOC excretion) compared to the maximum increases in the LNA C stock and rates by 3 %–15 % (3 ± 2 % for biomass, 6 ± 11 % for LDOC uptake, 15 ± 27 % for SDOC uptake, 8 ± 3 % for respiration, 7 ± 6 % for BP, 3 ± 2 % for viral mortality, 7 ± 18 % for grazing, and 3 ± 2 % for RDOC excretion). In contrast to bacterial C stocks and rates that increased consistently throughout the growth season, microzooplankton grazing rates showed seasonally mixed responses for both HNA and LNA cases, with the maximum decreases of 8 ± 32 % for HNA bacteria and of 4 ± 32 % for LNA bacteria. Similarly, there were the maximum increases in NPP and POC sinking flux by 14 ± 15 % and 3 ± 22 % and the maximum decreases in NPP and POC sinking flux by 4 ± 11 % and 3 ± 13 %, respectively. SDOC exhibited the maximum increase by 2 ± 1 % early in the season but became depleted strongly as the season progressed. LDOC decreased consistently in response to the perturbed conditions, with the maximum decrease by 10 ± 43 %.

4 Discussion

4.1 Model skill assessment

Despite the important role that heterotrophic marine bacteria play in the ocean carbon cycle, the vast majority of mechanistic biogeochemical models neither include them as a model state variable nor explicitly simulate their physiological processes. Most models parameterize the bacterial remineralization of sinking organic matter with depth by fitting the power law functions or other similarly derived empirical approaches (Buesseler et al., 2020; Cael and Bisson, 2018). Cellular functions, taxa, and functional gene expression of other prokaryotes, such as cyanobacteria (Hellweger, 2010; Martín-Figueroa et al., 2000; Miller et al., 2013), or a diverse suite of microbial functional groups (Coles et al., 2017; Dutkiewicz et al., 2020) have been modeled so far. However, our study is the first to explicitly model heterotrophic bacterial groups of different physiological traits and to link their relationship with the key ecosystem functions.

Only a subset of the parameters was optimized in our model to simulate microbial and ecological patterns for each year, consistent with other data assimilation modeling studies (Friedrichs, 2001; Friedrichs et al., 2006, 2007; Luo et al., 2010, 2012). In general, optimization of this class of marine ecosystem model requires adjustment of a small number of independent parameters to achieve well-posed model solutions because of the highly cross-correlated nature of the parameters in the inherently nonlinear model equations (Fennel et al., 2001; Harmon and Challenor, 1997; Matear, 1996; Prunet et al., 1996a, b). In our study, most of the constrained parameters were directly associated with the bacterial processes, and there were overall better model–observation fits for the bacterial data types compared to other data types. These results provide confidence in the simulated bacterial C stocks and rates.

Optimization also sheds light on major unknown parameters in bacterial grazing processes involving g_HNA and g_LNA (the half-saturation densities of HNA and LNA bacteria in microzooplankton grazing, respectively). Microzooplankton grazing of the given bacterial group is simulated using a Holling type 2 density-dependent grazing function with a preferential prey selection on diatoms, cryptophytes, and the other bacterial group, in which a single microzooplankton maximum grazing rate is implemented for both bacterial groups for model simplicity purposes (Tables S2–S6 in the Supplement; Kim et al., 2021). Thus, it is the half-saturation density that determines the degree of preferential grazing by microzooplankton on the given bacterial group, the change in which may ultimately depend on the group-specific C biomass. Due to the lack of a priori knowledge on the relative magnitude of g_HNA and g_LNA, we assigned an identical initial parameter value (Table 1) to let the data assimilation scheme determine the values that best fit the overall observations. Compared to g_LNA, smaller optimized g_HNA values (Tables S2–S6) reflect preferential grazing of HNA cells by microzooplankton, consistent with previous speculations that grazers selectively remove larger and more active bacterial cells (del Giorgio et al., 1996; Gonzalez et al., 1990; Sherr et al., 1992), so HNA bacteria (Garzio et al., 2013). Together with the higher mean cell-specific grazing rates for HNA bacteria (Sect. 3.2), our results suggest preferential grazing of HNA cells by microzooplankton.

The portability index in our study reflects the extent to which a single model framework represented by its distinctly optimized parameters in the same model equations captures the observed variability in different years, given variable environmental forcing and the accompanying shift in plankton ecosystem structure. The model parameter set optimized for 2012–13 was the most portable, while the model parameter set optimized for 2011–12 was the least portable (Table 3), in which the most (n=7 out of total 11) and the least numbers (n=5 out of total 11) of parameters were constrained, respectively (Tables S3 and S4). The other two years exhibited intermediate levels of model portability, with similar portability index values characterized by the same number of the constrained parameters (n=6 out of total 10 for 2010–11 and n=6 out of total 12 for 2013–14; Tables S2 and S5). In other words, it is the number of the constrained parameters that matters the most in driving high model portability, suggesting a connection between overfitting and the portability of the optimized parameter sets. Also, varying degrees of portability across the 4 study years rendered it difficult to choose one particular year's model solution to perform the climate change experiments (Sects. 3.4 and 4.4), consistent with the characteristics of the original WAP-1D-VAR v1.0 model. Instead, better model skill was found by utilizing the parameters from assimilating the climatological observations (i.e., the climatological model).

4.2 Bacterial carbon stocks and flows

The fact that cell-specific BP, respiration, and SDOC uptake rates of HNA bacteria were significantly higher than those of LNA bacteria (Sect. 3.2) is mainly because of the way the parameter optimization was conducted (Text S3). The higher initial parameter values assigned for HNA bacterial growth, RDOC excretion, mortality, and respiration rates (Table 1) might drive not only their faster cell-specific growth rates but also their higher DOC uptake rates to coexist with their LNA counterparts when the loss rates were relatively large for HNA bacteria. Though driven by the model assumptions, the important aspect of these results lies in the fact that the model can leverage such assumptions to examine the implications for the WAP food-web dynamics and biogeochemistry. As with phylogenetic groups (Fuchs et al., 2000; Teira et al., 2009; Yokokawa et al., 2004), cell-specific bacterial growth rates are expected to differ among distinct bacterial physiological groups, but there are limited studies focusing on group-specific cell activities (Gasol et al., 1999; del Giorgio et al., 1996; Günter et al., 2008; Longnecker et al., 2005; Moràn et al., 2011). Moràn et al. (2011) showed that HNA bacteria greatly outgrew LNA bacteria in Waquoit Bay estuary, with a cell-specific growth rate of up to 2.26 d⁻¹ for HNA cells versus < 0.5 d⁻¹ for LNA cells. Other studies have demonstrated that HNA bacteria might depend on phytoplankton substrates more than LNA bacteria (Li et al., 1995; Morán et al., 2007; Scharek and Latasa, 2007). The hypothesis that WAP bacteria might rely on SDOC when limited by LDOC availability has received indirect support previously (Ducklow et al., 2011; Kim and Ducklow, 2016; Luria et al., 2017), providing the basis for bacterial SDOC utilization in our model formulation.

The model also captured the rest of the ecosystem variables fairly well. The modeled nutrient stocks were above the detection limits, indicating no evidence of macronutrient limitations at the study site. The WAP typically exhibits strong interannual variability in physical forcing and ecological and biogeochemical processes (Ducklow et al., 2007), but the lack of strong interannual variability in the modeled microzooplankton C biomass is due to assimilating their climatological observations. One exception is krill C biomass that was modeled 3–8 times larger than the maximum value from the available field measurement in 2017–2018 (0.57 mmol C m⁻³; not shown). It should be noted that there were inconsistences in the nature of the assimilated data types, including a single-year observation of microzooplankton C biomass (versus each year-specific observation of other variables) and two unassimilated data types (e.g., krill C biomass). Also, there can be compensating errors in krill grazing rate and metabolism values given that krill are mobile laterally. These observational limitations make it challenging to construct a complete bacterial C budget without significant uncertainties. A more complete assimilation of zooplankton data should be the next effort to improve the model fits and minimize uncertainties in the bacterial variables. Another source of uncertainty in our study is that the model forcing does not seem to have sufficient information to capture small-scale and high-frequency sources of variability (e.g., local circulation and tidal flow near Palmer Station), resulting in relatively low standard deviation values of the modeled bacterial and ecosystem variables than those of the observations (e.g., Figs. 3 and S2). By contrast, our model adequately captures seasonal variations in modeled ecosystem dynamics likely because such high frequency processes do not strongly rectify in the seasonal cycles in the WAP ecosystem.

4.3 Bacterial physiological and taxonomic association with ecosystem functions

The positive associations of the observed fHNA with the modeled NPP and POC sinking flux suggest a relatively strong resource control on these actively growing HNA cells compared to slow-growing LNA cells. This is consistent with previous studies showing increased HNA growth rates in response to enhanced phytoplankton-derived organic substrate (Morán et al., 2011) and more abundant HNA cells in areas or periods in which bacterial assemblages were predominantly controlled by resources rather than grazing (Morán et al., 2007). It has been hypothesized that due to minimal inputs of terrestrial organic matter, bacteria must ultimately rely on in situ NPP as an organic matter source in the WAP (Ducklow et al., 2012a), supporting the importance of resource control on these actively growing bacterial populations.

In our study, modes 3, 5, and 7, characterized by copiotrophic taxa with large genomes and more 16S rRNA gene copies (Bowman et al., 2017), were associated with high values of the modeled NPP, POC sinking flux, and BCD, while modes 4 and 6, characterized by taxa associated with more oligotrophic conditions, were associated with low values of the modeled NPP, POC sinking flux, and BCD. Dokdonia sp. MED134, a common bacterial species of the modes associated with high NPP, POC sinking flux, and BCD, is a proteorhodopsin-containing marine flavobacterium that grows faster with light (Gómez-Consarnau et al., 2007; Kimura et al., 2011) and in conditions under which resources are abundant (Gómez-Consarnau et al., 2007). Given the coastal WAP being primarily light-limited (Ducklow et al., 2012b), the correspondence of D. Dokdonia MED134 to high values of the modeled NPP suggests light-enhanced growth rates and cell yields from sufficient irradiance. By contrast, mode 4, dominated by Planktomarina temperata RCA23, is a slow-growing bacterium that specializes in using complex organic substrates (Giebel et al., 2013). These attributes are consistent with the high occurrence of mode 4 during the periods of low values of the modeled NPP and POC sinking flux. Candidatus Pelagibacter, abundant in mode 6, is generally known as an oligotrophic specialist with a low DOC requirement but often appears during Antarctic phytoplankton blooms (Delmont et al., 2014; Luria et al., 2014), the characteristics of which support its occurrence during the periods of high values of the modeled NPP. In summary, our study provides a novel numerical framework combining the dynamics of different ecosystem functions and microbial physiology and taxonomy. Certain modes represent distinct WAP ecosystem states, and the mode–state associations are reasonably explained from microbial perspectives. However, we did not investigate a seasonal succession and development in mode itself or the mode association of the key WAP ecosystem states. Future investigations should focus on including a few dominant or seasonally distinct modes in the data assimilation process in order to fully resolve the seasonality of the mode–ecosystem state associations along the WAP.

4.4 Climate change experiments

The WAP has experienced significant atmospheric and ocean warming and resulting changes in marine ecological processes, and further climate change is projected for the next several decades. The magnitudes of the perturbations used in the climate change experiments (+0.5 ^∘C and +1.0 ^∘C compared to observed temperature fields and −5 % and −10 % compared to observed sea-ice fields) are within the range of the long-term changes in temperature and sea-ice duration along the WAP continental shelf. The temperature of the ACC water that has direct access to the WAP shelf has shown a large increase after the 1980s, equivalent to a uniform warming of the upper 300 m layer by 0.7 ^∘C (Ducklow et al., 2012b). The trend in the annual ice season duration is −1.5 d yr⁻¹ over 1979–80 to 2017–18 field season (Henley et al., 2019). The degree of melting (5 %–10 %) chosen for the climate change experiments is translated into the shortening of the ice season duration by 1–3 d (not shown), falling within the range of the trend in Henley et al. (2019).

Under combined warming and melting conditions, we expected that increased NPP and phytoplankton accumulations early in the season would result in a significant build-up of all DOC pools. However, this was the case only for SDOC, and bacteria were soon LDOC-limited due to their preferential LDOC uptake for their primary C source. Nonetheless, the growth of bacteria and increased bacterial rates under LDOC limitation were still possible because bacteria depended on SDOC to meet the rest of their C demand, resulting in the strong depletion of the SDOC pool later in the season (Fig. 7b). In other words, bacteria were more likely resource-limited, in particular by the labile DOC pool, and SDOC subsequently played an increasingly important role. This change was particularly important in HNA bacteria, as shown by the relatively large increase in HNA bacterial C demand via SDOC compared to LNA bacteria. Temperature is often regarded as a major factor regulating physiological rates by changing the rate of enzymatic reactions (Kirchman et al., 2009; White et al., 1991). In our study, the modeled C stock and rates of HNA bacteria increased under the warming alone conditions (Fig. S11) but equally or more than under the melting alone conditions (i.e., increased photosynthesis and resource availability; Fig. S12). This suggests that temperature per se is not necessarily a more important limiting factor for bacterial growth, at least for HNA bacteria, than resource availability (Ducklow et al., 2012a), and warming may rather enhance HNA bacterial utilization of the already increased organic matter from the increased phytoplankton productivity. Also, future climate may impact the (re)distribution of bacterial taxonomic groups, with a potential shift to more abundant HNA cells in the WAP bacterial communities owing to their preferential SDOC utilization.

The major limitation of our climate change experiments is the short duration of the simulations. An ideal set of climate change simulations should be performed for longer-term periods, as well as continuously across many years. However, our study could not accommodate these requirements because of the limited observations and existing data gaps in each year. Despite these challenges, we were able to validate the capacity of the climatological model to partly reproduce the already observed, climate-driven trends of some ecosystem variables along the WAP. Under each year's forcing fields, the climatological model parameter set reproduced the interannual variability fairly well compared to the observed interannual variability, except for only a few cases (e.g., overestimated BP and HNA biomass in 2011–12, underestimated PP in 2012–13 and 2013–14; Table S7). The period of 2011–12 was characterized by a negative temperature anomaly (−0.13 ± 0.83 ^∘C versus 0.03 ± 0.84 ^∘C for the 4-year climatology) and a positive sea-ice anomaly (24 ± 38 % versus 21 ± 29 % for the 4-year climatology), with lower temperature and higher sea-ice concentrations than the other 3 years (all p < 0.05, two-sample t test). This coldest year had the lowest values of BP, HNA biomass, and PP observations (Table S7), consistent with increases in the modeled BP, HNA biomass, and PP under the combined warming and melting conditions. A combination of low HNA biomass, low PP, and low POC flux was also modeled in 2011–12, being largely responsible for driving the positive association of the observed fHNA with the modeled NPP and POC sinking across years (Sect. 4.3). Sea ice did not retreat until mid-December in 2011–12 (Fig. S1), and as a result of subsequently low light levels PP was modeled to be low. The low modeled PP drove both low HNA biomass and low particle sinking flux, reinforcing the strong resource control on these fast-growing bacterial populations and the conventional “high PP-high export” paradigm along the WAP.

Finally, our climate change simulations share similar results with those performed using the WAP-1D-VAR v1.0 model with one bacterial compartment (Kim et al., 2021). In the original WAP-1D-VAR v1.0 model, combined warming and reduced sea-ice conditions also increased NPP, net community production, POC sinking flux, bulk bacterial productivity and biomass, and SDOC, in contrast to LDOC that was strongly limited early in the season. This potential shift to a more productive and efficient export system state is partially in agreement with the speculations suggested by previous studies that warming may induce more recycling-favorable and microbial-dominated food webs (Moline et al., 2004; Sailley et al., 2013). Despite the increased productivity and plankton accumulations, LDOC may become strongly depleted, and, therefore, bacteria may need to depend more on SDOC to meet a significant part of their C demand (i.e., an increasingly important role of SDOC for bulk bacterial communities). Most of these results convey the same story as our experiments, thereby adding confidence to the results of the climate change experiments in our study. Yet, it should be noted that the increased complexity of bacterial dynamics in our study's bacteria-oriented model adds two important contributions to the original WAP ecosystem model including (1) the dominance of HNA bacteria over LNA bacteria in the warming WAP waters and (2) bacterial taxonomic (i.e., mode) and physiological (i.e., fHNA) traits being a significant indicator of the key WAP ecosystem functions.

5 Conclusions

Heterotrophic microbial diversity has seldom been considered in detail in the formulation and analysis of marine pelagic ecosystem models, reflecting in part the lack of suitable field data for model evaluation. Utilizing genomic products to prescribe the taxonomic aspects of bacterial dynamics, our study demonstrates the association of bacterial abundance with different physiological states, bacterial community structure, and key ecosystem functions. The modeling approach in our study enables the observations in different bacterial populations to constrain the group-specific processes and model parameters that have been poorly understood. These include the partitioning of BP specific to HNA and LNA groups, the partitioning of the bacterial uptake of DOC pools with different lability, and the half-saturation density of each bacterial group in microzooplankton grazing. The model also serves as an effective numerical platform to explore the WAP microbial response to changing climate conditions, in which ocean warming and melting sea ice would induce a potential shift to the dominance of HNA bacteria in more productive waters due to their increasing dependence on SDOC.