Mixotrophic growth of a ubiquitous marine diatom

Abstract

Diatoms are major players in the global carbon cycle, and their metabolism is affected by ocean conditions. Understanding the impact of changing inorganic nutrients in the oceans on diatoms is crucial, given the changes in global carbon dioxide levels. Here, we present a genome-scale metabolic model (iMK1961) for Cylindrotheca closterium, an in silico resource to understand uncharacterized metabolic functions in this ubiquitous diatom. iMK1961 represents the largest diatom metabolic model to date, comprising 1961 open reading frames and 6718 reactions. With iMK1961, we identified the metabolic response signature to cope with drastic changes in growth conditions. Comparing model predictions with Tara Oceans transcriptomics data unraveled C. closterium's metabolism in situ. Unexpectedly, the diatom only grows photoautotrophically in 21% of the sunlit ocean samples, while the majority of the samples indicate a mixotrophic (71%) or, in some cases, even a heterotrophic (8%) lifestyle in the light. Our findings highlight C. closterium's metabolic flexibility and its potential role in global carbon cycling.

Fig. 1.. Metabolic network reconstruction and features of iMK1961.

(A) Reconstruction pipeline of the GEM of C. clostridium (see Supplementary Text and Materials and Methods). (B) Features of iMK1961. The donut chart represents the proportions of the gene (cyan blue) and non-gene–associated reactions (light cyan blue) in iMK1961. The bar plot illustrates the number of genes, metabolites, and reactions in iMK1961. (C) Distribution of reactions among different functional categories. These categories are defined on the basis of KEGG pathways. Reaction distribution shows that lipid metabolism constitutes the largest portion of the total reactions in the network. (D) Reaction distribution per compartment. Among the seven compartments, the cytosol contains most of the reactions.

Fig. 2.. Comparison of iMK1961 with previously published diatom GEMs.

(A) Comparative representation of genes, metabolites, and reactions in different diatom GEMs. C. closterium (iMK1961) and F. cylindrus GEM (iML830) include the highest and lowest number of reactions and metabolites. Apart from this plot, iThaps987 was not considered in this analysis because of its inconsistent reaction and metabolite identifiers. (B) The Multiple Correspondence Analysis (MCA) scatter plot shows dissimilarities between different diatom GEMs in terms of reactions and metabolites [n = 8803 (iMK1961), n = 7660 (iTps1432), n = 5909 (iLB1027_lipid), and n = 3005 (iML830)]. The scatter plot illustrates the first two components, which explained 61.85 and 22.70% of the variances. We compared the entire metabolic content (reactions and metabolites) of iMK1961 with previously published diatom GEMs using MCA (80). iMK1961 and iTps1432 (T. pseudonana) are placed closer to each other as they exclusively share 1274 reactions and 389 metabolites. F. cylindrus (iML830) is distinct from the other three GEMs because this model comprises the lowest number of reactions and metabolites. iMK1961 has 906 unique reactions and 429 unique metabolites, which are not present in any previous diatom models. These unique reactions are mostly distributed in pathways related to transport, lipid metabolism, amino acid metabolism, carbohydrate metabolism, cofactors and vitamins metabolism, nucleotide metabolism, energy metabolism, TCA cycle, and glycolysis/gluconeogenesis (fig. S2). (C) Venn diagrams represent the distinct and shared reactions and metabolites among different diatom GEMs. C. closterium GEM (iMK1961) includes the highest number of unique reactions and metabolites. It exclusively shared the maximum number of reactions with iTps1432 followed by iLB1027_lipid and iML830. (D) Comparison between diatom GEMs based on reaction distribution among different pathways. Reactions were defined in 10 functional categories based on KEGG metabolic pathways.

Fig. 3.. Validation of iMK1961.

(A) Heatmap represents experimentally confirmed predictions of metabolic phenotypes on various nutrient conditions composed of different carbon, nitrogen sources in photoautotrophic and heterotrophic, and mixotrophic conditions. Nitrate was used as a nitrogen source to stimulate growth in all three trophic conditions on different carbon sources. Similarly, to stimulate growth on different nitrogen sources, CO₂, succinate, and CO₂ + succinate were used in photoautotrophic, heterotrophic, and mixotrophic conditions, respectively. White cells in the heatmap represent no growth data available based on the definition of photoautotrophy, heterotrophy, and mixotrophy. (B) Comparison between in silico and experimental growth rates. Experimental data were used to constrain the inputs of the model. The biomass-producing reaction was used as an objective function during each simulation.

Fig. 4.. Model predicted differential metabolic fluxes under elevated concentrations of nutrients.

(A) Heatmap represents differential flux-carrying subsystems under elevated CO₂, HCO₃, NO₂, NO₂ + NO₃, PO₄, and Si conditions. It shows subsystems separately with increased and decreased fluxes due to change in nutrient uptake fluxes from low to high levels. The numbers in parentheses denote the number of differential flux-carrying subsystems for each nutrient. The details of altered subsystems under varying growth conditions can be seen in table S19. We compared the metabolic flux distributions through the network when the model was simulated using the low and high concentrations of nutrients. Low and high uptake fluxes of nutrients to the model were used to mimic the low and high concentrations that were obtained from metadata of Tara Oceans global ocean microbiome data (table S11). The low and high concentrations of nutrients were defined by comparing with average concentration values for each nutrient (see Materials and Methods). (B) This Venn diagram represents shared and unique subsystems of C. closterium with increased flux distributions under the elevated concentration of CO₂ and HCO₃. The shared and unique subsystems in other pairs of conditions can be seen in figs. S5 to S7 and tables S17 and S18.

Fig. 5.. Intersection between differential model–predicted fluxes and gene expressions.

(A) Bar plot represents subsystems that carry differential (increased or decreased) fluxes and are associated with differentially expressed genes under elevated concentrations of nutrients, such as CO₂, HCO₃, NO₂, NO₂ + NO₃, PO₄, and Si. (B) As an example, under elevated CO₂, 27 subsystems illustrated differential metabolic fluxes and differentially expressed genes. This bar plot represents the number of reactions involved in these subsystems. Affected subsystems under additional conditions are shown in figs. S9 to S13. (C) Model-predicted storage of TAG and carbohydrate in terms of accumulation fluxes of TAG and β-1,3 glucans (i.e., monomers of chrysolaminarin), respectively, under low and high CO₂. Each box plot denotes all allowable possible fluxes that were estimated using a sampling method (31) while simulating the model (see Materials and Methods). Significance (P value <0.001) was computed using two-side Student’s t test.

Fig. 6.. Trophic modes of C. closterium in marine environments.

(A) Predicted and experimental growth rate of C. closterium under photoautotrophic, heterotrophic, and mixotrophic conditions. Extracellular succinate was used as an organic carbon source to stimulate growth under heterotrophic and mixotrophic conditions. Box plots represent all feasible biomass flux values that were determined by simulating iMK1961 using a uniform random sampling method (31) (see Materials and Methods) and experimentally measured growth rates. To simulate growth under trophic conditions, the model was constrained using experimentally measured uptake fluxes table S9. (B) CO₂ fixation in C. closterium during photoautotrophy, heterotrophy, and mixotrophy. CO₂ fixation was represented in terms of the uptake flux of CO₂ while simulating the growth in different conditions. (C) The predicted flux distributions and metatranscriptomics data helped to identify trophic modes of C. closterium in the Tara Oceans samples. Bar plots represent the number of samples where solely photoautotrophic, heterotrophic, or mixotrophic growth of C. closterium was detected. Unique sets of active reactions, which overlapped with expressed genes in transcriptomics data for each trophic mode, were used to categorize samples in three different modes (table S21 and figs. S14 to S16). (D) The global distribution of various trophic modes of C. closterium was identified using iMK1961 and metatranscriptomics data from Tara Oceans samples. Dot size represents the predicted growth rate of C. closterium under the corresponding trophic condition. (E) To identify significantly differentially abundant marine prokaryotes between different trophic modes, linear discriminant analysis (LDA) effect size (79) was deployed on Tara Oceans metagenomic data. Bar plot represents 10 most differentially abundant prokaryotes in each trophic mode [P value <0.05 (Kruskal-Wallis test and pairwise Wilcoxon test); LDA score (log₁₀) > 2.0] (see Materials and Methods). A complete list of differentially abundant prokaryotes can be seen in table S26.