Received Date: 04-Apr-2014 Accept Date: 13-May-2014
Accepted Article
1 2
Genome scale metabolic reconstruction and constraints-based modelling of the Antarctic bacterium Pseudoalteromonas haloplanktis TAC125. 1
3 4 5 6 7 8
Marco Fondi1,2§, Isabel Maida1, Elena Perrin1, Alessandra Mellera1,2, Stefano Mocali3, Ermenegilda Parrilli4, Maria Luisa Tutino4, Pietro Liò5, Renato Fani1,2
9
10
1
11
Piano 6, 50019 Sesto Fiorentino (Firenze)
12
2
13
Fiorentino (Firenze)
14
3
15
(CRA-ABP), Piazza d’Azeglio 30, 50121 Firenze, Italy
16
4
17
Cintia, I-80126, Naples, Italy
18
5
19
Kingdom§Corresponding author
20
Dr. Marco Fondi
21
Laboratory of Microbial and Molecular Evolution, Department of Biology, University of Florence, Via Madonna del
22
Piano 6, 50019 Sesto Fiorentino (Firenze)
23
Tel +390554574736
24
Email: [email protected]
Laboratory of Microbial and Molecular Evolution, Department of Biology, University of Florence, Via Madonna del
ComBo, Florence Computational Biology group, University of Florence, Via Madonna del Piano 6, 50019 Sesto
Consiglio per la Ricerca e la Sperimentazione in Agricoltura, Centro di Ricerca per l’Agrobiologia e la Pedologia
Department of Chemical Sciences, University of Naples Federico II, Complesso Universitario M. S. Angelo, Via
Computer Laboratory, Cambridge University, William Gates Building 15, JJ Thomson Avenue, Cambridge, United
This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/1462-2920.12513 This article is protected by copyright. All rights reserved.
Summary
27 28 29
The Antarctic strain Pseudoalteromonas haloplanktis TAC125 is one of the model organisms of cold-adapted bacteria and is currently exploited as a new alternative expression host for numerous biotechnological applications.
30 31 32 33 34 35 36 37 38
Here, we investigated several metabolic features of this strain through in silico modelling and functional integration of –omics data. A genome-scale metabolic model of P. haloplanktis TAC125 was reconstructed, encompassing information on 721 genes, 1133 metabolites and 1322 reactions. The predictive potential of this model was validated against a set of experimentally determined growth rates and a large dataset of growth phenotypic data. Furthermore, evidence synthesis from proteomics, phenomics, physiology and metabolic modeling data revealed possible drawbacks of cold-dependent changes in gene expression on the overall metabolic network of P. haloplanktis TAC125. These included, for example, variations in its central metabolism, amino acids degradation and fatty acids biosynthesis.
39 40 41 42 43
The genome scale metabolic model described here is the first one reconstructed so far for an Antarctic microbial strain. It allowed a system-level investigation of variations in cellular metabolic fluxes following a temperature downshift. It represents a valuable platform for further investigations on P. haloplanktis TAC125 cellular functional states and for the design of more focused strategies for its possible biotechnological exploitation.
Accepted Article
25 26
44
This article is protected by copyright. All rights reserved.
Accepted Article
45 46
Introduction
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Antarctica is one of the most extreme environments on Earth due to the presence of stably low temperatures. Many (psychrophilic) microorganisms inhabit this environment, being responsible for most of the biomass of this ecological niche and playing a key role in maintaining proper ecosystem functions (Wilkins et al., 2013). Such cold-adapted microorganisms have begun to attract the attention of the scientific community from both fundamental and application viewpoints (Cavicchioli et al., 2002; D'Amico et al., 2006). Indeed, psychrophiles have evolved peculiar features (Feller and Gerday, 2003; Giordano et al., 2012) to face the profound impact of cold shock on several cellular processes (e.g. growth rates, ribosomal synthesis, cytoplasmic membrane composition, etc.) (Phadtare, 2004). Overall, survival at low temperature depends on the ability of microbes to sense changes in temperature and to transduce these signals to the genome, ultimately bringing changes in the regulation of cold-adaptation related genes (Shivaji and Prakash, 2010). These include: i) genes for fatty acid desaturases (Sato and Murata, 1980), ii) genes involved in replication (Jones et al., 1992; Graumann and Marahiel, 1999), transcription (Sledjeski et al., 1996; Inaba et al., 2003) and translation (Jones et al., 1987; Xia et al., 2003), iii) genes encoding cold shock proteins (Brandi et al., 1994; Atlung and Ingmer, 1997), as well as iv) many genes coding for still uncharacterized proteins (Kawamoto et al., 2007; Ting et al., 2010). Moreover, the increase of enzymes concentration (Willem et al., 1999) and the expression of specific cold-adapted enzymatic isoforms (Hoyoux et al., 2004) represent two additional examples of how psychrophilic enzymatic machineries may counteract the growth limiting effects of cold temperatures.
66 67 68 69 70 71 72 73 74 75 76 77
Recently, the spreading of –omics technologies has allowed the system-level investigation of the mechanisms involved in microbial cold adaptation (Seo et al., 2004; Goodchild et al., 2005; Bakermans et al., 2007; Zheng et al., 2007; Piette et al., 2010). For example, a transcriptomics analysis on Shewanella oneidensis MR-1 showed that more than 70% of its genes involved in energy metabolism were down-regulated following a temperature downshift (Gao et al., 2006). A proteomic study on Sphingopyxis alaskensis revealed that a large fraction of genes involved in metabolic processes (e.g. energy production/conversion, carbohydrate, amino acids, nucleotides and cofactors transport/metabolism) are less abundant at lower temperatures (Ting et al., 2010). Finally, from a metabolic viewpoint, a large-scale metabolomics analysis on the cold adaptation of Mesorhizobium sp. strain N33 revealed a key role of (unsaturated) fatty acids biosynthetic process. From such examples, it emerges that, despite some strategies may be shared, different microorganisms use different strategies to cope with cold environments (Kawamoto et al., 2007).
78 79 80 81 82 83 84 85 86
The Antarctic marine bacterium Pseudoalteromonas haloplanktis TAC125 (PhTAC125) has been isolated from sea water sampled along the Antarctic ice-shell, a permanently cold environment. PhTAC125 is capable of growing in a wide temperature range (4–25°C) and its lowest observed doubling time was detected at 20°C (Medigue et al., 2005). Several exceptional genomic and metabolic features were derived from the genome sequence of this bacterium, showing adaptation to periodic situations of nutrient abundance (Medigue et al., 2005). Indeed PhTAC125 is considered to be one of the model organisms of cold-adapted bacteria and has been suggested as an alternative host for the soluble overproduction of heterologous proteins, given its capability to grow fast at low temperatures (Duilio et al., 2004; Wilmes et al., 2010; Rippa et al., 2012; Corchero et al., 2013). This article is protected by copyright. All rights reserved.
Furthermore, bacteria belonging to the genus Pseudoalteromonas are known to possess an inhibitory activity against human pathogens belonging to the Burkholdeia cepacia complex (Bcc) and to be able to produce antibiofilm molecules (Papa et al., 2013; Papaleo et al., 2013), revealing the interesting biotechnological potential and metabolic biodiversity of these microorganisms. The increasing interest in PhTAC125 has led to the accumulation of different data types for this bacterium in the last few years, including its complete genome sequence (Medigue et al., 2005), its proteome (Piette et al., 2010; Piette et al., 2011) and detailed growth phenotypes (Wilmes et al., 2010; Giuliani, 2011). Thus, it is now possible to integrate such different data sources and perform a system-level investigation of PhTAC125 metabolism.
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
Genome scale metabolic modeling represents a valuable tool in this context. Indeed, this in silico approach can be adopted to quantitatively simulate chemical reactions fluxes within the cell, including metabolic adjustments in response to external perturbations (e.g. temperature downshift). Genome annotations are usually transformed into models by defining the boundaries of the system, a biomass assembly reaction, and exchange fluxes with the environment (Durot et al., 2009; Thiele and Palsson, 2010). Constraint-based modelling methods (e.g. Flux Balance Analysis, FBA) can then be used to compute the resulting balance of all the active cellular reactions in the cell and to simulate the maximal growth of a cell in a given environmental condition (Varma and Palsson, 1994; Schilling et al., 2000). In the past decade, this approach has been successfully applied for studying large-scale metabolic networks in microbes, with the aim of guiding rational engineering of biological systems for applications in industrial and medical biotechnology (Milne et al., 2009). Interestingly, different data types can be mapped onto metabolic models in order to elucidate more thoroughly the metabolism of a cell and its response to environmental factors. This is usually done by including functional characterization and accurate quantification of all the main cellular information levels of gene products, mRNA, proteins and metabolites, as well as their interaction (Zhang et al., 2010). In recent years, –omics-derived data have been used to refine, validate and/or integrate metabolic models, including transcriptomics (Colijn et al., 2009; Jensen and Papin, 2011), proteomics (Gille et al., 2010), fluxomics (Chen et al., 2011; Feng and Zhao, 2013) and phenomics (Fang et al., 2011).
115 116 117 118 119 120 121 122 123
In at least two study-cases, FBA and expression data have been merged to explore the effect of temperature downshift at the system level (Navid and Almaas, 2012; Tong et al., 2013). In particular, Tong et al. (2013) performed robustness analysis on the core metabolic model of Thermoanaerobacter tengcongensis to study the dynamic changes of the metabolic network following the perturbation of the culture temperature and collecting the bacterial growth rates and differential proteomes. Given the overall agreement between in silico simulations and observed phenotypes, this approach was shown to provide a reliable platform to systematically evaluate the mechanisms of bacterial metabolism and relevant switches in the presence of a temperature downshift.
124 125 126 127 128 129
In this study, a genome-scale reconstruction of P. haloplanktis TAC125’s metabolism was performed based on its genome annotation. The predictive capability of the model [named iMF721 according to the current naming convention (Reed et al., 2003)] was successfully validated comparing constraints-based modeling outcomes with experimentally determined growth rates and large scale growth phenotype data (Phenotype Microarray). The iMF721 model was then used to globally investigate possible metabolic adjustments of P. haloplanktis TAC125 during growth at
Accepted Article
87 88 89 90 91 92 93 94 95
This article is protected by copyright. All rights reserved.
low temperature by means of robustness analysis and functional integration of protein abundance data into the reconstructed network.
132 133 134
Finally, the metabolic reconstruction reported herein represents a reliable platform for the future design of experimental strategies aimed at the exploitation of the biotechnological potential of PhTAC125.
Accepted Article
130 131
135 136 137
Experimental procedures
138
Metabolic model reconstruction
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
An initial draft metabolic reconstruction of the strain PhTAC125 was constructed using RAST annotation system with default parameters (Overbeek et al., 2014) and then downloaded from the ModelSEED database (Henry et al., 2010). This reconstruction was then thoroughly inspected following the main steps listed in (Thiele and Palsson, 2010) and refined by integrating information from PhTAC125 original genome annotation (Medigue et al., 2005) and from different functional databases, including KEGG (Kanehisa, 2002) BRENDA (Scheer et al., 2011) and MetCyc (Caspi et al., 2006). At this stage, models of related organisms (e.g. Escherichia coli, Shewanella oneidensis MR-1) were used as reference in a comparative genomics workflow for identifying potentially missing reactions. BLAST (Altschul et al., 1997) searches (adopting the Bidirectional Best Hit criterion) on the PhTAC125 genome were carried out to confirm/exclude the inclusion of further reactions to the original model. The list of the PhTAC125’s cellular transporters was obtained probing the Transporter Classification Data Base [TCDB (Saier et al., 2006)]. At the end of this iterative, manual refinement procedure, the model gained information on 60 additional genes and (about) 200 reactions were added to the reconstruction (see Table S1). Also, to properly reconcile the model with experimental data, some of the reactions initially included by the adopted automatic reconstruction method were removed from the model during manual refinement.
155 156 157 158 159 160 161
To date, no detailed information on the biomass composition of PhTAC125 is available in scientific literature, except for specific constituents [e.g. the structure of its lipo-oligosaccharide fraction (Corsaro et al., 2001), RNA and DNA composition (Medigue et al., 2005)]. Accordingly, missing information concerning the biomass assembly reaction of PhTAC125 was derived from the closely related gamma-proteobacterium Shewanella oneidensis MR-1 and E. coli, for which detailed metabolic models were already available (Feist et al., 2007) (Pinchuk et al., 2010) (Orth et al., 2011). The description of the biomass composition can be found in Table S2.
162 163 164 165
The final PhTAC125 model was named iMF721 according to the current nomenclature standard (Reed et al., 2003) and was successfully validated with the SBML Validator tool available at http://sbml.org/Facilities/Validator/index.jsp. The SBML-formatted version of iMF721 is available as Document S1.
166 167
Metabolic modeling This article is protected by copyright. All rights reserved.
The Flux Balance Analysis (FBA) method was employed to simulate flux distribution in different conditions. Briefly, FBA is a constraint-based method relying on the representation of the biochemical system under investigation in the form of stoichiometric matrix S (m×n), where m is the number of metabolites and n the number of reactions. FBA is based on the assumption of the cellular pseudo-steady state, according to which the net sum of all the production and consumption rates of each internal metabolite within a cell is considered to be zero. Under this assumption, the system can be described by the set of linear equations:
Accepted Article
168 169 170 171 172 173 174
175 176 177 178
in which Xi is the concentration of metabolite i, Sij is the stoichiometric coefficient of the ith metabolite in the jth reaction, vj is the flux of the jth reaction, N the entire set of metabolites and M the entire set of reactions.
179 180
Upper and lower bounds of flux through each reaction act as further constraints and are expressed as:
181 182 183 184
where lb and ub are the lower and upper limits for reaction j, respectively. Finally, FBA exploits linear programming to determine a feasible steady state flux vector that optimizes a given objective function (e.g. biomass production).
185 186 187 188
The reconstructed model was analysed using COBRAToolbox-2.0 (Schellenberger et al., 2011) in MATLAB® R2009b (Mathworks Inc.). Gurobi 5.6 (www.gurobi.com) and GLPK 4.32 (http://www.gnu.org/software/glpk/) solvers were used for computational simulations presented herein.
189
Proteomic data and robustness analysis
190 191 192 193
Up- and down-regulation patterns of protein expression used in this work were obtained from previously published studies (Piette et al., 2010; Piette et al., 2011) in which the proteome of cells of PhTAC125 grown at 4°C and 18°C were analysed using two-dimensional differential in-gel electrophoresis.
194 195 196 197 198 199 200 201 202 203 204
These data were used to compare (by means of robustness analysis) the experimentally determined changes in protein abundance (following a temperature downshift) with the in silico inferred impact of various reaction flux variations on the cell growth rate. Robustness analysis consists in the calculation of suboptimal cellular growth (using FBA) when the reaction flux of a given reaction is varied around the optimal value. In this context, we perturbed the flux through each reaction whose corresponding genes i) were included in the model and ii) showed a significant change in expression during proteomics experiments. Results of this analysis can be visualized as a plot of the reaction flux (x-axis) versus the cellular growth rates (y-axis). We used the ad hoc implemented function in COBRA Toolbox to perform robustness analysis. This approach overall resembles the one recently used to investigate the response of Thermoanaerobacter tengcongensis adaptation to high temperatures (Tong et al., 2013). This article is protected by copyright. All rights reserved.
Proteomic data integration and fluxes visualization
206 207 208 209
Up- and down-regulation ratios (and corresponding p-values) of protein expression were mapped onto the PhTAC125 metabolic model using MADE (Metabolic Adjustment by Differential Expression) (Jensen and Papin, 2011). The visualization of the changes in reaction fluxes in the two conditions was performed using iPath 2.0 (Yamada et al., 2011).
Accepted Article
205
210 211 212 213
Phenotype Microarray
214 215 216 217 218 219 220 221 222 223 224 225
Cells to be tested were previously plated on TYP medium agar plates and incubated overnight at 25°C. Cells were swabbed from the plates after overnight growth and suspended in appropriate medium containing Inoculating fluid (BIOLOG) and Schatz salts (NaCl 10g/l, KH2PO4 1g/l, NH4NO3 1g/l, MgSO4 7H2O 0.2g/l, FeSO4 7H2O 10mg/l, CaCl2 2H2O 10mg/l) as additive solution until the 85% transmittance suspension of cells was obtained on a Biolog turbidimeter. In order to inoculate the microplates PM1 and PM2, 1% tetrazolium violet (vol/vol) (Dye Mix A, Biolog) was added to the suspension and 100 μL of such mixture were then inoculated in each well. Plates were incubated at 15°C for 1 week (167h) with readings taken manually three times a day using a Synergy HT (Biotek) system. The cellular growth was determined on the absorbance (OD 590nm) values of the kinetic curve of dye formation. In particular OD=0.45 was selected as cut-off value for cellular growth within each well of the PM microplates as tetrazolium colour development was visually detected only for OD>0.45.
226
This article is protected by copyright. All rights reserved.
Accepted Article
227 228
Results and discussion
229
Characteristics of the reconstructed metabolic model of PhTAC125
230 231 232 233 234 235
P. haloplanktis TAC125 was firstly isolated from an Antarctic coastal seawater sample. It is considered one of the model organisms of cold-adapted bacteria (Feller and Gerday, 2003) and, given its capability to growth at reduced temperatures in respect to other model organisms (e.g. E. coli), it is being exploited as a new alternative expression host for a number of heterologous proteins. The metabolism of this versatile representative of the marine bacterioplankton was here investigated by means of metabolic network reconstruction and constraints-based modelling.
236 237 238 239 240 241 242 243 244 245 246
The metabolic network of PhTAC125 was initially obtained from its genome annotation and integrated with additional functional information as described in Experimental procedures. The reconstructed genome-scale metabolic model (named iMF721) encompasses information on 721 ORFs (20.7% of the PhTAC125 protein encoding genes), 1133 metabolites and 1322 reactions (Table 1). The model includes non-gene-associated reactions accounting for i) the biomass assembly reaction (which also takes into consideration non-growth-associated ATP costs), ii) 48 reactions which filled gaps in the metabolic network (19 added by the “AUTOCOMPLETION” function of the RAST annotation system, and 29 added during the manual evaluation of the model), iii) 85 exchange reactions allowing the simulation of external conditions (e.g. nutrients exchange) and iv) 17 spontaneous reactions. Additionally, during the gap-filling process sink and demand reactions were added to the model when necessary.
247 248 249 250 251 252 253
Importantly, our model embeds almost all (93 out of 96, 97%) of the metabolism-related protein inventory expressed by P. haloplanktis TAC125 with a complex amino acid mixture as the only carbon and nitrogen source, as assessed by proteomics analysis (Wilmes et al., 2011). Also, the model falls well within the range of currently available models of (more or less) closely related microorganisms (Oberhardt et al., 2008; Puchalka et al., 2008; Flynn et al., 2012). On this basis, we conclude that the iMF721 model should, in principle, be able to provide a comprehensive picture of the metabolic features of PhTAC125.
254
iMF721 model validation
255 256 257 258 259 260 261
PhTAC125 seems to be well adapted to grow on rich media and in all media described for this bacterium, amino acids were used as carbon and nitrogen source (Wilmes et al., 2010). Also, the capability of this strain to grow on defined medium containing amino acids as the sole carbon source (including L-leucine, L-alanine, L-aspartate and L-glutamate) has been previously shown (Giuliani, 2011). To quantitatively assess the model’s accuracy in predicting growth rates, we simulated growth phenotypes for minimal medium supplemented with these different amino acids and compared the in vivo growth data with in silico prediction.
262 263 264 265
Accordingly, an in silico minimal growth medium was defined using exchange reactions present in the model and biomass optimization was selected as the model objective function (O.F.). More in detail, lower bounds of exchange reactions accounting for all the salts present in Schatz medium (Papa et al., 2007) were set to -1000 mmol/g*h-1, in order to mimic non-limiting conditions. Each of This article is protected by copyright. All rights reserved.
the aforementioned amino acids was then chosen as the unique carbon source of this in silico medium. The PhTAC125 enzymatic capacity for these compounds was calculated as the ratio of the growth rate to the biomass yield in batch experiments (Varma and Palsson, 1994) and was set to 0.7, 3.6, 2.5 and 3.4 (mmol/g*h-1) for leucine, alanine, aspartate and glutamate, respectively. The predicted growth rates were compared to those experimentally determined for PhTAC125 (Figure 1), revealing an overall agreement between experimentally determined growth rates and in silico predictions.
273 274 275 276 277 278 279 280 281 282 283 284 285 286
Furthermore, we used Biolog Phenotype Microarray (PM) data (obtained at 15°C) to evaluate and iteratively refine the iMF721 model. This is typically achieved by (qualitatively) comparing the estimated flux value across biomass assembly reaction of the model with the activity directly measured during phenotype microarray experiment. Of the 192 carbon sources tested with PM microplates, 64 (~33%) were accounted for by the iMF721 model and thus could be used to directly test model predictions. In silico growth on these substrates was simulated by setting each of them as sole carbon source and its uptake rate to the arbitrary value of 1 mmol/g*h-1 (under aerobic conditions). Simulation results (either “growth” or “no growth”) were compared with in vivo determined phenotypes. Inconsistencies between simulation results and PM data allowed the identification of metabolic gaps in the model and/or missing transport reactions. These included, for example, the gluconate:H+ symporter (encoded by PSHAb0479), the pyruvate transporter (putatively) encoded by PSHAa0587 and the ATP:D-fructose 6-phosphotransferase (encoded by PSHAb0209); these genes were missing in the initial draft reconstruction and thus precluded the model from using some of the tested carbon sources.
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
After this iterative refinement procedure, iMF721 growth phenotypes predictions were compared again with PM results, revealing that in 84% of the cases (54 out of 64) the outcomes of in silico simulations correctly matched growth phenotypes assessed by in vivo experiments (Table 2, Table S3). Discrepancies in experimental and in silico growth phenotypes may be due to several factors, including incomplete/incorrect homology-based gene annotation or regulatory mechanisms not currently accounted for by the iMF721 reconstruction. However, all the observed incongruences derived from the fact that the model and PM outcomes were “growth” and “no growth”, respectively. In these cases, the wrong incorporation of transport reaction(s) in the iMF721 model is the most likely explanation. It must be noted that, for at least two of the tested compounds (i.e. leucine and mannose), the model predictions are probably correct since it has been previously reported that PhTAC125 is able to utilize leucine and mannose as single carbon sources (Papa et al., 2006; Giuliani, 2011). However, it is also possible that such reactions were not detectable within the incubation time (167h) at 15°C due to the slow metabolism of PhTAC125 in presence of leucine as single c-source (Figure 1). Moreover the cut-off value considered in this work (“growth” if OD>0.450) could have led to an underestimation of the overall agreement between Biolog outcomes and in silico predictions as, for example, malate or succinate which resulted as “no growth” because of their OD values
Accepted Article
1 2
Genome scale metabolic reconstruction and constraints-based modelling of the Antarctic bacterium Pseudoalteromonas haloplanktis TAC125. 1
3 4 5 6 7 8
Marco Fondi1,2§, Isabel Maida1, Elena Perrin1, Alessandra Mellera1,2, Stefano Mocali3, Ermenegilda Parrilli4, Maria Luisa Tutino4, Pietro Liò5, Renato Fani1,2
9
10
1
11
Piano 6, 50019 Sesto Fiorentino (Firenze)
12
2
13
Fiorentino (Firenze)
14
3
15
(CRA-ABP), Piazza d’Azeglio 30, 50121 Firenze, Italy
16
4
17
Cintia, I-80126, Naples, Italy
18
5
19
Kingdom§Corresponding author
20
Dr. Marco Fondi
21
Laboratory of Microbial and Molecular Evolution, Department of Biology, University of Florence, Via Madonna del
22
Piano 6, 50019 Sesto Fiorentino (Firenze)
23
Tel +390554574736
24
Email: [email protected]
Laboratory of Microbial and Molecular Evolution, Department of Biology, University of Florence, Via Madonna del
ComBo, Florence Computational Biology group, University of Florence, Via Madonna del Piano 6, 50019 Sesto
Consiglio per la Ricerca e la Sperimentazione in Agricoltura, Centro di Ricerca per l’Agrobiologia e la Pedologia
Department of Chemical Sciences, University of Naples Federico II, Complesso Universitario M. S. Angelo, Via
Computer Laboratory, Cambridge University, William Gates Building 15, JJ Thomson Avenue, Cambridge, United
This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/1462-2920.12513 This article is protected by copyright. All rights reserved.
Summary
27 28 29
The Antarctic strain Pseudoalteromonas haloplanktis TAC125 is one of the model organisms of cold-adapted bacteria and is currently exploited as a new alternative expression host for numerous biotechnological applications.
30 31 32 33 34 35 36 37 38
Here, we investigated several metabolic features of this strain through in silico modelling and functional integration of –omics data. A genome-scale metabolic model of P. haloplanktis TAC125 was reconstructed, encompassing information on 721 genes, 1133 metabolites and 1322 reactions. The predictive potential of this model was validated against a set of experimentally determined growth rates and a large dataset of growth phenotypic data. Furthermore, evidence synthesis from proteomics, phenomics, physiology and metabolic modeling data revealed possible drawbacks of cold-dependent changes in gene expression on the overall metabolic network of P. haloplanktis TAC125. These included, for example, variations in its central metabolism, amino acids degradation and fatty acids biosynthesis.
39 40 41 42 43
The genome scale metabolic model described here is the first one reconstructed so far for an Antarctic microbial strain. It allowed a system-level investigation of variations in cellular metabolic fluxes following a temperature downshift. It represents a valuable platform for further investigations on P. haloplanktis TAC125 cellular functional states and for the design of more focused strategies for its possible biotechnological exploitation.
Accepted Article
25 26
44
This article is protected by copyright. All rights reserved.
Accepted Article
45 46
Introduction
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Antarctica is one of the most extreme environments on Earth due to the presence of stably low temperatures. Many (psychrophilic) microorganisms inhabit this environment, being responsible for most of the biomass of this ecological niche and playing a key role in maintaining proper ecosystem functions (Wilkins et al., 2013). Such cold-adapted microorganisms have begun to attract the attention of the scientific community from both fundamental and application viewpoints (Cavicchioli et al., 2002; D'Amico et al., 2006). Indeed, psychrophiles have evolved peculiar features (Feller and Gerday, 2003; Giordano et al., 2012) to face the profound impact of cold shock on several cellular processes (e.g. growth rates, ribosomal synthesis, cytoplasmic membrane composition, etc.) (Phadtare, 2004). Overall, survival at low temperature depends on the ability of microbes to sense changes in temperature and to transduce these signals to the genome, ultimately bringing changes in the regulation of cold-adaptation related genes (Shivaji and Prakash, 2010). These include: i) genes for fatty acid desaturases (Sato and Murata, 1980), ii) genes involved in replication (Jones et al., 1992; Graumann and Marahiel, 1999), transcription (Sledjeski et al., 1996; Inaba et al., 2003) and translation (Jones et al., 1987; Xia et al., 2003), iii) genes encoding cold shock proteins (Brandi et al., 1994; Atlung and Ingmer, 1997), as well as iv) many genes coding for still uncharacterized proteins (Kawamoto et al., 2007; Ting et al., 2010). Moreover, the increase of enzymes concentration (Willem et al., 1999) and the expression of specific cold-adapted enzymatic isoforms (Hoyoux et al., 2004) represent two additional examples of how psychrophilic enzymatic machineries may counteract the growth limiting effects of cold temperatures.
66 67 68 69 70 71 72 73 74 75 76 77
Recently, the spreading of –omics technologies has allowed the system-level investigation of the mechanisms involved in microbial cold adaptation (Seo et al., 2004; Goodchild et al., 2005; Bakermans et al., 2007; Zheng et al., 2007; Piette et al., 2010). For example, a transcriptomics analysis on Shewanella oneidensis MR-1 showed that more than 70% of its genes involved in energy metabolism were down-regulated following a temperature downshift (Gao et al., 2006). A proteomic study on Sphingopyxis alaskensis revealed that a large fraction of genes involved in metabolic processes (e.g. energy production/conversion, carbohydrate, amino acids, nucleotides and cofactors transport/metabolism) are less abundant at lower temperatures (Ting et al., 2010). Finally, from a metabolic viewpoint, a large-scale metabolomics analysis on the cold adaptation of Mesorhizobium sp. strain N33 revealed a key role of (unsaturated) fatty acids biosynthetic process. From such examples, it emerges that, despite some strategies may be shared, different microorganisms use different strategies to cope with cold environments (Kawamoto et al., 2007).
78 79 80 81 82 83 84 85 86
The Antarctic marine bacterium Pseudoalteromonas haloplanktis TAC125 (PhTAC125) has been isolated from sea water sampled along the Antarctic ice-shell, a permanently cold environment. PhTAC125 is capable of growing in a wide temperature range (4–25°C) and its lowest observed doubling time was detected at 20°C (Medigue et al., 2005). Several exceptional genomic and metabolic features were derived from the genome sequence of this bacterium, showing adaptation to periodic situations of nutrient abundance (Medigue et al., 2005). Indeed PhTAC125 is considered to be one of the model organisms of cold-adapted bacteria and has been suggested as an alternative host for the soluble overproduction of heterologous proteins, given its capability to grow fast at low temperatures (Duilio et al., 2004; Wilmes et al., 2010; Rippa et al., 2012; Corchero et al., 2013). This article is protected by copyright. All rights reserved.
Furthermore, bacteria belonging to the genus Pseudoalteromonas are known to possess an inhibitory activity against human pathogens belonging to the Burkholdeia cepacia complex (Bcc) and to be able to produce antibiofilm molecules (Papa et al., 2013; Papaleo et al., 2013), revealing the interesting biotechnological potential and metabolic biodiversity of these microorganisms. The increasing interest in PhTAC125 has led to the accumulation of different data types for this bacterium in the last few years, including its complete genome sequence (Medigue et al., 2005), its proteome (Piette et al., 2010; Piette et al., 2011) and detailed growth phenotypes (Wilmes et al., 2010; Giuliani, 2011). Thus, it is now possible to integrate such different data sources and perform a system-level investigation of PhTAC125 metabolism.
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
Genome scale metabolic modeling represents a valuable tool in this context. Indeed, this in silico approach can be adopted to quantitatively simulate chemical reactions fluxes within the cell, including metabolic adjustments in response to external perturbations (e.g. temperature downshift). Genome annotations are usually transformed into models by defining the boundaries of the system, a biomass assembly reaction, and exchange fluxes with the environment (Durot et al., 2009; Thiele and Palsson, 2010). Constraint-based modelling methods (e.g. Flux Balance Analysis, FBA) can then be used to compute the resulting balance of all the active cellular reactions in the cell and to simulate the maximal growth of a cell in a given environmental condition (Varma and Palsson, 1994; Schilling et al., 2000). In the past decade, this approach has been successfully applied for studying large-scale metabolic networks in microbes, with the aim of guiding rational engineering of biological systems for applications in industrial and medical biotechnology (Milne et al., 2009). Interestingly, different data types can be mapped onto metabolic models in order to elucidate more thoroughly the metabolism of a cell and its response to environmental factors. This is usually done by including functional characterization and accurate quantification of all the main cellular information levels of gene products, mRNA, proteins and metabolites, as well as their interaction (Zhang et al., 2010). In recent years, –omics-derived data have been used to refine, validate and/or integrate metabolic models, including transcriptomics (Colijn et al., 2009; Jensen and Papin, 2011), proteomics (Gille et al., 2010), fluxomics (Chen et al., 2011; Feng and Zhao, 2013) and phenomics (Fang et al., 2011).
115 116 117 118 119 120 121 122 123
In at least two study-cases, FBA and expression data have been merged to explore the effect of temperature downshift at the system level (Navid and Almaas, 2012; Tong et al., 2013). In particular, Tong et al. (2013) performed robustness analysis on the core metabolic model of Thermoanaerobacter tengcongensis to study the dynamic changes of the metabolic network following the perturbation of the culture temperature and collecting the bacterial growth rates and differential proteomes. Given the overall agreement between in silico simulations and observed phenotypes, this approach was shown to provide a reliable platform to systematically evaluate the mechanisms of bacterial metabolism and relevant switches in the presence of a temperature downshift.
124 125 126 127 128 129
In this study, a genome-scale reconstruction of P. haloplanktis TAC125’s metabolism was performed based on its genome annotation. The predictive capability of the model [named iMF721 according to the current naming convention (Reed et al., 2003)] was successfully validated comparing constraints-based modeling outcomes with experimentally determined growth rates and large scale growth phenotype data (Phenotype Microarray). The iMF721 model was then used to globally investigate possible metabolic adjustments of P. haloplanktis TAC125 during growth at
Accepted Article
87 88 89 90 91 92 93 94 95
This article is protected by copyright. All rights reserved.
low temperature by means of robustness analysis and functional integration of protein abundance data into the reconstructed network.
132 133 134
Finally, the metabolic reconstruction reported herein represents a reliable platform for the future design of experimental strategies aimed at the exploitation of the biotechnological potential of PhTAC125.
Accepted Article
130 131
135 136 137
Experimental procedures
138
Metabolic model reconstruction
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
An initial draft metabolic reconstruction of the strain PhTAC125 was constructed using RAST annotation system with default parameters (Overbeek et al., 2014) and then downloaded from the ModelSEED database (Henry et al., 2010). This reconstruction was then thoroughly inspected following the main steps listed in (Thiele and Palsson, 2010) and refined by integrating information from PhTAC125 original genome annotation (Medigue et al., 2005) and from different functional databases, including KEGG (Kanehisa, 2002) BRENDA (Scheer et al., 2011) and MetCyc (Caspi et al., 2006). At this stage, models of related organisms (e.g. Escherichia coli, Shewanella oneidensis MR-1) were used as reference in a comparative genomics workflow for identifying potentially missing reactions. BLAST (Altschul et al., 1997) searches (adopting the Bidirectional Best Hit criterion) on the PhTAC125 genome were carried out to confirm/exclude the inclusion of further reactions to the original model. The list of the PhTAC125’s cellular transporters was obtained probing the Transporter Classification Data Base [TCDB (Saier et al., 2006)]. At the end of this iterative, manual refinement procedure, the model gained information on 60 additional genes and (about) 200 reactions were added to the reconstruction (see Table S1). Also, to properly reconcile the model with experimental data, some of the reactions initially included by the adopted automatic reconstruction method were removed from the model during manual refinement.
155 156 157 158 159 160 161
To date, no detailed information on the biomass composition of PhTAC125 is available in scientific literature, except for specific constituents [e.g. the structure of its lipo-oligosaccharide fraction (Corsaro et al., 2001), RNA and DNA composition (Medigue et al., 2005)]. Accordingly, missing information concerning the biomass assembly reaction of PhTAC125 was derived from the closely related gamma-proteobacterium Shewanella oneidensis MR-1 and E. coli, for which detailed metabolic models were already available (Feist et al., 2007) (Pinchuk et al., 2010) (Orth et al., 2011). The description of the biomass composition can be found in Table S2.
162 163 164 165
The final PhTAC125 model was named iMF721 according to the current nomenclature standard (Reed et al., 2003) and was successfully validated with the SBML Validator tool available at http://sbml.org/Facilities/Validator/index.jsp. The SBML-formatted version of iMF721 is available as Document S1.
166 167
Metabolic modeling This article is protected by copyright. All rights reserved.
The Flux Balance Analysis (FBA) method was employed to simulate flux distribution in different conditions. Briefly, FBA is a constraint-based method relying on the representation of the biochemical system under investigation in the form of stoichiometric matrix S (m×n), where m is the number of metabolites and n the number of reactions. FBA is based on the assumption of the cellular pseudo-steady state, according to which the net sum of all the production and consumption rates of each internal metabolite within a cell is considered to be zero. Under this assumption, the system can be described by the set of linear equations:
Accepted Article
168 169 170 171 172 173 174
175 176 177 178
in which Xi is the concentration of metabolite i, Sij is the stoichiometric coefficient of the ith metabolite in the jth reaction, vj is the flux of the jth reaction, N the entire set of metabolites and M the entire set of reactions.
179 180
Upper and lower bounds of flux through each reaction act as further constraints and are expressed as:
181 182 183 184
where lb and ub are the lower and upper limits for reaction j, respectively. Finally, FBA exploits linear programming to determine a feasible steady state flux vector that optimizes a given objective function (e.g. biomass production).
185 186 187 188
The reconstructed model was analysed using COBRAToolbox-2.0 (Schellenberger et al., 2011) in MATLAB® R2009b (Mathworks Inc.). Gurobi 5.6 (www.gurobi.com) and GLPK 4.32 (http://www.gnu.org/software/glpk/) solvers were used for computational simulations presented herein.
189
Proteomic data and robustness analysis
190 191 192 193
Up- and down-regulation patterns of protein expression used in this work were obtained from previously published studies (Piette et al., 2010; Piette et al., 2011) in which the proteome of cells of PhTAC125 grown at 4°C and 18°C were analysed using two-dimensional differential in-gel electrophoresis.
194 195 196 197 198 199 200 201 202 203 204
These data were used to compare (by means of robustness analysis) the experimentally determined changes in protein abundance (following a temperature downshift) with the in silico inferred impact of various reaction flux variations on the cell growth rate. Robustness analysis consists in the calculation of suboptimal cellular growth (using FBA) when the reaction flux of a given reaction is varied around the optimal value. In this context, we perturbed the flux through each reaction whose corresponding genes i) were included in the model and ii) showed a significant change in expression during proteomics experiments. Results of this analysis can be visualized as a plot of the reaction flux (x-axis) versus the cellular growth rates (y-axis). We used the ad hoc implemented function in COBRA Toolbox to perform robustness analysis. This approach overall resembles the one recently used to investigate the response of Thermoanaerobacter tengcongensis adaptation to high temperatures (Tong et al., 2013). This article is protected by copyright. All rights reserved.
Proteomic data integration and fluxes visualization
206 207 208 209
Up- and down-regulation ratios (and corresponding p-values) of protein expression were mapped onto the PhTAC125 metabolic model using MADE (Metabolic Adjustment by Differential Expression) (Jensen and Papin, 2011). The visualization of the changes in reaction fluxes in the two conditions was performed using iPath 2.0 (Yamada et al., 2011).
Accepted Article
205
210 211 212 213
Phenotype Microarray
214 215 216 217 218 219 220 221 222 223 224 225
Cells to be tested were previously plated on TYP medium agar plates and incubated overnight at 25°C. Cells were swabbed from the plates after overnight growth and suspended in appropriate medium containing Inoculating fluid (BIOLOG) and Schatz salts (NaCl 10g/l, KH2PO4 1g/l, NH4NO3 1g/l, MgSO4 7H2O 0.2g/l, FeSO4 7H2O 10mg/l, CaCl2 2H2O 10mg/l) as additive solution until the 85% transmittance suspension of cells was obtained on a Biolog turbidimeter. In order to inoculate the microplates PM1 and PM2, 1% tetrazolium violet (vol/vol) (Dye Mix A, Biolog) was added to the suspension and 100 μL of such mixture were then inoculated in each well. Plates were incubated at 15°C for 1 week (167h) with readings taken manually three times a day using a Synergy HT (Biotek) system. The cellular growth was determined on the absorbance (OD 590nm) values of the kinetic curve of dye formation. In particular OD=0.45 was selected as cut-off value for cellular growth within each well of the PM microplates as tetrazolium colour development was visually detected only for OD>0.45.
226
This article is protected by copyright. All rights reserved.
Accepted Article
227 228
Results and discussion
229
Characteristics of the reconstructed metabolic model of PhTAC125
230 231 232 233 234 235
P. haloplanktis TAC125 was firstly isolated from an Antarctic coastal seawater sample. It is considered one of the model organisms of cold-adapted bacteria (Feller and Gerday, 2003) and, given its capability to growth at reduced temperatures in respect to other model organisms (e.g. E. coli), it is being exploited as a new alternative expression host for a number of heterologous proteins. The metabolism of this versatile representative of the marine bacterioplankton was here investigated by means of metabolic network reconstruction and constraints-based modelling.
236 237 238 239 240 241 242 243 244 245 246
The metabolic network of PhTAC125 was initially obtained from its genome annotation and integrated with additional functional information as described in Experimental procedures. The reconstructed genome-scale metabolic model (named iMF721) encompasses information on 721 ORFs (20.7% of the PhTAC125 protein encoding genes), 1133 metabolites and 1322 reactions (Table 1). The model includes non-gene-associated reactions accounting for i) the biomass assembly reaction (which also takes into consideration non-growth-associated ATP costs), ii) 48 reactions which filled gaps in the metabolic network (19 added by the “AUTOCOMPLETION” function of the RAST annotation system, and 29 added during the manual evaluation of the model), iii) 85 exchange reactions allowing the simulation of external conditions (e.g. nutrients exchange) and iv) 17 spontaneous reactions. Additionally, during the gap-filling process sink and demand reactions were added to the model when necessary.
247 248 249 250 251 252 253
Importantly, our model embeds almost all (93 out of 96, 97%) of the metabolism-related protein inventory expressed by P. haloplanktis TAC125 with a complex amino acid mixture as the only carbon and nitrogen source, as assessed by proteomics analysis (Wilmes et al., 2011). Also, the model falls well within the range of currently available models of (more or less) closely related microorganisms (Oberhardt et al., 2008; Puchalka et al., 2008; Flynn et al., 2012). On this basis, we conclude that the iMF721 model should, in principle, be able to provide a comprehensive picture of the metabolic features of PhTAC125.
254
iMF721 model validation
255 256 257 258 259 260 261
PhTAC125 seems to be well adapted to grow on rich media and in all media described for this bacterium, amino acids were used as carbon and nitrogen source (Wilmes et al., 2010). Also, the capability of this strain to grow on defined medium containing amino acids as the sole carbon source (including L-leucine, L-alanine, L-aspartate and L-glutamate) has been previously shown (Giuliani, 2011). To quantitatively assess the model’s accuracy in predicting growth rates, we simulated growth phenotypes for minimal medium supplemented with these different amino acids and compared the in vivo growth data with in silico prediction.
262 263 264 265
Accordingly, an in silico minimal growth medium was defined using exchange reactions present in the model and biomass optimization was selected as the model objective function (O.F.). More in detail, lower bounds of exchange reactions accounting for all the salts present in Schatz medium (Papa et al., 2007) were set to -1000 mmol/g*h-1, in order to mimic non-limiting conditions. Each of This article is protected by copyright. All rights reserved.
the aforementioned amino acids was then chosen as the unique carbon source of this in silico medium. The PhTAC125 enzymatic capacity for these compounds was calculated as the ratio of the growth rate to the biomass yield in batch experiments (Varma and Palsson, 1994) and was set to 0.7, 3.6, 2.5 and 3.4 (mmol/g*h-1) for leucine, alanine, aspartate and glutamate, respectively. The predicted growth rates were compared to those experimentally determined for PhTAC125 (Figure 1), revealing an overall agreement between experimentally determined growth rates and in silico predictions.
273 274 275 276 277 278 279 280 281 282 283 284 285 286
Furthermore, we used Biolog Phenotype Microarray (PM) data (obtained at 15°C) to evaluate and iteratively refine the iMF721 model. This is typically achieved by (qualitatively) comparing the estimated flux value across biomass assembly reaction of the model with the activity directly measured during phenotype microarray experiment. Of the 192 carbon sources tested with PM microplates, 64 (~33%) were accounted for by the iMF721 model and thus could be used to directly test model predictions. In silico growth on these substrates was simulated by setting each of them as sole carbon source and its uptake rate to the arbitrary value of 1 mmol/g*h-1 (under aerobic conditions). Simulation results (either “growth” or “no growth”) were compared with in vivo determined phenotypes. Inconsistencies between simulation results and PM data allowed the identification of metabolic gaps in the model and/or missing transport reactions. These included, for example, the gluconate:H+ symporter (encoded by PSHAb0479), the pyruvate transporter (putatively) encoded by PSHAa0587 and the ATP:D-fructose 6-phosphotransferase (encoded by PSHAb0209); these genes were missing in the initial draft reconstruction and thus precluded the model from using some of the tested carbon sources.
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
After this iterative refinement procedure, iMF721 growth phenotypes predictions were compared again with PM results, revealing that in 84% of the cases (54 out of 64) the outcomes of in silico simulations correctly matched growth phenotypes assessed by in vivo experiments (Table 2, Table S3). Discrepancies in experimental and in silico growth phenotypes may be due to several factors, including incomplete/incorrect homology-based gene annotation or regulatory mechanisms not currently accounted for by the iMF721 reconstruction. However, all the observed incongruences derived from the fact that the model and PM outcomes were “growth” and “no growth”, respectively. In these cases, the wrong incorporation of transport reaction(s) in the iMF721 model is the most likely explanation. It must be noted that, for at least two of the tested compounds (i.e. leucine and mannose), the model predictions are probably correct since it has been previously reported that PhTAC125 is able to utilize leucine and mannose as single carbon sources (Papa et al., 2006; Giuliani, 2011). However, it is also possible that such reactions were not detectable within the incubation time (167h) at 15°C due to the slow metabolism of PhTAC125 in presence of leucine as single c-source (Figure 1). Moreover the cut-off value considered in this work (“growth” if OD>0.450) could have led to an underestimation of the overall agreement between Biolog outcomes and in silico predictions as, for example, malate or succinate which resulted as “no growth” because of their OD values
