Modeling Selective Pressures on Phytoplankton in the Global Ocean

Our view of marine microbes is transforming, as culture-independent methods facilitate rapid characterization of microbial diversity. It is difficult to assimilate this information into our understanding of marine microbe ecology and evolution, because their distributions, traits, and genomes are shaped by forces that are complex and dynamic. Here we incorporate diverse forces—physical, biogeochemical, ecological, and mutational—into a global ocean model to study selective pressures on a simple trait in a widely distributed lineage of picophytoplankton: the nitrogen use abilities of Synechococcus and Prochlorococcus cyanobacteria. Some Prochlorococcus ecotypes have lost the ability to use nitrate, whereas their close relatives, marine Synechococcus, typically retain it. We impose mutations for the loss of nitrogen use abilities in modeled picophytoplankton, and ask: in which parts of the ocean are mutants most disadvantaged by losing the ability to use nitrate, and in which parts are they least disadvantaged? Our model predicts that this selective disadvantage is smallest for picophytoplankton that live in tropical regions where Prochlorococcus are abundant in the real ocean. Conversely, the selective disadvantage of losing the ability to use nitrate is larger for modeled picophytoplankton that live at higher latitudes, where Synechococcus are abundant. In regions where we expect Prochlorococcus and Synechococcus populations to cycle seasonally in the real ocean, we find that model ecotypes with seasonal population dynamics similar to Prochlorococcus are less disadvantaged by losing the ability to use nitrate than model ecotypes with seasonal population dynamics similar to Synechococcus. The model predictions for the selective advantage associated with nitrate use are broadly consistent with the distribution of this ability among marine picocyanobacteria, and at finer scales, can provide insights into interactions between temporally varying ocean processes and selective pressures that may be difficult or impossible to study by other means. More generally, and perhaps more importantly, this study introduces an approach for testing hypotheses about the processes that underlie genetic variation among marine microbes, embedded in the dynamic physical, chemical, and biological forces that generate and shape this diversity.


Global ocean ecosystem model
The ocean ecosystem model used in this study has been described previously by Follows et al. [1]. Changes in the concentration of the ith nutrient (N i ), jth phytoplankton group (P j ), kth zooplankton group (Z k ) and pools of particulate organic matter (POM i ) and dissolved organic matter (DOM i ) are modeled with the following equations: The parameters used in these equations are defined below (Table S1) and by Follows et al. [1] and Dutkiewicz et al. [2]. Inorganic nitrogen, phosphorus, iron and silicon nutrients (N i ) are taken up by phytoplankton. The rates of uptake are determined by the abundance (P j , expressed in terms of phosphorus biomass), growth rate ( j μ ) and elemental ratios (R ij ) of the phytoplankton. The inorganic nutrients are supplied by different processes that are represented here by S i and can include remineralization, supplies from external sources, and abiotic transformations. Phytoplankton growth rates ( j μ ) are determined by irradiance, nutrient availability, and temperature (see below).
Phytoplankton sink through the water column according to a sinking rate, P j w , and are lost by mortality P j m and grazing by zooplankton. The rate of grazing by the kth zooplankton on the jth phytoplankton is represented by g jk , which in turn is given by the . Here g mjk is the maximum rate at which the kth zooplankton can graze the jth phytoplankton per unit of phytoplankton biomass, η jk is the palatability of the jth phytoplankton to the kth zooplankton, A k is the sum of phytoplankton weighted by their palatability, The constants u and κ represent advection and diffusion (respectively) in the model ocean and were provided by the ECCO ("Estimating the Circulation and Climate of the Ocean") state estimates [3]. These equations were integrated in a 1º X 1º resolution ocean circulation model with 23 vertical layers [4]. The advection scheme used for these simulations was the 'superbee' second order flux limiter [5]. We also applied an additional flux limiter to prevent spurious negative concentrations that occur due to the enforcement that mixing occur along isopycnal surfaces.
These simulations were performed with 33 different types of phytoplankton, and 2 different types of zooplankton. The physiological properties and nutrient requirements of different phytoplankton groups were assigned randomly, subject to some simple taxonomic and allometric constraints, using the approach of Follows et al. [1]. In each integration, we initialized 15 'small' phytoplankton types and 18 'large' phytoplankton types. Small phytoplankton have slower maximum growth rates than large phytoplankton, but are able to grow at lower levels of light and nutrient availability (Fig.  S1). Initial concentrations of nutrients and grazers were determined using previous simulations [1]. An ensemble of 10 integrations was performed, using different randomly generated phytoplankton communities.

Phytoplankton growth and nutrient use
, where A, B, C and T O regulate the shape of the dependence of growth on temperature, and τ 1 and τ 2 are used to normalize T j γ to a maximum value (see Fig. S1A). The response of growth to local irradiance (I) is represented by the function ) exp( )) exp( 1 where k PAR controls the saturation of growth as a function of irradiance, k inhib controls the inhibition of growth by high irradiance, and F m is a factor that is used to normalize I j γ to a maximum value of 1 (Fig. S1B). The degree of nutrient limitation is determined by the ambient nutrient concentration that supports the lowest rate of growth, or ) For iron, phosphorus and silica (diatoms only), , where ij N k is the concentration of the ith nutrient at which the jth phytoplankton reaches half its maximum growth rate (Fig. S1C). Nitrogen is available in three inorganic forms, NO 3 , NO 2 and NH 4 . The use of NO 3 and NO 2 are repressed when NH 4 is available. This is represented in the model using the expression During the course of the simulations, we introduce phytoplankton mutants (see below) that are unable to use NO 3 , but retain the ability to use NO 2 and NH 4 . For these mutants, we use the expression Other mutants lose the ability to use both NO 3 and NO 2 , but retain the ability to use NH 4 . For these mutants, we represent limitation by nitrogen as

Mutation
In the present study, we introduce mutations for the loss (and gain) of specific nitrogen use abilities. We apply these mutations to each of the 15 'small' phytoplankton types, beginning after the third year of the simulation. At this point the model has settled into a relatively steady annually repeating cycle of phytoplankton abundances. Initially, all biomass in a specific phytoplankton type can use NO 3 , NO 2 and NH 4 as nitrogen sources. We refer to this biomass as comprising the 'parent' population. The process of mutation (see Fig. 1B) is simulated by transferring biomass between the parent population and an additional tracer that represents the biomass of 'mutant' phytoplankton that have identical properties to their parent, except for a specific difference that corresponds to the mutated phenotype.
The changes in the concentration of mutant (P M ) and parent (P P ) biomass are modeled as follows: where the growth of a specific parent population by cell division is described by , as in Eq. SE2 and described in Section 2 of the Methods. A fraction, Λ f , of cell divisions result in a mutation (such as the loss of the ability to use nitrate) producing mutant biomass (P M ) at a rate of where log 2 (e) adjusts the growth rate to cell divisions per day. We also allow a fraction, Γ f , of cell divisions of mutant biomass to result in mutation back to the parent phenotype. This results in the production of parent biomass at a rate of L P and L M represent the change in the concentration of parent and mutant phytoplankton (respectively) per unit biomass from sinking, mortality and grazing (as in Eq. SE2), or We use this approach to produce three types of mutants for each of the 15 different groups of small phytoplankton. These are (i) NO 3 loss mutants, which lose the ability to use nitrate, (ii) NO 3 /NO 2 loss mutants, which lose the ability to use both nitrate and nitrite, and (iii) null mutants, which retain the ability to use ammonium, nitrite and nitrate.
For each of the three mutant types, the proportion of replications of the parent population resulting in mutation (Λ f ) was set to a value of Λ f = 10 -8 . Observed mutation rates of microbes with DNA genomes are often approximately 0.003 per genome per replication [6,7]. For picocyanobacteria with genomes of size 1.75 Mbp -3 Mbp, we would therefore expect mutation rates of approximately 10 -9 per base pair per replication. Typically, mutations at a number of locations (bases) in the genome can lead to the loss of function of a particular gene or trait. This is called the 'target size' of a mutation. We do not know the target size for the loss of nitrate or nitrite use abilities in picocyanobacteria. We chose a value of Λ f that is ten-fold greater than the expected mutation rate (Λ f = 10 -8 = 10 -9 × 10), corresponding to a target size of approximately 10 bases. It is possible that target sizes for the loss of nitrogen use traits are substantially larger than this value. If so, we can think of the mutants in the model as representing cells with a mutation at a location that is part of a specific subset (i.e., 10 specific bases) of the locations in the mutational target for a phenotype. For instance, these might represent a specific group of mutations that can cause (i) the loss of function of nitrate reductase (NO 3 loss mutant phenotype) (ii) the loss of function of nitrite reductase or (NO 3 /NO 2 loss mutant phenotype) and (iii) no change in nitrogen use abilities (null mutant phenotype). Given the considerable uncertainty concerning mutation rates (per base pair) and the target size for nitrogen use targets, we performed additional simulations using a greater mutation rate (Λ f = 10 -6 see Results, Section 2. ii), to make sure our observations were robust to variation in this parameter.
We set the rate of mutations that return mutants to the parent phenotype to a value of Γ f = 10 -12 . It is not clear how often nitrate and nitrite use abilities might be newly acquired by cyanobacteria that lack them, though we anticipate it would be much smaller than the rate at which they are lost. We therefore chose a value that was several orders of magnitude smaller than Λ f . At this value, Γ f is so small (and values of P M remain so small) that is expected to have a negligible affect on the model solutions.

Supply of inorganic nutrients
The supply of inorganic nutrients is represented in equation SE1 by the term i N S , which includes different processes for the inorganic nutrients of different elements. The equations describing the supply of inorganic phosphorus, nitrogen (NO 3  (respectively) forms of the ith element. The term c scav is the rate of scavenging of free iron, whose concentration, Fe', is modeled following the approach of Parekh et al. [8].
The term F atmos represents the rate at which iron is deposited at the surface of the ocean in dust, and α represents the solubility of this iron. The oxidation of NO 2 to NO 3 is described by the term 3 NO ζ , and the oxidation of NH 4 to NO 2 is described by the term 2

NO
ζ . In addition, we assume that the oxidation of NH 4 to NO 2 is photoinhibited, such that it does not occur if the ambient light is above 10 μE m -2 s -1 .

Ambient nitrogen concentrations
The annual average concentrations of nitrate, nitrite and ammonium in surface waters (upper 10 m) of the modeled ocean are illustrated in Fig. S2 A, C and E. Annual average concentrations of nitrate, nitrite and ammonium along Atlantic Meridional Transect [9] 13 (AMT 13) in the modeled ocean are illustrated in Fig. S2 D, E and F.

Sensitivity of mutant accumulation to model parameters i. Time dependence of mutant accumulation
In the main text, our analyses concentrate on the abundance of mutants during the fifth year of a model integration (Figs 2 and 3). However, we expect the abundance of mutants in the model to change over time, as they are produced by parent populations. In Fig. S3, the distribution of phytoplankton in surface waters of the global ocean is illustrated for the tenth year of the integration. The distributions of small (parent) and large phytoplankton are similar in the fifth and tenth years. In the tenth year of the integration, mutants have accumulated to greater abundances than in the fifth year, but similar patterns are observed for the relative abundance of NO 3 /NO 2 loss mutants and null mutants in different ocean regions. That is, while NO 3 /NO 2 loss mutants reach abundances similar to null mutants in topical oligotrophic regions, null mutants are substantially more abundant than NO 3 /NO 2 loss mutants at higher latitudes. This means that the salient observations reported in the main text are not an artifact of the time at which the distribution of the mutants was studied.

ii. Mutation rate
We do not know the mutation rates for the loss of NO 3 use abilities or NO 3 /NO 2 use abilities in marine cyanobacteria (see Methods, Section 3). We therefore wanted to study the distribution of mutants in the model using a different value for the mutation rate, and performed a simulation where the proportion of cell divisions resulting in mutations (Λ f ) was set to Λ f = 10 -6 , instead of Λ f = 10 -8 . In the simulation where Λ f = 10 -6 ( Fig. S4) the abundance of NO 3 /NO 2 loss mutants and null mutants are approximately two orders of magnitude greater than in the simulation where Λ f = 10 -8 (Fig. 2). However, changing the parameter Λ f to 10 -6 did not substantially change the relative abundance of NO 3 /NO 2 loss mutants and null mutants in different ocean regions. This suggests that our observations are robust to variation in the mutation rate over the range of values that were studied.
iii. Oxidation rates of NH 4 and NO 2 We used a model that parameterizes important biogeochemical transformations of nitrogen (see Methods, section 4). There is considerable uncertainty surrounding the rates at which some of these transformations occur in the environment. We therefore performed simulations using different values for two such parameters. Specifically, we used smaller values for the rate of (light-inhibited) oxidation of ammonium to nitrite ζ result in small changes in the distribution of ammonium and nitrite, including greater concentrations of ammonium in the lower euphotic zone, and near the surface in several regions with high productivity. However, changing the rates of oxidation of NH 4 and NO 2 did not substantially alter the observation that NO 3 /NO 2 loss mutants reach similar abundances to null mutants in topical oligotrophic regions, but at higher latitudes, null mutants reach much higher abundances than NO 3 /NO 2 loss mutants.

iv. Randomly generated phytoplankton communities
We performed integrations of the model with 10 randomly generated phytoplankton communities. In the main text, results are presented for one of these communities. In Fig.  S7, the distribution of phytoplankton in global surface waters is illustrated for the remaining nine phytoplankton communities. There is some variation among these integrations, for example in the dominance of particular regions by large versus small phytoplankton. However, some important features of the distribution of NO 3 /NO 2 loss mutants and null mutants are common to each of these integrations. In each case, NO 3 /NO 2 loss mutants accumulate to similar abundances to null mutants in tropical oligotrophic regions, but at higher latitudes, null mutants occur at substantially greater abundances that NO 3 /NO 2 loss mutants, indicating that NO 3 /NO 2 loss mutants are disadvantaged in these regions.