Modelling lifestyles transitions by integrating transcriptomics and growth data.

1
Introduction
In the last few years, several studies have highlighted a significant association between cancer and symbiotic bacteria. Despite our limited knowledge of the actual host-microbiota interactions, it is clear that microbes can promote or prevent cancer development (Dzutsev et al., 2017; Fang et al., 2021; Shang and Liu, 2018; Whisner and Athena Aktipis, 2019). Commensal microbiota exerts a huge influence on the host biology, modulating its metabolism, hematopoiesis, immunity and inflammation. Accordingly, alterations in its composition can result in dysbiosis and chronic inflammation conditions, increasing the risk for neoplasms formation. Examples of carcinogenesis triggered by specific infective agents, such as Helicobacter pylori (Polk and Peek, 2010), human papillomavirus (HPV) (Crosbie et al., 2013) and Salmonella enterica (Di Domenico et al., 2017), are currently well known and are estimated to contribute roughly to 20 % of worldwide cancer burden.
Although the right mechanisms linking altered microbiome to cancer beginning and progression are still to be fully comprehended, the prevalence of some bacterial species in tumor samples led to the identification of putative “drivers”, i.e. bacteria promoting the cancer development. Of these, Fusobacterium nucleatum (F. nucleatum) is emerging as one of the species with the best characterized interaction mechanisms with cancer. A gram-negative, obligate anaerobic bacteria present as a constituent of the oral microbiome (Kolenbrander, 2000; Kolenbrander et al., 2006; Segata et al., 2012) F. nucleatum has been associated with different pathologies, including gastrointestinal abscesses (Brook, 2013), acute appendicitis (Swidsinski et al., 2011), intra-amniotic infections (Han et al., 2010), and different types of cancer (including colorectal, pancreatic and oral cancer) (Gallimidi et al., 2015; Mitsuhashi et al., 2015; Warren et al., 2013). In detail, significantly higher F. nucleatum levels were found in colon tissue of colorectal cancer (CRC) patients, compared with healthy people, with F. nucleatum colonizing on colorectal tissues in a half of CRC cases (Castellarin et al., 2012; Kostic et al., 2012). The in situ association with CRC (Castellarin et al., 2012; Kostic et al., 2012) promotes tumor progression directly, by inducing the activation of different oncogenes, and indirectly, protecting the cancer from the host’s immunity and promoting chemoresistance. In all these cases, the metabolism of the host-associated microbiota has been shown to play a key role for the progression of the pathological state and/or for affecting the efficacy of cancer drugs (Huang and Mao, 2022; Ren et al., 2021).
Recently, a considerable effort was made to characterize the interactions occurring between F. nucleatum and cancer, using different in vitro and in vivo models, and to find treatments (Bullman et al., 2017; Chen, et al., 2022; De Carvalho et al., 2019; Rubinstein et al., 2013). However, considering the complexity of the intra-tumoral environment and the number of multi-scale interactions occurring between cancer, bacteria and immune system, the nature of F. nucleatum -associated tumor is far from being fully characterized. In this scenario, understanding the details of F. nucleatum physiology is strategically relevant. Among all the other cellular processes, metabolism has been shown to play a key role during F. nucleatum host infection/colonization, based on a survey of differentially expressed genes (Cochrane et al., 2020). In detail, the rewiring of specific metabolic pathways (together with non-metabolic ones) has been linked to the adaptation of invasive F. nucleatum cells to an intracellular lifestyle, where they continue to proliferate (Cochrane et al., 2020)
In this work, we aim to build a computational model of F. nucleatum integrating, physiological, gene expression and phenomics data, and to explore the phenotypic variability in relation to different nutritional environments. We reconstructed an experimentally validated genome-scale metabolic model (GEM) of F. nucleatum subsp. nucleatum ATCC 25586 (hereafter Fn-GEM) and used it to investigate metabolic network reprogramming within a conserved core network associated with distinct lifestyle-related conditions, including free-living, adhesion and invasion states. Overall, we provide here a comprehensive and multi-layer characterization of F. nucleatum, a clinically relevant microorganism associated with multiple human diseases. We further anticipate that our integrated modelling framework and the metabolic reconstruction will represent a robust platform for future investigation of F. nucleatum biology.

2
Materials and methods
2.1
Cell culturing and growth conditions
Growth of F. nucleatum subsp. nucleatum ATCC 25586 was performed on Yeast Casitone Fatty Acids (YCFA) agar or broth, a rich medium commonly used for the cultivation of anaerobic bacteria (Duncan et al., 2002; Hosomi et al., 2025; Wornell et al., 2022). The medium composition corresponds to DSMZ Medium 1611, as provided by the German Collection of Microorganisms and Cell Cultures (DSMZ). Cultures were incubated at 37°C in a Whitley DG250 workstation under anaerobic conditions (80 % nitrogen, 10 % hydrogen, and 10 % carbon dioxide). To perform growth curves, F. nucleatum was grown on YCFA plates for 2 days, and then cultured in YCFA broth for 16 hours. F. nucleatum pre-inoculum was diluted 1:10, 1:20, 1:50 or 1:100 in fresh YCFA broth in triplicate and incubated with vigorous shaking. Positive acceleration phase, exponential phase and negative acceleration phase time points for every different starting dilution were determined by fitting a sigmoidal model on each curve, and the first derivate of the fit was used to estimate growth rates. For each condition, the time point corresponding to the maximum slope was taken as the onset of exponential phase. 30 ml of YCFA broth in 125 ml Erlenmeyer flasks were used to grow F. nucleatum in each independent experimental biological replicate that we performed.

2.2
RNA extraction, library preparation and sequencing
F. nucleatum cultures were sampled at three time points corresponding to different stages of active growth within the exponential phase (T1, T2 and T3), occurring approximately three, four and five hours after incubation start, respectively. Samples were incubated with Bacterial RNAprotect (Qiagen), vortexed and incubated at room temperature for 5 min. Samples were centrifuged at 4,500 x g, at 4°C for 15 min to obtain cell pellets. Supernatant was carefully discarded, and total RNA was extracted using the following procedure: briefly, bacterial pellets were resuspended in 1 ml of TRIzol (Thermo Fisher Scientific), transferred to Lysing matrix B tubes (MP Biomedicals), and processed twice in a FastPrep-FP120 bead-beating system at 6.5 m/s for 60 s. Samples were incubated on ice for 5 min between cycles. After centrifuging the samples at 16,000 x g for 10 min at 4°C, bacterial lysates were transferred to RNase-free Eppendorf tubes and Lysing matrix B was rinsed with one volume of absolute ethanol, which was then recovered and mixed with bacterial lysates. Samples were then added on Zymo-Spin IIC Columns and total RNA was purified using Direct-zol™ RNA miniprep kit (Zymo Research) with a 30 min DNAse treatment with Turbo DNase (Invitrogen) at 37°C with Ribolock RNase inhibitor (Thermo Scientific), followed by clean-up and concentration with RNA clean and concentrator™-5 kit (Zymo Research). Samples were quantified and analysed using an Agilent 2100 Bioanalyzer on a total RNA 6000 Nano kit chip. Samples were considered for further sequencing only if RIN was greater than 8.0. The samples were sequenced on a HiSeq 2000 lane which produced paired-end libraries averaging ∼26 million reads of 75 bp per sample (see Supplementary Table 1). For each time point, five biological replicates were analysed.

2.3
RNA-Seq data analysis
Sequence quality, GC content, K-mers length and artefacts presence have been checked through FastQC. Subsequently, reads have been mapped to the genome of F. nucleatum ATCC 25586 (NZ_CP028101.1) using Bowtie2 (Langmead et al., 2019), and gene expression was quantified using HTSeq 2.0. (Putri et al., 2022) with the “intersection-nonempty” mode to resolve reads spanning multiple genes. Raw read counts were used to perform differential expression analysis with DESeq2 (Love et al., 2014), separately contrasting samples at T1 vs T2, T2 vs T3 and T1 vs T3. Normalized values of gene expression were obtained computing for each gene the number of fragments per kilobase of transcript per million mapped fragments (FPKM).

2.4
Total cellular biomass collection
F. nucleatum was grown on YCFA plates for 2 days, then inoculated in YCFA broth for 16 hours. F. nucleatum overnight culture was diluted 1:10 in fresh YCFA broth and incubated with vigorous shaking for 4 hours. Five biological experimental replicates were performed for each cellular component biomass extraction and determination. To obtain F. nucleatum biomass, exponential phase cultures were centrifuged at 4,500 x g for 10 min at 4°C. Samples were washed once in PBS (phosphate-buffered saline), and dried at 45°C for 10 min in a SpeedVac system (Eppendorf). Exponential phase cultures for RNA biomass measurements were resuspended in bacterial RNA protect (Qiagen), vortexed, incubated at room temperature for 5 min prior to centrifugation and then dried in the SpeedVac system. Biomasses were weighed and kept at -80°C. Cell wall composition of biomass and fraction of inorganic material were assumed to be similar to that of Escherichia coli K12 MG1655 (Monk et al., 2017). Before macromolecule biomass determination, bacterial pellets were thawed and lysed by bead-beating using a FastPrep-FP120 system (6.0 m s⁻¹, three cycles of 30 s, with 30 s incubations on ice between cycles). The resulting bacterial lysate were used for the quantification of individual biomass components as described in the following sections.

2.5
Bacterial DNA biomass determination
Bacterial lysates were subjected to DNA extraction by adding one volume of phenol:chloroform:isoamyl alcohol (25:24:1) (v/v) reagent to the bacterial lysates. After vortexing for 20 s, samples were centrifuged at room temperature for 5 min at 16,000 x g. Aqueous phase was carefully transferred to a fresh tube where 0.5 x volume of sample of 7.5 M ammonium acetate (Sigma) and 2.5 x volume of sample of cold absolute ethanol were added. Samples were mixed and DNA precipitated overnight at -20°C. Samples were centrifuged at 4°C for 30 min at 16,000 x g to pellet the total DNA. Supernatant was carefully discarded, and samples were washed twice in cold 70 % ethanol. DNA was then dried at 45°C for 5 min in a SpeedVac system (Eppendorf) and resuspended in 100 μl of ultrapure water by pipetting. Quantification was performed using a Qubit 2.0 fluorometer (Invitrogen) with a Qubit dsDNA BR Assay Kit (Invitrogen).

2.6
Bacterial RNA biomass determination
Bacterial lysates obtained in the presence of TRIzol reagent were centrifuged at 16,000 x g for 10 min at 4°C, and transferred to RNase-free Eppendorf tubes. Lysing B matrix was rinsed with one volume of absolute ethanol, and centrifuged. Ethanol was then recovered, added to bacterial lysate, and mixed well. Samples were then added on Zymo-Spin IIC Columns and total RNA was purified using Direct-zol™ RNA miniprep kit (Zymo Research) with a 30 min DNAse treatment with Turbo DNase (Invitrogen) at 37°C with Ribolock RNase inhibitor (Thermo Scientific), followed by clean-up and concentration with RNA clean and concentrator™-5 kit (Zymo Research). Samples were quantified and analysed using a Qubit 2.0 fluorometer (Invitrogen) with a Qubit RNA BR Assay Kit (Invitrogen).

2.7
Bacterial protein biomass determination
Bacterial lysates were kept on ice and total protein content was assessed using Pierce™ BCA Protein Assay Kit (Thermo Scientific) following manufacturer instructions. Measurements were performed in triplicates.

2.8
Bacterial lipid biomass determination
Bacterial lysates underwent methanol-chloroform extraction as described by Izard and colleagues with minor modifications. Briefly, samples were transferred to a fresh tube. Methanol and chloroform were added to the sample with a ratio 2.5:1.25:1 (v/v), and used first to rinse Lysing matrix B, then recovered and finally added to the sample to remove all remaining bacterial lysate from the bead-beating tube. After a 10 min incubation at 4°C, lipids were separated from the water-soluble material by dilution of the extract with one volume of chloroform followed by one volume of water. Samples were then centrifuged for 15 min at 3000 x g. The chloroform layer was gently removed and dispensed in pre-weighed fresh tubes. Samples were allowed to evaporate overnight under a chemical hood and weights were recorded the following day.

2.9
Bacterial glycogen biomass determination
Glycogen present in bacterial lysates was digested to glucose using 2 U/ml of amyloglucosidase (Sigma-Aldrich) at 37°C for 30 min. Quantification of glucose was then carried out using the Glucose (HK) Assay Kit (Sigma-Aldrich). Measurements were performed in triplicates.

2.10
Biomass F. nucleatum determination
The general biomass reaction drafted through KBase was replaced by a curated one, specific for this strain, which was calculated in BOFdat (Lachance et al., 2019), based on to the experimental data obtained through the previously described wet-lab methods. This curated biomass reaction was incorporated into the reference Fn-GEM and was subsequently used unchanged in all downstream simulations (see Section 2.12).

2.11
Genome-scale metabolic reconstruction of F. nucleatum constraint-based modelling
A preliminary draft reconstruction of F. nucleatum metabolism was obtained using the KBase narrative interface (http://kbase.us). This draft model was then used for a bottom-up approach where each of the central metabolic processes (TCA cycle, Glycolysis, Pentose phosphate pathway, Fatty acids metabolism, Amino acids metabolism, etc.) was manually inspected and screened for missing GPRs in the initial draft reconstruction. Overall, we followed the main steps listed in the Thiele and Palsson protocol (Thiele and Palsson, 2010). Existing experimental data and online databases (Kanehisa et al., 2017; King et al., 2016; Moretti et al., 2016) were exploited during the curation phase. Pathways were especially cross-checked with KEGG (Kanehisa et al., 2007) and specific metabolic capabilities were added to the model on the base of literature references. Additional GPRs were added to the reconstruction on the basis of orthology relationships shared with E. coli, using one of the latest available GEM for this bacterium (Monk et al., 2017). Transporters were added based on data retrieved on TransportDB 2.0 (Elbourne et al., 2017).
All the reactions and the metabolites in the model were renamed after BiGG IDs (King et al., 2016). In case a given metabolite or reaction was not listed on BiGG we adopted alternative nomenclature systems, such as KEGG or Seed IDs (Devoid et al., 2013).
At this stage the reconstruction embedded a total of 452 genes, 1041 reactions and 938 metabolites. This reconstruction was subsequently formalized as a Fn-GEM (i.e., a mathematical representation of metabolism), and analysed through constrained-based metabolic modelling in COBRApy (Ebrahim et al., 2013) and COBRA (Schellenberger et al., 2011) environments.

2.12
Transcriptomic data mapping and flux analysis
Transcriptomic data (RPKM) from three experimental conditions (control, adhesion, and invasion), obtained for F. animalis strain 7_1 (Cochrane et al., 2020), were used as a functional input to interrogate metabolic network reprogramming within a conserved core network captured by the Fn-GEM. Accordingly, the aim of transcriptomic integration was not to reproduce strain- or species-specific expression patterns, but to contextualize conserved metabolic pathways underlying distinct physiological states (control, adhesion, and invasion). To ensure the biological validity of this cross-taxon integration, we quantitatively assessed the conservation of metabolic genes between the two subspecies prior to contextualization. Specifically, a comparative genomics analysis was performed. Reciprocal best-hit (RBH) BLASTp searches were conducted between the two complete proteomes using an e-value threshold of 1e−10. High confidence orthologs were defined by requiring ≥30 % sequence identity, ≥70 % query and subject coverage, and a minimum alignment length of 100 amino acids. Orthology relationships were further intersected with the gene content of the curated Fn-GEM to quantify the fraction of model-associated genes supported by F. animalis. Functional conservation was assessed by annotating orthologous genes using eggNOG-mapper v 2.1.12 (Cantalapiedra et al., 2021) and analysing their distribution across COG functional categories.
Prior to expression data integration, Fn-GEM was tested in a maximal growth medium consistent with the experimental setup from which the transcriptomic data were obtained. The rich medium was simulated by allowing unlimited uptake of all available carbon sources. Expression data were then mapped onto the Fn-GEM in each experimental condition using the Gene Inactivity Moderated by Metabolism and Expression (GIMME) algorithm (Becker and Palsson, 2008), generating three context-specific metabolic models. Genes with expression values below the selected threshold were treated as inactive, and their associated reactions were initially removed from model, preventing them from carrying flux. These reactions were reintroduced only when required to satisfy the imposed metabolic objective, in which case their inclusion contributed to the GIMME inconsistency score.
For each context-specific model generated, flux sampling analysis was performed using the OptGPSampler algorithm implemented in COBRApy (Megchelenbrink et al., 2014). A total of 1000 flux samples per condition were generated using a thinning factor of 500 and with a fixed random seed. The resulting output consisted of tables where each column represented a reaction in the model, and each row corresponded to a resampled flux value from the analysis. The data was subjected to two filtering steps: i) reactions with a total sum of sampled fluxes equal to zero were removed, and ii) reactions with a mean flux value below 1e-05 were also excluded.
To assess how transcriptomic constraints reshape the metabolic network, we compared the unconstrained Fn-GEM with each context-specific model in terms of network size, reaction retention, and flux-supported metabolic activity.
Subsequent analyses focused only on reactions shared among the three context-specific models, allowing for a direct comparison of flux distributions across conditions. Pairwise comparisons were performed following the physiological progression of the experiment: first between control and adhesion, and then between adhesion and invasion. For each selected reaction, the Shapiro Wilk test was applied to assess the normality of the flux distributions, confirming that data did not follow a normal distribution. Consequently, the non-parametric Mann Whitney U test was used to evaluate whether fluxes of each reaction differed significantly between the two conditions. In parallel, the absolute difference in the mean flux values between the two conditions under comparison was calculated for each reaction, to quantify the magnitude of change independently of direction. Reactions were then filtered by applying a dual threshold: i) a False Discovery Rate (FDR) < 0.05 to control for multiple testing and ensure statistical significance, and ii) an absolute mean flux difference > 1 to retain changes that are both statistically significant and biologically relevant. This cutoff was chosen to quantify the extent to which the mean flux of the same reaction differs between the two conditions under comparison (control vs adhesion and adhesion vs invasion), independently of direction. The robustness of the applied thresholds was supported by the distribution of mean absolute flux differences, which showed a strong enrichment of very small changes and a long-tailed distribution toward larger values (Supplementary Figure 1-A,B). Empirical cumulative distribution function (ECDF) analysis indicated that the selected cutoff lies beyond the region dominated by minor numerical variability (Supplementary Figure 1-C,D). Importantly, varying the mean flux difference threshold in both a more permissive (0.5) and more stringent (2) direction relative to the selected cutoff (1) resulted in only gradual and proportional changes in the number of significant reactions, indicating that the results are not driven by a specific or arbitrary threshold choice (Supplementary Figure 1-E,F). Finally, for each filtered reaction, a Raincloud plot was generated using a Python implementation (Allen et al., 2021). These plots combine a density plot, a boxplot and a scatterplot, allowing for the simultaneous visualization of flux distributions, mean values and variability for each condition. This approach also allowed a qualitative assessment of flow directionality based on the sign of the flux distributions (i.e., positive or negative), providing additional insight into the reaction behaviour under different conditions. The final subset of selected reactions was manually inspected using the KEGG (Kanehisa et al., 2017) and BioCyc (Karp et al., 2019) databases.

3
Results
3.1
F. nucleatum subsp. nucleatum ATCC 25586 physiology in a complex medium
To characterize the F. nucleatum physiological properties in a rich environment, we first determined its growth kinetics in YCFA at different dilution rates. Preliminary growth curves, obtained during a 48-hour experiment, show how F. nucleatum exponentially grows for the first 15–20 hours, depending on the dilution factor used, after which the culture enters the stationary phase (Supplementary Figure 2).
Looking at the variability between replicates, we found the dilution rate of 1/10 to be the most consistent, i.e., the variance across eight biological replicates was minimal compared with the other dilution rates; hence we decided to use this condition for all subsequent experiments. In addition, since it has been shown how bacterial growth during the exponential phase leads to optimal production of biomass, we focused on this phase to characterize the F. nucleatum physiology in a rich environment. First, we analysed the growth kinetics to identify three distinct time points within the exponential phase (T1, T2 and T3), corresponding to different stages of active growth. These time points occurred approximately three, four and five hours after incubation start, respectively (Fig. 1A). Hereafter these time points will be referenced as T1, T2 and T3. Then, we obtained biomass samples at T2 to characterize the relative abundance of their macromolecule constituents, i.e., lipids, nucleic acids, proteins and glycogen.
The Fig. 1B reports the estimated values of the biomass fractions for DNA, RNA, Proteins, Lipids and Glycogen. We were able to characterize the principal macromolecular components of F. nucleatum biomass, which together represented the bulk of the cellular dry weight: proteins (40 %), RNA (20 %), lipids (9 %), and glycogen (2.1 %).
Further, we characterized the F. nucleatum gene expression profiles at three distinct time points (T1, T2, and T3), each representing a different stage within the exponential growth phase (Fig. 1A). After sequencing, read mapping and counting, gene expression was quantified for each of the 1,992 F. nucleatum coding genes. The expression of 63 of them was neglectable in all samples (log2(FPKM)<1). Fig. 1C shows the clustering of samples and their corresponding expression values. With the only exception of sample 10 and sample 2 belonging to time point T2 and clustering with T1 samples, all the other samples cluster consistently, with overall gene expression at T2 being more similar to T3 than T1. Gene expression changes between the time points were assessed to identify differentially expressed genes (DEGs) (Supplementary Table 2). DEGs were defined as genes showing a false discovery rate (FDR) ≤ 0.05 and a log2 fold-change (|log2FC|) ≥ 0.585, corresponding to an absolute expression change of at least 50 %. Comparing T1 vs T2, there were 64 DEGs, of which 58 and 6 were up- and down-regulated, respectively. The other comparisons, T2 vs T3 and T1 vs T3, yielded respectively 102 (79 up- and 23 down-regulated) and 341 DEGs (220 up- and 121 down-regulated).
To gain a functional overview on the physiological states at the different time points and to guide the reconstruction of Fn-GEM (see below), we computed the enrichment of functional categories between up- and down-regulated genes (Fig. 1D). The enrichment was analysed separately of up- and down-regulated genes and a Fisher test was used to compute the statistical significance (with a threshold set to p-value < 0.05).
In the shift between T1 and T2 we identified two functional categories that were enriched in over-expressed genes, namely, Secondary Structure (Q) and Amino Acid metabolism and transport (E). Conversely, four categories were enriched in down-regulated genes, i.e., Inorganic ion transport and metabolism (P), Post-translational modification, protein turnover, chaperone functions (O), Lipid metabolism (I), Nucleotide metabolism and transport (F). In the contrast between T2 and T3 we observed three functional categories displaying a significant fraction of up- and down-regulated genes. Specifically, up-regulated genes were significantly enriched in genes belonging to the P category (Inorganic ion transport and metabolism) whereas down-regulated genes were mostly represented by genes belonging to categories J and E (Translation and Amino Acid metabolism and transport). Remarkably, this latter functional category is additionally over-represented among genes that significantly increase their expression in the transition between T2 and T3 (Fig. 1C). This, together with the observation that the expression of genes belonging to amino acids metabolism were also significantly, altered in the T1-T2 transition, underscores the relevance of amino acids metabolism in this bacterium, as already shown in other studies (Robinson et al., 2025; Sakanaka et al., 2022; Wu et al., 2023; Yu et al., 2022).

3.2
Reconstruction of the F. nucleatum subsp. nucleatum ATCC 25586 genome-scale metabolic model (Fn-GEM)
A bottom-up approach was used to reconstruct a genome-scale metabolic model of F. nucleatum subsp. nucleatum ATCC 25586 starting from a first draft reconstruction of the model built on KBase (Arkin et al., 2018) narrative interface by using the newly assembled NCBI genome GCA_003019295.1 (Todd et al., 2018) (Fig. 1B). This draft reconstruction was then intensively manually curated as described in Materials and Methods, for example by comparing the production rate of Fn-GEM by-products to those of the real strain (both grown in CDM medium), according to previously reported results (Buckel and Barker, 1974; Rogers et al., 1991). The final reconstruction comprised 1,041 reactions, 938 metabolites, and 452 genes. Following these controls, the model described most of the previously reported metabolic abilities, lastly leading to a more comprehensive reconstruction and the possibility of performing growth simulations of Fn-GEM in different nutritional settings. Afterwards, to assess whether the reconstructed network was capable to correctly predict growth phenotypes, we compared the model predictions against a set of different experimentally obtained growth rates. In-silico defined CDM-like media (exact compositions derived from (Jensen et al., 2020) and described in Supplementary Table 3) were simulated by setting all the nutrients uptake reactions to zero, except for those nutrients actually present in the medium. Overall, the model reproduced the observed growth rates across all tested media (0.06–0.10 h⁻¹). Although growth was slightly underestimated in CDM (∼35 %) and overestimated in CDMG, CDMF, and CDMFII (∼30 %), the overall agreement confirms that the reconstructed network reliably predicts the growth capacity of F. nucleatum (Fig. 1E). To quantitatively assess the agreement between simulated and experimental growth rates, we compared predicted and observed values across four media conditions using both root mean square error metrics (RMSE). RMSE was low (0.027 h⁻¹), and absolute deviations between predicted and observed growth rates were consistently small (≤ 0.035 h⁻¹). These results indicate that the model reproduces growth rates that are quantitatively close to the experimental measurements and correctly captures relative growth differences across media.

3.3
Context-specific modelling of a conserved core network in Fusobacterium
Some Fusobacterium species are known to exhibit a facultative intracellular lifestyle, being capable of invading and living inside various human cell types, including epithelial, endothelial, keratinocytes, and potentially immune cells (Gursoy et al., 2008; Strauss et al., 2011; Xu et al., 2007). This is likely accompanied by an overall reprogramming of many cellular features, including Fusobacterium metabolism (Cochrane et al., 2020). While the overall strategy adopted by invading Fusobacterium cells has been explored through gene expression profiles (Cochrane et al., 2020), the metabolic support this lifestyle switch have been poorly described up to now. Therefore, we employed Fn-GEM to simulate context-dependent metabolic network reprogramming associated with the transition from a free-living lifestyle to an invasive one, passing through an adhesion condition. To this aim, we leveraged publicly available F. animalis transcriptomic data generated under in the same three experimental conditions (control, adhesion and invasion) (Cochrane et al., 2020), a functional input to contextualize the Fn-GEM. For doing this, we leveraged on one of the main features of GEMs, i.e. the possibility, when the core metabolism is conserved (see below) to use them as scaffolds to bridge -omics data from different (although related) organisms. Accordingly, we first computed the degree of conservation of metabolic genes between the F. animalis and F. nucleatum. Reciprocal best-hit analysis identified 1,469 high-confidence orthologous protein pairs, corresponding to 72.5 % of the coding sequences of F. nucleatum strain ATCC 25586 (GCA_003019295.1) and 60.7 % of F. animalis strain 7_1 (GCA_000158275.2). Importantly, when restricting the analysis to genes explicitly represented in the curated Fn-GEM, 77.3 % of model-associated genes were supported by a direct ortholog in F. animalis. Functional annotation of orthology-supported model genes revealed a predominant association with core metabolic processes, including energy production and conversion, amino acid and carbohydrate metabolism, nucleotide metabolism, and cofactor metabolism (Supplementary Figure 3). These results indicate that the transcriptomic data robustly inform the central metabolic network captured by the model. Based on this evidence, gene expression data were integrated into the Fn-GEM using GIMME algorithm (Becker and Palsson, 2008), and condition-specific flux distributions were subsequently analysed.
Compared to the unconstrained Fn-GEM (1041 reactions), GIMME-based contextualization progressively reduced the network size, retaining 715 reactions in the control model (32 % pruned), 639 in the adhesion model (39 % pruned), and 609 in the invasion model (42 % pruned) (Supplementary Table 4). Accordingly, a quantitative analysis of flux distributions was required to identify reactions whose activity significantly changes during adaptation to the host environment at the species-level metabolic network. To this end, flux sampling was performed on each context-specific model to capture the range and variability of feasible steady state fluxes.
Following flux sampling, non-informative reactions were excluded by removing reaction with total flux sum equal zero or mean absolute flux < 1e-5. After this filtering step, 355 reactions were retained in the control condition, 325 in adhesion and 306 in invasion. Among these, 278 reactions were common to all three conditions and were therefore used for pairwise comparisons. Finally, applying combined statistical and magnitude-based criteria (FDR < 0.05 and |mean flux difference| > 1), 113 reactions were identified as significantly altered in the control vs adhesion comparison and 87 for adhesion vs invasion comparison (Supplementary Table 5). While a small number of reactions were found to be unique to either adhesion or invasion (i.e., not shared with control or with each other), none of these distinct reactions revealed was of biological interest (Supplementary Table 5, Fig. 2). This suggests that the lifestyle switch of Fusobacterium (from free-living to intracellular) mostly implies a fine-tuning of specific pathways, rather than the activation/deactivation of entire metabolic processes. Accordingly, a quantitative analysis is required to understand which reactions are significantly changing their activity to guarantee an efficient adaptation to the host environment. In the next sections, we address this point, initially distinguishing between a free living and a general “host-associated” Fusobacterium metabolism and later by analysing each transition (i.e., free living to adhesion and adhesion to invasion) separately. Finally, we checked whether using Fn-GEM to contextualize F. animalis affected the results reported above, i.e. the reaction sets shown to probably play a role in F. animalis lifestyle transitions. We found that, out of 113 reactions whose flux distributions were statistically altered in the control-adhesion transition, only 13 were (11 %) not assigned to a F. animalis gene and 15 out of 87 (17 %) in the adhesion-invasion transition. In other words, the 89 % and 83 % of the reactions showing a significant change in Fusobacterium control-adhesion and adhesion-invasion transitions, respectively, were reactions encoded by genes present both in F. nucleatum and F. animalis, thus supporting the use of Fn-GEM as a scaffold to explore the metabolic transitions of F. animalis. While the number of reactions missing a consistent gene annotation in F. animalis and still included in the Fn-GEM surely deserve further inspection, they represent a minor fraction of the core metabolism likely involved in adhesion to invasion transition and can be thought as the set of gap-filling reactions that are commonly included in metabolic reconstructions.

3.4
Overall network reprogramming during the transition from free-living to host-associated states
Here, we specifically address which metabolic pathways differentiate invading Fusobacterium metabolism during free-living and host-associated phases. Therefore, we analysed flux sampling outputs across the three conditions (control, adhesion and invasion), focusing on reactions displaying significantly altered flux distributions in control vs adhesion and control vs invasion comparisons. Broadly speaking, results of this analysis (Supplementary Table 5) indicate that metabolic reprogramming associated with adhesion is characterized by suppression of branched-chain amino acid (BCAA) degradation, and increased uptake of serine and methionine. In addition, adhesion is associated with a shift in carbon flux distribution toward mixed-acid fermentative pathways. This shift is characterized by increased pyruvate-to-formate conversion and secretion of fermentation end-products, consistent with a redox-balancing metabolic strategy. In contrast, the transition to the intracellular invasive phase involved a partial reversal of adhesion-specific features, including reactivation of BCAA biosynthesis, remodeling of one-carbon metabolism, and a broader redistribution of nitrogen fluxes. In the next section, we provide the details of the metabolic reprogramming occurring in F. animalis during life-style transitions.

3.5
Transition from free-living to host-adherent metabolic state
The Fusobacterium transition from a free-living state (control) to adhesion marks a profound metabolic reprogramming (Fig. 3; for more details, see Supplementary Table 5 and Supplementary Figure 4). Many of these changes involved amino acid transport and metabolism, redox balance, and carbon utilization, suggesting a coordinated adaptation to host cell adhesion. Amino acid metabolism and transport were among the most markedly altered processes during adhesion. In the control condition, several essential amino acids, including L-leucine, L-glutamine, valine, lysine, and threonine, were actively imported. However, once adhesion has occurred, these fluxes decreased dramatically or even reversed direction. For example, L-leucine transport via LEUt4 shifted from robust import to moderate export, while LEUt2r showed an opposite trend. This pattern suggests a probable compensatory role between the two transporters. In parallel, glutamine symport (GLNt4) was significantly decreased, mirroring a diminished cellular demand for key nitrogen donors. A strong reduction in valine uptake was also observed, as shown by the reduced flux through its proton symporter (VALt2r). Notably, the fluxes distribution in control condition for this reaction was wide, suggesting greater metabolic flexibility. In contrast, during adhesion fluxes were strictly constrained near zero, indicating restricted valine import (Fig. 3). Isoleucine catabolism was similarly affected, with a marked decrease in flux through ILETA. In addition, flux variability was greatly reduced in adhesion, indicating a strict suppression of isoleucine degradation. These combined trends suggest a partial repression of BCAA catabolism under adhesion, potentially reflecting a strategy to preserve biosynthetic intermediates and limit unneeded loss of nitrogen or carbon. In contrast, methionine and serine uptake increased significantly during adhesion. Specifically, methionine flux (METt2r) increased from negligible levels in control to one of the highest transport rates in adhesion, potentially linked to methylation processes or biosynthetic responses to stress. Serine uptake, estimated from mean fluxes through SERt4, was also doubled. This increase is consistent with its central role in nucleotide synthesis, lipid metabolism and one-carbon metabolism. Serine catabolism to pyruvate and ammonium (SERD_L) also increased, indicating its dual role of energy production and nitrogen supply. In addition, intracellular nitrogen processing appeared restructured. Interestingly, while GLUDxi showed decreased glutamate deamination, GLUDy reversed its direction by modestly regenerating glutamate from α-ketoglutarate and ammonium. These opposing flux patterns, evident from the shifts in mean fluxes, likely reflect a dynamic regulation of intracellular nitrogen balance. In this context, glutamate may function both as a key nitrogen donor for biosynthetic processes and as a substrate for energy metabolism, depending on the specific metabolic demands during adhesion. Flux through aspartate transaminase (ASPTA) not only decreased but also reversed direction, suggesting a shift in nitrogen assimilation and redistribution, while threonine synthesis via THRA also decreased, likely reflecting changes in amino acid biosynthesis. Tryptophan metabolism appeared to be more active during adhesion, as indicated by increased fluxes through tryptophanase (TRPAS2) and tryptophan synthase (TRPS2). This pattern suggests a bidirectional cycle between tryptophan and indole. Such cycling may play a role in adaptation to bacterial stress or intercellular signalling (Bansal et al., 2010; Lee and Lee, 2010). During adhesion, phenylalanine (PHEt2r) uptake was strongly reduced/suppressed, while tyrosine catabolism via TYRL showed only a modest decrease, indicating a selective conservation among aromatic amino acids.
In parallel, central carbon and energy metabolism also underwent a marked change after adhesion. An increase in flux through pyruvate formate lyase (PFL) suggests a greater involvement of mixed acid fermentation. This may reflect an alternative strategy for energy production or redox balancing during adhesion. Consistently, increased formate export (FORt2 and EX_for_e) was observed, indicating activation of fermentative pathways. Moreover, recent evidence shows that co-culture of F. nucleatum with tumor cells leads to elevated formate secretion, identified as a key oncometabolite in CRC progression (Ternes et al., 2022). Interestingly, although formate is a typical fermentation byproduct, its increased secretion could also reflect additional roles in intercellular communication or stress adaptation. Acetaldehyde dehydrogenase (ACALD) showed a complete reversal of directionality, switching from acetaldehyde consumption to its production upon adhesion. This reversal suggests a change in the way acetyl-CoA is metabolized. Notably, this change in metabolism was associated with an increase in flux through acetyl-CoA ACP transacylase (ACOATA), suggesting that more acetyl-CoA could be used for lipid biosynthesis. In parallel, the TCA cycle appeared to be partially activated. Citrate uptake increased significantly, as shown by the high fluxes in the transport reaction (CITt2r) and associated exchange (EX_cit_e). Key enzymatic steps such as aconitase (ACONT), isocitrate dehydrogenase (ICDHyr), and citrate lyase (CITL) also showed increased fluxes. These changes suggest that citrate was more actively used and recycled, probably to provide building blocks for biosynthesis or to keep metabolism flexible during adhesion. Fructose uptake increased during adhesion, as indicated by higher import fluxes through both extracellular exchange (EX_fru_e) and the specific transport reaction (FRUt4). Once inside the cell, fructose was rapidly metabolized, as shown by its conversion into glucose via aldose-ketose isomerase (XYLI2), suggesting that fructose represents a critical carbon source under this condition. Furthermore, notable changes were also detected in the interconversion of sugar phosphates. Several glucose-6-phosphate isomerases (such as G6PI2 and G6PI3) showed increased fluxes, indicating an intensification of the cycle between glycolysis and the pentose phosphate pathway. Despite these metabolic rearrangements, the fluxes through pyruvate kinase (PYK) remained substantially unchanged, suggesting a constant production capacity of ATP. Additionally, adhesion to host cells also induced specific changes in the redox and one-carbon metabolism of Fn. In detail, thioredoxin peroxidase (THIORDXi), a key enzyme in the defence against oxidative stress, showed an evident increase in fluxes. This enzyme catalyses the reduction of hydrogen peroxide through the thioredoxin system, suggesting an enhanced capacity for ROS detoxification (Zeller and Klug, 2006).

3.6
Metabolic network adaptation during intracellular colonization of human colon epithelial cells
The metabolic transition from adhesion to invasion in Fusobacterium was characterized by evident shifts in fluxes across several key metabolic processes (Fig. 4; for more details, see Supplementary 5 and Supplementary Figure 5). These changes probably reflect bacterial adaptation to the host intracellular environment. Specifically, the invasion condition revealed a marked reprogramming of amino acid metabolism, including transport and degradation pathways. Serine utilization is decreased, with lower fluxes through both the sodium symporter (SERt4) and the deaminase (SERD_L). This suggests a reduced requirement for serine in biosynthetic processes or cofactor-related functions. Notably, in invasion the mean flux through SERt4 decreased and returned to levels similar to those observed in the control condition, indicating a partial reversal of the serine uptake phenotype observed during adhesion. Tryptophan metabolism showed a tight regulation under different conditions, involving both its biosynthetic and degradative pathways. In detail, fluxes through TRPS2 and TRPAS2 were markedly decreased in the invasion condition, returning to levels like those of the control. The reduction in cycling may reflect a metabolic shift away from indole-mediated signalling and nitrogen or redox balancing. Instead, metabolism may favour a more stable metabolic state, likely oriented toward biomass maintenance or energy conservation during tissue invasion. Interestingly, the flux through isoleucine transaminase ILETA increased markedly in invasion, but with a negative sign, indicating a reversal of the direction of the reaction. Rather than supporting isoleucine degradation, this shift suggests increased utilization of 3-methyl-2-oxopentanoate and glutamate for isoleucine and α-ketoglutarate regeneration. This pattern is consistent with a possible anabolic requirement for BCAAs under invasive conditions. Glutamine uptake via GLNt4 increased during the invasion phase and was accompanied by higher activity of the NAD⁺-dependent glutamate dehydrogenase (GLUDxi). This could reflect increased glutamine demand in the tumor microenvironment (TME), in line with experimental evidence showing reduced extracellular glutamine levels in cancer cells upon co-culture with F. nucleatum (Ternes et al., 2022). The higher flux values ​​observed for GLUDxi indicated an enhanced conversion of glutamate to α-ketoglutarate, with ammonium release and reduction of NAD⁺ to NADH. Overall, fluxes had a broader distribution during invasion in respect to the other conditions (control and adhesion). This suggests a greater variability in how cells manage nitrogen and redox balance under this condition. A clear metabolic shift towards fermentative pathways was also observed during the invasion phase. Specifically, the flux through D-lactate dehydrogenase (LDH_D) increased significantly, indicating enhances oxidation of D-lactate to pyruvate and NADH. This pattern points to a reutilization of fermentation byproducts as substrates for central carbon metabolism or energy production. Such recycling may support biosynthetic demands or redox balancing under invasion condition. In contrast, the flux of pyruvate formate lyase (PFL), remained at high positive levels in invasion, confirming the continued conversion of pyruvate to acetyl-CoA and formate. However, the modest reduction in its flux compared to adhesion suggests a partial redirection of pyruvate away from formate production. Pyruvate may instead be incorporated into other acetyl-CoA-dependent biosynthesis or entry into the TCA cycle, consistent with the observed reutilization of lactate-derived pyruvate. This stage-dependent modulation could indicate that Fusobacterium relies more on formate during early adaptation and may redirect pyruvate once a stable niche is established, with the possibility that formate production increases again at later stages depending on growth dynamics. Whereas reactions involved in butyrate metabolism (BUTKr, PBUTT and BUTCT) showed a small but consistent decrease in flux during the invasion phase. The directionality of the flux indicates that PBUTT and BUTKr proceed in the direction of butyrate formation from butyryl-CoA, while BUTCT operates in reverse, suggesting that some butyrate was converted back to butyryl-CoA. This bidirectional behaviour highlights a dynamic balance between butyrate production and reutilization, likely related to CoA regeneration or metabolic flexibility during invasion. Furthermore, the early steps of the TCA cycle appeared less active during the invasion phase. Specifically, fluxes through citrate uptake (CITt2r), aconitase (ACONT) and NADP⁺-dependent isocitrate dehydrogenase (ICDHyr) were reduced. This pattern reflects decreased conversion of extracellular citrate into isocitrate and subsequently into α-ketoglutarate, along with a lower generation of NADPH and CO₂. This may be consistent with a shift toward fermentative metabolism during invasion, promoting redox-neutral or partially reductive pathways over fully oxidative respiration. Moreover, during the invasion phase, folate metabolism showed a marked reconfiguration. In detail, flux through the methylenetetrahydrofolate dehydrogenase (MTHFD) strongly increased and displayed a broader distribution in fluxes, suggesting enhanced NADPH production and flexible engagement of one-carbon metabolism. In contrast, the MTHFC reaction, which catalyses the downstream cyclohydrolase step, was dramatically reduced and became almost inactive. This indicates a possible bottleneck or redirection of folate intermediates under invasive conditions.

4
Discussion
F. nucleatum is one of the most alarming microbes to date, due to its capability to colonize different districts of the human body. Indeed, besides being one of the most important players in periodontitis-associated infections, it has been shown to be able to colonize on colorectal tissues in a remarkable amount of CRC cases (Amedei et al., 2021; Brennan and Garrett, 2019; Castellarin et al., 2012). In addition, evidence documents that F. nucleatum can adhere, invade and undergo a non-obligate intracellular life stage within human cells, including epithelial, endothelial and potentially immune ones (Gursoy et al., 2008; Strauss et al., 2011; Xu et al., 2007; Yang et al., 2017). Here, we reconstructed an experimentally validated GEM of F. nucleatum, providing a strain-specific resource that extends prior automated or partially curated reconstructions not designed for systematic validation under defined experimental conditions. We used this model to integrate existing transcriptomic data obtained under different conditions as a functional input to investigate context-dependent metabolic network reprogramming in invading Fusobacterium.
It is important to clarify that the transcriptomic data integrated into the model were derived from F. animalis strain 7_1 (Cochrane et al., 2020), whereas the reconstructed genome-scale metabolic model corresponds to F. nucleatum ATCC 25586. While metabolic responses may differ between taxa, our comparative genomics analyses show that the vast majority of genes represented in the model belong to a conserved metabolic core shared between F. nucleatum and F. animalis, that also includes the vast majority of those significantly involved in the life-style transition-mediated. Genome-scale metabolic models primarily capture conserved biochemical functions required for growth, energy balance, and redox homeostasis. Accordingly, the metabolic adaptations described here should be interpreted as reflecting responses encoded by a shared core metabolic network. Strain- or species-specific regulatory or virulence-associated differences may exist and warrant future investigation.
Data obtained from our analyses point towards a key role of amino acids metabolism regulation in the transition between these three stages. While amino acids are actively imported during free-living growth, consistent with the asaccharolytic nature of this bacterium, their overall utilization decreases during the adhesion phase and is selectively remodulated during the intracellular phase. Since F. nucleatum relies on amino acids metabolism to thrive (Robrish and Thompson, 1988; Rogers et al., 1998, 1991; Zilm et al., 2003), this scenario is compatible with a general downregulation of growth requirements during the adhesion phase, possibly leading to the redirection of metabolism for other processes. In obligate anaerobes, fermentative metabolism is intrinsically coupled to the generation and reoxidation of reducing equivalents to maintain intracellular redox balance (Buckel, 2021). Therefore, the modulation of amino acid catabolism may directly influence NADH production rates and the need for compensatory fermentative fluxes. Among them, we found an overall increase of the flux through three key enzymes that may lead to an increase in the pool of cellular acetyl-CoA (aconitase (ACONT), isocitrate dehydrogenase (ICDHyr), and citrate lyase (CITL)), paralleled by an increase in the activity of acetyl-CoA ACP transacylase (ACOATA), the entry point of lipids production. This overall fluxes’ redirection may indicate a general rewiring towards lipid biosynthesis, compatible with the synthesis of outer membrane vesicles (OMVs) by F. nucleatum during the adhesion to host cell. Indeed, F. nucleatum adhesion to human cells occurs by means of OMVs harbouring toxic bacterial components are continuously released during F. nucleatum growth (Chen et al., 2022; Liu et al., 2019; Zheng et al., 2024). OMVs produced by F. nucleatum have been shown to enter host cells, transfer virulence-associated factors and modulate inflammatory responses, including macrophage polarization and cytokine production (Chen et al., 2022; Liu et al., 2019). Thus, increased flux toward lipid biosynthesis during adhesion may mechanistically support membrane remodeling and OMV production, linking metabolic rewiring to virulence-associated host modulation (Chen et al., 2022; Liu et al., 2019). Invasion lifestyle marks the re-activation of the central metabolism oriented towards biomass production and/or energy conservation (clearly from a different set of input nutrients in respect to the complete medium of the control condition), consistently with the observation that intracellular F. nucleatum cells can multiply for several hours upon entrance in the host cell (Gursoy et al., 2008). Therefore, compared to the adhesion condition, the invasion is characterized by a substantial restoration of nitrogen assimilation, redox balance and fermentative metabolism, all adaptations essential for survival and proliferation within host cells. The observed increase in glutamine uptake (GLNt4) and the concomitant activation of NAD⁺-dependent glutamate dehydrogenase (GLUDxi) is emblematic of this shift. These reactions facilitate the conversion of glutamate to α-ketoglutarate, releasing ammonium and reducing NAD⁺ to NADH. In anaerobic metabolism, oxidative decarboxylation reactions generate reduced electron carriers that must be reoxidized through fermentative pathways to sustain metabolic flux (Buckel, 2021). Therefore, the coordinated increase in glutamine catabolism and fermentative reactions observed in invasion model may reflect a redox-balancing configuration typical of anaerobic bacterial adaptation. In particular, the reversal of flux through isoleucine transaminase (ILETA), from catabolic to anabolic, suggests that BCAAs are regenerated to meet biosynthetic needs, potentially for protein synthesis or membrane maintenance during intracellular adaptation. BCAAs have been shown to play a role in tumor-associated signaling pathways in the gut microbiota context. For instance, increased BCAA production by the anaerobic bacterium Clostridium symbiosum promotes colorectal cancer progression via activation of mTOR signaling (Ren et al., 2024). While this mechanism has not been directly demonstrated for F. nucleatum, our predicted modulation of BCAA metabolism suggests that alterations in local BCAA availability could represent a potential axis of host-microbe metabolic interaction within the tumor microenvironment. This metabolic configuration is particularly intriguing in the context of the TME. F. nucleatum does not utilize glucose as a primary energy source. Instead, is known to exploit amino acids and peptides as carbon and nitrogen sources, consistent with its asaccharolytic metabolism. This metabolic configuration limits direct competition with tumor cells for glucose, a substrate largely consumed by cancer cells via aerobic glycolysis (Yu et al., 2022). Interestingly, the metabolic profile observed during invasion not only reflects a reactivation of central metabolism for intracellular adaptation, but also suggests a potential contribution of Fusobacterium to the TME. Indeed, the increased fluxes through fermentative pathways (e.g., LDH_D and PFL), along with the dynamic utilization and recycling of butyrate-related metabolites (BUTKr, PBUTT, BUTCT), point to an active production of short-chain fatty acids and organic acids such as formate, lactate, and butyrate. These compounds, widely reported to accumulate in colorectal tumors, are known to modulate host immune responses, promote cancer stem cell self-renewal, and stimulate epithelial proliferation (Ternes et al., 2022). In addition, F. nucleatum has been shown to activate oncogenic pathways such as Wnt/β-catenin signaling via interaction of its adhesin FadA with host E-cadherin, thereby promoting epithelial proliferation and tumor progression (Yang et al., 2026). Consistently, a study investigating the metabolic activity of F. nucleatum in colorectal tumors reported a significant upregulation of pyruvate fermentation pathways, along with inosine, adenosine, and rhamnose biosynthesis (Zhou et al., 2025). These findings further support the idea that Fusobacterium adapts to the TME through fermentative metabolism and nucleotide-related biosynthetic pathways. In parallel, microbiota-derived formate has been shown to actively promote CRC progression by stimulating glycolysis and activating hypoxia-inducible factor 1-alpha (HIF1α) signalling in tumor cells, thereby accelerating malignant transformation (Ternes et al., 2022). From a bacterial perspective, the increased formate and lactate production observed during invasion could also reflect an adaptive response to oxidative stress, as previously demonstrated in F. nucleatum cultures grown under continuous conditions (Diaz et al., 2000), supporting the idea that fermentative metabolism plays a dual role in both energy/redox homeostasis and host interaction.
In the present framework, transcriptomic data integrated into the model describe discrete experimental conditions (control, adhesion, and invasion) and provide static snapshots of gene expression associated with each physiological state. Therefore, the contextualized models are intended to capture representative metabolic configurations characteristic of these conditions, rather than the full temporal dynamics of the transitions between them. Within the constraint-based modelling framework, this approach allows the identification of stable, condition-specific metabolic states. In contrast, dynamic aspects of metabolic adaptation would require time-resolved transcriptomic or multi-omics datasets, which represent a natural extension of our study.
Finally, while our flux-based analyses predict substantial reorganization of amino acid utilization, nitrogen assimilation, and fermentative metabolism across different physiological states of Fusobacterium associated with host interaction. Integrating targeted or untargeted metabolomics data would represent an important next step to experimentally validate predicted changes in metabolic fluxes and metabolite secretion, particularly in the context of host-microbe interactions and TME modulation.

5
Conclusion
In this study, we reconstructed and validated the GEM of F. nucleatum subsp. nucleatum ATCC 25586, capturing its core physiological features under controlled growth conditions. By integrating gene expression data from different conditions into Fn-GEM, we generated context-specific models that revealed how invading Fusobacterium (at least F. nucleatum and F. animalis) reprograms its metabolism to adapt to host-associated lifestyles.
Building on the present study, additional modelling approaches could provide further insight into Fn metabolism within the TME. Coupling the reconstructed bacterial GEM with host cell metabolic models would enable explicit investigation of host-microbe metabolic interactions, including nutrient competition, metabolite exchange, and redox coupling. In addition, embedding the model within a community-level context would allow exploration of metabolic cross-feeding and cooperative interactions with other gut microbes enriched in CRC. Further, the integration of modelling techniques as those used here with experimental platforms such as the colon-on-a-chip microfluidic devices (Mitrofanova et al., 2024) would permit to get to a system-level characterization of the complex interplay between the microbiota and the human intestinal cells. Taken together, these approaches represent promising strategies to better contextualize invading Fusobacterium metabolic activity within the complex and dynamic landscape of the TME.
Recent advances in genome-scale metabolic modelling have enabled mechanistic investigation of host-microbe metabolic interactions through integrated host-microbiome and community-based approaches that explicitly capture metabolite exchange and cross-feeding dynamics (Heinken et al., 2021, 2013; Lee et al., 2025; Molina Ortiz et al., 2025; Srinak et al., 2025). In this context, our species-specific, flux-resolved modelling of Fusobacterium provides a focused and complementary perspective on metabolic adaptation during host interaction and can be readily extended to multi-organism or host-integrated modelling strategies.
Overall, our study offers a robust systems-level resource to explore invading Fusobacterium’s metabolic strategies during host interaction and provides a solid foundation for future investigations aimed at targeting its metabolic vulnerabilities in disease contexts such as colorectal cancer, to prevent or contrast the pathologic role of F. nucleatum.

Data availability
RNA-seq raw reads of Fusobacterium nucleatum subsp. nucleatum strain ATCC 25586, generated from three time points within the exponential growth phase (T1, T2 and T3) with five biological replicates each, have been deposited under NCBI BioProject PRJNA1327068.
The reconstructed genome-scale metabolic model (GEM) of F. nucleatum has been deposited in Zenodo under the https://doi.org/10.5281/zenodo.17098853.

Funding
The research was founded with a grant from the regional contribution of “The Programma Attuativo Regionale (Toscana)” funded by FAS (now FSC), grant MICpROBIMM. Funding was awarded by R. Fani and A. Amedei.

CREDIT author statement
MG: Methodology; Validation; Formal analysis; Data curation; Visualization; Writing - original draft; Writing - review & editing.
EB: Methodology; Formal analysis; Investigation.
WV: Formal analysis.
LP: Methodology; Validation; Formal analysis.
EV: Review & editing.
IB: Review & editing.
GN: Review & editing.
MDS: Investigation; Resources.
HB: Investigation; Resources.
TDL: Investigation; Resources.
TS: Investigation.
EL: Investigation.
JK: Investigation.
AA: Funding acquisition.
RF: Project administration; Funding acquisition.
MF: Conceptualization; Methodology; Writing – original draft; Writing – review & editing; Supervision; Project administration.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Renato Fani and Amedeo Amedei reports financial support was provided by Tuscany Region. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.