Insights into metabolic osmoadaptation of the ectoines-producer bacterium Chromohalobacter salexigens through a high-quality genome scale metabolic model

Piubeli et al. Microb Cell Fact (2018) 17:2 https://doi.org/10.1186/s12934-017-0852-0 RESEARCH Insights into metabolic osmoadaptation of the ectoines-producer bacterium Chromohalobacter salexigens through a high-quality genome scale metabolic model Francine Piubeli1, Manuel Salvador1, Montserrat Argandoña1, Joaquín J. Nieto1, Vicente Bernal2,3, Jose M. Pastor2, Manuel Cánovas2 and Carmen Vargas1* Abstract Background: The halophilic bacterium Chromohalobacter salexigens is a natural producer of ectoines, compatible solutes with current and potential biotechnological applications. As production of ectoines is an osmoregulated process that draws away TCA intermediates, bacterial metabolism needs to be adapted to cope with salinity changes. To explore and use C. salexigens as cell factory for ectoine(s) production, a comprehensive knowledge at the systems level of its metabolism is essential. For this purpose, the construction of a robust and high-quality genome-based metabolic model of C. salexigens was approached. Results: We generated and validated a high quality genome-based C. salexigens metabolic model (iFP764). This comprised an exhaustive reconstruction process based on experimental information, analysis of genome sequence, manual re-annotation of metabolic genes, and in-depth refinement. The model included three compartments (periplasmic, cytoplasmic and external medium), and two salinity-specific biomass compositions, partially based on experimental results from C. salexigens. Using previous metabolic data as constraints, the metabolic model allowed us to simulate and analyse the metabolic osmoadaptation of C. salexigens under conditions for low and high production of ectoines. The iFP764 model was able to reproduce the major metabolic features of C. salexigens. Flux Balance Analysis (FBA) and Monte Carlo Random sampling analysis showed salinity-specific essential metabolic genes and different distribution of fluxes and variation in the patterns of correlation of reaction sets belonging to central C and N metabolism, in response to salinity. Some of them were related to bioenergetics or production of reducing equivalents, and probably related to demand for ectoines. Ectoines metabolic reactions were distributed according to its correlation in four modules. Interestingly, the four modules were independent both at low and high salinity conditions, as they did not correlate to each other, and they were not correlated with other subsystems. Conclusions: Our validated model is one of the most complete curated networks of halophilic bacteria. It is a powerful tool to simulate and explore C. salexigens metabolism at low and high salinity conditions, driving to low and high production of ectoines. In addition, it can be useful to optimize the metabolism of other halophilic bacteria for metabolite production. Keywords: Genome–scale metabolic model, Flux balance analysis, Chromohalobacter salexigens, Metabolic osmoadaptation © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/ publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Open Access Microbial Cell Factories *Correspondence: [email protected] 1 Department of Microbiology and Parasitology, Faculty of Pharmacy, University of Sevilla, C/Profesor García González 2, 41012 Sevilla, Spain Full list of author information is available at the end of the article Page 2 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 Background Chromohalobacter salexigens is a moderately halophilic gamma-proteobacterium of the family Halomonadaceae, phylogenetically related to the genus Halomonas [1]. It displays one of the broadest known growth salinity ranges [2], and is considered as a biological model to study prokaryotic osmoadaptation. C. salexigens has been proposed, among other natural producers, as a cell factory for the production of the biostabilizing compatible solutes ectoine and hydroxyectoine [3, 4]. The microorganism has several advantages as a cell factory: (i) it is easily cultured by using conventional media, (ii) it has broad metabolic versatility being capable of using different carbon sources [1], (iii) a high diversity of genetic tools are available for its study [5], and (iv) its genome has been sequenced [6]. To further explore and use C. salexigens as cell factory for ectoine(s), a comprehensive understanding of its metabolism is essential. In C. salexigens, production of ectoines is an osmoregulated process that burdens central metabolic routes by quantitatively drawing TCA cycle intermediaries. The ectoine biosynthetic pathway starts with the phosphorylation of L-aspartate and consumes oxaloacetate (OAA) and acetyl-CoA, the latter of which is produced by oxidative decarboxylation of pyruvate by pyruvate dehydrogenase. This use of intermediaries of the TCA cycle is of special relevance, since active anaplerotic pathways are needed in order to replenish the cycle [7]. Consequently, central metabolism is adapted to the requirements of this biosynthetic route. Genome-scale metabolic models (GEMs) are currently the only approach that enables the modeling and global analysis of an organism’s metabolic and transport network by a global analysis [8]. A genome-wide constraint-based model consists of a stoichiometric reconstruction of all reactions known to act in the metabolism of the organism, along with an accompanying set of constraints on the fluxes of each reaction in the system [9]. A major advantage of this approach is that the model does not require knowledge on the kinetics of the reactions. These models define the organism’s global metabolic space, network structural properties, potential flux distribution, and provide a framework with which to navigate through the metabolic wiring of the cell [10, 11]. Through various analysis techniques, constraint-based models can help to predict cellular phenotypes, given particular environmental conditions. Constraint-based analysis techniques (i.e. flux-balance analysis [FBA]), have been instrumental in elucidating metabolic features in a variety of organisms [9] and so far have been used for certain biotechnology endeavors [12, 13]. Here, we describe a highly detailed and curated metabolic reconstruction of Chromohalobacter salexigens DSM 3043, one of the mayor ectoine(s) natural producers. The reconstruction was built using the COBRA approach [10, 14] and validated using flux balance analysis [9]. The in silico metabolic network was further evaluated by comparing predicted growth rate capacities in different carbon sources. The model was used to gain insight into the physiology of C. salexigens metabolic stress response at conditions of low and high ectoines production. In addition to optimize C. salexigens strains for ectoines production, iFP764 may become an important tool to understand the metabolism of halophilic microorganisms, and may serve as a platform to test metabolic capabilities of halophilic bacteria at low and high salinity. Methods Metabolic reconstruction and refinement The metabolic reconstruction of C. salexigens was generated by using a bottom-up approach. This was based on experimental information reported in previous studies, analysis of genome sequence and manual re-annotation of metabolic genes. An initial core model was built up based on previous studies on metabolic osmoadaptation and accumulation of compatible solutes by C. salexigens, as well as a thorough literature revision. For this purpose, genes and reactions for the synthesis and degradation of compatible solutes and their connection with the central metabolism, as well as the main C and N pathways, were compiled. Sequentially, this core was expanded with all necessary routes for cell growth. Additionally, gap filling was performed to connect routes and to complete the model. We made sure that main pathways involved in salinity dependent-metabolic adaptations were included in the reconstruction. Thus proteomic and transcriptomic data from experiments performed at low and high salinity (0.6M and 2.5M of NaCl respectively) (unpublished results), and previous metabolic data from our laboratory [7], were considered. All pathways, reactions, metabolites, and genes included in the reconstruction were manually inspected. Orthologous sequences of all enzymes of interest were identified and compared using BLAST [15], analyses of domains (Conserved Domain Database [16], and SMART [17]) and Genomic context (STRING [18], and ARCHAEA [19]). Moreover, phylogenetic analyses (ClustalW [20] and MEGA 6 [21] were carried out to correctly assign different isoenzymes to their specific reactions. Reaction directionalities were inspected, and appropriate changes were made based on the BRENDA database [22]. Reactions without any gene association were added based on evidence from the literature, presence in databases, or the need to fill functional gaps to achieve a functional network (as in the case of putative transporters). Exchange reactions were also included, as they are reactions that represent the supply to, or removal of, metabolites from the extra-organism “space” [23]. The reaction rate through the non-growth- associated ATP maintenance reaction (ATPM) was set Page 3 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 to 7.6mmolg−1 (dry weight)h−1 to account for upkeep. Growth-associated ATP maintenance was included in the two-biomass equations used. Since both equations differ in overall macromolecular composition, the energy cost for these macromolecules synthesis differed concordantly. Thus, we used the Unknown Growth-asso- ciated ATP maintenance of 36.94mmolg−1 calculated by Feist etal. [28] and followed their procedure to calculate the use of ATP per mmol of protein, DNA and RNA polymerization at low and high salinity (23.04 and 8.79mmolg−1 respectively). In thatway, theGrowthassociatedATP maintenance wasset to 45.73mmolg−1 (dry weight) at high salinity and 59.98 mmolg−1 (dry weight) at low salinity. The molecular formulae and charges of the metabolites in the model were determined assuming a pH of 7.2. The structure and the name of reaction and metabolites were assigned by following an established protocol [23]. All data bases and bioinformatics programs used are listed in Table1. Formulation of biomass reaction Given the metabolic adaptation to the salinity of C. salexigens [7], the biomass composition defined in the model was adapted to simulate high and low salinity conditions. For that purpose, two different biomass reaction compositions (BIO) were estimated: “BIO_H” to simulate high salinity (2.5M NaCl), and “BIO_L” to simulate low salinity (0.6M NaCl). According to Thiele and Palsson [23], the detailed biomass composition needs to be experimentally determined. However, in some occasions, it may not be possible to obtain a detailed biomass composition of the organism of interest. In this case, the main alternative strategy is the utilization of the experimentally determined biomass composition of a non-related microorganism (i.e. Escherichia coli) [23]. A second, more accurate, strategy is to partially reformulate an experimentally-determined biomass composition by estimating the relative fraction of the biomass metabolites of the organism of interest that can be experimentally calculated. In this work, both biomass compositions were obtained by partially reformulating Escherichia coli biomass composition [25–27] by using previously published and experimental data from C. salexigens [7, 24]. Appropriate amino acid molar ratio for both biomass reactions were calculated assuming that all proteins encoded by C. salexigens genome were expressed equally [26, 28]. The proportion of total proteins in response to salinity was calculated based on previous results by our group [7] and the molar ratio of membrane phospholipids in response to salinity were previously determined by Vargas and co-workers [24]. The percentage of G+C in C. salexigens genome was experimentally estimated [1]. This was utilized to calculate the genome nucleotide composition, which was used to establish the deoxyribonucleotide molar ratio. The ribonucleotide composition was deduced from the composition of C. salexigens rRNA and tRNA in the proportions defined previously [26]. Generally, the relative composition of the nucleic acids does not change with osmolarity, as the overall DNA and RNA content related to dry weight does not seem to be affected by osmotic stress [27]. This assertion was assumed when establishing the relative composition of DNA and RNA. The proportions of the other compounds were assumed not to vary between E. coli and C. salexigens, and these coefficients were derived from values reported in the literature [25]. Ectoine and hydroxyectoine were included in the biomass composition, because they are growthassociated metabolites and are accumulated in the cells in response to salinity [3]. The coefficients for these compounds were determined based on previous results [7]. Units of each component of a biomass reaction were mmolgDW−1 (milimoles per gram cell dry weight). Flux through a biomass reaction has h−1 units and is equivalent to the exponential growth rate of the organism [23]. Detailed information on the biomass composition calculation is included in Additional file1: Tables S1 and S2. Conversion of the reconstruction into a mathematical model The reactions and their participating metabolites in the metabolic network were connected via the Table 1 Database sources used for C. salexigens iFP764 metabolic reconstruction Name Direction Archaea http://archaea.u-psud.fr/absynte/ BioCyc http://biocyc.org/ BioMet ToolBox 2.0 http://biomet-toolbox.org/index. php?page=about Blast http://blast.ncbi.nlm.nih.gov/Blast.cgi Brenda http://www.brenda-enzymes.org/ ClustalW http://www.ebi.ac.uk/Tools/msa/ clustalw2/ Conserved Domain Database http://www.ncbi.nlm.nih.gov/cdd EcoCyc http://ecocyc.org/ GeneOntology http://geneontology.org/ Kegg http://www.genome.jp/kegg/ Mega http://www.megasoftware.net/ MicrobesOnline http://www.microbesonline.org/ NCBI http://www.ncbi.nlm.nih.gov/ Pathway-tools http://brg.ai.sri.com/ptools/ Smart http://smart.embl-heidelberg.de/ String http://string-db.org TCDB http://www.tcdb.org/ TransportDB http://www.membranetransport.org/ index.html Page 4 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 stoichiometric matrix (S), which contains the stoichiometric coefficients for each metabolite (row) in each of the reactions of the metabolic model (column). The values in this matrix correspond to the stoichiometry of the metabolite in the reaction. A negative number represents consumption and a positive number represents production of the metabolite. This matrix is the base of the metabolic model that considers metabolite mass conservation and that the accumulation of metabolites (b) in the system is equivalent to the product of S matrix and the vector of reaction fluxes (v): S·v=b. When using FBA to solve the differential equations system, one of the requirement is that the organism is in steady-state growth so, b=0, resulting in a linear system of equations [10]. Then the system is bound to a large solution space which needs to be constrained by applying upper and lower bounds to each individual reaction flux (vi): vi lower≤vi≤vi upper. Thus, in this process, system boundaries were defined, converting the general reconstruction of C. salexigens into a specific model. The initial model differed in scope and boundaries from the final model, which was accomplished after multiple iterations of validation and refinement, and which was used to simulate phenotypic behaviour in a prospective manner. The upper and lower bounds were defined for each reaction based on the flux they could carry. Reversible reactions have an upper bound value of 1000 mmolgDW−1h−1 and a lower bound value of −1000 mmolgDW−1h−1, making them practically unconstrained, while irreversible reaction have a lower bound value of zero. By default, the biomass reaction was set as the objective to be maximized. All exchange reactions have zero lower bound value except for the carbon source (glucose, −10 mmolgDW−1h−1), and oxygen (−20 mmolgDW−1h−1), as simulations were run considering growth on a glucose mineral medium in aerobic conditions. Upper and lower bounds of exchange reactions for inorganic ions required by biomass reaction were settled as 1000 and −1000 mmolgDW−1h−1, respectively. The default value for lower bound of glucose uptake was based on the typical glucose uptake rates [25]. Biomass reaction was set as objective function and was maximized to calculate the feasible flux distribution. The obtained value for this reaction was equal to maximum cell growth rate [29]. The final reconstruction metabolic network was named iFP764 according to the current standard nomenclature [30]. The complete iFP764 model is included in the Additional file2. All model simulations were performed using COBRA Toolbox1.8 [31] in MATLAB (R2014b) (The MathWorks Inc., Natick, MA), linked to the GLPK or Gurobi linear program solver. Computational network analysis Minimal medium simulation The minimal medium M63 [32] was simulated in silico by allowing entrance of some simple salts and ions, water, and O2 (aerobic simulations) in the system. For that purpose, upper and lower bounds of exchange reactions for those metabolites exchange reactions were settled as 1000 and −1000 mmolgDW−1h−1, respectively. The composition of the M63 medium is described in Additional file1: Table S3. The carbon and nitrogen source was changed according to the objective of simulation as well as the secretion rate of pyruvate and acetate (Additional file1: Table S4). Gap analysis All blocked metabolites in the iPF764 model were identified by using the GapFind MILP algorithm [33] from the COBRA Toolbox. GapFind was run twice, once with each mass balance constraint option. Root no-production and no-consumption metabolites were identified from the model stoichiometric matrix (S) by searching for rows containing only negative or positive coefficients, respectively. Downstream no-production and upstream noconsumption gaps were identified by removing the root gaps from the GapFind outputs [34]. In silico validation of GEM model: prediction of C. salexigens growth phenotype To validate the model for the prediction of growth phenotypes of C. salexigens, the ability to utilize different compounds as the sole carbon sources for growth was determined. For this purpose, the lower bound of glucose exchange reaction was constrained to zero, the lower bound of each different carbon exchange reaction, one at a time, was set to −10 mmolgDW−1h−1 (a typical uptake rate for growth-supporting substrates), and growth was maximized by FBA. The target substrate was considered growth supporting if the predicted growth rate was above zero. The simulations were carried out in minimal medium using BIO_L as biomass equation (see Additional file1: Tables S2, S4). The results arising from the simulations were compared with experimental data [1] and used to validate the model. Gene essentiality To assess the robustness of the metabolic network to genetic perturbations (e.g., knock-out mutations), we carried out an in silico analysis of the essentiality of single genes and reactions. This enabled us to identify the most fragile nodes of the iFP764 network. Gene essentiality simulations were performed by systematically removing each gene from the network and in silico assessing (via FBA) the ability of the model to produce biomass Page 5 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 in minimal medium with a sole carbon source (at low and high salinity). In terms of the metabolic simulations implemented in this study, no-growth following perturbation (deletion of gene-associated reaction) indicates the complete disruption of robustness in C. salexigens and indicates that this gene/reaction is essential for growth. Simulation at low and high salinity conditions The effect of salinity over C. salexigens metabolism at low and high salinity (conditions for low and high ectoines production, respectively) was examined by flux balance analysis (FBA). Rates for uptake and secretion of metabolites in minimal medium M63 with 0.6M or 2.5M NaCl previously described [7], were used as constraints to simulate aerobic growth at both conditions. For this purpose, metabolites present in the medium were first allowed to enter and leave the network by setting the flux limits on the corresponding exchange reactions vi, min>−1000 mmolgDW−1h−1 and vi, max<1000 mmolgDW−1h−1. Second, the glucose and ammonium uptake rates, and the acetate and pyruvate secretion rates, were changed according to experimentally measured fluxes at those salinities, and a biomass reaction (“BIO_L” or “BIO_H”) was selected depending on the salinity simulation conditions. Calculation of fluxes distribution at high and low salinity by Monte Carlo Random sampling The steady-state flux space is a finite polytope containing all possible steady-state flux distributions [35]. Thus, computing the exact volume for large, n-dimensional polytopes is a hard problem [34]. However, it is possible to solve systems of low dimension to examine the accuracy by Monte Carlo approximations. The exact polytope volume calculation algorithm [36] was implemented in Matlab (Mathworks Inc, Massachusetts). The first step in the Monte Carlo calculation was to determine all the independent variables in the system. This essentially ‘projects’ the flux polytope into a fulldimensional object. The full-dimensional polytope was sampled uniformly over the independent variables, and the dependent variables were back-calculated based on the dependencies in the stoichiometric matrix, S. The result was a random flux distribution, V. The next step was to determine if the flux distribution was valid. Back-calculating the dependent flux values often results in invalid flux distributions (i.e. negative flux values or flux values that violate the Vmax constraints). A valid flux distribution was included in the sampled set. Once a set of approximately 20,000 valid points were collected, the process of determining the fractional volumes could start. These valid points were all part of the ‘base’ flux polytope. In this work, Monte Carlo Random Sampling was carried out using the BIO_L and BIO_H to verify the influence of the composition of the formula of the biomass on the internal flows. In addition, specificsalin- ity constraints based on previous experimental data [7] (Additional file1: Table S4) were imposed, and the fraction of points meeting the new criteria determined the reduction in volume. The set of all possible flux values for each reaction were represented as a distribution histogram. Each histogram showed one-dimensional information on its x axis, in terms of the extent of possible values for that particular flux. The y axis represents the “size” of space in the other r–1 dimensions resulting from slicing the metabolic solution space along a specific value of the flux through the indicated reaction [35]. The median of the values obtained was normalized by the salinityspecific glucose uptake rate, and was used to analyze the internal flows in all scenarios calculated. The variation in the flux distribution between both conditions (low and high salinity) was evaluated. For that purpose, flux vector values obtained in the sampling for each reaction in both conditions were compared and evaluated by a Wilcoxon rank sum test. Those reactions which showed significant different distribution (Z statistic over 50 and p value<0.01) were selected. Additionally, to see the extent of those differences, 2000 points from each reaction vector were randomly selected from the distribution values obtained after sampling in both conditions. Those random points were subtracted to each other between conditions. Then, the average from all computed differences for each reaction was considered as an estimate of flux distribution variation between both conditions. Correlation analysis of metabolic fluxes at high and low salinity To detect correlated reaction sets in the metabolic network within each condition, a Spearman Rank correlation matrix was performed taking as an input the 20,000 values of each reaction flux vector obtained during the Monte Carlo sampling. A Spearman Rank correlation value over 0.7 or below −0.7 was considered as the minimal threshold to select reactions that were correlated to each other, respectively. The resulting filtered adjacency matrix was visualized and analyzed in CYTOSCAPE 3.2. as a correlation network, were each node represented a metabolic reaction and each edge represented an observed correlation. Nodes were also grouped according to their membership to different pathways subsystems, which were defined in the metabolic reconstruction iFP764. Page 6 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 Results and discussion Construction of the high‑quality metabolic model iFP764 The iterative process to build up the metabolic reconstruction of the C. salexigens was based on a bottomup approach and exhaustive manual refinement of the reconstruction that systematizes the current knowledge about the metabolism of C. salexigens. A major value of a manual model-building effort is the careful compilation of genes and reactions by revision of the literature and experimental data, BLAST searches, phylogenetic analysis, database comparison, gaps filling, etc., to achieve a robust metabolic model with a highly reliable prediction capacity. All the workflow for the C. salexigens metabolic reconstruction followed standard procedures described in the literature [23], and is detailed in "Methods" and Fig.1. The final reconstruction metabolic network was named iFP764 according to the current standard nomenclature [30]. The most relevant steps are described in detail in (Additional file1). Compartmentalization In C. salexigens, the differential use of the periplasmic or cytoplasmic variants of the Entner–Doudoroff pathway contributes to modulating the production of reducing equivalents (NADPH), adapting it to the requirements for ectoine synthesis [7]. For this reason, the C. salexigens metabolic model iFP764 generated in this work included three distinct cellular compartments: the cytoplasm, the periplasm and the extracellular space. This compartmentalization allowed the inclusion of transport systems in both, the inner and outer membrane, and represents more accurately the metabolic machinery available in C. salexigens, including pathways for compatible solute metabolism, and associated central metabolism (Fig.2a, b). Apart from periplasmic oxidation of glucose to gluconate, the most important routes included in the model were those involved in the transport of d-glucono-1,5-lactone (gl), and 2-keto-d-gluco- nate (2ddglcn), two compounds that can be used by C. salexigens as sole carbon sources [7] (Table2, Fig.2b). Formulation of two biomass reactions to simulate low salinity and high salinity conditions Two specific biomass reactions, “BIO_H” and “BIO_L”, were formulated to simulate high and low salinity conditions, respectively, and incorporated into the reconstruction. The formulation of both reactions was partially based on previous C. salexigens experimental data, and on the formulation defined for Escherichia coli [28]. Ectoine and hydroxyectoine were included in both biomass reactions because they are growth-associated metabolites [3], and their specific contents are affected by the salinity of growth medium [7]. This provided two condition-specific models that will serve as platforms to test metabolic capabilities at low and high salinity. Detailed information on the biomass composition formulation can be found in “Methods” and Additional file1: Tables S1 and S2. Properties of the iFP764 reconstruction The genome of C. salexigens is composed of a circular chromosome of 3696kb encoding 3412 total genes [6] (Fig. 3a). The iFP764 reconstruction accounts for 764 genes, 1530 reactions and 1123 metabolites segregated into 3 compartments: cytoplasmic, periplasmic and extracellular. The metabolic reactions include 943 cytoplasmic reactions, 42 periplasmic reactions, 5 extracellular reactions and 33 multiple compartments reactions (Table2). The C. salexigens reconstructed metabolic network presented in this work covers 22.30% of the total genes present in the genome (Fig.3b). The coverage of the genome is similar to those of metabolic networks reconstructions published, such as those of Pseudomonas aeruginosa (19.3% coverage; [37]), Escherichia coli (32% coverage; [27]), Bacillus subtilis (20% coverage; [38]), Clostridium acetobutylicum (13.0% coverage; [39]) and the halophilic bacterium Halomonas smyrnensis AAD6T (25.3% coverage; [40]). All enzymes included in the iFP764 model were clustered into one of the 12 Orthologous Groups (COGs) functional categories of proteins [41] (Fig.3b). Comparison of iFP764 with the previous C. salexigens metabolic reconstruction iOA584 After completion of the iFP764 metabolic network reconstruction, it was compared with the previous semiautomatically-generated, reconstruction of C. salexigens (iOA584) [42]. Table2 details the number of reactions of iFP764 in each subsystem and the comparison with the previous reconstruction iOA584. Up to 169 genes that were found in the iOA584 model were not include in the iFP764 model, since these routes have not been described in C. salexigens, and the reactions catalyzed by these enzymes were not connected to any pathway. Whereas periplasmic and extracellular metabolic reactions were considered in iFP764, iOA584 did not include any transport reaction other than exchange reactions. In addition, the semi-automated reconstruction contains just single genes associated with reactions, but genes associated to enzymatic complexes or isoenzymes were not annotated, as in iFP764. This improvement is very relevant to implement strategies for in silico knockout analysis or strain design for metabolic engineering. A detailed comparison of the number of ORFs assigned from each COG functional category in both models is presented in Fig.4a. The largest increase in Page 7 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 coverage compared with the previous reconstruction corresponded to transport and metabolism of amino acids (61% increase, with addition of 67 ORFS). Additionally, complete synthesis and degradation pathways of the main compatible solutes of C. salexigens, which are necessary for osmoadaptation and drive metabolic adaptation to salinity, were included in iFP764 (see Additional files 1, 2), whereas some of them were found wrong-gene associated or not present in iOA584. Gap analysis The final reconstruction iFP764 included some incomplete metabolic pathways that contained gaps or blocked reactions, which generated 98 dead-end metabolites: 27 root no-production metabolites (metabolites that participates in consuming reactions but not in producing reactions) and 71 no-consumption metabolites (metabolites that are present in producing reactions but not in consuming reactions) (Additional file1: Table S5). In total, •Comparison with previous metabolic reconstruction • Gap and Dead End analysis •Test if the model produce biomass and predicts growth •Compare predicted physiological properties with known properties • Metabolic osmoadaptation analysis •Gene essentiality analysis •Flux Balance Analysis • Monte Carlo Random Sampling •S matrix and Boundaries definition • Load the reconstruction into Matlab (SBML) •Verify S matrix using COBRA Toolbox • Rational Bottom-up approach partially based on experimental data •Exhaustive review of information from databases, published metabolic reconstructions and literature •Manual re-annotation of metabolic genes • Phylogenetic studies for assignment of isoenzymes • Gap Analysis • Assignment of directionality • Compartmentalization • Determination of biomass composition at low and high salinity Conversion into a mathematical format Computational Network analysis Metabolic Reconstruction and Refinement Genome annotation Platform for Osmoadaptation metabolic analysis and strain design Fig. 1 Scheme of the C. salexigens metabolic network reconstruction and analysis process Page 8 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 b Glc-D[e] Glucose assimilation Glc-D[p] Glc-D[c] g6p[c] g3p[c] gl[p] Glcn[p] Glcn[e] 2-Ketgluc[p] 2-Ketgluc[e] Glcn[c] 2-Ketgluc[c] 2dhdglu6[c] gl[c] 6pgc[c] Gcd (Csal0170) Gcd_2 (Csal2831) Pgl Gnuk KguD (Csal0273) KguK (Csal0274) Glk Edd Eda Zwf 6pgl[c] f6p[c] fdp[c] Gluc6 PP_PfK (Csal1534) EM ED 2ddg6p[c] Asp-L[c] aspsa[c] 4pasp[c] AspK 24abtn[c] AsD DoeC (Csal2725) DoeD (Csal2724) EctB EctA EctC EctDE (Csal0542) DoeA (Csal2732) DoeB (Csal2731) NG4aL24dab[c]ect-L[c] hdect[c] NA4aL24dab[c] ECTtex ECTtra OHtex Compatible solute metabolism Hom-L[c] Homode phom[c] Homok thr-L[c] Tres gly[c] Sarc[c] dmgly[c]glyb[c] betald[c] chol[c] Treal Dimeoxi (Csal0990; Csal0991) saroxi Homosmet EutB (Csal2723) EutC (Csal2722) 2mdi6car[c] BetB BetA chltex GnclI (Csal0453; Csal1727) GnclI (Csal2779) Glc-D[e] Glc-D[c] Glc-D[p] Periplasmic space Periplasmic space g6p[c] tre6p[c] OtsA tre[c] OtsB TreF (Csal0235; Csal1261) (Csal0935; Csal3175) (Csal0137) a Fig. 2 Examples of relevant curated pathways in C. salexigens. a Pathways related with compatible solutes metabolism. b Entner–Doudoroff (ED) and partial glycolysis (EM) pathways of central metabolism. Names of reactions and associated genes are shown in red. A discontinuous arrow means a new reaction or a re-annotated gen-associated reaction in iFP764, with respect to iOA548. A red-truncated arrow represents a geneassociated reaction non-included in iFP764, respect to iOA548. a glc-d[c]: d-glucose; f6p[c]: beta-d-fructose-6-phosphate; tre6p[c]: alpha, alphatrehalose 6-phosphate; tre[c]: alpha, alpha, trehalose; chol[c]: Choline; dmgly[c]: N,N-dimethylglycine; sarc[c]: Sarcosine; gly[c]: Glycine; thr-l[c]: l-Threonine; phom[c]: O-phospho-l-homoserine; hom-l[c]: l-homoserine; asp-L[c]: l-aspartate; aspsa[c]: l-aspartate-4-semialdehyde; 4pasp[c]: 4-phospho-l-aspartate; 24abtn[c]: l-2-4-diamino-butanoate; ect-l[c]: ectoine; 2mdi6car[c]: 2-ethyl-4,5-dihydropyrimidine-6-carboxylate; hdect[c]: 5-Hydroxyectoine; nA4aL24dab[c]: N-alpha-acetyl-l-2,4-diaminobutyrate; nG4aL24dab[c]: N-4-Acetyl-l-2,4-diaminobutanoate. b glc-d[c]: d-glucose; f6p[c]: beta-d-fructose-6-phosphate; fdp[c]: beta-d-fructose-1,6-bisphosphate; g6p[c]: d-glucose 6-phosphate; g3p[c]: glyceraldehyde-3-phosphate; gl[c]: d-glucono-1,5-lactone; 6pgl[c]: d-glucono-1,5-lactone 6-phosphate; glcn[c]: d-gluconic-acid; 6pgc[c]: 6-phospho-d-gluconate; 2ddg6p[c]: 2-Dehydro-3-deoxy-6-phospho-d-gluconate; 2-ketgluc[c]: 2-keto-d-gluconic acid; 2dhdglu6[c]: 6-phospho-2-dehydro-d-gluconate Page 9 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 11.17% of the metabolites in iFP764 were blocked under all conditions due to gaps. The major root no-production metabolites were those implicated in the production of those, like tRNAs, serve for non-metabolic functions once they are produced. This problem is also present in the E. coli iJO1366 metabolic reconstruction [25]. On the other hand, out of the 71 no-consumption metabolites detected mainly belonged to pathways related to tRNA (9.2%), cell envelopes, including cell wall and cell membranes (15.6%), cofactors and prosthetic groups synthesis (8.6%) or carbon metabolism (7.1%). With the exception of the metabolites found in the routes of carbon metabolism, the rest of dead ends metabolites were found in pathways insufficiently described in C. salexigens. An important flaw in the automatically generated iOA584 was the high number of gaps and dead ends metabolites, as we found 272 no-production metabolites and 645 noconsumption metabolites. The effort in the exhaustive refinement process leading to the iFP764 metabolic network resulted in an improvement of connectivity compared with iOA584. A graphic representation of both networks, performed using CYTOSCAPE, evidenced a connectivity failure in the iOA584, as several disconnected reactions and dead-end metabolites from the central network were found. This was highly solved up in iFP764 (Fig.4b). Biomass production Even the most complete models are not definitive, if they contain gaps or missing information [25]. However, filling in the gaps as much as possible is important to make the model more realistic, improving the predictive capability, and this degree of completion is especially needed when the model does not predict growth under some specific condition (i.e., given sets of nutrients and secretions) when submitted to FBA simulation [43]. This occurs due to the fact that the reactions producing or consuming dead-end metabolites can never carry flux under any condition, so one fundamental requirement for a reconstructed metabolic network is its ability to produce flux through the biomass reaction (i.e., to predict growth) [44]. When the iOA584 model was analyzed by FBA with the COBRA Tool box, it could not produce biomass, and consequently it could not simulate cell growth, in contrast to iFP764. This is likely due to the numerous flaws present in the automated network that are absent (or, at least, present in a much lower number) in the curated iFP764 model. Among these, incorrect reaction stoichiometry, repetitions, dead-end metabolites, problems with proton balances, and severe troubles in connectivity, were found in iOA584. In summary, in addition to the enlargement, this comparison revealed great differences between the two models, and a large number of drawbacks associated to semi-automatically-generated networks. Model validation Flux balance analysis (FBA) can be used with a constraint-based model to predict metabolic flux distributions, growth rates, and substrate uptake and production secretion rates. Thus, FBA analyses of the model with a set objective function (i.e. biomass reaction that is equal to growth) to generate maximal amounts of biomass precursors (i.e., simulated optimal growth) from substrates available in growth media can generate results that are quantitatively consistent with experimental data [25, 33]. In this context, the utilization of different carbon sources by the iFP764 network was simulated and compared to previous reports [1, 7]. The growth in minimal medium M63 was simulated by fixing the set values for the exchange reactions for the import of metabolites present in the simulated growth medium. Each substrate was taken as the sole carbon and energy source in the simulation and the uptake rates were Table 2 Properties of metabolic reconstructions of C. salexigens * Reactions can occur in or between multiple compartments and metabolites can be present in more than one compartment Features In silico metabolic model iFP764 iOA584 [42] No. of total reactions 1530 1387 No. of total ORFs associated 1190 874 No. ORFs no associated 340 513 No. of total metabolic reactions 1023 880 Multiple compartments* 33 0 Cytoplasmatic 943 880 Periplasmic 42 0 Extracellular 5 0 No. of total transport reactions 507 507 Exchange 178 507 Extracellular to periplasmic 149 0 Periplasmic to cytoplasmic 137 0 Periplasmic to extracellular 6 0 Cytoplasmic to periplasmic 34 0 Extracellular to cytoplasmic 3 0 No. of total metabolites 1123 1411 Cytoplasmatic 790 920 Periplasmic 176 0 Extracellular 157 491 Unique functional proteins 764 579 Single gene 458 579 Genes involved in complexes 155 0 Instances of isozymes 151 0 Page 16 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 Cs (CITSY)], and with additional reactions of the pentose phosphate pathway. The third cluster of reactions (both of them belonging to the ED pathway, but included in the pentose phosphate subsystem) showed correlation at high salinity with Pdh (Glycolysis/Gluconeogenesis subsystem). The ED-EMP cycle has been suggested to be an essential distributor of carbon through the central metabolic metabolism and an enhanced supply of reducing equivalents [55]. Therefore, the above results could pave the way for future experimental design, in order to figure out how this cycle is functioning in C. salexigens, and if it is partially coupled with the TCA cycle at high salinity, where high ectoines demand requires high levels of reducing power and carbon. Regarding ectoines metabolism, the correlation analysis distributed this subgroup into four cyclic modules. The first module included reactions for the synthesis and degradation of hydroxyectoine (EutC, EutB and EctDE), the second module corresponded to the futile cycle for ectoine synthesis and degradation (EctC, EctA, DoeA, DoeB and EctA) [56], the third module was responsible for the l-aspartate semialdehyde conversion into diamino butyric acid, or viceversa, with an interconversion coupled reaction of glutamate to 2-oxoglutarate (EctB and GlucoseGlucono1,5-lactoneGluconate Glucose Glucose-6P Gluconate Glucono1,5-lactone 6P-Gluconate 6P-Glucono lactone KDPG G3P Pyruvate Fructose-1,6-P2 F6P DHAP PEP Acetyl-CoA citrate isocitrate 2-oxoglutarate succinate fumarate malate OAA glyoxylate Glutamine Glutamate L-Aspartate Aspartyl-P Asp-β-semialdehide DABA NγAcDABA NαAcDABA Ectoine Hydroxiectoine 2mdi6car Ru-5P Xu-5P r5P S7p G3P Acetate Acetyl-P Lactate NH4 0.92 0.92 0.063 0.821 0.912 1.002 0.018 0.489 0.15 0.085 0.008 0.008 0.128 0.489 0.02 0.03 0.03 0.02 0.569 0.537 0.016 0.138 0.128 0.2356 0.166 0.166 0.026 0.629 0.02 0.071 0.062 0.89 0.177 0.018 0.01 0.03 0.014 0.045 0.045 27.46 27.46 0.062 0.062 0.003 0.003 0.007 0.072 0.072 0.5337 0.03 0.859 0.028 0.821 1.002 0.113 Ru-1,5 diP 3PG CO2 1,3 bi-PG 0.0174 0.017 0.75 0.75 0.356 0.069 0.025 e4p 0.085 1,3 bi-PG 0.75 3PG 0.75 GlucoseGlucono1,5-lactone Gluconate Glucose Glucose-6P Gluconate Glucono1,5-lactone 6P-Gluconate 6P-Glucono lactone KDPG G3P Pyruvate Fructose-1,6-P2 F6P DHAP PEP Acetyl-CoA citrate isocitrate 2-oxoglutarate succinate fumarate malate OAA glyoxylate Glutamine Glutamate L-Aspartate Aspartyl-P Asp-β-semialdehide DABA NγAcDABA NαAcDABA Ectoine Hydroxiectoine 2mdi6car Ru-5P r5P G3P Acetate Acetyl-P Lactate NH4 0.841 0.841 0.155 0.625 0.625 0.881 1.587 0.068 0.142 0.724 0.027 0.027 0.801 0.142 0.0497 0.026 0.026 0.417 0.686 0.546 0.139 0.2463 0.813 0.755 0.755 0.339 1.535 0.043 0.585 0.574 1.811 0.919 0.056 0.09 0.102 0.029 0.071 0.071 329 329 0.04 0.04 0.006 0.006 0.012 0.0545 0.0545 1.555 0.062 0.751 1.587 0.07 Ru-1,5 diP 3PG CO2 1,3 bi-PG 0.046 0.736 0.736 0.031 0.031 0.033 0.064 0.106 Xu-5P S7p e4p 1,3 bi-PG 3PG 0.736 0.736 0.78 0.398 0.398 0.398 Low salinity(0.6 M NaCl)High salinity(2.5 M NaCl) ab EX_glc-D(e) 1.0 EX_glc-D(e) 1.0 Fig. 7 Overview of central and ectoines metabolism fluxes distribution at low and high ectoines production conditions (only the first scenario is shown). Median values of fluxes calculated by Monte Carlo Random sampling for each reaction are normalized relative to the glucose consumption rate of each specific salinity-condition a Low salinity (0.6 M NaCl). b High salinity (2.5 M NaCl). Ru-5P d-ribulose-5-phosphate, Xu-5P d-xylulose- 5-phosphate, r5P d-ribose-5-phosphate, S7p sedoheptulose-7-phosphate, e4p d-erythrose-4-phosphate, G3P glyceraldehyde-3-phosphate, F6P fructose-6-phosphate, 1,3 bi-PG 3-phospho-d-glyceroyl-phosphate, 3PG 3-phospho-d-glycerate, PEP phosphoenolpyruvate, OAA oxaloacetate, DABA l-2-4-diamino-butanoate; NαAcDABA N-alpha-acetyl-l-2,4-diaminobutyrate, 2mdi6car 2-methyl-4,5-dihydropyrimidine-6-carboxylate Page 17 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 TCA Gdh Gly/Ser/Trh PP Glycolisis/ gluconeog Purine/pirimidine Val/Leu/Ile Tyr/Trp/Phe Trh/Lys His AsD DoeC Met Cys Folate Ectoines metabolism DoeD EctD Ala/Asp EctA DoeB DoeA EctC EctB Glutamate Metabolism AlternateC Trehalose metabolism Ppc Arg/Pro EutC EutB Pyruvate metabolism TCA Gdh Gly/Ser/Trh PP PP C fix Glycolisis/ Gluconeog Purine/pirimidine Val/Leu/Ile Tyr/Trp/Phe Trh/Lys His AsD DoeC Met Cys Folate Folate Ectoines metabolism DoeD EctD Ala/Asp Arg/Pro EctA DoeB DoeA EctC EctB Glutamate metabolism Alternate C Trehalose metabolism LacDH Arg/Pro EutC Ms EutB OtsA OtsB TreF Ict Oxidative phosphorilation Pyruvate metabolism OtsA OtsBTreF a b Aspk Aspk Arg/Pro Folate PP C fix Oxidative phosphorilation Ms Anaplerotic reactions Fig. 8 Fluxes correlation network of compatible solute and main carbon and nitrogen metabolism pathways at low (a) and high (b) salinity. Each node represents a metabolic reaction and each edge represents an observed correlation (with Spearman Rank correlation value larger than 0.7). The adjacency matrix was visualized in CYTOSCAPE (3.2). Nodes were grouped according to their membership to different pathway subsystems and some reactions/enzymes was specified. Aspk aspartate kinase, AsD aspartate semialdehyde dehydrogenase, C fix carbon fixation, Ict isocitrate lyase, Ms malate synthase, PP pentoses phosphate, Ppc phosphoenolpyruvate carboxylase Page 18 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 DoeD). The last module included reactions responsible of the synthesis of L-aspartate to L-aspartate semialdehyde, or viceversa, [DoeC, AsD, aspartate kinase (ASPK)]. Interestingly, all the four modules were independent both at low and high salinity conditions, as they did not correlate to each other, and they were not correlated with other subsystems (Fig.8). Regarding other compatible solutes, trehalose synthesis and degradation reactions (OtsA, OtsB and TreF) formed a separate module including internal reactions of central C subsystem, particularly those involved in the conversion of glucose-1-P into UDP-glucose (UTPGLPHO, UTP-glucose 1-P-uridylyl transferase) and glucose-6P (PGM, phosphoglucomutase). Both products are precursors of trehalose-6-phosphate. As observed for ectoines metabolism, this correlation pattern of the trehalose metabolism module did not vary with salinity (Fig.8). The predicted modular distribution of correlated compatible solutes metabolism sets of reactions is very intriguing, and deserves experimental work to better understand C. salexigens osmoadaptive mechanisms, and especially ectoines production. Most of them were cyclic (not only the ectoine futile cycle), independent from other subsystems, and with no interdependence pattern variation in response to salinity. It is plausible that reactions included in each module respond to a similar regulation, and/or that a coordinated functioning of cyclic modules is necessary to adapt solute compatible metabolism, and therefore metabolic osmoadaptation, to external conditions. Our results are only an example (“the tip of the iceberg”) of the potential of the joint rational use of computational tools, together with experimental data, in order to explore metabolism in halophilic microorganisms. Conclusions The information resulting from this work showed that the iFP764 model is a powerful tool to simulate and explore C. salexigens metabolism at different salinity conditions. This halophilic bacterium has one of the broadest salinity range of growth, and is considered as a model to study osmoadaptation. Notwithstanding the existing limitations in reconstructing and validating a C. salexigens specific network, the iFP764 model is one of the most complete and curated networks of halophilic bacteria, and accounts for all existing relevant experimental information. Thus, it may become an important tool to understand the metabolism of halophilic microorganisms, helping to drive the generation of new experimental hypotheses. In addition, C. salexigens is considered as a good cell factory for ectoines, whose production is salinity-dependent. Therefore, this model could be a robust and reliable tool to be used for the design of new C. salexigens strains optimized for the production of ectoines and other biotechnologically interesting compounds, avoiding the expenditure of resources. Authors’ contributions MA and CV conceived and supervised the study. FP, MS, and VB performed the model construction and simulations. FP, MS, MA, VB, JMP, JJN, MC and CV were involved in the analysis and discussion of results. FP, MA and CV drafted the manuscript. All authors read and approved the final manuscript. Author details 1 Department of Microbiology and Parasitology, Faculty of Pharmacy, University of Sevilla, C/Profesor García González 2, 41012 Sevilla, Spain. 2 Department of Biochemistry and Molecular Biology B and Immunology, Faculty of Chemistry, Campus Regional de Excelencia Internacional “Campus Mare Nostrum”, University of Murcia, 30100 Murcia, Spain. 3 Centro de Tecnología de Repsol, REPSOL S.A. Calle Agustín de Betancourt, s/n. 28935, Móstoles, Madrid, Spain. Acknowledgements We thank Dr. L.N. Csonka for language editing of the manuscript. Competing interests The authors declare that they have no competing interests. Availability of data and materials All data generated or analyzed during this study are included in this published article and its additional files. Consent for publication Not applicable. Ethics approval and consent to participate Not applicable. Funding This work has been partly funded by Projects BIO2014-54411-C2-1-R and BIO2015-63949-R (MINECO/FEDER, UE), Junta de Andalucía (Spain; Project P08-CVI-03724), and Spanish National Network on Extremophilic Microorganisms (BIO2015-71815-REDT). Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Additional files Additional file 1: Table S1. Biomass composition at high salinity (´´BIO_ H``). Table S2. Biomass composition at low salinity (´´BIO_L``); Table S3. Composition of M63 minimal medium and the exchange reactions formulated to simulate the uptake of metabolites from medium; Table S4. In silico (computational) constraints used for simulation at high or low salinity; Description of the bottom-up building and exhaustive manual refinement of the C. salexigens metabolic reconstruction iFP764; Table S5. The total number of dead ends metabolites found in the iFP764 model: root no-production metabolites and root no-consumption metabolites; Table S6. In silico prediction of utilization of various metabolites as carbon sources in C. salexigens; Figure S1. Selected histograms of possible flux values obtained in the first scenario by Monte Carlo sampling at low and high salinity relative to salinity-specific glucose consumption rate. Additional file 2. Genome-scale metabolic reconstruction of C. salexigens. Details of iPF764 containing all genes associated-reactions and metabolites. Additional file 3: Table S7. Essential genes at low and high salinity and their associated COGs. Page 19 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 Received: 27 July 2017 Accepted: 20 December 2017 References 1. Arahal DR, García MT, Vargas C, Cánovas D, Nieto JJ, Ventosa A. Chromohalobacter salexigens sp. nov., a moderately halophilic species that includes Halomonas elongata DSM 3043 and ATCC 33174. Int J Syst Evol Microbiol. 2001;51:1457–62. 2. Vargas C, Argandoña M, Reina-Bueno M, Rodríguez-Moya J, Fernández- Aunión C, Nieto JJ. Unravelling the adaptation responses to osmotic and temperature stress in Chromohalobacter salexigens, a bacterium with broad salinity tolerance. Saline Syst. 2008;4:14. 3. Fallet C, Rohe P, Franco-Lara E. Process optimization of the integrated synthesis and secretion of ectoine and hydroxyectoine under hyper/ hypo-osmotic stress. Biotechnol Bioeng. 2010;107(1):124–33. 4. Rodríguez-Moya J, Argandoña M, Iglesias-Guerra F, Nieto JJ, Vargas C. Temperature- and salinity-decoupled overproduction of hydroxyectoine by Chromohalobacter salexigens. Appl Environ Microbiol. 2013;79(3):1018–23. 5. Argandoña M, Vargas C, Reina-Bueno M, Rodríguez-Moya J, Salvador M, Nieto JJ. An extended suite of genetic tools for use in bacteria of the Halomonadaceae: an overview. In: Lorence A, editor. Recombinant gene expression: reviews and protocols. Methods molecular biology, vol. 824. Ciudad: SpringerLink; 2012. p. 167–201. 6. Copeland A, O’Connor K, Lucas S, Lapidus A, Berry KW, Detter JC, Del Rio TG, Hammon N, Dalin E, Tice H, Pitluck S, Bruce D, Goodwin L, Han C, Tapia R, Saunders E, Schmutz J, Brettin T, Larimer F, Land M, Hauser L. Complete genome sequence of the halophilic and highly halotolerant Chromohalobacter salexigens type strain (1H11T). Stand Genomic Sci. 2011;5(3):379–88. 7. Pastor JM, Bernal V, Salvador M, Argandoña M, Vargas C, Csonka L, Sevilla A, Iborra JL, Nieto JJ, Cánovas M. Role of central metabolism in the osmoadaptation of the halophilic bacterium Chromohalobacter salexigens. J Biol Chem. 2013;288(24):17769–81. 8. O’Brien EJ, Monk JM, Palsson BO. Using genome-scale models to predict biological capabilities. Cell. 2015;161(5):971–87. 9. Lewis NE, Nagarajan H, Palsson BO. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat Rev Microbiol. 2012;10(4):291. 10. Price ND, Reed JL, Palsson BO. Genome scale models of microbian cells: evaluating the consequences of constraints. Nat Rev Microbiol. 2004;2:886–97. 11. Papin JA, Price ND, Wiback SJ, Fell DA, Palsson BO. Metabolic pathways in the post-genome era. Trends Biochem Sci. 2003;28:250–8. 12. Lee KH, Park JH, Kim TY, Kim HU, Lee SY. Systems metabolic engineering of Escherichia coli for l-threonine production. Mol Syst Biol. 2007;3:149. 13. Pharkya P, Burgard AP, Maranas CD. Exploring the overproduction of amino acids using the bilevel optimization framework OptKnock. Biotechnol Bioeng. 2003;84:887–99. 14. Reed JL, Famili I, Thiele I, Palsson BO. Towards multidimensional genome annotation. Nat Rev Genet. 2006;7(2):130–41. 15. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10. 16. Marchler-Bauer A, Panchenko AR, Shoemaker BA, Thiessen PA, Geer LY, Bryant SH. CDD: a database of conserved domain alignments with links to domain three-dimensional structure. Nucleic Acids Res. 2002;30(1):281–3. 17. Schultz J, Milpetz F, Bork P, Ponting CP. SMART, a simple modular architecture research tool: identification of signaling domains. Proc Natl Acad Sci. 1998;95:5857–64. 18. Snel B, Lehmann G, Bork P, Huynen MA. STRING: a web-server to retrieve and display the repeatedly occurring neighborhood of a gene. Nucleic Acids Res. 2000;28(18):3442–4. 19. Despalins A, Marsit S, Oberto J. Absynte: a web tool to analyze the evolution of orthologous archaeal and bacterial gene clusters. Bioinformatics. 2011;27(20):2905–6. 20. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1990;22(22):4673–80. 21. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9. 22. Schomburg I, Chang A, Schomburg D. BRENDA, enzyme data and metabolic information. Nucleic Acids Res. 2002;30:47–9. 23. Thiele I, Palsson BO. A protocol for generating a high-quality genomescale metabolic reconstruction. Nat Protoc. 2010;5(1):93–121. 24. Vargas C, Kallimanis A, Koukkou AI, Calderon MI, Canovas D, Iglesias- Guerra F, Drainas C, Ventosa A, Nieto JJ. Contribution of chemical changes in membrane lipids to the osmoadaptation of the halophilic bacterium Chromohalobacter salexigens. Syst Appl Microbiol. 2005;28(7):571–81. 25. Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, Palsson BO. A comprehensive genome-scale reconstruction of Escherichia coli metabolism. Mol Syst Biol. 2011;7:535. 26. Neidhardt FC, Van Bogelen RA. Heat shock response. Escherichia coli and Salmonella typhimurium; cellular and molecular biology. In: Neidhardt FC, Ingraham JL, Low KB, Magasanik B, Schaechter M, Umbarger HE, editors. American society for microbiology. Washington: American Society for Microbiology; 1987. 27. Cayley S, Record MT. Roles of cytoplasmic osmolytes, water, and crowding in the response of Escherichia coli to osmotic stress: biophysical basis of osmoprotection by glycine betaine. Biochemistry. 2003;42(43):12596–609. 28. Feist AM, Henry CS, Reed JL, Krummenacker M, Joyce AR, Karp PD, Broadbelt LJ, Hatzimanikatis V, Palsson BO. A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Mol Syst Biol. 2007;3:121. 29. Feist AM, Palsson BO. The biomass objective function. Curr Opin Microbiol. 2010;3:344–9. 30. Reed JL, Patel TR, Chen KH, Joyce AR, Applebee MK, Herring CD, Bui OT, Knight EM, Fong SS, Palsson BO. Systems approach to refining genome annotation. Proc Natl Acad Sci USA. 2006;103:17480–4. 31. Becker SA, Feist AM, Mo ML, Hannum G, Palsson BO, Herrgard MJ. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox. Nat Protoc. 2007;2(3):727–38. 32. Csonka LN. A third l-proline permease in Salmonella typhymurium which functions in media of elevated osmotic strength. J Bacteriol. 1982;151:1433–43. 33. Kumar SV, Dasika MS, Maranas CD. Optimization based automated curation of metabolic reconstructions. BMC Bioinform. 2007;8:212. 34. Orth JD, Palsson BO. Systematizing the generation of missing metabolic knowledge. Biotechnol Bioeng. 2010;107(3):403–12. 35. Wiback SJ, Famili I, Greengerg HJ, Palsson BO. Monte Carlo sampling can be used to determine the size and shape of the steady-state flux space. J Theor Biol. 2004;228:437–47. 36. Lawrence J. Polytope volume computation. Math Comput. 1991;195:259–71. 37. Bartell JA, Blazier AS, Yen P, Thøgersen JC, Jelsbak L, Goldberg JB, Papin JA. Reconstruction of the metabolic network of Pseudomonas aeruginosa to interrogate virulence factor synthesis. Nat Commun. 2017;8:14631. 38. Oh Y-K, Palsson BØ, Park SM, Schilling CH, Mahadevan R. Genomescale reconstruction of metabolic network in Bacillus subtilis based on high-throughput phenotyping and gene essentiality data. J Biol Chem. 2007;282:28791–9. 39. McAnulty MJ, Yen JY, Freedman BG, Senger RS. Genome-scale modeling using flux ratio constraints to enable metabolic engineering of clostridial metabolism in silico. BMC Syst Biol. 2012;6:42. 40. Diken E, Ozer T, Arikan M, Emrence Z, Oner ET, Ustek D, Arga KY. Genomic analysis reveals the biotechnological and industrial potential of levan producing halophilic extremophile, Halomonas smyrnensis AAD6T. SpringerPlus. 2015;4:393. 41. Tatusov RL, Koonin EV, Lipman DJ. A genomic perspective on protein families. Science. 1997;278(5338):631–7. 42. Ates O, Oner ET, Arga KY. Genome-scale reconstruction of metabolic network for a halophilic extremophile, Chromohalobacter salexigens DSM3043. BMC Syst Biol. 2011;5:12. 43. Ponce-de-León M, Montero F, Peretó J. Solving gap metabolites and blocked reactions in genome-scale models: application to the metabolic network of Blattabacterium cuenoti. BMC Syst Biol. 2013;7:114. Page 20 of 20 Piubeli et al. Microb Cell Fact (2018) 17:2 • We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal • We provide round the clock customer support • Convenient online submission • Thorough peer review • Inclusion in PubMed and all major indexing services • Maximum visibility for your research Submit your manuscript at www.biomedcentral.com/submit Submit your next manuscript to BioMed Central and we will help you at every step: 44. Larocque M, Chénard T, Najmanovich R. A curated C. difficile strain 630 metabolic network: prediction of essential targets and inhibitors. BMC Syst Biol. 2014;8:117. 45. Schatschneider S, Persicke M, Watt SA, Hublik G, Puhler A, Niehaus K, Vorholter FJ. Establishment, in silico analysis, and experimental verification of a large-scale metabolic network of the xanthan producing Xanthomonas campestris pv. campestris strain B100. J Biotechnol. 2013;167:123–34. 46. Kitano H. Towards a theory of biological robustness. Mol Syst Biol. 2007;3:137. 47. Ye YN, Hua ZG, Huang J, Rao N. Guo FB1. CEG: a database of essential gene clusters. BMC Genom. 2013;14:769. 48. Oren A. Bioenergetic aspects of halophilism. Microbiol Mol Biol Rev. 1999;63:334–48. 49. Puchałka J, Oberhardt MA, Godinho M, Bielecka A, Regenhardt D, Timmis KN, Papin JA, Martins dos Santos VA. Genome-scale reconstruction and analysis of the Pseudomonas putida KT2440 metabolic network facilitates applications in biotechnology. PLoS Comput Biol. 2008;4(10):e1000210. 50. Pfeiffer T, Schuster S. Game-theoretical approaches to studying the evolution of biochemical systems. Trends Biochem Sci. 2005;30:20–5. 51. Schuster S, Pfeiffer T, Fell DA. Is maximization of molar yield in metabolic networks favored by evolution? J Theor Biol. 2008;252:497–504. 52. Metris A, George S, Baranyi S. Modelling osmotic stress by Flux Balance Analysis at the genomic scale. Int J Food Microbiol. 2012;152(3):123–8. 53. Covert MW, Palsson BO. Transcriptional regulation in constraints-based metabolic models of Escherichia coli. J Biol Chem. 2002;277(31):28058–64. 54. Papin JA, Reed JL, Palsson BO. Hierarchical thinking in network biology: the unbiased modularization of biochemical networks. Trends Biochem Sci. 2004;29(12):641–7. 55. Nikkel PI, Chavarría M, Fuhrer T, Sauer U, de Lorenzo V. Pseudomonas putida KT2440 strain metabolizes glucose through a cycle formed by enzymes of the Entner–Doudoroff, Embden–Meyerhof–Parnas and Pentose phosphate pathways. J Biol Chem. 2015;43:25920–32. 56. Schwibbert K, Marin-Sanguino A, Bagyan I, Heidrich G, Lentzen G, Seitz H, Rampp M, Schuster SC, Klenk HP, Pfeiffer F, Oesterhelt D, Kunte HJ. A blueprint of ectoine metabolism from the genome of the industrial producer Halomonas elongata DSM 2581T. Environ Microbiol. 2011;13(8):1973–94.