Gempipe: a tool for drafting, curating, and analyzing pan and multi-strain genome-scale metabolic models

TL;DR

Gempipe is a new tool that improves the creation and analysis of genome-scale metabolic models for multiple bacterial strains, combining reference-based and reference-free methods.

Contribution

Gempipe introduces a hybrid method for multi-strain GSMM reconstruction, enabling accurate pan-GSMM creation from single-strain references.

Findings

Gempipe improves accuracy in multi-strain GSMM reconstruction compared to existing tools.

It enables the creation of pan-GSMMs even when only a single-strain reference is available.

Metabolic diversity in Limosilactobacillus reuteri was explored, revealing metabolically coherent strain clusters.

Abstract

Genome-scale metabolic models (GSMMs) can mechanistically explain phenotypic differences among closely related bacterial strains. However, high-throughput multi-strain reconstructions of GSMMs are still challenging: reference-based methods inherit curated information while missing new contents; alternatively (universe-based), reference-free methods could cover strain-specific reactions, but they disregard curated information. Ideally, references should be curated pan-GSMMs for species (or genus), but their reconstruction is extremely demanding, making them still rare in the literature. Here, Gempipe is presented, a computational tool streamlining the multi-strain reconstruction and analysis of GSMMs, going through the production of a pan-GSMM. Its reconstruction method is hybrid; as an optional reference, GSMM is automatically expanded with extra reactions taken from a reference-free reconstruction. Gempipe also downloads, filters, and annotates genomes; performs in-depth gene recovery; annotates models’ contents; and predicts strain-specific capabilities. The companion programming interface includes functions ranging from the (pan-)GSMMs’ curation to the multi-strain analysis. Gempipe was validated using multi-strain data sets, showing improved accuracy when compared with state-of-the-art tools. Moreover, metabolic diversities within Limosilactobacillus reuteri were explored, grouping strains into metabolically coherent clusters and systematically predicting health-related metabolites’ biosynthesis.

Available genome-scale metabolic model (GSMM) reconstruction tools present major limitations in the context of multi-strain modeling. Gempipe surpasses these limitations by implementing a novel, hybrid reconstruction strategy. Not only does it produce more accurate strain-specific GSMMs, but it also produces pan-GSMMs when the only available reference is a manually curated model for a single strain, which is currently the most common case. With the vast availability of genome sequences, the high-throughput, multi-strain GSMM reconstruction and analysis approach provided by Gempipe will facilitate large-scale studies of exploration and bioprospecting of strain-level bacterial metabolic diversity, moving a step forward in strains’ screening and rational selection.

Figures

Overview of Gempipe. “gempipe recon” creates a draft pan-GSMM and the presence-absence matrix (PAM), containing the information relative to the protein clustering; “gempipe derive” derives strain-specific GSMMs, starting from the PAM and the pan-GSMM; the latter should be manually curated beforehand, for example, by using dedicated functions of the Gempipe API. When using “gempipe autopilot,” the manual curation is substituted by an automatic gap-filling of the draft pan-GSMM. Finally, the Gempipe API can be used once again to perform multi-strain analyses. When strains of different species are inputted in the same run, species-specific GSMMs are also produced, defined by the set of reactions always present in all strain-specific GSMMs belonging to the same species.

Comparison of the general reconstruction metrics. “gempipe_rf” indicates Gempipe ran without reference. (A) Content comparison. “G”: number of genes; “R”: number of reactions excluding exchange reactions; “uM”: number of unique metabolites, i.e., not considering their compartment; “exr”: number of exchange reactions. Bar height corresponds to the mean between strains, while error bars represent standard deviations. Hatched area represents the contents in common with the reference. (B) Similarity between tools based on reaction content. Cells report the mean Jaccard index along the strains, computed for the reaction IDs.

Comparison between experimental and simulated Biolog PM growth assays. Bar height corresponds to the mean between strains, while error bars represent standard deviations. “gempipe_rf” indicates Gempipe ran without reference. The “_pan” suffix indicates the use of a manually-curated pan-GSMM as the reference model, instead of a strain-specific one. Data set composition is the following: “01_klebsiella,” 37 strains belonging to the Klebsiella pneumoniae species complex; “02_ralstonia,” 11 strains belonging to the Ralstonia solanacearum species complex; “03_pseudomonas,” 36 strains of Pseudomonas chlororaphis. (A) Outcome of the comparison considering single substrates. TP, true positive; TN, true negative; FP, false positive; FN, false negative. (B) Overall metrics for the substrate utilization prediction.

Comparison between experimental and simulated gene essentiality assays. Bar height corresponds to the mean between strains, while error bars represent standard deviations. “gempipe_rf” indicates Gempipe ran without reference. Data set composition is the following: “04_streptococcus,” 17 strains of Streptococcus pneumoniae grown on rich THY medium; “05_pseudomonas,” 9 strains of Pseudomonas aeruginosa grown on rich LB medium; “06_pseudomonas,” the same 9 strains of P. aeruginosa but grown on minimal M9 medium. (A) Comparison between modeled genes (background gray bars, with gray point representing individual strains) and modeled genes for which an experimental outcome is available (colored bar, with black points representing individual strains). (B) Outcome of the comparison considering single genes. TP, true positive; TN, true negative; FP, false positive; FN, false negative. (C) Overall metrics for the gene essentiality prediction.

Relative number of modeled reactions with no GPR (orphans). Exchange, sink, and demand reactions are excluded, as well as reactions containing the substring “diffusion” in their name. Bar height corresponds to the mean between strains, while error bars represent standard deviations. The hatched area represents orphans in common with the reference.

Multi-strain analysis using the Gempipe API functions. (A) Phylometabolic tree built using presence/absence data for reactions, auxotrophies, and growth on alternative C substrates. Six clusters of metabolically coherent strains were extracted from the tree and shown alongside the subspecies attribute. Only features not constant across strains are represented. (B) Potential strain-specific production of health-related metabolites, vitamin B12 (“adeadocbl_c”), reuterin (“3hppnl_c”), and histamine (“hista_c”), and presence of urease (“UREA”). The subspecies attribute is reported.

Full text

INTRODUCTION

Different strains of the same bacterial species can exhibit marked differences at the phenotypic level, such as the ability to catabolize different substrates, the presence of specific auxotrophies, or the acquisition of biosynthetic pathways through lateral gene transfer (1, 2).

Genome-scale metabolic models (GSMMs) are systems-biology tools that describe the metabolic potential encoded by a genome. Assuming a steady state and specifying a biomass composition, their constraint-based simulations enable predictions of cellular growth under specific nutritive inputs (3). Therefore, the availability of a GSMM for each strain of interest enables in silico screenings of phenotypic characteristics, offering a faster and more cost-effective alternative to traditional experimental methods, generating hypotheses to be subsequently validated experimentally.

Unfortunately, the creation of high-quality GSMMs is time-expensive because manual curation is required (4). This bottleneck remains a key driver for the development of automated tools (5–9), some of which quickly gained traction due to their ability to produce simulation-ready GSMMs, using a reference-free, universe-based top-down approach (10).

Since the first pioneering work by Monk and colleagues on Escherichia coli (11), the general steps for the multi-strain reconstruction of GSMMs have remained mostly the same (12–18). Briefly, genomes are collected and filtered for quality; then, genes are predicted and clustered, creating orthologous gene families sometimes referred to as the pangenome; one of the strains includes a high-quality, manually curated GSMM that is used as reference; for each strain, a copy of the reference is made and all the genes that do not have an ortholog with the reference are subtracted, consequently removing associated metabolic reactions; finally, gap-filling is usually limited to minimal media as the starting genomes are quality-filtered and excessive gap-fillings could hide true strain-specificities.

The above reference-based method was formalized in 2019 by Norsigian et al. (19), where metabolic functions are inherited from the reference GSMM after orthologous genes are detected via a blastp (20) best reciprocal hits (BRH) alignment. The method was then implemented in Bactabolize (21), a recent tool published in 2023. Overall, its efficacy is dependent on the availability of a curated and phylogenetically close GSMM taken as reference.

However, the method has a key limitation: output GSMMs contain subsets of the reactions in the reference, excluding unmodeled strain-specific reactions. Indeed, the protocol (19) requires a manual curation of the output GSMMs, adding new reactions that were not originally in the reference, a step that was not automated in Bactabolize (21). For this reason, to fully capture strain-specific metabolic features, a curated pan-GSMM that encompasses the metabolic diversity of the entire species (or genus) should be provided as a reference instead of a strain-specific GSMM, as suggested in Bactabolize.

However, while curated strain-specific GSMMs are time-consuming to produce and therefore often lacking for non-model organisms, comprehensive pan-GSMMs are even more challenging to obtain. When manually curated, pan-GSMMs can require years of development (16, 22, 23), and indeed they are still rare in the literature (21): the few covered species include Klebsiella pneumoniae (23), Escherichia coli (24), Salmonella enterica (12), and Bacillus subtilis (17, 18).

Given the vast number of strain-specific genome sequences now available and the scarcity of comprehensive pan-GSMMs, there is a growing need for tools performing multi-strain reconstructions efficiently, even in the absence of a pan-GSMM. These tools should not only be capable of using a reference strain as a starting point, but also of autonomously integrating new, strain-specific reactions to fully capture the metabolic diversity across strains.

In this work, Gempipe is introduced, a novel package that, to the best of our knowledge, is the first to offer pan and multi-strain reconstruction of GSMMs by implementing a hybrid reconstruction method where an optional reference GSMM is automatically expanded with new contents transferred from an independent reference-free reconstruction. Gempipe also provides additional features: retrieval and quality-filtering of genomes; gene annotation; in-depth gene recovery; re-annotation of modeled contents; a companion application programming interface (API) helping the manual curation of (pan-)GSMMs; an “autopilot” mode skipping the (recommended) manual curation; various flux-balance analysis (FBA)-based predictions of strain-specific metabolic features; dedicated API functions for the multi-strain analysis. In this sense, Gempipe is not only a reconstruction tool, but also an analysis tool in the context of biodiversity exploration/bioprospecting.

MATERIALS AND METHODS

From genomes to gene clusters

Gempipe is composed of three command-line programs, “gempipe recon,” “gempipe derive,” and “gempipe autopilot,” along with a Python API (Fig. 1). “gempipe recon” reconstructs draft pan-GSMMs and supports four types of inputs: proteomes in Genbank format, proteomes in FASTA format, genome assemblies in FASTA format, and a list of NCBI Species Taxonomy IDs (taxids). When taxids are given in input, all the available assemblies for the indicated species are automatically downloaded from Genbank. Each genome/proteome is treated as a separate strain.

Overview of Gempipe. “gempipe recon” creates a draft pan-GSMM and the presence-absence matrix (PAM), containing the information relative to the protein clustering; “gempipe derive” derives strain-specific GSMMs, starting from the PAM and the pan-GSMM; the latter should be manually curated beforehand, for example, by using dedicated functions of the Gempipe API. When using “gempipe autopilot,” the manual curation is substituted by an automatic gap-filling of the draft pan-GSMM. Finally, the Gempipe API can be used once again to perform multi-strain analyses. When strains of different species are inputted in the same run, species-specific GSMMs are also produced, defined by the set of reactions always present in all strain-specific GSMMs belonging to the same species.

When genomes are inputted, they are first subjected to gene annotation running Prodigal v2.6.3+ (25) via Prokka v1.14.6+ (26), the latter used just for its gene naming feature. Next, BUSCO v5.4.0+ (27) is run on the proteomes, indicating the database of interest, which is automatically downloaded, to obtain the number of missing and fragmented expected single-copy orthologs. Subsequently, seqkit stats v2.2.0+ (28) is run to compute the number of contigs and N50 for each assembly. After these steps, strains that do not fulfill all the following user-specified thresholds are discarded from subsequent analysis: maximum number of missing (default 2%) or fragmented (default 100%) BUSCO orthologs, minimum N50 (default 50,000 [29]), and maximum number of contigs (default 200 [29]).

When a proteome is obtained for each strain, amino acid sequences are grouped into clusters based on high global sequence identity (90%). This is done by using CD-HIT v4.8.1+ (30) with parameters -M -g 1 -aL 0.70 -aS 0.70 -d 0 -c 0.90, obtaining a representative sequence for each cluster. Using the clustering information, an initial gene presence/absence matrix (PAM) is created, with the cluster IDs in rows, strains in columns, and IDs of strain-specific member genes in cells.

Next, a three-step gene recovery is applied to mitigate possible errors arising from genome assembly or gene calling. Three scenarios are addressed (see Supplementary Information 1.1.1 at https://zenodo.org/records/17799453https://doi.org/10.5281/zenodo.15544430): (i) a premature stop codon that breaks a protein sequence in two pieces; (ii) sequences located in genomic regions overlooked by the gene caller; (iii) sequences overlapping an annotated gene (31) and a previously overlooked region. The PAM is updated with the recovered sequences. When proteomes are inputted, they are directly used in subsequent analysis, and no gene recovery is performed.

Reference-free draft pan-GSMM generation

The representative sequences of clusters are processed to create a draft pan-GSMM using a reference-free reconstruction approach. This pan-GSMM is used to expand the optional reference GSMM with new strain-specific contents; alternatively, when no reference GSMM is provided, it is directly used in downstream analysis.

This reference-free reconstruction phase is based on the bacterial universes (Gram-positive or negative) provided by CarveMe v1.5.2 (5). Even if Gempipe and CarveMe share the same universes (and underlying gene database), the reference-free reconstruction algorithm differs: Gempipe is more conservative and accounts for different gene isoforms while preserving the original enzyme complexes as defined in BiGG (32), a database of manually curated GSMM, leading to enhanced gene-to-reaction association (GPR) rules (see Supplementary Information 1.1.2 at https://zenodo.org/records/17799453https://doi.org/10.5281/zenodo.15544430).

Expanded reference-based draft pan-GSMM generation

When an optional reference GSMM is provided (together with its associated proteome), it is used as the cornerstone of the reconstruction process. First, reference gene IDs are translated into cluster IDs, making the reference GSMM compatible with the PAM and the previously made reference-free draft pan-GSMM. This translation is based on orthologs determination via blastp (20) BRH alignments between each strain and the reference proteome, similarly to the approach suggested by reference 19 (see Supplementary Information 1.1.3 at https://doi.org/10.5281/zenodo.15544430).

The translated reference GSMM is then expanded with new gene clusters, reactions, and metabolites taken from the reference-free draft pan-GSMM used as a repository of new contents, leading to the production of an expanded reference-based draft pan-GSMM. The latter inherits from the reference key features, such as the non-growth associated maintenance energy (4) and the biomass equation (33), which would be otherwise inherited from the CarveMe universe (5). Moreover, during the expansion phase, curated information contained in the reference GSMM in terms of metabolites’ mass/charge and reactions’ balancing is respected (see Supplementary Information 1.1.4 at https://doi.org/10.5281/zenodo.15544430).

The final draft pan-GSMM is subjected to reannotation of its metabolites, reactions, genes, and Systems Biology Ontology terms (34). This facilitates prospective uses of the output GSMMs and can lead to better scores in the community standard test suite MEMOTE (35). This reannotation is mainly (but not exclusively) based on MetaNetX v4.4 (36) (see Supplementary Information 1.1.5 at https://doi.org/10.5281/zenodo.15544430).

Derivation of strain-specific GSMMs

Once the pan-GSMM has been sufficiently curated (the Gempipe API can be used for this task, see Supplementary Information 1.1.6 at https://doi.org/10.5281/zenodo.15544430), it is inputted into “gempipe derive” together with the PAM, producing a strain-specific GSMM for each strain. Briefly, for each strain, a copy of the pan-GSMM is made. If the strain has no genes in a cluster, the cluster is removed, potentially leading to the loss of associated reactions. Similarly, if all genes in a cluster have premature stop codons, the cluster is removed. Next, reactions are iterated while updating their GPR: each remaining cluster ID is replaced by the corresponding strain-specific genes.

Each strain-specific GSMM is then gap-filled using a user-provided recipe for a medium, preferably minimal, known or assumed to support the growth of all the input strains. More than one recipe can be provided, leading to multiple rounds of gap-filling. If no media file is provided, a generic minimal aerobic medium recipe is used, having glucose, ammonia, phosphate, and sulfate as sole C, N, P, and S sources, respectively. The COBRApy (37) gap-filling algorithm is applied, using the pan-GSMM as a source of reactions and a user-selectable minimum flux through the objective reaction. Moreover, the strain-specific gap-filling step can optionally be skipped, which is useful, for example, when auxotrophies have to be studied on minimal media (24). At this point, strain-specific GSMMs have the minimum requirements to be used in simulations.

Before proceeding with “gempipe derive,” the draft pan-model should be curated, for example, by using dedicated functions of the Gempipe API. As an alternative, strain-specific GSMMs can be seamlessly produced from genomes/proteomes by using “gempipe autopilot,” which skips the manual curation by applying an automated gap-filling to the draft pan-GSMM. This gap-filling is prioritized by using penalties derived from alignment metrics (see Supplementary Information 1.1.7 at https://doi.org/10.5281/zenodo.15544430).

Multi-strain predictions and analyses

Once the strain-specific GSMMs have been obtained, specific metabolic features can be predicted, including the capability to catabolize alternative C, N, P, or S substrates, the presence of auxotrophies for amino acids and vitamins, and the potential biosynthesis of specific metabolites (see Supplementary Information 1.1.8 at https://doi.org/10.5281/zenodo.15544430). These features, together with the presence of reactions in strains, are stored as binary feature tables (BFTs). These tables have strains in columns, binary features in rows, and 1 (feature presence) or 0 (absence) in cells.

The Gempipe API contains dedicated functions for multi-strain analysis, where any number of BFTs can be inputted. Briefly, BFTs are combined into a single table, and the pairwise similarity between strains is then calculated using the Jaccard index. This produces a distance matrix, which is further processed to create a dendrogram using Ward’s agglomerative clustering (11). The latter is referred to as a “phylometabolic tree,” where strains with similar metabolic potential are placed closely. Clusters of metabolically coherent strains can be extracted from a phylometabolic tree and their characteristic features identified (see Supplementary Information 1.1.9 at https://doi.org/10.5281/zenodo.15544430). Tutorials for multi-strain analyses using the Gempipe API are available in the Gempipe documentation.

RESULTS

Models’ contents and similarity between tools

To evaluate contents and validate substrate usage predictions of reconstructed strain-specific GSMMs, three datasets were used: “01_klebsiella,” “02_ralstonia,” and “03_pseudomonas,” composed of 37 strains belonging to the Klebsiella pneumoniae species complex (38), 11 strains belonging to the Ralstonia solanacearum species complex (39), and 36 strains of Pseudomonas chlororaphis, respectively (see Supplementary Information 1.2.1 at https://doi.org/10.5281/zenodo.15544430). Comparisons were made with current state-of-the-art reference-free reconstruction tools, namely CarveMe (5) and gapseq (7), as well as a recent reference-based tool focused on strain-specificity studies, Bactabolize (21) (see Supplementary Information 1.2.2, https://doi.org/10.5281/zenodo.15544430). Comparisons were focused on automation; therefore, manual curation was skipped, and “gempipe autopilot” was used.

Contents of reconstructed strain-specific GSMMs were compared against their relative manually curated reference (Fig. 2A; Fig. S1 at https://doi.org/10.5281/zenodo.15544430). Moreover, every tool was compared against each other using the mean Jaccard index of the reaction content (Fig. 2B) (see Supplementary Information 1.2.3 at https://doi.org/10.5281/zenodo.15544430).

Comparison of the general reconstruction metrics. “gempipe_rf” indicates Gempipe ran without reference. (A) Content comparison. “G”: number of genes; “R”: number of reactions excluding exchange reactions; “uM”: number of unique metabolites, i.e., not considering their compartment; “exr”: number of exchange reactions. Bar height corresponds to the mean between strains, while error bars represent standard deviations. Hatched area represents the contents in common with the reference. (B) Similarity between tools based on reaction content. Cells report the mean Jaccard index along the strains, computed for the reaction IDs.

Leveraging its hybrid reconstruction method (expanded reference-based), Gempipe had a generally better reference coverage than reference-free tools, comparable to purely reference-based tools (Bactabolize). On the other hand, output models went beyond the reference with a higher number of modeled contents, aligning with the reference-free tools. The mean Jaccard index confirmed a remarkable reference coverage for Gempipe, where the apparent lower performances compared to Bactabolize are only due to the addition of new reactions during the reference expansion phase. When run without a reference, Gempipe’s models were more similar to the ones produced by CarveMe, likely due to the shared BiGG-based database. However, the introduction of reactions and metabolites in Gempipe is clearly more conservative when using default parameters. gapseq models, despite having a consistently higher number of metabolites, were the most divergent from the reference. However, the conversion between SEED and BiGG IDs provided by MetaNetX (36) is not perfect, so the coverage metrics reported for gapseq could have been underestimated.

Phenotype prediction accuracy

To evaluate the ability to recapitulate phenotypic traits, publicly available binarized Biolog PM data were used as a benchmark (Fig. 3; Fig. S2; Table S1; Supplementary Information 1.2.4 at https://doi.org/10.5281/zenodo.15544430), where the kinetic signal was converted into a binary response “can grow”/“cannot grow.”

Comparison between experimental and simulated Biolog PM growth assays. Bar height corresponds to the mean between strains, while error bars represent standard deviations. “gempipe_rf” indicates Gempipe ran without reference. The “_pan” suffix indicates the use of a manually-curated pan-GSMM as the reference model, instead of a strain-specific one. Data set composition is the following: “01_klebsiella,” 37 strains belonging to the Klebsiella pneumoniae species complex; “02_ralstonia,” 11 strains belonging to the Ralstonia solanacearum species complex; “03_pseudomonas,” 36 strains of Pseudomonas chlororaphis. (A) Outcome of the comparison considering single substrates. TP, true positive; TN, true negative; FP, false positive; FN, false negative. (B) Overall metrics for the substrate utilization prediction.

In general, the accuracy of Gempipe was better or in line compared to the other tools (Fig. 3B), particularly when using its hybrid reconstruction mode. While CarveMe was close to Gempipe in terms of mean accuracy, it was observed that its internal gap-filling algorithm, designed to “enforce network connectivity” (5), tended to maximize the number of substrates for which growth is predicted; this led not only to fewer FNs (Fig. 3A), resulting in better recall (Fig. 3B), but also to more FPs, resulting in detrimental specificity. Despite using the same CarveMe’s assets (gene database and reaction universes), the different implementation of Gempipe better represented the no-growth phenotypes. Benefits are even more evident when modeling genera still without reference GSMMs deposited in BiGG (32), as they cannot be represented in CarveMe’s assets, like, for example, in “02_ralstonia.” Despite the lowering of the default identity threshold and the medium-specific gap-filling step in common with the other tools (see Supplementary Information 1.2.2 at https://doi.org/10.5281/zenodo.15544430), Bactabolize-generated reaction networks resulted excessively gapped in the data set “02_rastonia,” preventing any positive growth prediction (both TPs and FPs) and questioning the utility of Bactabolize-generated models for this single data set.

When a manually curated pan-GSMM was used as the reference model instead of a strain-specific one, the performances of both Gempipe and Bactabolize substantially improved (Fig. 3). Specifically, the purely reference-based reconstruction method implemented in Bactabolize, which subsets a single, coherent set of reference reactions, led to slightly better performances compared to Gempipe. The latter, instead, has to cope with the expansion of the reference, which can introduce spurious reactions affecting simulations. While investing time and resources in the manual curation of pan-models surely pays off in the long term (23), Gempipe provides the overall best-performing option when only a strain-specific curated model is available as a reference, which is the most common case, and simultaneously provides the first draft of the pan-model, thanks to its reference-expansion capabilities, on which the manual curation can subsequently be based.

In addition to substrate usage predictions, gene essentiality predictions were also evaluated (Fig. 4; Fig. S3; Table S2 at https://doi.org/10.5281/zenodo.15544430). In this case, three other different data sets were used, providing transposon insertion sequencing outcomes used as benchmark: “04_streptococcus,” “05_pseudomonas,” and “06_pseudomonas,” composed of 17 strains of Streptococcus pneumoniae grown on THY rich medium (40), 9 strains of Pseudomonas aeruginosa grown on LB rich medium, and the same 9 strains grown on M9 minimal medium (41), respectively (see Supplementary Information 1.2.1 at https://doi.org/10.5281/zenodo.15544430).

Comparison between experimental and simulated gene essentiality assays. Bar height corresponds to the mean between strains, while error bars represent standard deviations. “gempipe_rf” indicates Gempipe ran without reference. Data set composition is the following: “04_streptococcus,” 17 strains of Streptococcus pneumoniae grown on rich THY medium; “05_pseudomonas,” 9 strains of Pseudomonas aeruginosa grown on rich LB medium; “06_pseudomonas,” the same 9 strains of P. aeruginosa but grown on minimal M9 medium. (A) Comparison between modeled genes (background gray bars, with gray point representing individual strains) and modeled genes for which an experimental outcome is available (colored bar, with black points representing individual strains). (B) Outcome of the comparison considering single genes. TP, true positive; TN, true negative; FP, false positive; FN, false negative. (C) Overall metrics for the gene essentiality prediction.

In general, Gempipe modeled a higher number of genes compared to the other tools when leveraging its hybrid reconstruction mode (Fig. 4A). In the context of gene essentiality, reference-based approaches seemed to give better performances; however, while both Gempipe and Bactabolize led to higher precision, the number of modeled genes was much lower in Bactabolize, despite the lowering of the default identity threshold (see Supplementary Information 1.2.2 at https://doi.org/10.5281/zenodo.15544430). Again, the particular gap-filling strategy of CarveMe (see above) emerged: in this context, it hindered the detection of essential genes, leading to lower recall, with respect to the reference-free run of Gempipe, which uses the exact same CarveMe’s assets, but a more conservative gap-filling (Fig. 4B); this difficulty in detecting essential genes is even more accentuated when the genus is not accounted for in BiGG (32), nor consequently in CarveMe’s assets, like, for example, “04_streptococcus.” In this scenario, Gempipe gave a clear improvement in overall accuracy; otherwise, it was in line with the other tools (Fig. 4C). Finally, when reconstruction and simulation of the same strains were based on a minimal medium instead of a rich one, the number of essential genes increased as expected (compare “05_pseudomonas” with “06_pseudomonas”).

Remaining orphan reactions

The quality of reconstructions may also be evaluated by the number of modeled metabolic reactions (not exchanges, sinks, nor demands) which have not been associated with genes and, at the same time, have not been labeled as spontaneous. These reactions, also known as “orphan” reactions (42), can be left by internal gap fillers of automated reconstruction tools to improve network connectivity (5). Ideally, the number of orphans should be minimized by manually checking and reassociation to the corresponding genes: a high number of orphans may indicate insufficient curation and excessive reliance on gap-filling.

The presence of orphan reactions was compared (Fig. 5), and Gempipe reconstructions contained the lowest fraction in every data set. In Bactabolize reconstructions, orphans were copied from the reference, but the fraction in “02_ralstonia” and “03_pseudomonas” data sets was inflated due to the gap-filling. In Gempipe, orphans are also copied from the reference, but during the reference expansion phase, their GPRs are supplemented with missing gene clusters taken from the independent reference-free reconstruction. This resulted in a fraction of orphans lower than the reference in all three data sets, and even lower with respect to the reference-free tools. With respect to CarveMe, in particular, the difference seemed to be more accentuated when the reference was not part of the BiGG collection (32), as in the “02_ralstonia” data set. When Gempipe was run without a reference, the number of orphans reached its minimum; in this reconstruction mode, remaining orphans are a consequence of the two biomass-centered gap-filling steps: the first applied to the draft pan-GSMM, and the second to the strain-specific GSMMs.

Relative number of modeled reactions with no GPR (orphans). Exchange, sink, and demand reactions are excluded, as well as reactions containing the substring “diffusion” in their name. Bar height corresponds to the mean between strains, while error bars represent standard deviations. The hatched area represents orphans in common with the reference.

Metabolic biodiversity of Limosilactobacillus reuteri

Limosilactobacillus reuteri (Lr, formerly Lactobacillus reuteri) is a species of lactic acid bacteria (LAB) adapted to the gastrointestinal tract (GIT) of vertebrates and widely studied for its probiotic potential (43). Recently, six subspecies of Lr were formally proposed by Li et al. (43) and characterized both phylogenetically and phenotypically. Subspecies reflect the host: Lr subsp. reuteri is adapted to humans and herbivores; Lr subsp. kinnaridis to humans and poultry; Lr subsp. porcinus and Lr subsp. suis to pig; Lr subsp. murium and Lr subsp. rodentium to rodents (43).

Gempipe was used to explore the metabolic biodiversity of Lr at the strain level. A total of 1,056 Lr assemblies, of which 597 derived from metagenomes (MAGs), were retrieved. Contigs belonging to other species were removed from MAGs (see Supplementary Information 1.3.1 at https://doi.org/10.5281/zenodo.15544430). The 545 genomes that remained after taxonomy and quality filtering (see Table S3 at https://doi.org/10.5281/zenodo.15544430) were assigned to subspecies based on ANI thresholds reported in reference 43: 69 resulted classified as kinnaridis, 63 as reuteri, 52 as rodentium, 21 as suis, 2 as porcinus, 115 as murium, while the remaining were not classified (see Fig. S4 at https://doi.org/10.5281/zenodo.15544430).

A curated GSMM for Lr JCM1112 (44) was used as a reference in Gempipe (see Supplementary Information 1.3.2 at https://doi.org/10.5281/zenodo.15544430). Due to the hybrid reconstruction mode, it was expanded with new strain-specific contents, generating a draft pan-GSMM. From the latter, 545 GSMMs were derived. BFTs for metabolic reactions, auxotrophies, and growth on alternative C sources were used to build a phylometabolic tree (see Fig. S5 at https://doi.org/10.5281/zenodo.15544430). Since vitamin B12 is unrelated to the ecological niche/subspecies (45), the 24 reactions forming the B12 biosynthetic pathway were removed before generating the tree. A number of clusters equal to the subspecies (6) was extracted. Clusters were generally consistent with the subspecies: reuteri and porcinus were substantially contained in Cluster_3, kinnaridis in Cluster_5, suis in Cluster_1, murium and rodentium in Cluster_2 (Fig. 6A).

Multi-strain analysis using the Gempipe API functions. (A) Phylometabolic tree built using presence/absence data for reactions, auxotrophies, and growth on alternative C substrates. Six clusters of metabolically coherent strains were extracted from the tree and shown alongside the subspecies attribute. Only features not constant across strains are represented. (B) Potential strain-specific production of health-related metabolites, vitamin B12 (“adeadocbl_c”), reuterin (“3hppnl_c”), and histamine (“hista_c”), and presence of urease (“UREA”). The subspecies attribute is reported.

Consumption data for 49 C substrates characterizing the six subspecies were reported by Li et al*.* (43) (“general” data set). These data were derived by the same authors from experiments on 12 strains, two for each subspecies (“specific” data set) (see Table S4 at https://doi.org/10.5281/zenodo.15544430). The “specific” data set was compared with strain-specific simulations, which resulted in 93% accuracy (mean computed on 11 strains, as GCA_000712565.2 was quality-filtered). The “general” data set was compared with the feature’s relative frequency within clusters (see Supplementary Information 1.3.3 at https://doi.org/10.5281/zenodo.15544430); this comparison revealed interesting discrepancies between literature and modeled phenotypes.

For example, growth on galactose is reported as species capability (43); however, ~8% of strains (45), mainly belonging to murium and rodentium subspecies, were predicted to be unable to catabolize this substrate. Upon further investigation, these strains were consistently lacking the aldose 1-epimerase, possibly explaining the deficient phenotype (see Supplementary Information 2.1 at https://doi.org/10.5281/zenodo.15544430). Among compared substrates, D-xylose was particularly interesting for three reasons: (i) its catabolic route was not in the reference GSMM, but it was automatically included by Gempipe during the reference expansion phase; (ii) its utilization was predicted to have high variability across strains, ranging from ~9% in murium to ~88% in rodentium; (iii) its simulation accuracy with the “specific” data set was 100%. For instance, only ~38% of the 21 strains strictly classified as suis were predicted to grow on D-xylose, while this is reported to hold for the entire subspecies (43). This further corroborated the idea that phenotypic descriptions of species/subspecies provided in literature are actually not always valid, possibly because they are based on a number of strains (two in this case) too low for generalization, or due to the limited reliability of phenotypic testing (46). In this context, the GSMM-based analysis gave more comprehensive results than traditional phenotypic characterization.

In humans, Lr is found in many sites other than the GIT, including breast milk, skin, and urinary tract. Lr is reported as a probiotic, and some of the metabolites it produces have a proven health-related effect (47, 48). For example, reuterin (a mixture of different forms of 3-hydroxypropionaldehyde) is released by several strains and exerts antimicrobial activity against Gram-negative bacteria, making Lr effective against GIT infections (47, 49). Histamine is another strain-specific metabolite, acting as an intestinal immunomodulator and anti-inflammatory agent (47, 49). Vitamin B12 (cobalamin) is an essential vitamin in humans, introduced with the diet; only four B12-producing Lr strains were clearly identified as of 2018 (49). Apart from health-related metabolite production, other metabolic features are of interest to assess adaptation strategies in different hosts; one of them is the conversion of urea to ammonium and CO_2_ (urease), likely involved in acid resistance (50).

Gempipe enabled a systematic evaluation of these four key metabolic features in all the 545 filtered strains (Fig. 6B). The majority of rodent-associated strains (Lr subsp. murium and rodentium) were predicted to be capable of urease activity, probably needed for survival in rodent GIT (51, 52); this is in accordance with a recent genome-wide association study (53). Moreover, most but not all of the human- and poultry-associated strains were predicted as potential producers of reuterin and vitamin B12, with the underlying genes being indeed part of the same pdu-cbi-cob-hem gene cluster (51). Interestingly, histamine production was an exclusive feature of subsp. reuteri (43), with only a few reuteri strains lacking this trait.

DISCUSSION

In the present work, Gempipe was introduced, a multi-purpose package for pan and multi-strain genome-scale metabolic modeling. It adopts a hybrid reconstruction approach that lies between reference-free and reference-based methods, as an optional reference is expanded with new contents coming from a universal GSMM. Together with an internal clustering of strain-specific genes, an effective in-depth gene recovery (see Supplementary Information 2.2 at https://doi.org/10.5281/zenodo.15544430), and conservative generation of GPRs, the implemented approach proved to be effective for multi-strain reconstructions, with better or similar performances compared to current established reconstruction tools, when focusing on metabolic features without considering manual curation. A summary of the main features of Gempipe compared to the benchmarked tools is provided in Table S8, https://doi.org/10.5281/zenodo.15544430.

Gempipe represents a third option in the panorama of reconstruction tools: reference-based methods use a manually curated model (or sorted list of models, from the phylogenetically closest to the most curated) to be used as template, from which reactions are copied after ortholog genes are determined (39, 54, 55); instead, reference-free methods use a semi-curated universal model with a generic biomass equation, from which reactions are copied and gap-filled based on sequence homology (one-way alignment) (5, 7). From the first approach, Gempipe inherits the ortholog determination to translate reference genes into equivalent gene clusters; from the second, Gempipe inherits the homology-based insertion of new, strain-specific reactions, with GPRs already based on gene clusters.

In the case of multi-strain reconstruction, reference-based methods are usually applied (11–18). The main concern, in this context, is that the template must be representative of the metabolic diversity of the entire species (or genera) (21); otherwise, strain-specific reactions must be added afterward to each generated model (19). This is why a comprehensive (pan) GSMM is usually curated prior to running reference-based reconstruction tools (23). However, such pan-GSMMs representing a species (or genus) are complex to build and still rare in the literature (12, 18, 23, 24), with as little as four species reported (21). Indeed, in some GSMM-based biodiversity studies, strain-specific models have been used as a reference (15–17). This is a limiting approach, as generated models will just be a subset of a single strain. In this context, Gempipe was able to grasp strain-specificities better than purely reference-based methods like Bactabolize (21) when a strain-specific GSMM is used as reference instead of a pan-GSMM, which is a common case. Moreover, inheriting contents from a manually curated reference, the models generated by Gempipe provided more accurate predictions with respect to those built with reference-free methods like CarveMe (5) and gapseq (7). All this was achieved while minimizing the number of orphan reactions and maximizing the MEMOTE metrics (35) (see Supplementary Information 2.3 at https://doi.org/10.5281/zenodo.15544430).

Given the limited number of pan-GSMMs, Gempipe also emerges as a valuable tool to quickly build and curate pan-GSMMs to be used in biodiversity or bioprospecting studies. It must be noted, however, that the concept of pan-GSMM was also introduced in the context of metagenomics-derived models (56), which is remarkably different. In this context, indeed, pan-GSMMs are instrumental to cope with the incompleteness and contamination of MAGs, and strain-specificity (in terms of genes and reactions) is lost in favor of a consensus/mean reconstruction representing a species-level genome bin (56). Instead, in the GSMM-guided exploration of biodiversity, the context where Gempipe operates, strain-specificity must be retained and emphasized as it is needed for the subsequent generation of strain-specific GSMMs (19).

In the development of Gempipe, particular attention was put into how metabolic features are encoded in the model. In this regard, it must be noted that metrics commonly used to evaluate and compare reconstructions (accuracy, precision, recall, specificity) (5, 7, 21) do not take into consideration the faithfulness of the reaction network in representing the organism. Indeed, behind a true positive match with the experimental data (e.g., Biolog screenings), there could be cases of (i) wrong reaction mechanisms (e.g., erroneous transporter types); (ii) presence of metabolic reactions not supported by genes (orphans); (iii) reactions with seriously impaired GPRs, lacking many components of a protein complex. Provided that manual curation still remains essential for obtaining truthful GSMMs, users should be aware of the reconstruction principles that each tool follows to draft a model. Gempipe includes reactions only when a protein complex is fully supported by genes; otherwise, mismatches (false negatives) will lead the manual curation in closing the gaps. On the other side, tools providing (i) internal gap-filling procedures aimed to improve the network connectivity and (ii) too permissive mechanisms of GPR generation and reaction inclusion could lead not only to predict a higher number of growth-supporting substrates (higher recall and lower specificity) but, most importantly, to match the experimental data with a biologically inaccurate representation of the metabolism, which can be easily overlooked during the manual curation.

Finally, Gempipe is not meant to be just a reconstruction tool, but also an analysis tool in the context of biodiversity exploration. Indeed, while the Gempipe API includes functions to curate models, it also contains functions to analyze the deck of strain-specific GSMMs in output. For example, functions are included to cluster strains according to their metabolism and to visually compare metabolic clusters with respect to other attributes, such as niche metadata or the formal species or subspecies classification. This aids users to achieve goals including (i) the screening of strains for desired metabolic traits, (ii) the classification of strains according to their metabolic capabilities, and (iii) the definition of species and subspecies, and possibly other taxonomic ranks, in terms of their core metabolic potential. In this context, the case study on L. reuteri here reported provided insights also from a taxonomic point of view (see Supplementary Information 3.1 at https://doi.org/10.5281/zenodo.15544430).

Limitations of Gempipe are mostly due to the resources it relies on. The BiGG namespace (57) was adopted as it was convenient for two main reasons: (i) the availability of high-quality, manually curated reference models based on the same namespace (32); (ii) the human-readable IDs, particularly useful during manual curation, when GSMM-based metabolic maps have to be hand-drawn or interpreted (58, 59) (for example, D-glucose is “glc__D,” immediately recognizable with respect to its ModelSEED (60) equivalent “cpd00027”). While convenient, the BiGG database has imperfections due to its structure: it is not a coherent biochemical database, but rather a collection of models from which a biochemical database is derived. Depending on the model of origin, the same metabolite or reaction can be defined differently. Consequences are many: (i) the same reaction can be represented with different reversibility (e.g., “ILETA2”); (ii) the same metabolite can be represented with different IDs, leading to duplicate metabolites (e.g., “ind_c”/“indole_c”); (iii) the same reaction can use duplicate metabolites, leading to duplicate reactions (e.g., “TRPS2”/“TRPS2_1”); (iv) metabolites with the same ID can be represented with different chemical formula or charge (e.g., “fmn_c”). Therefore, the integration of a BiGG-based model with reactions coming from another BiGG-based model can lead to the introduction of unbalanced reactions or, even worse, stoichiometric inconsistencies (61). Gempipe tries to circumvent this issue by enabling users to superimpose particular metabolite charges and formulas or reaction balances during the reference expansion phase.

Another limitation, depending on BiGG (32), is the limited representation of bacterial diversity, mainly in terms of genes. Indeed, in its current version (v1.6), BiGG contains as little as 108 GSMMs, of which 88 are prokaryotic, and only 22 do not belong to the Escherichia or Shigella genera. The set of BiGG genes, used by Gempipe and CarveMe (5), is therefore biased on model species and does not cover much bacterial diversity, potentially leading to missed genes (and thus reactions) during the reference-free reconstruction phase. Gempipe tries to limit this problem by comparing the eggNOG-mapper (62) functional annotation of clusters’ representative sequences (see Supplementary Information 1.1.2 at https://doi.org/10.5281/zenodo.15544430). In future versions, the set of BiGG genes could be expanded by considering BiGG-compliant GSMMs stored in BioModels (63).

GSMMs are useful tools to catch and explore strain-level differences in metabolism. Secondary metabolite production, like antibiotic production, can in principle be described by GSMMs, as long as the underlying biosynthetic pathways can be dissected in terms of stoichiometric equations and involved enzyme complexes. However, to date, the coverage of secondary metabolism in GSMMs is generally poor, even in manually curated models. Indeed, aside from a few well-studied groups, such as polyketides and nonribosomal peptides, many biosynthetic pathways for strain-specific secondary metabolites present knowledge gaps even in large, general-purpose metabolic databases like MetaCyc (64) and KEGG (65). The representation of secondary metabolism is then even worse in GSMM reconstruction tools, especially in BiGG-based tools (see above). Therefore, the modeling of secondary metabolism still heavily relies on manual curation, which, however, may remain insufficient due to the current knowledge gaps. Another major challenge comes from the modeling framework itself. Plain FBA, which does not account for gene regulation, can have good predictive or descriptive potential only for growth-coupled metabolites under steady state. This means that classic GSMM/FBA often fail to compute fluxes in secondary metabolisms, as they are typically tightly gene-regulated and triggered during the stationary phase, in stress conditions, or by environmental signals. Interested readers are referred to reference 66 for one of the most comprehensive and updated reviews on the topic.

Further development on Gempipe will benefit from the introduction of new universes, such as one for yeasts (67), and from new API functions for the analysis of the deck of strain-specific GSMMs created.

In conclusion, Gempipe will facilitate metabolic biodiversity studies for a wide range of bacterial species, including those not having a dedicated pan-GSMM, which are currently the large majority.