Evidence of episodic positive selection in Corynebacterium diphtheriae complex of species and its implementations in identification of drug and vaccine targets

Samples and taxonomy

Episodic positive selection across species or other phylogenetic groups of the diphtheriae group was investigated using a branch-site test (Zhang, 2005). This test is more appropriate for inter-specific samples, because it assumes that the observed mutations have already been fixed by selection (Kryazhimskiy & Plotkin, 2008; Anisimova & Liberles, 2012; Kosiol & Anisimova, 2012), and one genome could represent a species. For this reason, we limited the samples to one per species, subspecies or biovar. For the foregrounds (target groups), we selected the six type strains from the C. diphtheria complex and other strains to represent biovars and lineages. C. diphtheriae biovars could not be used as foregrounds as they are united in a single clade (Sangal et al., 2014). As background, we selected 40 total genomes that included 30 representative (O’Leary et al., 2016), eight complete and three WGS genomes of Corynebacterium species, all of which were available in the Pathosystems Resource Integration Center (PATRIC) (Davis et al., 2020). These were annotated by RASTtk (Brettin et al., 2015) and downloaded from PATRIC (Table S1).

The taxonomy of the samples was verified using TYGS (Meier-Kolthoff & Göker, 2019). TYGS determines the closest related type strains using the MASH algorithm (Ondov et al., 2016) for entire genomes and BLASTn for 16S sequences. It calculates the pairwise distances of 16S and genome sequences using GBDP (Meier-Kolthoff et al., 2013), followed by inference of 16S and genome phylogenies based on the pairwise GBDP distances using FastME (Lefort, Desper & Gascuel, 2015), digital DNA-DNA hybridization (dDDH) using GGDC (Meier-Kolthoff et al., 2013), and deviation of G+C content. Genomes with >70% of dDDH and <1% of G+C content deviation are considered to be in the same species (Meier-Kolthoff & Göker, 2019).

Positive selection analysis

A genome-scale positive selection analysis was performed using the PosiGene pipeline (Sahm et al., 2017) on the Corynebacterium genomes. Ten foreground genomes were tested (Table S2) representing 10 target clades or subclades. Six clades represent the species from the C. diphtheria complex (C. belfantii, C. diphtheriae, C. pseudotuberculosis, C. rouxii, C. silvaticum and C. ulcerans), two subclades represent C. pseudotuberculosis (biovars equi and ovis), and two subclades representing C. ulcerans lineages (lineage 1 and 2).

In the module “create_catalog”, homologous genes were identified using BLASTp (Camacho et al., 2009) with the best-bidirectional hit criterion (Altenhoff & Dessimoz, 2009). In the module “alignments”, orthologous genes are identified, gene trees are built, and a species tree is built from the gene trees. In this module, the anchor species is the genome that the orthologous gene sequences are aligned to, and the reference species is the genome from which the gene names are extracted. C. diphtheriae NCTC 11397^T was selected as both the reference and anchor genome. In the first step, orthologous gene sequences were aligned to the anchor genome sequences using CLUSTALW (Larkin et al., 2007), with the parameters’ minimum identity and minimum pairs identity set to 40%. The aligned sequences were filtered by GBLOCKS for gaps and unreliable alignment columns (Jordan & Goldman, 2012). In the next step, a phylogeny was built for each gene using the parsimony method and jackknifing implemented in DNAPARS from the PHYLIP package (Felsenstein, 2005). In the third step, a species tree was built based on the gene trees consensus, using CONSENSE from the PHYLIP package. The species tree is required to test for positive selection along specific lineages (Yang & Nielsen, 2002). The consensus tree was visualized using FigTree v1.4.4 and rooted using C. kroppenstedtii DSM 44385 as an outgroup. This strain was identified as an outgroup based on another tree generated by the same pipeline including M. tuberculosis H37RV to find the correct Corynebacterium species to use as an outgroup (Table S1) in that tree. The M. tuberculosis tree was not included in the downstream analysis.

In the module “positive_selection”, a likelihood ratio test compares the non-synonymous to synonymous substitution rate ω = d_N/d_S in the foreground and the background. Here, d_N is the number of non-synonymous substitutions per non-synonymous site and d_S is the number of synonymous substitutions per synonymous site. The episodic positive selection model assumes ω > 1 in the foreground and ω = 0 or ω < 1 in the background, while the null model assumes ω = 0 or ω < 1 for foreground and background (Yang & Dos Reis, 2011). We considered positive selection if p < 0.05 for False Discovery Rate, as this correction is more suitable for genome wide experiments (Storey & Tibshirani, 2003).

Due to an assumption of no recombination by the branch-site models we used (Yang, 2005), recombination could cause false positive results (Anisimova, Nielsen & Yang, 2003). To avoid that artifact, the genes identified as positively selected by PosiGene were then tested for intragenic recombination using PhiPack (Bruen, Philippe & Bryant, 2006) that calculates Pairwise Homoplasy Index method (PHI) (Bruen, Philippe & Bryant, 2006), Neighbor Similarity Score (NSS) (Jakobsen & Easteal, 1996), and Maximum Chi-Square (Smith, 1992). In our analysis, we considered that recombination had occurred when q < 0.05 for PHI and at least NSS or Maximum Chi-Square (Hongo et al., 2015). Genes identified as recombinant were discarded for downstream analysis.

To minimize the false positive results caused by misalignments, frameshifts and ortholog prediction (Schneider et al., 2009; Markova-Raina & Petrov, 2011), we visually checked the alignments and checked the homology prediction by comparing the orthologs protein domains using PATRIC’s annotation of Local and Global Families (Davis et al., 2016) and the Conserved Domain Database (CDD) (Lu et al., 2020).

Gene annotation

Gene function was predicted using annotations from PATRIC and eggNOG-mapper (Huerta-Cepas et al., 2017), and InterProScan (Jones et al., 2014) was used to examine protein domains. Subcellular localization of the proteins was assessed using SufG+ v1.2.1 (Barinov et al., 2009). Protter (Omasits et al., 2014) was used to identify the position of the positively selected amino acids in relation to the cytoplasmic, transmembrane, and surface exposed portions of the proteins.

GIPSy v1.1.2 (Soares et al., 2016) was used to identify genomic islands (GIs) in C. diphtheriae NCTC 11397^T using C. glutamicum ATCC 1302 (NC_006958.1) as the non-pathogenic reference. The islands were predicted by genes with C+G content and codon usage deviation, presence of transposases, presence of specific genes, non-conservation in comparison to reference genome, and flanking tRNA genes. The islands were classified according to the proportion of specific genes as pathogenicity islands, metabolic islands, resistance islands and symbiotic islands. Finally, prophages were predicted using PHASTER (Arndt et al., 2016).

Prediction of drug and vaccine targets

Virulence genes and drug targets were predicted using the Pipeline Builder for Identification of Targets (PBIT) (Shende et al., 2016). PBIT predicted drug targets among the cytoplasmic proteins by a subtractive approach, identifying sequences of interest using BLASTp and specific databases in the following way. First, homologs to the human proteome, anti-targets and gut microbiota proteomes were filtering out to avoid cross-reactivity of drugs. Then the essential genes were identified using the Database of Essential Genes (Zhang, 2004) and virulence genes using the Virulence Factor Database (Chen et al., 2005). The druggability of the remaining candidates was predicted by similarity to experimentally validated druggable targets from the Therapeutic Target Database (Li et al., 2018).

Vaccine targets were predicted from the transmembrane, putative surface exposed and secreted proteins predicted using SurfG+ and Vaxign (He, Xiang & Mobley, 2010) with prediction of human MHC Class I and II epitopes.

C. diphtheriae

In C. diphtheriae (foreground Cd), we identified two genes encoding proteins predicted to be essential and involved in translation, amino acid transport and metabolism (Table 1 and Table S15). The gene ansB is in GI10 and encodes secreted L-asparaginase type II which is a high-affinity enzyme that catalyzes the conversion of L-asparagine to L-aspartate and ammonia. The E. coli, Dickeya dadantii and human homologs to this gene are used for leukemia treatment, where the consequent low L-asparagine levels in plasma leads to apoptosis of the leukemia cells (Lubkowski & Wlodawer, 2021). This gene was suggested as a candidate vaccine target due to its classification as a secreted protein and predicted epitopes (Table 2). The second protein, SSU ribosomal protein S3p (rpsC), is a 30S ribosomal subunit that binds to the initiator Met-tRNA (Burd & Dreyfuss, 1994). Possible reasons for the selective pressure on these genes could be the effects on L-aspartate uptake (ansB) and translation efficiency (rpsC). rpsC was suggested as a drug target but has homology to a protein in the human gut microbiota proteome, requiring compounds that can selectively inhibit it.

C. pseudotuberculosis

When C. pseudotuberculosis was the foreground (Cp), eight positively selected genes were identified. Five of the genes were tagged as essential, and three as virulence factors. They are involved in translation, coenzyme transport and metabolism, inorganic ion transport, defense from foreign DNA, nucleotide transport and metabolism and adhesion (Table 1 and Table S15). Among the essential genes, an ABC transporter permease protein (mntC) plays a role in the transport of Mn²⁺ and Zn²⁺ (Claverys, 2001) and was classified as a virulence factor. The Dihydropteroate synthase 2 (folP2) is nonfunctional according to PATRIC annotation but functional genes with this annotation are essential for the de novo synthesis of folate (Bertacine Dias et al., 2018). The HNH endonucleases degrade foreign DNA, and can also be involved in DNA repair, replication and recombination (Wu, Lin & Yuan, 2020). Peptide chain release factor 1 (prfA) recognizes the stop codons UAA and UAG, promoting the end of translation (Scolnick et al., 1968). The Putative phosphoglycerate mutase (pgmB) is capable of interconverting 2- and 3-phosphoglycerate in glycolysis (Rigden, 2008), although this particular gene is annotated as putative.

The other three identified genes (add, spaE and the putative oxidoreductase) are not characterized as essential, so may not be suitable drug targets. Adenosine deaminase (add) catalyzes the hydrolytic deamination of adenosine into inosine (Chang et al., 1991). spaE, from the operon spaDEF, encodes the minor pilin SpaE in C. diphtheriae (Mandlik et al., 2008) and is in GI5. It’s ortholog in C. pseudotuberculosis also encodes the minor pilin and was first described as spaB from the operon spaABC (Trost et al., 2010). The putative oxidoreductase is a flavoenzyme with a “FAD-binding domain, ferredoxin reductase-type” (IPR017927), but its specific reaction is unknown.

Why would these particular genes be under positive selection? One could hypothesize that there would be more efficient manganese uptake (mntC), tissue adhesion on a new host range (spaE), improved efficiency for defense against foreign DNA (HNH endonuclease), translation (prfA), and metabolism of nucleotides (add) and carbohydrates (pgmB).

The vaccine targets (mntC and spaE) were indicated due to either their membrane location, their predicted epitopes, and that they might have surface exposed sites that are under positive selection. There were four drug targets (Table 2). pgmB was predicted as a target by PBIT and is this same gene is a drug target in helminth parasites (Timson, 2016). folP2 was also predicted and is a well know target of sulfa and imidazole derivatives in human pathogens such as Staphylococcus aureus, M. tuberculosis, Bacillus anthracis, Streptococcus pneumoniae, Burkholderia cenocepacia and Yersinia pestis (Bertacine Dias et al., 2018). Although it was annotated as non-functional, the evidence of positive selection in this protein suggests an unknown adaptation due to specific amino acids that could be targeted. prfA has homology to human and human gut microbiota proteome and has no predicted druggability, but it is a known target of the drug Apidaecin in gram negative bacteria (Matsumoto et al., 2017). The HNH endonuclease had no predicted druggability.

C. pseudotuberculosis biovar equi

When C. pseudotuberculosis biovar equi was used as the foreground (Cpequi), two genes were identified as being under positive selection, and were also characterized as essential. These two genes (mapB and tyrS) are both involved in translation (Table 1 and Table S15). Methionine aminopeptidase (mapB) cleaves the initiator methionine from newly synthesized polypeptides (Helgren et al., 2016; Pillalamarri et al., 2021). Tyrosyl-tRNA synthetase (tyrS) attaches the amino acid tyrosine to the appropriate tRNA (Hughes et al., 2020; Othman et al., 2021). These two genes could be under positive selection as it could affect translation efficiency. Both genes were predicted as homologs to human gut microbiota and have been identified as drug targets in other studies (Table 2). tyrS was predicted as a virulence factor and an ortholog from Pseudomonas aeruginosa was found to be targeted by four drug-like compounds (Hughes et al., 2020), and a S. aureus ortholog to this gene was targeted by new pyrazolone and dipyrazolotriazine derivatives (Othman et al., 2021). mapB from M. tuberculosis and S. pneumoniae were shown to be selectively targeted despite homology to human protein (Krátký et al., 2012; Arya et al., 2015).

C. rouxii

A single hypothetical protein (ERS451417_00470) was identified when C. rouxii was used as the foreground (Cr). It had no predicted domains, but it could be a vaccine target candidate (Table 2) due to its transmembrane location and the surface exposed sites under positive selection.

C. ulcerans

A single essential gene (glyA) was identified when C. ulcerans was used as the foreground (Cul). glyA encodes a serine hydroxymethyltransferase enzyme (Table 1 and Table S15). This gene is known to participate in the one-carbon metabolism of serine/glycine interconversion and also in the folate/methionine cycle (Batool et al., 2020), which could explain its being under selective pressure. This gene was also identified as a potential drug target (Table 2), but it does have homology to human and human gut microbiota proteome. It has been shown to play a key role in lysostaphin resistance in Staphylococcus aureus (Batool et al., 2020).

C. ulcerans lineage 1

Three genes were identified when the C. ulcerans lineage 1 genome was used as the foreground (Cul1). Two of them were essential, involved in carbohydrate transport and metabolism (Table 1 and Table S15). The Phosphoenolpyruvate-dihydroxyacetone phosphotransferase, dihydroxyacetone binding subunit DhaK (dhaK) gene is in GI34. This enzyme phosphorylates ketones and short chain aldoses using adenosine triphosphate (ATP) (Peiro et al., 2019). The second gene is annotated in PATRIC as “Similar to citrate lyase beta chain, 3” (citE) and is probably one of the catalytic subunits of citrate lyase, the enzyme that catalyzes the cleavage of citrate to acetate and oxaloacetate during citrate fermentation (Schneider, Dimroth & Bott, 2000). Both proteins were predicted as virulence factors by PBIT. The third gene encodes a hypothetical protein (ERS451417_00635) with no predicted domain or cellular localization.

The two genes with predicted function (dhaK and citE) appear to be related to metabolism inside the host. In Listeria monocytogenes, Phosphoenolpyruvate-dihydroxyacetone phosphotransferase (DhaK and other subunits) is required to utilize carbon sources for amino acid synthesis inside murine macrophages (Eylert et al., 2008). In Enterococcus faecalis, mutants of citrate fermentation genes (citE and others) were less pathogenic for the model Galleria mellonella (Martino et al., 2018). These same two genes are possible drug targets (Table 2). citE was predicted as a virulence factor and drug target candidate, while dhaK was predicted as a virulence factor and suggested despite the homology to a protein in the gut microbiota homology.

C. ulcerans lineage 2

Two essential genes (Table 1 and Table S15) were identified when the C. ulcerans lineage 2 was used as the foreground (Cul2). The DNA polymerase III epsilon subunit (dnaQ) has a domain with 3′–5′ exonuclease proofreading activity (Raia, Delarue & Sauguet, 2019). The other gene, cobK, encodes Precorrin-6A reductase which is involved in part I of the cobalamin cofactor (vitamin B12) biosynthesis pathway (Kipkorir et al., 2021). The mutations seen in these genes could provide the organism with a more efficient means of DNA replication (dnaQ) and biosynthesis of the essential cofactor cobalamin (cobK). Neither of these genes had any predictable druggability.

Probable adaptations across groups

It is reported that most of the genes identified as being under positive selection are exposed on the surface and are involved in host colonization, and resistance to phage and antibiotics (Petersen et al., 2007; Anisimova & Liberles, 2012). Those under positive selection that are not surface exposed have been shown to be involved in metabolism (Petersen et al., 2007; Rao, Sivakumar & Jayakumar, 2019) or gene regulation (Zhang et al., 2011). Considering the function of the identified genes, most of the probable adaptations appear to be related to metabolism. A notable exception is in C. pseudotuberculosis, where the pilin SpaE could have enhanced adhesion to different host species tissues. Although the specific adaptations are not clear, an amino acid fixed by positive selection is an attractive target for a therapeutic molecule, as a non-synonymous mutation that could avoid interaction would decrease fitness. These genes could be used for reverse vaccinology and in silico drug targeting methods.