Genome-scale model of Rothia mucilaginosa predicts gene essentialities and reveals metabolic capabilities
ABSTRACT
Cystic fibrosis (CF), an inherited genetic disorder caused by mutations in the cystic fibrosis transmembrane conductance regulator gene, results in sticky and thick mucosal fluids. This environment facilitates the colonization of various microorganisms, some of which can cause acute and chronic lung infections, while others may positively impact the disease. Rothia mucilaginosa, an oral commensal, is relatively abundant in the lungs of CF patients. Recent studies have unveiled its anti-inflammatory properties using in vitro three-dimensional lung epithelial cell cultures and in vivo mouse models relevant to chronic lung diseases. Apart from this, R. mucilaginosa has been associated with severe infections. However, its metabolic capabilities and genotype-phenotype relationships remain largely unknown. To gain insights into its cellular metabolism and genetic content, we developed the first manually curated genome-scale metabolic model, iRM23NL. Through growth kinetics and high-throughput phenotypic microarray testings, we defined its complete catabolic phenome. Subsequently, we assessed the model’s effectiveness in accurately predicting growth behaviors and utilizing multiple substrates. We used constraint-based modeling techniques to formulate novel hypotheses that could expedite the development of antimicrobial strategies. More specifically, we detected putative essential genes and assessed their effect on metabolism under varying nutritional conditions. These predictions could offer novel potential antimicrobial targets without laborious large-scale screening of knockouts and mutant transposon libraries. Overall, iRM23NL demonstrates a solid capability to predict cellular phenotypes and holds immense potential as a valuable resource for accurate predictions in advancing antimicrobial therapies. Moreover, it can guide metabolic engineering to tailor R. mucilaginosa’s metabolism for desired performance.
IMPORTANCE
Cystic fibrosis (CF) is a genetic disorder characterized by thick mucosal secretions, leading to chronic lung infections. Rothia mucilaginosa is a common bacterium found in various parts of the human body, acting as a normal part of the flora. In people with weakened immune systems, it can become an opportunistic pathogen, while it is prevalent and active in CF airways. Recent studies have highlighted its anti-inflammatory properties in the lower pulmonary system, indicating the intricate relationship between microbes and human health. Herein, we have developed the first manually curated metabolic model of R. mucilaginosa. Our study examined the previously unknown relationships between the bacterium’s genotype and phenotype and identified essential genes that impact the metabolism under various conditions. With this, we opt for paving the way for developing new strategies in antimicrobial therapy and metabolic engineering, leading to enhanced therapeutic outcomes in cystic fibrosis and related conditions.
INTRODUCTION
Rothia mucilaginosa is a Gram-positive, encapsulated, non-motile, and non-spore-forming bacterium of the Micrococcaceae family (1, 2). While it is mainly aerobic, it may also grow anaerobically as it can switch to fermentation or other non-oxygen-involving pathways. R. mucilaginosa is a common commensal of the normal oral, upper and lower respiratory tract, and part of the skin florae in humans (1, 3–6). This means it coexists harmlessly within the host and may even provide benefits. Nonetheless, it can also act as an opportunistic pathogen, particularly in individuals with weakened immune systems, as an etiological agent of serious infections such as endocarditis, sepsis, and meningitis (7). Janek et al. highlighted the high prevalence of R. mucilaginosa within the nasal microbiome (8). Moreover, they report its susceptibility to certain staphylococcal bacteriocins, indicating its major competition with the nasal staphylococci and the substantial impact of bacteriocins in shaping the nasal microbiota. In 2020, Uranga et al. revealed that R. mucilaginosa produces the strongest Fe3+-binding archetypal siderophore known, called enterobactin (9). This attribute contributes to its competition with oral microbiota (the cariogenic S. mutans, A. timonensis, and Streptococcus sp.) as well as four methicillin-resistant strains of S. aureus (MRSA). Enterobactin is a type of siderophore produced by bacteria to scavenge, chelate, and transport ferric irons from their surrounding environment. These are essential for bacteria when iron is scarce as they facilitate their acquisition necessary for their growth and metabolic processes.
Prior metagenomic sequencing analyses have unveiled the prevalence of R. mucilaginosa at high abundances and its enhanced metabolic activity in the lungs of cystic fibrosis (CF) patients (10, 11). CF is caused by the hereditary mutation of the cystic fibrosis transmembrane conductance regulator (CFTR) gene that disrupts the transepithelial movement of ions, leading to an aberrant accumulation of thick and sticky mucus within the airways. The impaired immune clearance creates a hypoxic environment (12) promoting the polymicrobial colonization of opportunistic microbes together with fungi and viruses, ultimately resulting in persistent and recurring infections (13). Guss et al. and Bittar et al. declared R. mucilaginosa as an emerging CF bacterium (14, 15), while Lim et al. provided evidence supporting that R. mucilaginosa is a frequently encountered and metabolically active inhabitant of CF airways (16). Additionally, a study from 2018 shows that the opportunistic pathogen Pseudomonas aeruginosa, which frequently causes infections in CF patients, builds essential primary metabolites, like glutamate, by utilizing compounds produced by R. mucilaginosa (17). This symbiotic interaction implies that P. aeruginosa benefits from its neighboring microbes, which contributes to its pathogenesis in the CF lungs. On the other hand, Rigauts et al. revealed the anti-inflammatory activity of R. mucilaginosa in the lower respiratory tract, which could impact the seriousness of chronic lung diseases (18).
In systems biology, genome-scale metabolic models (GEMs) represent comprehensive reconstructions of organisms’ metabolic networks. They are built using genomic sequences and comprise all known biochemical reactions and associated genes. These models provide systems-level insights into cellular metabolism, allowing researchers to simulate and analyze the flow of metabolites through these networks (19). The interactions among reactions and metabolites in a metabolic model are mathematically represented with a stoichiometric matrix (20). In the past years, an array of in silico methods has been developed to analyze GEMs and derive valuable hypotheses. Flux balance analysis (FBA) is such a powerful computational technique that operates on the principle of achieving a steady state by optimizing the flux (rate) of metabolites through reactions while accounting for various constraints such as stoichiometry, thermodynamics, and uptake/secretion boundaries (21). Applying flux balance analysis on a GEM provides insights into the intricate biological system interactions. This analytical approach facilitates the prediction of cellular phenotypes and identification of promising drug targets and contributes to optimizing biotechnological processes (22). Moreover, such models can guide genetic engineering by suggesting genetic modifications that could enhance desired product formation or cellular behavior. Further applications include ameliorating culture media by incorporating components that increase bacterial growth rates. So far, GEMs have been an invaluable resource in the systems biology field that helped untangle the metabolism of various organisms and especially of high-threat pathogens (23, 24). As described above, R. mucilaginosa has gained great interest in the context of polymicrobial CF environments. However, its metabolic capabilities and genotype-phenotype relationships in isolated monoculture settings remain largely unexplored.
Here, we present the first manually curated and high-quality GEM of R. mucilaginosa, iRM23NL, striving to understand its metabolism and unique phenotypes under diverse conditions. Our simulation-ready network accounts for thousands of reactions and is available in a standardized format following the community guidelines (25). Through growth kinetic experiments and high-throughput phenotypic microarray assays, we validated iRM23NL’s accuracy in predicting growth and substrate utilization patterns. We refined the reconstruction by comparing the in vitro results to in silico simulations, resulting in novel metabolic reactions and genes. To our knowledge, this is the first study presenting high-throughput nutrient utilization and comprehensive growth data for R. mucilaginosa. Finally, we employed FBA to formulate novel gene essentiality hypotheses that could expedite the development of antimicrobial strategies. Figure 1 summarizes the experimental and computational work presented here.
Fig 1
RESULTS
Reconstruction of a high-quality metabolic model for R. mucilaginosa DSM20746
The pipeline we previously developed (26) was used to build the first high-quality and manually curated GEM of R. mucilaginosa DSM20746 (ATCC 25296). An initial draft metabolic model was derived with CarveMe (27) and is based on the Biochemical, Genetical, and Genomical (BiGG) identifiers (28). The translated sequence with over 1,700 proteins and the Gram-positive-specific template were employed. This enabled us to build a more precise reconstruction considering information on the peptidoglycan layer for the biomass objective function (BOF). The draft network contained 1,015 reactions (141 pseudo-reactions), 788 metabolites, and 265 genes (Fig. 2). In the first gap-filling stage (Draft_2), we expanded the list of reactions based on the annotated genome and growth kinetics data in diverse growth environments. For this, we extensively indexed organism-specific literature and databases and included additional enzymatic reactions together with 50 new gene-protein-reaction associations (GPRs). Subsequently, high-throughput nutrient utilization assays and model validation incorporated further reactions and their associated metabolic genes. Non-metabolic genes, which take part in other cellular processes, e.g., signaling pathways or transcription, were not considered. In total, 95 reactions, together with their associated metabolites, were newly added into the model, along with 121 novel GPRs, increasing the genetic coverage. Over 60% of the transport reactions have a GPR assigned, while 63% of the total enzymatic reactions have at least one gene assigned. Moreover, missing exchange reactions were added to all extracellular metabolites to represent the exchange of substrates between the extracellular environment and the model. The strain-specific BioCyc31 database was further employed to correct the reversibility of biochemical reactions, while duplicated reactions and metabolites were eliminated. In all cases, when no organism-specific information was available, we leveraged data from closely related species based on our phylogenomic analysis (Fig. 3). According to the calculated average nucleotide identity (ANI) matrix, R. mucilaginosa exhibits a similarity to six out of the 12 tested Rothia genomes. More specifically, it shares a greater resemblance with R. aeria and R. dentocariosa underscoring a closer evolutionary relationship between these species.
Fig 2
Fig 3
R. mucilaginosa is primarily aerobic, efficiently generating ATP through oxic respiration; however, in low-oxygen or oxygen-absent conditions, it shifts to anaerobic metabolism to produce energy. This metabolic adaptability enables R. mucilaginosa to adapt in microaerophilic environments like the oxygen-restricted conditions in CF lungs (16). Our draft model lacked the ability to demonstrate anaerobic growth. Therefore, we investigated the metabolic cascade and systematically incorporated missing enzymes to ensure that the model can simulate growth even in the absence of oxygen by identifying and integrating alternative pathways. This refinement included the incorporation of enzymatic reactions, such as the superoxide dismutase (SPODM) and catalase (CAT) that are responsible for the breakdown of radical reactive oxygen species (ROS) and shielding the cell against oxidative damage (Fig. 4 Panel A). Such scavenging enzymes play an integral role in counteracting the harmful effects of ROS during anaerobic respiration (31). However, during this process, we found no associated GPRs for CAT within the organism-specific BioCyc database. Additional scavenging enzymes like glutathione and thioredoxin reductases essential for maintaining the redox balance (32) were already present in the initial draft model (GTHOr, GTHRDabc2pp, and TRDR). Altogether, the final model, iRM23NL, contains 1,162 reactions (619 gene-associated; 65 catalyzed by enzyme complexes, 70 catalyzed by isozymes, and 484 by simple gene association), 171 exchange and sink reactions, 874 metabolites (558 in cytoplasm, 148 in periplasm, and 168 in the extracellular space), and 372 genes (Fig. 2). The model’s metabolic coverage is at 3.12%, which indicates a high level of modeling detail regarding reactions, enzymes, and their associated genes (33). Additionally, we enriched the model elements with numerous database cross-references (34), while appropriate and precise Systems Biology Ontology (SBO) terms were assigned to each model entity using the SBOannotator package (35). The presence of no energy-generating cycles (EGCs) was ensured and controlled after each curation stage, and the mass- and charge-imbalances were corrected. With this, the final Metabolic Model Testing (memote) (36) score of iRM23NL is 89%, while with highly specific SBO terms, the score drops by 2%. The final curated model, iRM23NL, is available as a supplementary file in Systems Biology Markup Language (SBML) Level 3 Version 1 (37) and JavaScript Object Notation (JSON) formats with the flux balance constraints (fbc) and groups plugins available.
Fig 4
The first validation step of iRM23NL aimed to evaluate its ability to correctly simulate biomass production across diverse environmental conditions and growth media formulations. To elucidate the bacterium’s optimal conditions and metabolic preferences, we experimentally tested five commonly used media, including three general nutrient media; brain heart infusion (BHI) and Luria-Bertani (LB), and tryptic soy broth (TSB), and two defined media; M9 minimal medium (M9) pure and Roswell Park Memorial Institute (RPMI) (Fig. 4 Panel B). The BHI medium was used as a baseline for the in vitro experiments since it is a known and well-established environment for the growth of R. mucilaginosa and enabled us to compare the bacterium’s growth characteristics to the newly tested media. For the in silico simulations, we applied FBA and added additional constraints to the linear programming problem defined in Equation 5. In more detail, we specified the flux constraints such that only extracellular metabolites defined in the medium of interest could flow freely through the system (unconstrained, finite fluxes), while the remaining fluxes were constrained to zero. We compared the in vitro to the in silico observed growth using the FCOD as a qualitative measure of growth (see Materials and Methods). Furthermore, we compared the OD at the start and the end of the experiment, considering a statistically significant difference between these measurements as an indication of growth. Our metabolic network, iRM23NL, simulated positive fluxes through the biomass reaction for all tested media except for the M9 pure medium, where a zero flux was observed. These findings align with the experimentally observed data. More specifically, there is no statistically significant difference in OD between the initial and final time-points in M9 pure medium (P-value = 0.1202 and FCOD < 1.4) indicating no significant growth. Conversely, in the remaining examined media, statistically significant growth was observed (P-value = 0.00006–0.00142 and FCOD > 1.4) indicating significant growth in these settings. The highest aerobic growth rate was predicted in TSB [1.6 mmol/(gDW · h)], while the lowest biomass production flux was recorded for the M9 pure medium containing only essential salts. However, the RPMI medium followed as the second-highest in supporting bacterial in vitro cellular growth, offering a defined medium suitable for R. mucilaginosa’s cultivation. Although R. mucilaginosa increased its biomass after 24 h, it slightly declined after 48 h. On the other hand, the simulated network resulted in a contrary outcome compared to the expected experimental effect. More specifically, iRM23NL simulated a lower flux through biomass [0.44 mmol/(gDW · h)] with RPMI when compared to LB. It is important to note here that in order to simulate growth in RPMI medium, six metal ions [cobalt (Co2+), copper (Cu2+), manganese (Mn2+), zinc (Zn2+), ferric iron (Fe3+), and ferrous iron (Fe2+)] were supplemented. These compounds were missing from the providers’ medium formulation. Our findings underscored R. mucilaginosa’s adaptability to various nutritional environments, growing best in nutrient-rich conditions while revealing specific growth requirements beyond minimal settings.
We further employed iRM23NL to examine whether it could generate biomass within the human nasal environment and the CF lungs. For this purpose, we performed in silico simulations using the synthetic cystic fibrosis sputum medium (SCFM) (38) and synthetic nasal medium (SNM) (39) (Fig. 4 Panel C). Our computational model successfully simulated positive growth in both media, with a growth rate of 0.43 mmol/(gDW · h) in SNM and 0.45 mmol/(gDW · h) in SCFM. These results align with the documented metabolic activity of R. mucilaginosa in CF lungs and its frequent isolation from the human nasal cavity. Notably, the observed growth rates closely resembled the flux rate predicted for biomass production in RPMI medium. Additionally, we confirmed that iRM23NL accurately represented R. mucilaginosa’s capacity for facultative anaerobic respiration. In more detail, when the oxygen uptake was turned off iRM23NL could successfully exhibit growth using alternative metabolic pathways across all tested nutritional media. When the oxygen level was decreased, the model predicted up to 68% reduction in biomass yield compared to aerobic conditions. Consequently, the remarkably lower anaerobic rates in all tested media mimic R. mucilaginosa’s inherent facultative anaerobic capabilities.
Nutrient utilization profile of R. mucilaginosa and predictive performance of iRM23NL
We experimentally characterized the metabolic phenotype of R. mucilaginosa DSM20746 using four 96-well Biolog PM microplates (Fig. 5). These high-throughput assays serve as proxies for bacterial growth by measuring cellular respiration across several conditions. Active respiration in the minimal medium is detected by the reduction of tetrazolium dye over time, indicating the utilization of the provided sole source (40). We cultivated our strain in a minimal medium supplemented with various sources, and growth was monitored over 48 h to identify suitable nutrients for the bacterium (as described in Materials and Methods). The derived OD measurements were normalized according to the average growth over replicates per plate and converted to qualitative data representing non-growth (NG) or growth (G). In total, we tested the uptake and utilization of 379 distinct carbon, nitrogen, phosphorus, and sulfur substrates. R. mucilaginosa demonstrated the ability to utilize 61 of 190 tested carbon substrates, including carboxylates, saccharides, and amino acids, while 10 of 95 were found to be viable nitrogen sources (Fig. 5 Panel B). Furthermore, out of 59 tested phosphorus sources, R. mucilaginosa exhibited a loss of metabolic activity for 28 compounds, resulting in a non-viable phenotype, while only 71.4% of all analyzed sulfur substrates supported positive growth. More specifically, 6 inorganic phosphorus (IP), 14 organic phosphorus (OP), 2 cyclic nucleoside monophosphates (cNMPs), and 9 nucleoside monophosphates (NMPs) were successfully utilized as sole phosphorus sources (Fig. 5, Panel C). The experimentally defined nutrient utilization phenome of R. mucilaginosa can be found in Fig. S1. An overview of all experimentally tested substrates, along with the assay results, can be found in Table S3. We independently confirmed the Biolog nutrient utilization data by testing the ability of DSM20746 to grow on minimal media in the presence of 10 compounds (see Materials and Methods, Fig. S2).
Fig 5
Additionally, we evaluated the predictive performance of our metabolic model by using various C-, N-, P-, and S-containing substrates. All compounds from the high-throughput phenotypic data were mapped to BiGG (28) identifiers and subsequently to iRM23NL. In total, 286 could be successfully mapped to the BiGG database. From these, 126 existed as extracellular metabolites in iRM23NL and were considered for further analysis. Model simulations were performed under aerobic conditions with the minimal medium defined in Table S2 and FBA (see Materials and Methods). An extracellular reaction was enabled for each tested substrate to force the model to use its transporters. Discrepancies between the Biolog data and the model simulations were utilized as basis for hypotheses to further improve and refine the network reconstruction. We resolved most inconsistencies via extensive literature mining and iterative gap analysis. For this, we used the organism- and strain-specific BioCyc (33) database. Throughout this process, we encountered different scenarios regarding incorrect model predictions. These included compounds present in all compartments, including the extracellular space, as well as substrates defined within the intracellular space and periplasm, with no transporter defined toward the extracellular space. If the experimental results indicated utilization of an undefined compound, we searched BioCyc (33) to find strain-specific and gene-based missing transporters or enzymatic reactions. When no organism-specific evidence was available, we sought supporting data from genomically identical species (Fig. 3). For instance, the compound 3-sulfino-l-alanine (3sala) was initially absent from any compartment in the preliminary draft model. Since no strain-specific information was available, we conducted a homology-based search using Basic Local Alignment Search Tool (BLAST) (41) to find genes with high similarity (similarity threshold: > 80%) in related species. Subsequently, we identified cysteine desulfurase (SULFCYS) along with three associated transport reactions (proton-mediated; SULFCYSpp, diffusion; SULFCYStex, and ABC transport; SULFCYSabc) that displayed a similarity over 80% with R. dentocariosa. These components were consequently incorporated into iRM23NL, resulting in the expected positive utilization phenotype. Generally, false negative or false positive predictions arise from missing or erroneous involvement of transporters, respectively. We resolved false positives by removing transport reactions lacking supporting gene evidence or adjusting their reversibility to facilitate export solely. More specifically, initial model predictions indicated that iRM23NL could not sustain growth when supplied with either l-cysteate (Lcyst) or AMP (amp) as sole sources, while Biolog assays indicated the opposite. To rectify this, we introduced the corresponding irreversible transporters (LCYStex and AMPt) and enabled their in silico utilization of these compounds. Moreover, several metabolites (e.g., phosphoenolpyruvate; pep, trimetaphosphate; tmp, hypotaurine; hyptaur, and inorganic triphosphate; pppi) which were absent from the initial draft model but exhibited positive growth in utilization assays were subsequently incorporated into the final network, leading to additional true positives predictions. All in all, over 50 transport reactions were added into the network, while 37 wrongly added enzymatic functions were removed. We also incorporated novel GPRs encoding over 60 biochemical reactions. Nevertheless, we identified approximately 20 instances where the resolution of inconsistencies necessitated the inclusion of metabolic reactions lacking associated gene evidence. For instance, to enable the utilization of l-aspartate, we introduced a transporter via diffusion from extracellular to periplasm (ASPtex), for which no associated GPR was available. Similar scenarios arose for other compounds, e.g., d-galactose, d-glucuronate, and acetate. These instances underscore knowledge gaps in the metabolism of DSM20746 that require in-depth investigation. In total, 14 carbon and nitrogen sources failed to promote growth in iRM23NL. Surprisingly, all of these sources had corresponding transport reactions iRM23NL but still remained ineffective (e.g., l-fucose, l-arabinose, and l-rhamnose) and nitrogen (l-tyrosin). We could not find further information on their transport or metabolic mechanism either in the genome annotation or the literature.
In summary, the final prediction accuracy of nutrient assimilation and utilization achieved by iRM23NL was 77% for carbon sources (MCC for PM1 = 0.52 and PM2A = 0.58), 94.4% for nitrogen sources (MCC = 0.82), 97% for phosphorus and sulfur sources (ACC = 100%; MCC = 1.0 and ACC = 92.3%; MCC = 0.82, respectively) (Fig. 6). Our model’s performance was notably increased by 40% post-comprehensive curation compared to the initial draft model. Our refinement reduced false positive predictions by 17, leaving only 3 unresolved mismatches. The most remarkable improvement was in nitrogen and phosphorus sources predictions. The high predictive accuracy indicates that core metabolic pathways and multiple catabolic routes of DSM20746 have been accurately reconstructed within iRM23NL. Consequently, the network can predict the catabolism of numerous common compounds, such as sugars and amino acids.
Fig 6
Gene essentiality predictions using iRM23NL
Given the increased percentage of gene-associated reactions (Fig. 2, Panel C) and the high predictive accuracy of the metabolic reconstruction, we employed iRM23NL further to predict exploitable single gene knockouts. For this purpose, we systematically removed each biochemical reaction from the network and optimized iRM23NL to produce biomass using FBA. To mitigate the inherent variability of the optimization algorithms, we repeated our FBA simulation 100 times. Additionally, we employed parsimonious enzyme usage flux balance analysis (pFBA), which involves solving two sequential linear optimization problems to determine the flux distribution of the optimal solution while minimizing the total sum of flux. Then, we compared the predicted growth rates before and after introducing the simulated gene deletion. The FCgr between the knocked-out and wild-type growth rates was employed as a proxy for the gene’s essentiality. We proceeded with condition-specific in silico single gene deletions. For this purpose, we utilized a minimal and a nutrient-rich medium (M9 supplemented with glucose and LB) as well as two growth media that mimic the intra-human nasal passages and the lungs of CF patients [SNM (39) and SCFM (38)]. Generally, when subjected to nutrient-limited conditions, the model predicted a higher number of genes as essential for growth, while the count of essential genes remained consistent among oxic and anoxic conditions (Fig. 7, Panel A). In total, 4 metabolic genes exhibited a partially essential effect across all tested media. This indicates that these genes promote cellular fitness, and their deletion partially impairs the bacterium’s capacity to generate biomass. These genes are the Trka family potassium uptake protein (WP_005506372.1), ribulose-phosphate 3-epimerase (WP_005507411.1), glucose-6-phosphate isomerase (WP_005508482.1), and transaldolase (WP_005509117.1). The majority of essential genes are involved in nucleotide metabolism, peptidoglycan biosynthesis, or the energy metabolism. These over-represented subsystems among the identified essential genes suggest their importance in supporting the bacterium’s respiration (Fig. S4). Nevertheless, in nutrient-poor conditions (M9 medium) genes from the biosynthesis of leucine (2-isopropylmalate synthase; WP_005508679.1 and 3-isopropylmalate dehydratase; WP_005507445.1), valine (ketol-acid reductoisomerase; WP_005508646.1 and dihydroxy-acid dehydratase; WP_005509229.1), and chorismate (shikimate kinase; WP_005508729.1 and 3-dehydroquinate dehydratase; WP_005504658.1) were found to be critical for the organism’s survival. Tables S4 and S5 list in detail the predicted essential genes, each corresponding to specific approaches employed in this study.
Fig 7
Subsequently, we conducted a protein sequence homology analysis with BLAST (42) against the human proteome to identify potential targets that could be exploited in future therapeutic strategies. For this, only genes highlighted as essential in both laboratory and synthetically defined media were considered (Fig. 7, Panel B). Overall, 35 essential genes were common in LB and M9, of which 20 common genes reported homologous counterparts in the human genome. Further analysis revealed that among these genes, five genes exhibited over 50% sequence similarity with homologous proteins although none resulted in over 80% similarity. Similarly, when iRM23NL was simulated with SCFM and SNM in both aerobic and anaerobic conditions, 45 shared genes were predicted to be essential. Homology analysis against the human genome yielded 25 genes with exhibited homology in the human genomes, with 7 demonstrating over 50% sequence similarity. For instance, genes encoding proteins such as phosphopyruvate hydratase (WP_005506838.1), CTP synthase (WP_044141843.1), and adenylosuccinate synthase (WP_005509175.1) consistently exhibited human counterparts with similarity exceeding 50% across all tested growth media and oxygen levels. Among the essential genes shared between both LB and M9, 15 of them did not have any homologous hits. The same was observed for 20 common essential genes in SCFM and SNM. Some examples of these genes include orotate phosphoribosyltransferase (WP_005507935.1), type I pantothenate kinase (WP_005505041.1), dihydroneopterin aldolase (WP_005507619.1), and pantetheine-phosphate adenylyltransferase (WP_005508106.1). A more detailed comparison can be found in Table S6.
Our in silico transposon mutant analysis using iRM23NL could serve as a basis for several research and practical applications from rational and condition-specific drug target development to biotechnological applications and metabolic engineering.
DISCUSSION
The metabolic phenome of R. mucilaginosa, a bacterium with both beneficial and pathogenic behavior, remains still largely unexplored. Investigating its metabolic traits is of great importance as it holds the potential to unveil unique capabilities, including substrate utilization, byproduct production, and contributions to host-microbe interactions. R. mucilaginosa is a versatile microbe found in humans’ oral, respiratory, and skin flora, where it coexists harmoniously. However, in immunocompromised individuals, R. mucilaginosa can act as an opportunistic pathogen, causing severe infections. Our study focuses on the metabolic aspects of R. mucilaginosa, particularly its behavior in isolated cultures. In 2019, a 17-species bacterial community model was reconstructed to simulate the polymicrobial community of the CF airways (43). This model accurately predicted the abundance of specific bacteria within patients’ CF lung communities by linking metabolomics and 16S rRNA gene sequencing data. However, studying a bacterium’s metabolism and genotype-phenotype relationships in monoculture provides a more controlled knowledge base. This allows for the precise manipulation of variables, enhancing our understanding of its individual traits, genetic makeup, metabolic pathways, and responses to stimuli (22, 23). Moreover, one can elucidate the bacterium’s unique contributions to nutrient uptake, substrate production, and growth dynamics, crucial for understanding its role in a broader ecosystem. Monoculture studies identify key genes and pathways, revealing how the bacterium functions autonomously. Such analysis serves as a valuable reference, differentiating inherent characteristics from those influenced by external interactions. To this end, we empirically analyzed the metabolic phenome of R. mucilaginosa DSM20746 and developed the first high-quality strain-specific GEM of R. mucilaginosa, called iRM23NL. We considered literature and database organism-specific evidence to manually gap-fill the model and include highly relevant biochemical reactions. Phylogenetic analysis of further Rothia species provided insights into the relationship and genetic diversity between these species and was utilized to extend the metabolic network’s completeness. Our model is simulation-ready, follows strictly community standards (25), and exhibits a high content quality memote score.
R. mucilaginosa is primarily aerobic and can perform oxic respiration by efficiently generating energy in the form of adenosine triphosphate (ATP) (1). However, when oxygen is limited or absent, R. mucilaginosa switches to anaerobic metabolism, which may involve fermentation or other alternative pathways to generate energy. As already mentioned, R. mucilaginosa has been previously found to be metabolically active in CF lungs where the oxygen levels are notably restricted (16). This indicates that the bacterium undergoes metabolic shift and can survive in microaerophilic environments. Various ROS products emerge as byproducts in the bacterial response to the fluctuating oxygen levels (34). In more detail, the cascade of ROS is initiated by the formation of O2− upon univalent oxygen reduction within the electron transport chain (ETC). Extreme oxygen fluctuations may be lethal and can ultimately damage cellular structure. The detoxifying pathway includes the enzymes superoxide dismutase (SOD), catalase, and peroxidase that break down lethal radicals to water and oxygen enabling the cell to neutralize the oxidative stress (44) (see Fig. 4). However, the exact anaerobic respiration mechanism of R. mucilaginosa must be thoroughly examined in experimental settings.
Since R. mucilaginosa’s metabolic behavior and adaptability are mainly yet unknown, we started by testing its growth behavior in various nutrient media. Exploring how bacteria react to various growth conditions within the human body is pivotal for understanding diseases and developing effective treatments. Moreover, they are essential for evaluating their evolution and adaptation to different environmental conditions, leading to new ecological niches in which the bacterium could be metabolically active. We ultimately validated iRM23NL using our growth kinetics data in various growth media. Overall, iRM23NL’s predictions were in line with the experimental observations. R. mucilaginosa demonstrated higher experimental growth in nutrient-rich media. The model successfully simulated growth for most media, while no biomass production was achieved in the M9 pure medium. When comparing LB to RPMI, the simulated growth rate was higher in LB, while the empirical growth in RPMI was twice as high as that in LB. This can be attributed to the fact that computer models cannot mimic the entire experimental settings and lack kinetic parameters. As of September 2023, bacteria like S. aureus, B. subtilis, and E. coli have been extensively researched for decades, with hundreds of thousands of PubMed (45) entries since the early 1990s. In contrast, R. mucilaginosa’s scientific prominence only began in the 21st century, with only 423 publications to date, indicating significant knowledge gaps crucial for metabolic reconstructions. More specialized BOF would enhance the predictive power and would reflect a more organism-specific metabolism. Therefore, this scarcity underscores the urgent need for further research efforts to uncover the hidden facets of R. mucilaginosa’s metabolism and its significance. Notably, to simulate in silico growth in RPMI and SCFM media, six metal ions needed to be supplemented. These metals have also been confirmed as essential for the in silico growth of S. aureus in RPMI (41). According to the model’s predictions RPMI, supplementation with manganese, zinc, and molybdate was required. Transition metals could be highly toxic; however, in controlled levels are important in the survival of all living organisms (46). For instance, they are involved in redox catalysis, needed for energy production through respiration, and in non-redox catalysis, necessary for many biosynthetic and metabolic processes. Additionally, transition metals are required for the activity of many enzymes, including those involved in genomic replication and repair and nitrogen fixation. However, since these compounds were absent from the providers’ medium formulation for RPMI, we speculate that the provided medium definition may not be exact. In all cases, the suggested metal co-factor promiscuity in R. mucilaginosa by iRM23NL needs to be examined to shed light on whether the bacterium could survive in the absence of one of the suggested metals.
Moreover, we experimentally characterized the strain’s ability to assimilate and utilize substrates using high-throughput phenotypic microarray assays. The utilization of various nitrogen sources did not result in active respiration, indicating that the bacterial genome lacks genes encoding for respective transporters. We used the phenotypic results to validate and extend our metabolic reconstruction, iRM23NL. Inconsistencies between the model and the phenotypic microarray results served as a basis for further model refinement. We enriched the model with missing transport reactions and their respective GPRs by referring to the organism- and strain-specific BioCyc (47) registry and the General Feature Format (GFF) annotation file. All in all, characterizing and determining the repertoire of nutrient sources a strain can use or assimilate is a critical factor of pathogenesis. It provides valuable insights into how pathogens adapt to host environments and evade host defenses. Our transporter-augmented model reflects a high accuracy degree with the experimental data regarding using carbon, nitrogen, phosphorus, and sulfur sources. Discrepancies between computational and empirical results highlight areas of current uncertainty knowledge regarding the metabolism of R. mucilaginosa. They could be attributed to non-metabolic factors that fall beyond the metabolic models’ scope, including regulatory processes, gene expression, and signaling pathways. However, targeted experiments are needed to fill the remaining network gaps and reveal novel enzymatic processes.
Considering the predictive precision of our metabolic reconstruction, we utilized iRM23NL to derive novel hypotheses. We examined the effects of condition-specific single gene knockouts on the bacterial capacity to produce biomass. Gene essentiality analysis is inherently contingent upon specific conditions. In the context of constraint-based metabolic modeling, a plethora of constraints are established, with the availability of nutrients, often the growth medium, being the most prevalent. By altering the availability of these nutrients, the environmental conditions are modified, consequently exerting a profound influence on the metabolic state and growth of an organism (48, 49). However, the true strength and versatility of GEMs lie in their ability to rapidly generate condition-specific hypotheses on a large scale, circumventing the need for labor-intensive and expensive screenings that may not always yield direct success. Various models, spanning organisms like A. baumannii, E. coli, S. cerevisiae, P. falciparum, and P. aeruginosa, demonstrated predictive accuracies ranging from 72% to 93% (50–60). Additionally, gene essentiality analysis has been instrumental in identifying potential drug targets for diseases such as cancer and viral infections, aligning well with both in vitro (61–63) and in vivo (64) data. Therefore, we utlized our GEM and created a high-throughput in silico-derived transposon mutant library considering two standard growth media, LB and M9, along with two growth media formulated to mimic the environment within the human body, SNM and SCFM. In this regard, we identified putative essential and partially essential genes and assessed their potential vulnerability under varying nutritional environments. With this, we opted for detecting candidate genes that could be considered in future antimicrobial and anti-inflammatory strategies in immunocompromised and CF patients. With this, we opted for identifying candidate genes for future research that hold promise for experimental validation. Determining which essential genes have human counterparts is of great importance for antibiotic drug development, as it helps assess potential side effects and cross-species effects on human genes targeted by antibiotics. Moreover, it provides insights into the molecular mechanisms of host-pathogen interactions, explaining how pathogens manipulate host cells and evade the immune system. Utilizing our GEM offers promising venues for future targeted engineering strategies without the need for laborious large-scale screening of knockouts and mutant libraries. This methodology would facilitate the rapid design of metabolic gene knockout strains by eliminating the associated reaction(s) from the model. Finally, CF lungs represent a highly dynamic environment (65, 66). However, GEMs are adaptable and can be tailored to reflect the metabolic capabilities of bacteria across diverse environmental conditions.
The main objective in our endeavor to combat R. mucilaginosa as an opportunistic pathogen causing infections (7) is identifying essential genes, particularly those without human counterparts. Determining these essential genes is crucial as we aim to neutralize the pathogen without harming the host. Simultaneously, we are exploring R. mucilaginosa as an agent with anti-inflammatory properties (18). In this context, we opt for promoting Rothia’s growth, focusing on modulating the environmental conditions that have been reported to do so. Once the key pathways involved in the beneficial functions of R. mucilaginosa are known, our gene essentiality predictions can be exploited to boost activation of these pathways. Nonetheless, being aware of R. mucilaginosa-specific essential genes is crucial to avoid inadvertently targeting them during therapeutic treatments, ensuring both the bacterium’s growth and its anti-inflammatory activities. With this dual perspective, we indicate R. mucilaginosa’s therapeutic variety, including developing strategies to combat the bacterium, when it is detrimental while increasing cell biomass production when its anti-inflammatory properties are beneficial. The latter could benefit human health in the context of cystic fibrosis. However, these model-driven hypotheses need to be extensively validated via in vitro and in vivo experiments.
Altogether, creating a genome-scale metabolic network for R. mucilaginosa reveals insights that would have been resource-intensive to acquire using traditional wet-lab means. Understanding the metabolic complexities of R. mucilaginosa is essential for expanding our basic understanding of bacterium’s microbiology and would benefit various practical applications. In medicine, it could facilitate the development of strategies to deal with caused infections, while in biotechnology, it would allow us to use its metabolic abilities for bioprocessing and bioengineering purposes. Hence, our high-quality metabolic network, iRM23NL, could provide a systematic and detailed framework for analyzing R. mucilaginosa’s metabolism, yielding valuable insights with broad-reaching impacts.
MATERIALS AND METHODS
Experimental settings
Bacterial strain and growth conditions
The R. mucilaginosa DSM20746 (ATCC 25296) used for the experimental work in this study is a type strain, and it was purchased from the American Type Culture Collection (ATCC, US). To create an inoculum, the bacterium was streaked onto nutrient agar (NA, Neogen, Heywood, UK) plates from a cryopreserved glycerol stock stored at −80°C using a sterile loop. Subsequently, the plates were incubated at 37°C for 48 h to form colonies (pure cultures). It is important to note that each biological replicate was conducted using pure cultures derived from the initial frozen stock (no sub-culturing). This ensures maintaining the genetic and phenotypic characteristics of the strain without introducing any potential mutations or adaptations.
Growth kinetics protocol
R. mucilaginosa overnight liquid cultures were prepared by adding bacterial colonies from pure cultures to 5 mL BHI (Neogen, Heywood, UK) and were put at 37 mL in a shaking incubator for 24 h. The initial OD was assessed and, if necessary, adjusted via up-concentration or dilution to achieve OD590nm = 0.25. Then, the bacterial suspension was subjected to centrifugation at 10,000 RPM for 5 min, and the resulting pellet was re-suspended in the medium of interest at a dilution of 1:10. Ultimately, the inoculated growth media were transferred to a sterile 96 well-plate, including three technical replicates for each tested condition together with their corresponding control conditions (sterile growth media). The outer wells were filled with milliQ water (MQ) to prevent evaporation. The respective OD590nm was measured aerobically at three distinct time points (0, 24, and 48 h) using an EnVision microplate reader (Perkin Elmer, Waltham, MA, USA). The microplates were incubated at 37°C during the interim periods between measurements. The final growth curves were generated for three biological replicates (n = 3) for the following growth media: BHI (baseline medium), LB (Neogen, Heywood, UK), M9 pure, RPMI medium (RPMI-1640 Sigma-Aldrich), and TSB (Neogen, Heywood, UK). In the M9 pure medium, only salts were considered. For detailed information regarding the constitution of M9, see Table S1. The rest of the media were prepared according to the providers’ instructions.
The raw data were normalized by subtracting the blank values from the measured ODs and were summarized by calculating the arithmetic mean across all replicates. To interpret the growth of bacterial cells in all tested media and compare their growth characteristics, we employed the FCOD ratio, which is defined as follows:
(1)
In this context, we define FCOD below 1.4 as no growth, while FCOD ratios greater than 1.4 indicate a growth increase over time. This FCOD threshold was chosen by analyzing experimental data and growth curves. A value of 1.4 was selected, considering the range of calculated growth rates. Statistical tests, as described below, validated the threshold’s reliability in accurately discerning growth from no growth in the bacterial cultures.
Phenotypic microarray screenings
DSM20746 was tested for utilizing multiple carbon, nitrogen, phosphorus, and sulfur sources. Biolog Phenotype Microarrays (PM, Hayward, CA, USA) were employed to test the utilization of 190 carbon (PM1 and PM2A), 95 nitrogen (PM3B), 59 phosphorus (PM4A), and 35 sulfur sources (PM4A). These assays use a tetrazolium redox dye to enable a colorimetric detection of active cell respiration across different nutrient sources (40). Normal cell respiration is indicated by the formation of a purple color as a result of the reduction of the colorless dye during incubation.
The PM plates were prepared following the manufacturer’s protocol for Gram-positive bacteria. Table 1 lists the assay set up for of PM plates. However, modifications were made during the cell suspension preparation. The strain was grown on nutrient agar plates without undergoing sub-culturing. Using an inoculation loop, individual colonies were picked and suspended in an inoculating fluid (IF-0) at an absorbance of 0.0915 at 590 nm. Per the established protocol, 81% of transmittance (T) should be achieved. Given our measurement of OD, the subsequent conversion of transmittance to absorbance was carried out employing the following formula:
(2)
TABLE 1
| For 1× PM | |
|---|---|
| IF-0a GN/GP (1.2×) | 10.0 |
| Dye mix (100×) | 0.12 |
| PM additive (12x) | 1.0 |
| 81%T cell suspension | 0.88 |
| Total volume | 12.0 |
a
Volumes are expressed in μL. The provided volume quantities are adequate for inoculating the specified number of plates in this study, using 100 μL/well with an additional excess.
In each well of a 96 well-plate, we introduced 100 μL of cell suspension, followed by a 48 h incubation period at 37°C. Bacterial growth was measured by the OD at 590 nm using a victor Nivo Multimode microplate reader. Each PM plate was tested in duplicate.
The subsequent analysis of the acquired data included calculating the arithmetic mean across all technical and biological replicates for all measured n time points. Background noise was also removed, and the data were normalized by subtracting the blank values from the actual measurements. The area under curve (AUC) was used to distinguish between growth (AUC ≥ 50) and no growth (AUC < 50). The computation of the AUCs was carried out by leveraging the linear trapezoidal rule that expresses the interpolation between data points:
(3)
where ti is the respective measured time point and . More specifically, the trapezoidal rule is iteratively applied to adjacent data points defined along the curve whose summation resulted in the final AUC value. Hence for n measured data points, the final AUC value is defined as follows:
(4)
Finally, we repeated this across the spectrum of tested compounds within the microarray plates.
Independent confirmatory testings of Biolog data
To independently confirm the Biolog data, we applied the growth kinetics protocol described above to 10 compounds. Although the base inoculating fluid (IF) used for the metabolic PM plates is proprietary, it is considered to reflect a minimal medium composed mainly of salts and buffers (40, 67). Hence, we used the M9 pure medium supplemented with different substrates to perform the independent tests (Fig. S2). The following compounds were examined: α-d-glucose, d-mannose, adonitol, l-ornithine, l-methionine, salicin, succinate, l-alanine, l-malate, and l-histidine. We also included negative controls of substrates with the Biolog inoculation fluid zero (IF-0). To ensure accuracy, triplicates were carried out for each tested compound. The M9 pure medium and the exact concentrations of added substrates are described in Table S1. All bacterial cell suspensions were prepared in 1:10 dilutions, and the ODs590nm were measured at 0, 24, and 48 h using an EnVision microplate reader (Perkin Elmer, Waltham, MA, US) and the associated software package.
We computed the arithmetic mean across the three replicates from the collected data set for each measured time point. Additionally, we performed a background correction to mitigate the influence of background noise or unwanted signal interference present in the measured ODs.
Statistical hypothesis analysis
We conducted statistical tests to evaluate the chosen threshold and potential statistically significant differences between measurements at the initial and final time-points, thereby indicating the significant growth or no growth. Specifically, we employed the Student’s t-test for each experimental condition, taking into account the data from the three biological replicates. The null hypothesis is formulated as following: there is no significant difference between the measured OD values in starting and end time-points. Prior to hypothesis testing, we checked the correctness of associated assumptions. More specifically, we assessed data normality through the Shapiro–Wilk test and verified the homogeneity of variances using the Levene’s test.
Computational framework and modeling methodology
Phylogenomic analysis
We supported the gap-filling process using evidence of closely related species within the Rothia genus. Employing ANIclustermap v.1.1.0 (68), we conducted a comprehensive genomic comparison involving R. mucilaginosa DSM20746 and 12 distinct Rothia species: R. koreensis, R. kristinae, R. santali, R. halotolerans, R. aeria, R. dentocariosa, R. terrae, R. amarae, R. nasimurium, R. mucilaginosa, R. aerolata, R. nasisuis, and R. endophytica (see Fig. 3). In brief, ANIclustermap creates an all-vs-all genome ANI clustermap and groups microbial genomes based on their genetic similarity. ANI is a pairwise measure to classify bacterial genomes according to their genetic similarity. It is defined as the genetic similarity across all orthologous genes shared between any two genomes (69, 70). It serves as a powerful tool for distinguishing strains of the same species or closely related species.
Draft model reconstruction and curation
The proteome of R. mucilaginosa DSM20746 (GCF_000175615.1) served as the basis for reconstructing a draft metabolic network. The DSM20746 (ATCC 25296) represents a type strain obtained from the throat, and its genetic and proteomic sequences were retrieved from National Centre for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov). The genome sequence was annotated using the NCBI Prokaryotic Genome Annotation Pipeline (PGAP) (71). An initial draft model was built using CarveMe 1.5.1 (27). CarveMe uses mixed-integer linear programming (MILP) to convert a universal model into an organism-specific one by deleting metabolites and reactions with low occurrence scores within the specific organism of interest. The universal BOF might yield incorrect gene essentiality predictions for biosynthesis pathways that rely on precursors unique to Gram-positive bacteria due to the absence of specific membrane and cell wall information. Hence, we chose the specialized Gram-positive template instead of the universal one to build our model more accurately. This Gram-positive template incorporates cell wall and membrane components specific to Gram-positive bacteria in contrast to the universal biomass reaction defined in CarveMe (27). Specifically, the BOF developed for iRM23NL includes essential macromolecules such as nucleotides and amino acids, co-enzymes, and inorganic ions. Moreover, it encompasses cell wall components like lipoteichoic acids, a peptidoglycan unit, and glycerol teichoic acids. The growth-associated (GAM) energy requirements are integrated into the biomass reaction (labeled as Growth), while non-growth associated maintenance (NGAM) is explicitly considered in the model, expressed by the reaction ATPM.
We conducted an extensive two-staged iterative gap-filling to address incomplete or missing information in the metabolic model. Gaps or missing reactions can arise for various reasons, such as incomplete genome annotations or undiscovered enzymatic activities. For this purpose, we leveraged information from both the bibliome and biochemical databases, including BioCyc (47). Thus, we ensured that the model could support the growth and viability of the organism under specific conditions.
We applied our previously published pipeline (26) to curate further the model based on community standards. The pipeline consists of eight steps, from which five (step 3–step 4) are related to model curation and ensure a high quality of the final model. In Summary, ModelPolisher (34) and SBOannotator (35) were employed to enrich the model with multiple cross-references, while the mass- and charge-unbalanced reactions were fixed. Further annotations integrated into the model encompassed: Evidence and Conclusion Ontology (ECO) terms representing the confidence level and the assertion method (biological qualifier: BQB_IS_DESCRIBED_BY), KEGG (30) subsystems as groups:member (biological qualifier: BQB_OCCURS_IN), and gene annotations. The latter was done by mapping the gene locus tags to the old tags using the GenBank GFF (72). Finally, we checked the presence of potential EGCs that could bias the final predictions (73). To manipulate the model structure, we employed the libSBML library (74).
The SBML Validator from libSBML (74) was used to assure a correct syntax of the model, while the quality control was carried out employing memote (36). However, it is worth noting that, as we discussed in our previous publication, memote considers only the parent nodes of the SBO directed acyclic graph excluding their respective children (26). Hence, memote was used carefully and not as an absolute quality indicator.
Linear programming: formulation of assumptions and constraints
FBA is used to determine the flux distribution through optimization of the objective function, typically the maximization of biomass production rate, under steady-state conditions (21). To address the under-determined nature of the system, constraints are imposed to define an allowable solution space that aligns with cellular functions. These constraints, encompassing mass balance, thermodynamics, and capacity, contribute to the FBA maximization problem. The linear programming problem used to obtain growth rates is described as follows:
(5)
where is the vector of fluxes within the network, S is the stoichiometric matrix, Z is the linear objective function, is the vector of coefficients, and I represents an index set containing the indices of all irreversible reactions. The dimensionality of vector matches the number of reactions, denoted as n in the system, and is consistent with the n columns in the matrix S.
The unit mmol/(gDW · h) is utilized to denote the predicted growth rates since the biomass consistency (rate at which biomass is produced per unit of dry weight in the model) of iRM23NL approaches 1 mmol/(gDW · h) (36, 75, 76). Consequently, direct comparisons can be made between experimentally observed and predicted growth rates. These conversions maintain their validity under the assumptions of constant volume and steady-state inherent in constraint-based modeling (75).
Bacterial growth analysis and nutrient utilization assays
Bacterial cell growth within various media and multiple substrate utilization evaluations were determined by solving Equation 5. The medium and the nutrient source of interest defined additional constraints. To achieve this objective, adjustments were made to the upper and lower limits of exchange reactions, as appropriate. We set specific uptake rates for key components within the growth medium as follows: the uptake rate of transition metals was set at 5.0 mmol/(gDW · h), the uptake rate of oxygen under aerobic conditions was established at 20.0 mmol/(gDW · h), and the rest media components equal to 10.0 mmol/(gDW · h). As previously mentioned, the M9 pure medium was used for the substrate utilization assays. Only substrates present in the metabolic network as intra- or extracellular metabolites were considered for the in silico validation. The results from the experimental and the in silico growth tests were categorized into “growth” (G) or “non-growth” (NG). Here, “growth” indicates the network’s ability to generate biomass and, therefore, a positive growth rate. The model’s overall prediction performance was assessed using the following statistical parameters: overall agreement (ACC):
(6)
and Matthews Correlation Coefficient (MCC):
(7)
where true negative (TN) and true positive (TP) represent accurate predictions, and false negative (FN) and false positive (FP) indicate incorrect predictions. Inconsistencies were resolved via iterative manual network gap-filling. For all FBA simulations, we employed the Constraints-Based Reconstruction and Analysis for Python (COBRApy) (77) package. All growth media definitions are listed in Table S2.
Gene lethality analysis
The in silico single-gene knockouts were performed as described in our previous study using FBA (26). To address the degeneracy issue of optimization, we additionally ran our FBA simulations in a total of 100 independent runs. Furthermore, we utilized pFBA, a method that allows us to ascertain the flux distribution of the optimal solution while concurrently minimizing the overall flux sum (78). The results were categorized as either essential FCgr = 0, inessential (FCgr = 1), or partially essential (0 < FCgr < 1), where FCgr denotes the FC bacterial growth rate before and after deletion (27). Shared essential genes between FBA and pFBA, as well as all tested conditions, were further aligned against the human genome using BLAST (42).
ACKNOWLEDGMENTS
N.L. is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—EXC 2124-390838134, and by the Cluster of Excellence ‘Controlling Microbes to Fight Infections’ (CMFI). A.D. is supported by the German Center for Infection Research (DZIF, doi: 10.13039/100009139) within the Federal Ministry of Education and Research (BMBF, German Centers for Health Research of the Federal Ministry of Education and Research), grant no. 8020708703. The authors acknowledge the support of the Open Access Publishing Fund of the University of Tübingen (https://uni-tuebingen.de/en/216529). The authors also thank Charlotte Rigauts and Anouk Van Hauwermeire for assistance with the Biolog assays.
SUPPLEMENTAL MATERIAL
Figure S1 - spectrum.04006-23-s0001.pdf
Experimentally-derived catabolic phenome of R. mucilaginosa DSM20746.
- Download
- 3.75 MB
Figure S2 - spectrum.04006-23-s0002.pdf
Growth curves of the independent confirmatory tests for validating the Biolog PM results.
- Download
- 186.86 KB
Figure S3 - spectrum.04006-23-s0003.pdf
Detailed comparative analysis of gene essentiality in silico predictions using iRM23NL.
- Download
- 139.86 KB
Figure S4 - spectrum.04006-23-s0004.pdf
Distribution of essential genes in metabolic subsystems.
- Download
- 65.19 KB
Table S1 - spectrum.04006-23-s0005.xlsx
Reconstitution of 10X M9 medium and amount of substrates supplemented.
- Download
- 14.79 KB
Table S2 - spectrum.04006-23-s0006.xlsx
Detailed definition of growth media used for in silico simulations.
- Download
- 15.55 KB
Table S3 - spectrum.04006-23-s0007.xlsx
Summary of high-throughput Biolog PM results.
- Download
- 24.74 KB
Table S4 - spectrum.04006-23-s0008.xlsx
List of essential genes predicted by FBA in different nutritional media.
- Download
- 29.22 KB
Table S5 - spectrum.04006-23-s0009.xlsx
List of essential genes predicted by pFBA in different nutritional media.
- Download
- 25.59 KB
Table S6 - spectrum.04006-23-s0010.xlsx
Comparative analysis of predicted essential genes for R. mucilaginosa in different nutrient environments.
- Download
- 12.21 KB
OPEN PEER REVIEW - reviewer-comments.pdf
An accounting of the reviewer comments and feedback.
- Download
- 553.68 KB
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.
REFERENCES
1.
Bergan T, Kocur M. 1982. Stomatococcus mucilaginosus gen.nov., sp.nov., ep. rev., a member of the family Micrococcaceae. Int J Syst Bacteriol 32:374–377.
2.
Collins MD, Hutson RA, Båverud V, Falsen E. 2000. Characterization of a Rothia-like organism from a mouse: description of Rothia nasimurium sp. nov. and reclassification of Stomatococcus mucilaginosus as Rothia mucilaginosa comb. nov. Int J Syst Evol Microbiol 50:1247–1251.
3.
Olsen I, Preza D, Aas JA, Paster BJ. 2009. Cultivated and not-yet-cultivated bacteria in oral biofilms. Microb Ecol Health Dis 21:65–71.
4.
Guglielmetti S, Taverniti V, Minuzzo M, Arioli S, Stuknyte M, Karp M, Mora D. 2010. Oral bacteria as potential probiotics for the pharyngeal mucosa. Appl Environ Microbiol 76:3948–3958.
5.
Wos-Oxley ML, Plumeier I, von Eiff C, Taudien S, Platzer M, Vilchez-Vargas R, Becker K, Pieper DH. 2010. A poke into the diversity and associations within human anterior nare microbial communities. ISME J 4:839–851.
6.
Wilbert SA, Mark Welch JL, Borisy GG. 2020. Spatial ecology of the human tongue dorsum microbiome. Cell Rep 30:4003–4015.
7.
Dastager SG, Srinivasan Krishnamurthi NR, Dharne M. 2014. The family Micrococcaceae. In The prokaryotes
8.
Janek D, Zipperer A, Kulik A, Krismer B, Peschel A. 2016. High frequency and diversity of antimicrobial activities produced by nasal Staphylococcus strains against bacterial competitors. PLoS Pathog 12:e1005812.
9.
Uranga CC, Arroyo P, Duggan BM, Gerwick WH, Edlund A. 2020. Commensal oral Rothia mucilaginosa produces enterobactin, a metal-chelating siderophore. mSystems 5:e00161-20.
10.
Tunney MM, Field TR, Moriarty TF, Patrick S, Doering G, Muhlebach MS, Wolfgang MC, Boucher R, Gilpin DF, McDowell A, Elborn JS. 2008. Detection of anaerobic bacteria in high numbers in sputum from patients with cystic fibrosis. Am J Respir Crit Care Med 177:995–1001.
11.
Lim YW, Schmieder R, Haynes M, Willner D, Furlan M, Youle M, Abbott K, Edwards R, Evangelista J, Conrad D, Rohwer F. 2013. Metagenomics and metatranscriptomics: windows on CF-associated viral and microbial communities. J Cyst Fibros 12:154–164.
12.
Worlitzsch D, Tarran R, Ulrich M, Schwab U, Cekici A, Meyer KC, Birrer P, Bellon G, Berger J, Weiss T, Botzenhart K, Yankaskas JR, Randell S, Boucher RC, Döring G. 2002. Effects of reduced mucus oxygen concentration in airway Pseudomonas infections of cystic fibrosis patients. J Clin Invest 109:317–325.
13.
Filkins LM, O’Toole GA. 2015. Cystic fibrosis lung infections: polymicrobial, complex, and hard to treat. PLoS Pathog 11:e1005258.
14.
Guss AM, Roeselers G, Newton ILG, Young CR, Klepac-Ceraj V, Lory S, Cavanaugh CM. 2011. Phylogenetic and metabolic diversity of bacteria associated with cystic fibrosis. ISME J 5:20–29.
15.
Bittar F, Richet H, Dubus J-C, Reynaud-Gaubert M, Stremler N, Sarles J, Raoult D, Rolain J-M. 2008. Molecular detection of multiple emerging pathogens in sputa from cystic fibrosis patients. PLoS One 3:e2908.
16.
Lim YW, Schmieder R, Haynes M, Furlan M, Matthews TD, Whiteson K, Poole SJ, Hayes CS, Low DA, Maughan H, Edwards R, Conrad D, Rohwer F. 2013. Mechanistic model of Rothia mucilaginosa adaptation toward persistence in the CF lung, based on a genome reconstructed from metagenomic data. PLoS One 8:e64285.
17.
Gao B, Gallagher T, Zhang Y, Elbadawi-Sidhu M, Lai Z, Fiehn O, Whiteson KL. 2018. Tracking polymicrobial metabolism in cystic fibrosis airways: Pseudomonas aeruginosa metabolism and physiology are influenced by Rothia mucilaginosa-derived metabolites. mSphere 3:10–1128.
18.
Rigauts C, Aizawa J, Taylor SL, Rogers GB, Govaerts M, Cos P, Ostyn L, Sims S, Vandeplassche E, Sze M, Dondelinger Y, Vereecke L, Van Acker H, Simpson JL, Burr L, Willems A, Tunney MM, Cigana C, Bragonzi A, Coenye T, Crabbé A. 2022. Rothia mucilaginosa is an anti-inflammatory bacterium in the respiratory tract of patients with chronic lung disease. Eur Respir J 59:2101293.
19.
Passi A, Tibocha-Bonilla JD, Kumar M, Tec-Campos D, Zengler K, Zuniga C. 2021. Genome-scale metabolic modeling enables in-depth understanding of big data. Metabolites 12:14.
20.
Palsson BØ. 2006. Systems biology: properties of reconstructed networks. Cambridge university press.
21.
Orth JD, Thiele I, Palsson BØ. 2010. What is flux balance analysis? Nat Biotechnol 28:245–248.
22.
Oberhardt MA, Palsson BØ, Papin JA. 2009. Applications of genome-scale metabolic reconstructions. Mol Syst Biol 5:320.
23.
Gu C, Kim GB, Kim WJ, Kim HU, Lee SY. 2019. Current status and applications of genome-scale metabolic models. Genome Biol 20:121.
24.
Sertbas M, Ulgen KO. 2020. Genome-scale metabolic modeling for unraveling molecular mechanisms of high threat pathogens. Front Cell Dev Biol 8:566702.
25.
Carey MA, Dräger A, Beber ME, Papin JA, Yurkovich JT. 2020. Community standards to facilitate development and address challenges in metabolic modeling. Mol Syst Biol 16:e9235.
26.
Leonidou N, Xia Y, Friedrich L, Schütz MS, Dräger A. 2023. Exploring the metabolic profiling of A. baumannii for antimicrobial development using genome-scale modeling. bioRxiv.
27.
Machado D, Andrejev S, Tramontano M, Patil KR. 2018. Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Res 46:7542–7553.
28.
Norsigian CJ, Pusarla N, McConn JL, Yurkovich JT, Dräger A, Palsson BO, King Z. 2020. BiGG models 2020: multi-strain genome-scale models and expansion across the phylogenetic tree. Nucleic Acids Res 48:D402–D406.
29.
Lex A, Gehlenborg N, Strobelt H, Vuillemot R, Pfister H. 2014. UpSet: visualization of intersecting sets. IEEE Trans Vis Comput Graph 20:1983–1992.
30.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. 2021. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res 49:D545–D551.
31.
Zhao X, Drlica K. 2014. Reactive oxygen species and the bacterial response to lethal stress. Curr Opin Microbiol 21:1–6.
32.
Korge P, Calmettes G, Weiss JN. 2015. Increased reactive oxygen species production during reductive stress: the roles of mitochondrial glutathione and thioredoxin reductases. Biochim Biophys Acta 1847:514–525.
33.
Monk J, Nogales J, Palsson BO. 2014. Optimizing genome-scale network reconstructions. Nat Biotechnol 32:447–452.
34.
Römer M, Eichner J, Dräger A, Wrzodek C, Wrzodek F, Zell A. 2016. ZBIT bioinformatics toolbox: a web-platform for systems biology and expression data analysis. PLoS One 11:e0149263.
35.
Leonidou N, Fritze E, Renz A, Dräger A. 2023. SBOannotator: a python tool for the automated assignment of systems biology ontology terms. Bioinformatics 39:btad437.
36.
Lieven C, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, Bartell JA, Blank LM, Chauhan S, Correia K, et al. 2020. MEMOTE for standardized genome-scale metabolic model testing. Nat Biotechnol 38:272–276.
37.
Hucka M, Bergmann FT, Dräger A, Hoops S, Keating SM, Le Novère N, Myers CJ, Olivier BG, Sahle S, Schaff JC, Smith LP, Waltemath D, Wilkinson DJ. 2018. The systems biology markup language (SBML): language specification for level 3 version 1 core. J Integr Bioinform 15:20170080.
38.
Turner KH, Wessel AK, Palmer GC, Murray JL, Whiteley M. 2015. Essential genome of Pseudomonas aeruginosa in cystic fibrosis sputum. Proc Natl Acad Sci U S A 112:4110–4115.
39.
Krismer B, Liebeke M, Janek D, Nega M, Rautenberg M, Hornig G, Unger C, Weidenmaier C, Lalk M, Peschel A. 2014. Nutrient limitation governs Staphylococcus aureus metabolism and niche adaptation in the human nose. PLoS Pathog 10:e1003862.
40.
Bochner BR, Gadzinski P, Panomitros E. 2001. Phenotype microarrays for high-throughput phenotypic testing and assay of gene function. Genome Res 11:1246–1255.
41.
Seif Y, Monk JM, Mih N, Tsunemoto H, Poudel S, Zuniga C, Broddrick J, Zengler K, Palsson BO. 2019. A computational knowledge-base elucidates the response of Staphylococcus aureus to different media types. PLoS Comput Biol 15:e1006644.
42.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. J Mol Biol 215:403–410.
43.
Henson MA, Orazi G, Phalak P, O’Toole GA. 2019. Metabolic modeling of cystic fibrosis airway communities predicts mechanisms of pathogen dominance. mSystems 4:e00026-19.
44.
Baron S. 1996. Medical microbiology
45.
Wheeler DL, Chappey C, Lash AE, Leipe DD, Madden TL, Schuler GD, Tatusova TA, Rapp BA. 2000. Database resources of the national center for biotechnology information. Nucleic Acids Res 28:10–14.
46.
Hood MI, Skaar EP. 2012. Nutritional immunity: transition metals at the pathogen-host interface. Nat Rev Microbiol 10:525–537.
47.
Karp PD, Billington R, Caspi R, Fulcher CA, Latendresse M, Kothari A, Keseler IM, Krummenacker M, Midford PE, Ong Q, Ong WK, Paley SM, Subhraveti P. 2019. The BioCyc collection of microbial genomes and metabolic pathways. Brief Bioinform 20:1085–1093.
48.
Gasperotti A, Brameyer S, Fabiani F, Jung K. 2020. Phenotypic heterogeneity of microbial populations under nutrient limitation. Curr Opin Biotechnol 62:160–167.
49.
Camacho A, Rochera C, Picazo A. 2022. Effect of experimentally increased nutrient availability on the structure, metabolic activities, and potential microbial functions of a maritime Antarctic microbial mat. Front Microbiol 13:900158.
50.
Presta L, Bosi E, Mansouri L, Dijkshoorn L, Fani R, Fondi M. 2017. Constraint-based modeling identifies new putative targets to fight colistin-resistant A. baumannii infections. Sci Rep 7:3706.
51.
Norsigian CJ, Kavvas E, Seif Y, Palsson BO, Monk JM. 2018. iCN718, an updated and improved genome-scale metabolic network reconstruction of Acinetobacter baumannii AYE. Front Genet 9:121.
52.
Zhao J, Zhu Y, Han J, Lin Y-W, Aichem M, Wang J, Chen K, Velkov T, Schreiber F, Li J. 2020. Genome-scale metabolic modeling reveals metabolic alterations of multidrug-resistant Acinetobacter baumannii in a murine bloodstream infection model. Microorganisms 8:11.
53.
Zhu Y, Zhao J, Maifiah MHM, Velkov T, Schreiber F, Li J. 2019. Metabolic responses to polymyxin treatment in Acinetobacter baumannii ATCC 19606: integrating transcriptomics and metabolomics with genome-scale metabolic modeling. mSystems 4:e00157-18.
54.
Monk JM, Lloyd CJ, Brunk E, Mih N, Sastry A, King Z, Takeuchi R, Nomura W, Zhang Z, Mori H, Feist AM, Palsson BO. 2017. iML1515, a knowledgebase that computes Escherichia coli traits. Nat Biotechnol 35:904–908.
55.
Joyce AR, Reed JL, White A, Edwards R, Osterman A, Baba T, Mori H, Lesely SA, Palsson BØ, Agarwalla S. 2006. Experimental and computational assessment of conditionally essential genes in Escherichia coli. J Bacteriol 188:8259–8271.
56.
Abdel-Haleem AM, Hefzi H, Mineta K, Gao X, Gojobori T, Palsson BO, Lewis NE, Jamshidi N. 2018. Functional interrogation of Plasmodium genus metabolism identifies species- and stage-specific differences in nutrient essentiality and drug targeting. PLoS Comput Biol 14:e1005895.
57.
Oh Y-K, Palsson BO, Park SM, Schilling CH, Mahadevan R. 2007. Genome-scale reconstruction of metabolic network in Bacillus subtilis based on high-throughput phenotyping and gene essentiality data. J Biol Chem 282:28791–28799.
58.
Nanda P, Patra P, Das M, Ghosh A. 2020. Reconstruction and analysis of genome-scale metabolic model of weak Crabtree positive yeast Lachancea kluyveri. Sci Rep 10:16314.
59.
Oberhardt MA, Puchałka J, Fryer KE, Martins dos Santos VAP, Papin JA. 2008. Genome-scale metabolic network analysis of the opportunistic pathogen Pseudomonas aeruginosa PAO1. J Bacteriol 190:2790–2803.
60.
Harrison R, Papp B, Pál C, Oliver SG, Delneri D. 2007. Plasticity of genetic interactions in metabolic networks of yeast. Proc Natl Acad Sci U S A 104:2307–2312.
61.
Aller S, Scott A, Sarkar-Tyson M, Soyer OS. 2018. Integrated human-virus metabolic stoichiometric modelling predicts host-based antiviral targets against Chikungunya, Dengue and Zika viruses. J R Soc Interface 15:146.
62.
Bidkhori G, Benfeitas R, Elmas E, Kararoudi MN, Arif M, Uhlen M, Nielsen J, Mardinoglu A. 2018. Metabolic network-based identification and prioritization of anticancer targets based on expression data in hepatocellular carcinoma. Front Physiol 9:916.
63.
Agren R, Mardinoglu A, Asplund A, Kampf C, Uhlen M, Nielsen J. 2014. Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling. Mol Syst Biol 10:721.
64.
Renz A, Hohner M, Breitenbach M, Josephs-Spaulding J, Dürrwald J, Best L, Jami R, Marinos G, Cabreiro F, Dräger A, Schindler M, Kaleta C. 2023. Metabolic modeling elucidates phenformin and atpenin A5 as broad-spectrum antiviral drugs. Biol Life sci.
65.
Zhao J, Schloss PD, Kalikin LM, Carmody LA, Foster BK, Petrosino JF, Cavalcoli JD, VanDevanter DR, Murray S, Li JZ, Young VB, LiPuma JJ. 2012. Decade-long bacterial community dynamics in cystic fibrosis airways. Proc Natl Acad Sci U S A 109:5809–5814.
66.
Greenwald MA, Wolfgang MC. 2022. The changing landscape of the cystic fibrosis lung environment: from the perspective of Pseudomonas aeruginosa. Curr Opin Pharmacol 65:102262.
67.
Zhou L, Lei X-H, Bochner BR, Wanner BL. 2003. Phenotype microarray analysis of Escherichia coli K-12 mutants with deletions of all two-component systems. J Bacteriol 185:4956–4972.
68.
Shimoyama Y. 2022. ANIclustermap: a tool for drawing ANI clustermap between all-vs-all microbial genomes
69.
Konstantinidis KT, Tiedje JM. 2005. Genomic insights that advance the species definition for prokaryotes. Proc Natl Acad Sci U S A 102:2567–2572.
70.
Goris J, Konstantinidis KT, Klappenbach JA, Coenye T, Vandamme P, Tiedje JM. 2007. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int J Syst Evol Microbiol 57:81–91.
71.
Tatusova T, DiCuccio M, Badretdin A, Chetvernin V, Nawrocki EP, Zaslavsky L, Lomsadze A, Pruitt KD, Borodovsky M, Ostell J. 2016. NCBI prokaryotic genome annotation pipeline. Nucleic Acids Res 44:6614–6624.
72.
Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW. 2016. GenBank. Nucleic Acids Res 44:D67–D72.
73.
Fritzemeier CJ, Hartleb D, Szappanos B, Papp B, Lercher MJ. 2017. Erroneous energy-generating cycles in published genome scale metabolic networks: identification and removal. PLoS Comput Biol 13:e1005494.
74.
Bornstein BJ, Keating SM, Jouraku A, Hucka M. 2008. LibSBML: an API library for SBML. Bioinformatics 24:880–881.
75.
Gottstein W, Olivier BG, Bruggeman FJ, Teusink B. 2016. Constraint-based stoichiometric modelling from single organisms to microbial communities. J R Soc Interface 13:124.
76.
Chan SHJ, Cai J, Wang L, Simons-Senftle MN, Maranas CD. 2017. Standardizing biomass reactions and ensuring complete mass balance in genome-scale metabolic models. Bioinformatics 33:3603–3609.
77.
Ebrahim A, Lerman JA, Palsson BO, Hyduke DR. 2013. COBRApy: constraints-based reconstruction and analysis for python. BMC Syst Biol 7:1–6.
78.
Lewis NE, Hixson KK, Conrad TM, Lerman JA, Charusanti P, Polpitiya AD, Adkins JN, Schramm G, Purvine SO, Lopez-Ferrer D, Weitz KK, Eils R, König R, Smith RD, Palsson BØ. 2010. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol 6:390.
79.
Malik-Sheriff RS, Glont M, Nguyen TVN, Tiwari K, Roberts MG, Xavier A, Vu MT, Men J, Maire M, Kananathan S, Fairbanks EL, Meyer JP, Arankalle C, Varusai TM, Knight-Schrijver V, Li L, Dueñas-Roca C, Dass G, Keating SM, Park YM, Buso N, Rodriguez N, Hucka M, Hermjakob H. 2020. BioModels—15 years of sharing computational models in life science. Nucleic Acids Res 48:D407–D415.
80.
Bergmann FT, Adams R, Moodie S, Cooper J, Glont M, Golebiewski M, Hucka M, Laibe C, Miller AK, Nickerson DP, Olivier BG, Rodriguez N, Sauro HM, Scharm M, Soiland-Reyes S, Waltemath D, Yvon F, Le Novère N. 2014. COMBINE Archive and OMEX format: One file to share all information to reproduce a modeling project. BMC Bioinformatics 15:369.
81.
Neal ML, König M, Nickerson D, Mısırlı G, Kalbasi R, Dräger A, Atalag K, Chelliah V, Cooling MT, Cook DL, et al. 2019. Harmonizing semantic annotations for computational models in biology. Brief Bioinform 20:540–550.
Information & Contributors
Information
Published In
Microbiology Spectrum
Volume 12 • Number 6 • 4 June 2024
eLocator: e04006-23
Editor: Angela Re, Politecnico di Torino, Torino, Italy
PubMed: 38652457
Copyright
© 2024 Leonidou et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.
History
Received: 28 November 2023
Accepted: 20 March 2024
Published online: 23 April 2024
Peer Review History
Download review history as PDF.
Keywords
Data Availability
Supplementary data are available along with this article. Additionally, iRM23NL is available at the BioModels Database (79) as an SBML Level 3 Version 1 (37) file wrapped in an Open Modelling EXchange format (OMEX) (80) archive file together with a metadata file (81). Access the model at https://identifiers.org/biomodels.db/MODEL2310240001.
Contributors
Editor
Angela Re
Editor
Politecnico di Torino, Torino, Italy
Metrics & Citations
Metrics
Note:
- For recently published articles, the TOTAL download count will appear as zero until a new month starts.
- There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.
- Citation counts come from the Crossref Cited by service.
Citations
Citation
2024.
Genome-scale model of Rothia mucilaginosa predicts gene essentialities and reveals metabolic capabilities. Microbiol Spectr
12:e04006-23.
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. For an editable text file, please select Medlars format which will download as a .txt file. Simply select your manager software from the list below and click Download.