Redefining genomic view of Clostridioides difficile through pangenome analysis and identification of drug targets from its core genome

Introduction

Clostridioides difficile, earlier known as Clostridium difficile, is a toxin-producing Gram-positive, anaerobic bacteria (1). During infection, it releases toxins that disrupt the intestinal epithelia, resulting in a variety of diseases ranging from mild, self-limiting diarrhoea to the fatal pseudomembranous colitis (PMC) (2-5). During the past two decades, there has been dramatic increase in the incidence and severity of C. difficile infection (CDI) (6,7). CDI is usually followed by the antibiotic treatment that impairs the protective gut microflora (8).

C. difficile was first identified from microbial flora of faeces of healthy newborn infants and was considered that it has no deleterious effects in human (9). But later on, it was identified as the cause of antibiotic-associated PMC (10). C. difficile has been reported to have genetic heterogeneity because of its wide ecological adaptability. Hence, in the past 20 years, significant changes in CDI epidemiology have been reported (11). The differences in the severity of the infection, presence of pathogen at multiple sites (human, animal and environment) and their genetic differences have revived interest in the genomic comparison of C. difficile.

Several strains of C. difficile have been isolated and sequenced from different ecological niches; these sequences were procured for genomic comparison. The availability of genome sequence and further development in genomics and related sciences has provided a platform to understand the functions of various proteins encoded by its genome. The concept of pangenome can be applied to identify different genomes: core genome, that is, genes present in all strains of the dataset; dispensable or accessory genome, which are genes present in few strains of the dataset; and strain specific or unique genome, which are genes present in only one strain and absent in others (12). Core gene(s) can be utilized to identify the drug targets and design broad-spectrum antibiotics for pathogenic species, whereas accessory and unique genes are supposed to give them advantage in survival, pathogenicity or habitat adaptation (13). The analysis of gene functions reveals their incorporation in different genes and their proteins, which are functional in various metabolic processes that help pathogens to survive in the different ecological niches (14).

There are various bioinformatics approaches to investigate the drug target from the genome such as ligand-based interaction fingerprint, proteochemometrics modelling, linear interaction energy modelling and many more (15). Here we have used the core genome of C. difficile to identify the drug target by an integrative approach using sequence and interactome analysis. Conventional methods for drug discovery are very costly and time consuming; however, using computer analysis at initial stages can reduce the cost and time. In the present study, pangenome analysis was done and its core genome has been used to identify the drug targets. This method has been developed for the first time to identify the drug targets from the core genome of C. difficile, which can be used in the future for other pathogens too.

Methodology

Interactome analysis

To search the key proteins (from selected proteins), the choke point analysis using pathway tool is performed to find out the proteins which were specific in the metabolic network and whose function cannot be replaced by any other protein (31). A choke point protein is compulsory for any pathway and is essential for pathogen survival. Targeting choke point protein affects the metabolic pathway, which results in death of the pathogen. Further analysis was done using interactome studies, in which the input of all choke point proteins was given to the STRING database and network was created at high confidence level (32). The metabolic functional interaction was created using various methods such as gene fusion, neighbourhood, co-occurrence, co-expression and text mining. The interactome was downloaded from STRING in xml format and its analysis was done using cytoscape (33).

In cytoscape, various critical network parameters such as clustering coefficient, characteristic path length and network centralization were calculated for each node of the network. Clustering coefficient C_n for node n was calculated using:

C_n = 2e_n /k(k_n – 1)

where e_n is the number of connected pairs between all neighbours of n and k_n is the number of neighbours. Characteristic path length is the distance between nodes. Network centralization is the measure of network association around the central node; node having value close to 1 is central to network and value near to 0 shows decentralization. The values for clustering coefficient and characteristic path length are calculated for a node as well as after deleting the node. The difference in these two values shows the impact of node in the network (33). The complete methodology is shown in Figure 1.

Results and discussion

The 97 complete genomes and associated proteomes of C. difficile available till the present study were downloaded from the NCBI. Their information such as accession number, name of the strain, country from where isolated and genome statistics is provided in Supplementary File 1. The pangenome analysis of these 97 strains reveals 6,286 gene families (pangenome), out of which 1,481 are core genes (present in all species). Thus the core genes form 23.5% of the total genome, which signifies a high genetic diversity among different strains. The same feature is represented in core-pan genome plot (Fig. 2). As the genomes are added, the size of pangenome increases, whereas the size of core genome declines. The curve of pangenome (yellow colour) is still progressing, indicating the likelihood of addition of more genes, that is, global gene repertoire is likely to change in the near future and its pangenome is almost open.

The power law regression model equation,

f(x) = a.X^b

where f(x) is the pangenome size, X is the number of genomes used, and a and b are fitting parameters used to find the openness and closeness of pangenome, has been used in the present study (17). In our study with C. difficile genes, the values of f(x) = 6,286, X = 97, a = 2687.17 and b = 0.185459 indicated that the pangenome is open at b > 0; otherwise, the pangenome would be considered to be closed. The value of b = 0.185459 indicates that the pangenome is open but soon may be closed with increase in genome data. C. difficile has 186,308 accessory genes which are present in a few strains and has 976 unique genes which are present in specific strains of C. difficile. The conservation level of C. difficile does not seem to be very high. The genome level analysis reveals high genetic variability; this may be due to its existence in different niches.

On functional analysis, it is observed that unique genes (shown in blue) and accessory genes (shown in red) are mostly responsible for metabolism and transporter function (Fig. 3). This shows that both of them are more diversified in C. difficile. Recently in 2021, Kulecka et al. also reported the variability in the metabolism genes in recurrent CDI cases (34). In addition, the core genes like PolC-type deoxyribonucleic acid (DNA) polymerase III, exonuclease subunit C, cell wall–binding protein Cwp20, sensor histidine kinase KdpD and alanine-tRNA ligase (shown in green) are the genes for cellular processes such as cell division, cell cycle and its control, cell motility, cell wall/membrane biogenesis, transcription, translation, ribosomal structure and biogenesis, and are mostly conserved. This signifies that genes involved in important cellular mechanism are conserved, while genes required to adapt to the new ecological niches are variable.

The pan-phylogeny-based phylogenetic tree is shown in Figure 4. It is observed that strains are divided into different clades based on genome similarity. On observing each clade, it is noticed that similarity is mainly based on the country from where the samples are isolated. It suggests that C. difficile adapts to different environmental conditions by expression of relative proportions of the different gene products. For example, strains Cd9, Cd12, MT5121, Cd23 and W0023a are all assembled in one clade and are isolated from the USA. Similarly strains CD-10-00484, 10-00078, DSM 102860, DSM 102978 and DSM 29745 isolated from Germany show genome similarity and are assembled in one clade.

With high genetic variability and drug resistance for CDI, it is very necessary to design a drug that targets the core genes of the pathogen, as core genes are present in all the strains of the pathogen and are essential for the survival of the pathogen. Therefore, targeting core gene will surely help to overcome CDI. We have used an integrative approach based on sequence and interactome analysis to find the drug target against C. difficile.

From the core genome, genes that are homologous to humans are excluded in the first step, as it may adversely affect the host metabolism. A total of 1,130 proteins are found to be non-homologous to human (Supplementary File 2). On further screening, essential proteins that are vital for the survival of pathogens are searched using DEG. Among them, 370 proteins were found to be essential and crucial for C. difficile survival (Supplementary File 2). Essential proteins were further screened for their virulence, as these factors are responsible for pathogenesis. From 370 essential proteins, 130 proteins were found to be virulence-associated factors (Supplementary File 2).

All 130 proteins were checked for their physicochemical properties. Proteins having more than 100 amino acids, less molecular weight, low pI, negative GRAVY value, high aliphatic index and membrane or cytoplasmic localization were further considered (26-30). All these are the physicochemical properties required for the potent drug target. A total of 94 proteins were obtained after all physicochemical checks (Supplementary File 2); they are further used for choke point analysis.

On choke point analysis, only 39 proteins involved in the unique metabolic pathways were identified (Supplementary File 2). For these 39 proteins, interactome is created using STRING as shown in Figure 5.

On interactome analysis with cytoscape, eight potent drug targets were found. Their interactome analysis results are shown in Table I. Out of the eight drug targets, three targets, namely, aminodeoxychorismate synthase (PabB), D-alanyl-D-alanine carboxypeptidase (DD-CPase) and undecaprenyl diphospho-muramoyl pentapeptide beta-N-acetylglucosaminyl transferase (MurG transferase) were identified which have been previously reported (35).

PabB is involved in folate synthesis; its inhibition affects DNA and protein synthesis adversely. It is reported that it is targeted by the antibiotics 6-fluoroshikimic acid and atrop-abyssomycin C (35). Bifunctional diguanylate cyclase/phosphodiesterase (cyclic-diGMP) is a messenger protein that regulates motility, virulence and biofilm formation attributed to pathogenicity (36,37). Inhibition of cyclic-diGMP affects many processes of the pathogen that result in the death of the pathogen. Another identified target is sporulation transcription factor Spo0A, which is a key factor for entry into sporulation in stress conditions and biofilm formation (38).

Histidinol-phosphate transaminase (HisC) is a transferase that is mainly involved in transferring nitrogenous group. It is involved in synthesis and metabolism of many amino acids (39). 3-deoxy-7 phosphoheptulonate synthase (DAHP synthase) is involved in the shikimate pathway and is also responsible for the synthesis of aromatic amino acids, such as tyrosine, phenylalanine and tryptophan which are essential for bacterial metabolism (40). Another identified drug target, MurG transferase, is involved in the peptidoglycan biosynthesis and is reported to be the drug target for many pathogens such as Neisseria meningitidis (41), Acinetobacter baumannii (42) and Mycobacterium tuberculosis (43).

Another identified potent target, DD-CPase is a reported drug target and is involved in peptidoglycan biosynthesis and remodelling and is inhibited by β-lactam antibiotics (44). C-di-GMP phosphodiesterase (PdcA) is a potent drug target that regulates bacterial pathogenesis as well as is involved in surface adherence and biofilm development (45).

Conclusion

CDI is a challenging situation worldwide. The unavailability of specific antibiotic and emergence of antibiotic resistance against C. difficile is a matter of concern. Diversity in C. difficile genome drives the genomic comparison of the pathogen. Pangenome analysis reveals the open pangenome of C. difficile that may be soon closed. The diversity is due to its adaptability in different niches and different hosts. Due to such genomic diversity, a drug target can be designed only from its core gene, whose inhibition affects all strains of C. difficile. From sequence and interactome analysis of core genes, eight potent drug targets are reported. Out of these, three of the targets – PabB, DD-CPase and MurG transferase – are also reported as drug target for other pathogens, whereas bifunctional cyclic-diGMP, Spo0A, HisC, DAHP synthase and PdcA are newly reported targets. This indicates that the method originated in the present study has a high rate of success and saves considerable time and money. The same method can be used for other pathogens also.

Limitations of the study

This computational method uses multiple genomes with complete genome assembly, therefore this method cannot be employed for pathogens whose majority of the strains have not been sequenced. Pangenome analysis can only be performed if multiple sequenced strains are available for the organism. Pangenomics can easily be done only when proper software or tool is available; it is very complex to handle the genome manually. Another limitation with this study is the use of in silico method for drug target identification that saves considerable time and money, but its accuracy is still questionable. We have used the core genes which are conserved and involved in essential processes for drug target identification, but the drug targets from novel genes which the bacteria inherits from its adaptation in new niches cannot be considered in this type of study.

Acknowledgements

The authors acknowledge the facilities of the Department of Biotechnology, Ministry of Science and Technology, Government of India, New Delhi (DBT), under the Bioinformatics Sub Centre as well as M.Sc. Biotechnology programme used in the present work.

Disclosures

Conflict of interest: The authors confirm that they have no conflict of interest.

Financial support: This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.