Isotope labelling and phylogenetic analysis of a natural marine bacterioplankton population at sea

Mediterranean seawater was collected during August 2017 (station N1200, 32.45° N, 34.37 °E) from 11 depths by Niskin bottles and divided into triplicate 250 ml polycarbonate bottles. Two bottles from each depth were labelled with 1 mM sodium bicarbonate-13C and 1 mM ammonium-15N chloride (Sigma-Aldrich), and all three bottles (two labelled and one control) were incubated at the original depth and station at sea for 3.5 h around mid-day. The stable isotopes were chosen to enable direct comparison of C and N uptake in single cells, and the short incubation time was chosen to minimize isotope dilution and potential recycling and transfer of 13C and 15N between community members25. After incubation, bottles were brought back on board and the incubations were stopped by fixing with 2× electron-microscopy-grade glutaraldehyde (2.5% final concentration) and stored at 4 °C until sorting analysis. Cell sorting, NanoSIMS analyses and the calculation of uptake rates were performed as described in Roth-Rosenberg et al.26.

DNA collection and extraction from seawater

Samples for DNA were collected on 0.22 µm Sterivex filters (Millipore). Excess water was removed using a syringe, 1 ml lysis buffer (40 mM EDTA, 50 mM Tris pH 8.3, and 0.75 M sucrose) was added and both ends of the filter were closed with parafilm. Samples were kept at −80 °C until extraction. DNA was extracted by using a semi-automated protocol including manual chemical cell lysis before automated steps using the QIAamp DNA Mini Protocol: DNA Purification from Blood or Body Fluids (Spin Protocol, starting from step 6, at the BioRap unit, Faculty of Medicine, Technion). The manual protocol began with thawing the samples, then the storage buffer was removed using a syringe and 170 µl lysis buffer added to the filters. Thirty microlitres of Lysozyme (20 mg ml−1) were added to the filters and incubated at 37 °C for 30 min. After incubation, 20 µl proteinase K and 200 µl buffer AL (from the Qiagen kit) were added to the tube for 1 h at 56 °C (with agitation). The supernatant was transferred to a new tube, and DNA was extracted using the QIAcube automated system. All DNA samples were eluted in 100 μl DNA-free distilled water.

ITS PCR amplification

PCR amplification of the ITS was carried out with specific primers for Prochlorococcus CS1_16S_1247F (5′-ACACTGACGACATGGTTCTACACGTACTACAATGCTACGG) and Cs2_ITS_Ar (5′-TACGGTAGCAGAGACTTGGTCTGGACCTCACCCTTATCAGGG)21,22. The first PCR was performed in triplicate in a total volume of 25 μl containing 0.5 ng of template, 12.5 μl of MyTaq Red Mix (Bioline) and 0.5 μl of 10 μM of each primer. The amplification conditions comprised steps at 95 °C for 5 min, 28/25 (16 S/ITS) cycles at 95 °C for 30 s, 50 °C for 30 s and 72 °C for 1 min followed by one step of 5 min at 72 °C. All PCR products were validated on a 1% agarose gel, and triplicates were pooled. Subsequently, a second PCR amplification was performed to prepare libraries. These were pooled and after a quality control sequenced (2 × 250 paired-end reads) using an Illumina MiSeq sequencer. Library preparation and pooling were performed at the DNA Services facility, Research Resources Center, University of Illinois at Chicago. MiSeq sequencing was performed at the W.M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign.

ITS sequence processing

Paired-end reads were analysed using the Dada2 pipeline46. The quality of the sequences per sample was examined using the Dada2 ‘plotQualityProfile’ command. Quality filtering was performed using the Dada2 ‘filterAndTrim’ command with parameters for quality filtering truncLen=c(290,260), maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, trimLeft=c(20,20). Following error estimation and dereplication, the Dada2 algorithm was used to correct sequences. Merging of the forward and reverse reads was done with minimum overlap of 4 bp. Detection and removal of suspected chimaeras was done with command ‘removeBimeraDenovo’. In total, 388,417 sequences in 484 amplicon sequence variants were counted. The amplicon sequence variants were aligned in MEGA6 (ref. 47), and the first ~295 nucleotides, corresponding to the 16S gene, were trimmed. The ITS sequences were then classified using BLASTn against a custom database of ITS sequences from cultured Prochlorococcus and Synechococcus strains as well as from uncultured HL and LL clades.

Individual-based model

PlanktonIndividuals.jl (v0.1.9) was used to run the individual-based simulations48. Briefly, the cells fix inorganic carbon through photosynthesis and nitrogen, phosphorus and DOC from the water column into intracellular quotas and grow until division or grazing. Cell division is modelled as a probabilistic function of cell size. Grazing is represented by a quadratic probabilistic function of cell population. Cells consume nutrient resources, which are represented as Eulerian, density-based tracers. A full documentation of state variables and model equations are available online at https://juliaocean.github.io/PlanktonIndividuals.jl/dev/. Equations related to mixotrophy are shown below as an addition to the online documentation.

$$V_}} = V_}}^}} \cdot }}\left( }}\left( }^}} - q_}}}}^}} - q_}^}}}}} \right)} \right) \cdot \frac}}}} + K_}}^}}}}$$

(1)

$$f_}} = \frac}}}} + V_}}}}$$

(2)

$$V_}} = 0,\,}\,f_}} < f_}}^}}$$

(3)

where VDOC is the cell-specific DOC uptake rate (mol C cell−1 s−1), \(V_}}^}}\) is the maximum cell-specific DOC uptake rate (mol C cell−1 s−1), \(q_}^}}\) is the maximum cell carbon quota (mol C cell−1), \(q_}^}}\) is the minimum cell carbon quota (mol C cell−1). The maximum and minimum functions here is used to keep qC between \(q_}^}}\) and \(q_}^}}\). \(K_}}^}}\) is the half-saturation constant for DOC uptake (mol C m−3). fPS is the fraction of fixed C originating from photosynthesis (PS, mol C cell−1 s−1). DOC uptake stops when fPS is smaller than \(f_}}^}}\)(minimum fraction of fixed C originating form photosynthesis, 1% by default) according to laboratory studies of Prochlorococcus that showed that they cannot survive long exposure to darkness (beyond several days) even when supplied with organic carbon sources13. (1 − fPS) is also shown in Fig. 3 as the contribution of DOC uptake.

We set up two separate simulations; each of them has a population of either an obligate photo-autotroph or a mixotroph that also consumes DOC. The initial conditions and parameters (Supplementary Table 3) are the same for the two simulations except the ability of mixotrophy. The simulations were run with a timestep of 1 min for 360 simulated days to achieve a steady state. We run the two simulations for multiple times in order to get the range of the stochastic processes.

Evaluation of autotrophic growth rates

We evaluated the carbon-specific, daily-averaged carbon fixation rate, ℙ as a function of light intensity (I, µE), following Platt et al.33:

$$ = \frac}_0^ }}}}}}}} P_}^}}\left( }}I/P_}^}}}} \right)e^}}I/P_}^}}}\Delta t$$

(4)

Here, \(P_}^}}\), αChl and βChl are empirically determined coefficients representing the chlorophyll-a-specific carbon fixation rate (mol C (mol Chl)−1 s−1), the initial slope of the photosynthesis–light relationship and photo-inhibition effects at high photon fluxes, respectively. We impose empirically determined values for \(P_}^}}\), αChl and βChl from the published study of Moore and Chisholm24. The natural Prochlorococcus community comprises HL and LL ecotypes, which have different values of \(P_}^}}\), αChl and βChl, and the community growth rate is expected to be between that of HL extremes and LL extremes. Therefore, we use photo-physiological parameters for an HL-adapted ecotype (MIT9215), acclimated at 70 µmol photons m−2 s−1 and an LL-adapted ecotype (MIT9211), acclimated 9 µmol photons m−2 s−1. The models with these values are shown as the different lines in Fig. 2b,d. I is the hourly PAR, estimated by scaling the observed noon value at each depth with a diurnal variation evaluated from astronomical formulae based on geographic location and time of year37,38.

\(\frac}}}}}}}\) is the molar chlorophyll-a to carbon ratio, which is modelled as a function of growth rate and light intensity using the Inomura34 model (equation 17 therein) where parameters were calibrated with laboratory data from Healey49. In addition, the maximum growth rate (\(\mu _}}^I\)) based on macromolecular allocation is also estimated using the Inomura model (equation 30 therein). An initial guess of the growth rate and the empirically informed light intensity are used to estimate \(\frac}}}}}}}\), which is then used to evaluate the light-limited, photoautotrophic growth rate

$$_}^}} = \min \left( - K_},\mu _}}^I} \right)$$

(5)

from which the \(\frac}}}}}}}\) is again updated. The light-limited growth rate is used to re-evaluate the \(\frac}}}}}}}\). Repeating this sequence until the values converge, \(_}^}}\) and \(\frac}}}}}}}\) are solved iteratively.

The nitrogen-specific uptake rate of fixed nitrogen (day−1) is modelled as

$$_}}} = _}^}}\frac}}}\frac}}}}}$$

(6)

where values of the maximum uptake rate, \(_}^}}\), and half-saturation, KN, are determined from empirical allometric scalings35, along with a nitrogen cell quota QN from Bertilsson et al.39.

The P-limited growth rate, or the phosphorus-specific uptake rate of phosphate (day−1), is modelled as

$$_} = _}^}}\frac}}}\frac_}^}}_}^ + K_}}}$$

(7)

where values of the maximum uptake rate, \(_}^}}\). and half-saturation, KP, are determined from empirical allometric scalings35, along with a nitrogen cell quota QP from Bertilsson et al.39.

Iron uptake is modelled as a linear function of cell surface area (SA), with rate constant (\(k_}}^}}\)) following Lis et al.36.

$$_}} = k_}}^}} \cdot }\frac}}}}}$$

(8)

The potential light-, nitrogen-, phosphorus- and iron-limited growth rates (\(_},_},_},_}}\)) were evaluated at each depth in the water column and the minimum is the local modelled photo-autotrophic growth rate estimate, assuming no mixotrophy (Fig. 2b,d, blue lines). Parameters used in this evaluation are listed in Supplementary Table 2.

An important premise of this study is that heterotrophy is providing for the shortfall in carbon under very low light conditions, but not nitrogen. It is known that Prochlorococcus can assimilate amino acids9 and therefore the stoichiometry of the heterotrophic contribution might alter the interpretations. However, it is also known that Prochlorococcus can exude amino acids40, which might cancel out the effects on the stoichiometry of Prochlorococcus.

For the estimates of phototrophic growth rate from local environmental conditions (Fig. 2) we employed photo-physiological parameters from laboratory cultures of Prochlorococcus24. For the purposes of this study, we have assumed that the photosynthetic rates predicted are net primary production, which means that autotrophic respiration has been accounted for in the measurement. However, the incubations in that study were of relatively short timescale (45 min), which might suggest they are perhaps more representative of gross primary production. If this is the case, our estimates of photo-autotrophic would be even lower after accounting for autotrophic respiration, and thus would demand a higher contribution from heterotrophic carbon uptake. In this regard, our estimates might be considered a lower bound for organic carbon assimilation.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.