A curated C. difficile strain 630 metabolic network: prediction of essential targets and inhibitors

Creation of the network

The genome of C. difficile strain 630 is composed of a circular chromosome of 4,290,252 bp coding for 3968 open reading frames (ORFs) as well as a plasmid of 7881 bp coding for 11 ORFs [25]. The C. difficile strain 630 draft reconstructed metabolic network presented here covers 20.3% of the ORFs present in the chromosomal genome of the bacteria with 806 ORFs. These 806 genes code for proteins catalyzing 769 metabolic, 145 transport and 117 exchange reactions. A total of 592 unique metabolites (703 in total not considering extracellular or intracellular state) are involved in the 1031 reactions in the network. The coverage of the genome is similar to those of previously published reconstructed metabolic networks such as B. subtilis with 20% [23] and higher than the most recent network of C. acetobutylicum with 13.0% [26]. Most reactions have at least one gene association (77.9%). Reactions without any gene association were added based on the existence of evidence from the literature such as in the case of Stickland reactions [27]-[30], presence in databases such as xanthine amido hydrase or to fill functional gaps to obtain a functional network as in the case of putative transporters for end-products of fermentation.

The final version of the network is available in 3 different formats: 1. An excel file that shows on different spreadsheets the reactions, metabolites, genes, and compartments that comprise the network and the definitions of the network based on the standard described in the RAVEN toolbox [31]. This file is meant to be easily readable by humans; 2. A tab-separated format of the network amenable to analysis in the R Statistical Computing environment (www.r-project.org) using Sybil [32]; and lastly, 3. A SBML Level 2 formatted network [33] that can be used with tools such as Matlab or other SBML compliant software. While a naming convention has been suggested recently for metabolic reconstructions [16], we feel that a naming convention that does not allude to the name of the organism is insufficient. Therefore, in the present work, the C. difficile strain 630 metabolic network reconstruction is called i MLTC806cdf as per the suggested convention with the added cdf suffix denoting the KEGG [34] three-letter organism ID representing C. difficile strain 630. The SBML version of the model has been deposited to the BioModels database [35] and assigned the identifier MODEL1409240004 .

Validation

Essential metabolites

In order to validate the network, experiments involving the removal or addition of certain metabolites from the media were reproduced in silico. Each essential amino acid (cysteine, leucine, Isoleucine, proline, tryptophan and valine) was confirmed essential in the network as their removal prevented the production of biomass in any medium. None of the three essential vitamins are essential in the network (Additional file 1: Table S1). The essentiality of two of these, biotin and pyridoxine, is due to their implication in the regulation of processes which could not be simulated in the metabolic network. Panthothenate (which is important in lipid metabolism) was in the past determined to be essential for a number of C. difficile strains tested [37]. More recently, a new ketopantoate reductase (KPR) gene panG was discovered in a number of pathogenic bacteria and found to have a homolog in C. difficile strain 630 [40]. Therefore, whereas panthothenate is commonly thought to be essential, this essentiality is strain specific and absent in C. difficile strain 630. A ?panG mutant in Francisella tularensis did not have any differences compared to wild type infections in a mouse model for pneumonic tularemia [40]. This is likely due to the fact that panthothenate (vitamin B5) is widely available in food and most bacteria are able to import panthothenate through a sodium co-transport mechanism [40].

Non-essential metabolites

Removal of non-essential metabolites did not have an important effect on growth. The non-essential amino acid methionine is known to enhance growth of the bacteria and is used in the minimal medium to increase growth rate. Interestingly, the removal of methionine from the minimal medium leads to a slight reduction in biomass production (less than 1%), a small but qualitatively correct effect. We simulated the removal of arginine and histidine (both non-essential amino acids) from the rich medium and in both cases this lead to the qualitatively correct result of a decrease in biomass production in agreement with the experimental evidence [36]. Most (8 out of 11) non-essential amino acids that when removed from a complete media do not affect growth rate experimentally also have no effect in silico when removed from the complex medium. Furthermore, the addition of certain non-essential amino acids (8 out of 13) or any nucleosides to the minimal medium leads to an augmentation of biomass production in agreement with experimental data (Additional file 1: Table S1).

Carbon sources

The utilisation of different carbon sources in the absence of glucose in the network was simulated and compared to experimental data. Such data is not always specific to C. difficile strain 630 and small differences among strains do exist [36] (Additional file 1: Table S2). The removal of glucose produced the largest decrease in biomass production (~33%) but had a smaller effect in the complex medium (~6%). For 15 experimentally tested carbon sources (out of 20 carbon sources tested in silico) we obtain 100% agreement between the predicted utilization of alternative carbon sources by C. difficile, including the impossibility to use lactose as a carbon source. Furthermore, we predict that C. difficile strain 630 should not be able to use rhamnose or myo-inositol (for which there is some evidence of usage in other strains but no experimental information for strain 630) while malate, glycerol and chorismate would lead to increased growth rates.

Comparison with existing metabolic network reconstructions

We also compared i MLTC806cdf to a curated metabolic network reconstruction of the closely related bacterium Clostridium acetobutylicum. Three curated metabolic networks exist for this organism [26],[42],[43] all focusing on metabolic engineering of the bacteria to maximize the production of butanol. The latest network [26], called i CAC490, was chosen for the comparative analysis as it was the largest, the most recent and the only one available in a commonly accepted usable format (SBML). The analysis of reactions shared between the two reconstructions was possible due to the extensive use of KEGG identifiers in both networks. In the case of transport reactions, two transporters were considered similar if they transported the same molecule with one transporter from one network potentially matching multiple transporters in the other network. We did not differentiate between phosphoenolpyruvate (PEP):carbohydrate phosphotransferase system (PTS) [44], ion channels [45] or ATP driven transporters [46] as long as the transported molecules were the same in the two networks. i MLTC806cdf contains 174 more reactions than the reconstructed network of C. acetobutylicum and 416 in common with it (representing 45.5% of our reactions). i MLTC806cdf has 62 less unique metabolites and 389 metabolites in common (representing 66% of our unique metabolites). The presence of 186 (28% of the metabolome) dead-end metabolites (a complete list can be found in Additional file 1: Table S3) suggests that the C. acetobutylicum network still has incomplete pathways and gaps that could affect biomass production. A detailed comparison is presented in Table 2 or in the form of a Venn diagram in Additional file 1: Figure S2.

Single gene deletions

We performed an in silico gene deletion study using both Synthetic Accessibility (SA) and Flux Balance Analysis (FBA) on i MLTC806cdf in order to identify potential essential genes that may lead to the discovery of novel therapeutic targets. This analysis removed reactions that were catalyzed by each gene alone (or in pairs, next section) or by a complex that involved that gene product and then measured the capacity of the network to produce biomass (either a flux in biomass production for FBA or S_net for SA, see Methods). The complex medium (as describe in Table 1) was the one used for the gene deletion studies since it reproduces the high concentration and diversity of nutriments found in the intestinal tract. Furthermore, the medium used is the same as the simulation of the Bacillus subtilis metabolic network [23], which in turn is an approximation of the one used to perform the experimental validation of lethality of single gene deletions in that organism [47].

Essential genes were compared to experimental results for the gram-positive bacteria Bacillus subtilis[47], which is the closest relative of Clostridium difficile with experimental essentiality data for all of its genes. Since the simulation involved the deletion of genes via deletion of metabolic reactions, functional homologs (genes responsible for reactions that share the same EC number and catalyze similar reactions) were used for the comparison. Among the 76 genes with a predicted effect on biomass production in C. difficile, 46 have homologs that are essential in B. subtilis, 8 did not have any functional homolog and 22 had a homolog that was not essential in B. subtilis (Table 3). We were able to compare 618 C. difficile genes (76.8% of the genes in i MLTC806cdf) for which we could detect a B. subtilis functional homolog with an overall prediction accuracy of 89.2% (Additional file 1: Table S4). A similar rate of 89.0% was obtained based on the comparison of 525 genes using sequence homology (E value <1e^?5, sequence identity above 30% and alignment overlap over 80% of the C. difficile sequence, Additional file 1: Table S4). While an accuracy of around 89% is extremely high, it is important to keep in mind that despite being closely related, differences are expected between the B. subtilis and C. difficile.

We utilized i CAC490 to perform the prediction of essential genes in C. acetobutylicum and compared the results to those above for C. difficile. Essential genes in the C. acetobutylicum network were identified in a similar manner than in i MLTC806cdf, using Sybil [32] within the R environment for statistical computing. Of the 658 C. difficile proteins with C. acetobutylicum sequence homologs, 368 are present in i CAC490. Based on these 368 proteins, we obtain an agreement of 72.3% between i MLTC806cdf i CAC490 (Additional file 1: Table S5). It is interesting to note that we obtain a higher agreement between our computational results for C. difficile and experimental results for B. subtilis than with computational results for i CAC490 representing a more closely related species.

The difficulty of performing genetic manipulations in C. difficile is notorious and severely restricts our ability to compare our results with experimental information. Despite the lack of extensive information on experimentally verified essential genes, the little evidence that exists, supports our predicted essential roles for a number of genes: metK [49], guaA [50] and ntpA-B-C-D-E-F-K [51]. Other known essential genes such as secA1-A2 [52], metG and gyrA-B [53], trpS [49] and gldA [54] are not present in the network and are involved in non-metabolic processes. Lastly, the gene prdF has been mutated and was shown to be non-essential [55], in agreement with its predicted non-essential role in i MLTC806cdf.

We compared the list of 76 genes predicted to be essential in C. difficile using i MLTC806cdf to the genes in the Database of Essential Genes (DEG) [56]. Interestingly, and serving as further validation of i MLTC806cdf, a total of 61 of these genes are present in DEG, i.e., these genes have homologs known to be essential in other species. The remaining 15 predicted essential genes that are not present in DEG (Table 3) include ntpA,D,E,F,K (all subunits V-type ATP synthase) as well as metF (involved in amino acid synthesis) and xpt (guanine synthesis) among others (Additional file 1: Table S6).

The inhibition of the product of essential genes that are upregulated during CDI may require a smaller drug dose to generate an effective response, thus decreasing side effects. We utilized transcriptomics data associated to the differential expression of genes during infection [57] to annotate predicted essential genes in view of their potential use as therapeutic targets. Eight predicted essential genes are downregulated in vivo while 6 are upregulated during infection. Some of the genes that are upregulated during infection and predicted to be essential such as fabD, serA or CD2549, are of additional interest as they could not only affect growth, but also colonisation and pathogenesis processes [57] (Table 3).

Detection of potential human cross-reactivity targets

Local structural similarities

There is a possibility that potential cross-reactivity targets may perform different functions (EC numbers) and have little sequence similarity yet still have sufficient 3D atomic binding-site similarities to be inhibited by a drug developed against a C. difficile target. We proceeded to analyse binding-site similarities to detect potential cross-reactivity targets for the 24 C. difficile proteins where sequence or functional similarity search did not detect any human homolog. Binding-site similarities are measured in terms of detected geometrically and chemically equivalent atoms in common between two binding-sites [58],[59]. We found an average of 36 atoms in common between binding-sites (average p-value 0.038 and z-score 3.05) for 21 out of 24 models with members of Pfam families [60] that include human proteins. In three cases, ispF, CD2549 and dapH, no significant level of binding-site similarity was found to a Pfam family that contains human homologs (Additional file 1: Table S7).

It is hard to judge if the found matches are significant or not considering that no threshold for binding-site similarity can be uniquely defined above which cross-reactivity is certain [59]. However, in 17 cases out of 24 cases, the top-scoring detected binding-site similarities for each case represent binding-sites in proteins that bind ligands that are similar to at least one of the substrates of the reaction catalyzed by the modelled C. difficile protein. In seven of those cases, the top matching Pfam family that contains human homologs binds a similar ligand (Additional file 1: Table S6). Taking as an example the case of the enzyme encoded by the asd gene, we detect 39 atoms (Z-score 3.92, p-value 0.012) in common to a glyceraldehyde-phosphate dehydrogenase from spinach (PDB ID 2PKQ) bound to NADPH, a member of Pfam family PF00044 that has human homologs (Figure 1). Five out of the top 7 most similar binding-sites also bind NADPH or NADP, all from different Pfam families. The superimposition of these diverse binding-sites based on their similarities to the asd gene product binding-site leads to an extremely good superposition of their respective bound-ligands (Inset Figure 1). This suggests that the detected similarities are biologically significant. The quality of the resulting superimpositions together with the detection of similarities across families that bind similar ligands to those that bind the C. difficile targets reinforces the confidence in the biological significance of our predictions. The quality of the alignment of the NADP molecules across different families via the detected similarities suggests that these capture the molecular determinants responsible for binding. As such determinants are conserved across protein families, there is a possibility that these are also conserved within families and thus present in the human homolog. In most cases, the similarities that were detected actually represent commonly used cofactors or other ubiquitously used ligands, such as NADP above or ATP. These results do not necessarily mean that a drug targeting the C. difficile protein will bind the human homolog belonging to the detected Pfam families, but these should be used as potential cross-reactivity targets in the rational design of inhibitors against the C. difficile protein in question. Furthermore, given that the detected similarities focus on common cofactors and ubiquitous molecules such as ATP, the results also suggest that targeting the sub-pockets of less ubiquitously used substrates may reduce the chance of cross-reactivity.

Double gene deletions

We performed double gene deletions to identify potential polypharmacological targets and to target reactions that are catalyzed by isoenzymes. Based on FBA analysis, 203 gene pairs involving 69 unique genes that had small or no effect in single gene deletion were deleterious when removed in pairs. An additional 3 essential gene pairs involving 6 new unique genes were found using SA. Eight gene pairs were considered essential in both SA and FBA (Additional file 1: Table S8).

Some double mutants show a synergistic effect, defined as an effect greater that an additional 1% reduction in biomass production in the double mutant than the sum of effects of each single mutant. The 39 synergistic double mutants were analysed in more detail (Additional file 1: Table S9). Thirteen of these synergistic gene combinations resulted in total abolition of biomass production in FBA or prevented the biosynthesis of at least one element of the biomass in SA, 11 had an effect between 10% and 20% and the remaining 15 had a small effect on biomass production (between 5 and 10%). Twelve of the essential pairs of genes are isoenzymes that catalyze the same reactions. Sixteen gene pairs represent enzymes involved in pathways with the same functional category while the remaining 11 gene pairs affect different pathways (Figure 2). All essential pairs identified by both SA and FBA are isoenzymes whose removal results in a total arrest of biomass production.

Isoenzymes usually result in a higher biomass loss than relatively distant pairs (Figure 2). For the 12 isoenzymes, the deletion of the two genes in a pair is required to remove a reaction that is catalysed by both. For the remaining essential gene pairs, the reactions associated with both genes are used in parallel in the wild type, a case of metabolic plasticity [61], or the reactions from only one of the genes is used while the reaction from the other member of the pair can act as a backup, a case of metabolic redundancy [61]. Depending on the category in which an essential gene pair falls, different strategies may be required in order to target the pair [61]. From the 39 essential genes pairs, 23 represent cases of plasticity of the network and 4 cases of redundancy in the network (Additional file 1: Table S9).

Distribution of predicted essential genes across pathways

In order to see which parts of metabolism are more enriched with essential genes, we classified reactions into 8 functional pathway classifications (Additional file 1: Table S10). A hypergeometric test [62] for over and under-representation was performed to identify pathways enrichment in essential genes. This analysis confirmed that linear pathways, like lipid synthesis, tend to have more essential genes (25 out of 71 single gene deletions, overrepresentation p-value 1.81e^?10) due to the lack of alternative ways leading to the production of biomass metabolites (Figure 3, Additional file 1: Table S11).

Creation and curation of the draft metabolic network

We created a collection of metabolic and transport reactions associated with C. difficile strain 630 from the KEGG [34], MetaCyc [78] and TransportDB [79] databases. Reactions involving polymers (glycogen, peptides and others) were ignored to avoid the spontaneous creation of matter due to the presence of molecules of undefined length. Also, some of these polymers, like glycogen, are only used for energy storage [80] and would not have any impact in simulations that optimize biomass production. Reactions involved in spore formation, the conjugation process, RNA, peptide or DNA modification, cell repair, and other non-metabolic enzymatic reactions were not added to the network since these reactions are not directly contributing to the production of biomass constituents. Finally, the reconstructed network concentrates solely on metabolism without including signalling, gene regulation and post-translational modification of proteins even if these can clearly affect metabolism.

The initial draft network is little more than a collection of reactions. In that state it cannot be used for any sensible application such as the simulation of biomass production. This is due to the numerous inaccuracies present in the databases. These inaccuracies consist of, but are not limited to, the presence of generic and dead-end metabolites, missing or erroneous pathways, missing genes and unbalanced reactions (Additional file 1: Table S12). Two cases exemplify some of these inaccuracies. First, noting that aerobic pathways should not be present at all in the C. difficile, the inclusion of incomplete versions of both aerobic and anaerobic pathways for the biosynthesis of vitamin B12 is problematic (Additional file 1: Figure S3). Second, the omission of most reactions in the Stickland and amino acid fermentation pathways (Additional file 1: Figure S4), important sources of energy for the bacteria, which had to be completed based on literature. Both of these corrections a numerous others in of the same nature were necessary to create a functional network (Figure 4).

Most transport reactions were based on TransportDB [79] although putative transporters not present in the databases had to be added for molecules known to be exported or imported by the bacteria. An exchange reaction (reaction that simulates interaction with the media via the appearance or disappearance of the given metabolites in the network) was created for each metabolite with an extracellular version. These exchange reactions are set to allow the export of a metabolite unless it is part of the tested growth medium in which case import is also possible.

The curation not only involved the addition and suppression of reactions from the initial draft, many characteristics of each reaction such as the directionality, the presence of complexes or the assignment of gene-reaction associations and their inclusion as part of a pathway had to be analyzed manually. The directionality of reactions was based on information obtained from the MetaCyc database when available. Reactions found exclusively in KEGG were kept bidirectional unless leafing to the production of a highly energetic compounds (ATP, NAD+, NADP+, etc.) that were not known to be produced in such a manner. Examples of reactions producing highly energetic compounds are ATP synthases, amino acid fermentation and glycolysis.

Protein-protein interactions are highly important in the network since they can modify the essentiality of a gene based on knowledge of the involvement of its protein product as part of a protein complex. The possibility that a protein complex is responsible for catalysis was investigated for every reaction that had more than one gene associated to the reaction and for any gene whose protein is identified as a subunit of a complex. We used information from TransportDB, UniProt, Brenda and literature data (either for the C. difficile protein of interests or a homolog of same function that could indicate a similar interaction between genes products). In every case the STRING version 9.1 [81] score was calculated to evaluate the confidence score attributed to each complex (Additional file 1: Table S13).

We used KEGG pathway identifiers since they are generic and allow maintaining a small number of pathways in the network. Pathways containing less than five reactions were manually merged to the closest relevant pathway. Reactions that were not associated to any pathway in KEGG (32% of the reactions in KEGG) were linked to existing pathways. The pathway assigned to the reaction was the one in which the less linked metabolite of the reaction was most often present in the network (Additional file 1: Table S10).

Lipid biosynthesis and cell membrane composition

A problem with lipid biosynthesis in metabolic networks is the fact that a variety of different fatty acids exist in cells with different lengths and saturation states and can be used in numerous cellular processes, both metabolic and non-metabolic. Also, the lipid composition of the cellular membrane is a mix of various phospholipids and glycerolipids of varying lengths and composition and this composition varies depending on the growth conditions [82]. A complete definition of lipid metabolism is likely impossible to define at this moment given the lack of experimental information specific to C. difficile membranes. Furthermore, the resulting network would be dominated by lipid reactions with the same few genes repeated for every possible length and saturation state for every lipid type. Therefore, having an exhaustive definition of lipid metabolism would not bring any additional relevant information on metabolism. In most reconstructed metabolic networks, lipids are used almost exclusively in membrane formation. Following [42], in order to reduce the complexity of lipid metabolism while keeping it as close to the real bacteria as possible, the fatty acid composition of all lipids is held at a constant 16:0 (carbon chain length: number of double bonds), which is the most abundant fatty acid observed in C. acetobutylicum[83]. The last step in simplifying lipid biosynthesis was to combine the elongation of fatty acids from acetyl (2:0) to palmitate (16:0) into a single reaction and beta-oxidation of the palmitate back into acetyl as another single reaction.

Biomass

In order to simulate bacterial growth, the biomass (ensemble of macromolecules necessary for cellular growth and division) composed of DNA, RNA, cell wall, proteins, solute pool and lipids was defined based on the following elements. Nucleic acid composition of DNA used in the network was calculated based on the nucleotide content of the genome and plasmid of Clostridium difficile strain 630, RNA composition is based on the content of the transcriptome using the UCSD genome browser [84] and protein composition using the proteome. Lipid, cell wall and solute pool composition as well as the overall biomass composition were taken from the C. acetobutylicum network [42] due to a lack of literature specific to C. difficile. A detailed composition of the biomass is available in supplementary Additional file 1: Table S14.

Simulation of growth and gene essentiality

Two methods were used to simulate bacterial growth and determine gene essentiality: flux balance analysis (FBA) [21] and synthetic availability (SA) [22].

Flux Balance Analysis (FBA)

Flux balance analysis [21] is a constraint-based modeling method commonly used in the study of genome-scale metabolic networks. FBA firsts creates a stoichiometry matrix (S) from the network where each row represents a metabolite and each column a reaction. The values in this matrix correspond to the stoichiometry of the metabolite in the reaction with a negative number representing consumption and a positive number representing production of the metabolite. A system of linear equations is produced by multiplying S with a column vector v representing the fluxes through each reaction. FBA creates a steady-state distribution of fluxes where the product of the previous multiplication must equal zero S???v?=?0. Since the resulting system of linear equations is undetermined, FBA uses linear programming to maximize a particular objective function Z, in our case biomass as the representation of growth, by maximizing the multiplication of a row vector c containing the weight of each reaction on Z with the column vector v used previously (maximize Z?=?c???v). Values in v are constrained by lower and upper bounds representing various factors like enzyme directionality of the reaction, capacity, uptake, secretion rates, etc. FBA finds a distribution of fluxes in the network respecting the constraints on v and maximizing Z at the same time. In our studies, the Sybil package version 1.1.11 [32] available for the free R environment for statistical computing (version 2.15.2) was used in order to run FBA simulations. Other tools also exist that can use FBA with different interfaces like the COBRA package [85] that runs on the proprietary Matlab computing environment.

Synthetic Accessibility (SA)

Synthetic accessibility [22] is a parameter-free method to predict the essentiality of genes through their deletions in metabolic networks. SA uses network topology to calculate the number of reaction steps needed to produce the outputs (biomass) of the network from the inputs metabolites available in the growth medium. SA works by examining all reactions that use only input metabolites and marks those reactions and their products as 'accessible'. In an iterative manner, successive iterations search for new reactions that have all required substrates marked as accessible until no new reaction can be added. Each metabolite j have a synthetic accessibility value S _j representing the number of iterations before this metabolite became accessible and the Synthetic Accessibility of the network as a whole, S _net , is the sum of S _j of each of the output metabolites. Increases in S _net resulting from gene deletions are predicted as being deleterious. Synthetic accessibility is a simpler method than FBA and gives comparable results on the prediction of essential genes demonstrating that the topology of the network is the principal factor influencing essentiality. We use our own implementation of the algorithm.

Local structure similarity

The identification of local structure similarities is separated into three steps: the creation of a three dimensional model of the protein, the identification of the probable binding site, the comparison of the binding site to a dataset of known binding sites. The models used where created using I-Tasser [86], a tool which builds 3D models based on multiple-threading alignment and iterative fragment based simulations. Detailed methodology for I-Tasser can be found elsewhere [86]. We used the Isocleft Finder (Kurbatova et al., [58]) web-interface to compare the largest cavity of each of the models (representing the binding site in 83% of cases [87]), identified using their own implementation of the SURFNET algorithm [88], to a non-redundant dataset of 7339 binding-sites of unique combinations of protein families bound to distinct ligands [58] using Isocleft [58],[59]. IsoCleft is a graph-matching based method for the detection of structural and chemical similarities between pairs of protein cavities. In each case, we selected the most similar match found to a protein in the IsoCleft Finder non-redundant dataset that contains human homologs. It is important to note however that the representative of such family in the non-redundant dataset may itself not be necessarily a human protein.