Identifying potential drug targets in the kinomes of two monogenean species

Kinomes were identified in the predicted proteins from the transcriptomes of R. viridisi and S. longicornis (Caña-Bozada et al., 2022) by using an integrated bioinformatics pipeline, similar to Giuliani et al. (2018), with modifications (Fig. 1A). To that end, a preliminary group of kinases was identified and classified using Kinannote 1.0 software (Goldberg et al., 2013), with the parameter “-m” of metazoan. Incomplete sequences detected in the proteins of S. longicornis were processed with the CAP3 program (Huang & Madan, 1991) to obtain longer sequences. Then, another preliminary group of kinases was identified and classified by retrieving R. viridisi and S. longicornis kinases with orthology to Schistosoma japonicum (Giuliani et al., 2018), Schistosoma mansoni (Andrade et al., 2011), Taenia asiatica, and Taenia saginata (Wang et al., 2016) using Orthologous MAtrix (OMA) Standalone 2.5.0 software (Altenhoff et al., 2019). To improve the annotation, the kinase sequences identified from Kinannote and OMA were merged, and their classifications were compared to detect possible misannotations. These sequences were aligned against the NCBI nonredundant protein database version 5 (Sayers et al., 2022) using BLASTp 2.12.0 (e-value < 0.05) (Camacho et al., 2009) to eliminate non-kinase sequences and possible contamination of bacteria, viruses, fungi, and fish. To avoid overrepresentation of genes, the longest isoform per gene was extracted using the Trinity script “get_longest_isoform_seq_per_trinity_gene.pl”.

Fig. 1.

In addition, the proteins were classified by aligning the kinase sequences against the database of Kinbase (Manning et al., 2022) using BLASTp. Sequences with an e-value > 10e⁻⁵ were removed. The sequences obtained from Kinannote, OMA, and BLASTp were used to classify the putative kinases by group, family, and subfamily, according to KinBase, which is a database containing information on human kinases and their homologs in other eukaryotes (Fig. 1B). Then, this classification was confirmed by phylogenetic analysis, as described below.

Phylogenetic analysis was performed for each major group of kinases using their conserved catalytic domain (PROSITE domain type PS50011), which was predicted using PS-SCAN 2.0 (de Castro et al., 2006). The catalytic domain sequences were aligned using MAFFT 7.31 (Katoh & Standley, 2013) with the parameters “--localpair” and “maxiterate 1000”, and then the gaps were removed with Trimal (Capella-Gutiérrez et al., 2009) using the parameter “-gappyout”. The alignment was refined in MUSCLE 3.8.31 (Edgar, 2004). The analysis included reference kinases from Homo sapiens, Drosophila melanogaster, Saccharomyces cerevisiae, and Caenorhabditis elegans obtained from the KinBase database.

The tree was constructed with IQ-TREE 1.6.12 (Nguyen et al., 2015), using the Shimodaira–Hasegawa-like approximate likelihood ratio test (SH-aLRT; 1000 replicates). The best evolutionary models were selected according to the Bayesian Information Criterion (BIC) with the ModelFinder program (Kalyaanamoorthy et al., 2017), which is integrated into IQ-TREE. The trees were visualized and annotated in FigTree 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/). Nodes with support values ≥ 80 % were considered functional predictions with support.

Putative kinases were functionally annotated using Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000). The PANNZER2 server (Törönen et al., 2018) was used to retrieve GO terms. The KEGG Automatic Annotation Server (KAAS 2.0) (Moriya et al., 2007) and BlastKOALA 3.0 (Kanehisa et al., 2016) were used to retrieve KEGG Orthology (KO) terms and then map to the KEGG terms using the KO database as reference (https://www.genome.jp/keggbin/get_htext?ko00001.keg; accessed May 6, 2021).

Monogenean drug targets (MDTs) were selected from predicted kinomes according to homology (Fig. 1C). To that end, the R. viridisi and S. longicornis kinases with orthology to the Caenorhabditis elegans lethal RNAi phenotype, which represent potential targets for parasite control (Gahoi et al., 2019), were retrieved using BLASTp (e-value < 10e⁻⁵). Then, the retrieved monogenean sequences were aligned against the ChEMBL database (Gaulton et al., 2017) to identify potential kinase drug targets using BLASTp (e-value < 10e⁻³⁰). Finally, the retained sequences were submitted to the Orthovenn2 web platform (Xu et al., 2019) to identify kinases shared between the monogeneans, the fish hosts C. viridis (NCBI SRA: SRP165941) and O. niloticus (NCBI Assembly: GCA_001858045.2), and humans. Potential target kinases of monogeneans were considered MDTs if they had orthologs in neither fish nor humans but did in the lethal phenotype of C. elegans. Subsequent analyses included only MDTs that we deemed priorities, which were selected as follows. Particular attention was paid to KEGG pathways involving only one or two kinases. If those kinases corresponded to an MDT by manual inspection, they were assigned priority (Fig. 1D). The rationale here was based on the assumption that an essential biological process is easier to inhibit if it depends on only a few proteins instead of several (Stroehlein et al., 2018). In addition, priority MDTs were also selected through a network analysis of protein–protein interactions (PPIs) performed for each monogenean species. To that end, each kinome was submitted to the STRING server (Snel et al., 2000) using the PPI database for S. mansoni and the parameters “high confidence cutoff of 0.7” and “zero node addition”. The interaction network was exported, visualized, and analyzed within the Cytoscape 3.9.1 platform (Shannon et al., 2003). The cytoNCA 2.1.6 plugin (Tang et al., 2015) in Cytoscape was used to calculate the betweenness centrality, which is a value representing the participation of nodes in multiple pathways and that presents a positive correlation with its degree (Rai et al., 2021). MDTs with better centrality indices were assigned priority.

Priority MDTs were mapped against the Protein Data Bank (PDB) database (www.rcsb.org) (Burley et al., 2017), using BLASTp (e-value < 10e⁻⁵), to identify those with known 3D structures. The molecular structures of MDTs were generated using the SWISS-MODEL server (Schwede et al., 2003) and downloaded in PDB file format. For each of the five target kinases selected, only the best model predicted by SWISS-MODEL was used in posterior analyses. The best predicted structures were selected based on the sequence identity and the evidence of the template used (Crystal structure of apo murine Nf-kappaB inducing kinase (NIK), Crystal structures of human p70S6K1-T389A (form II), Crystal structures of human p70S6K1-T389E Discovery of a Novel, Potent and Selective Inhibitor of 3-Phosphoinositide Dependent Kinase (PDK1), and Structure of the N-terminal domain of SAD). The predicted structure for each candidate MDT was used to perform molecular docking against all US Food and Drug Administration (FDA)-approved drugs from the ZINC20 database (Sterling & Irwin, 2015). The drug database included the 3D structure in Mol2 file format of approximately 2000 drugs. The prepare_ligand4.py script of AutoDock tools (Morris et al., 2009) was used to translate each PDB file, of both MDT and drugs, into Protein Data Bank Partial Charge & Atom Type (PDBQT) format. The 3D structures of the modeled proteins along with model information are available in Supplementary File S1. For each of the combinations (receptor-ligand) of the five receptors and 2000 ligands, we explore the whole receptor space by generating a custom script to perform the molecular docking with AutoDock Vina 1.1.2 (Trott & Olson, 2010) in a 3-dimensional grid of 50x50x50 Angstroms and exhaustiveness of 16 (as recommended in the AutoDock Vina Manual, https://vina.scripps.edu/manual/).

Putative kinases of R. viridisi and S. longicornis were identified and classified using an integrative bioinformatic pipeline. The curated kinome of R. viridisi is composed of 159 ePKs and 1 aPK (Table 1), distributed across 64 families and 69 subfamilies. The curated kinome of S. longicornis comprises 191 ePKs and 2 aPKs (Table 1), distributed into 74 families and 80 subfamilies. More information about the results of each step in the bioinformatic pipeline is given in Supplementary Table S1 and the sequences in Fasta format are available in Supplementary File S2.

The kinase group having the most members was CMGC, followed by AGC, whereas the groups least represented were RGC, PKL, and aPKs (Table 1). The most represented families were CDK and CAMKL (Supplementary Table S1).

The classification of kinases described above was well supported (bootstrap ≥ 80 %) by the phylogenetic analysis, except for 13 kinases of R. viridisi and 18 of S. longicornis (Supplementary Figs. S1 – S9). The 31 kinases that had low bootstrap support were grouped with kinases from subfamilies or families that differed from the classification obtained using Kinannote, OMA, and BLASTp; thus, those kinases were only classified to the level supported by phylogenetic analysis (major group or family; Supplementary Table S1).

Of the 160 kinases in R. viridisi, 106 were mapped to 1109 GO terms, whereas of the 193 S. longicornis kinases, 122 were mapped to 1125 GO terms (Supplementary Table S2). The GO term distribution was similar in both species. For the two species of monogeneans, the most represented GO terms were ATP binding (GO:0005524), protein kinase activity (GO:0004672), and protein serine/threonine kinase activity (GO:0004674) within of the Molecular Function category; protein phosphorylation (GO:0006468), peptidyl-serine modification (GO:0018209), and intracellular signal transduction (GO:0035556) within the Biological Process category; and integral component of membrane (GO:0016020), nucleus (GO:0005634), and cytoplasm (GO:0005737) within the Cellular Component category (Supplementary Fig. S10A).

Kinases from the two monogenean species mapped to 224 KEGG pathways (Supplementary Table S3). According to the KEGG database, the kinases mainly belonged to Brite Hierarchies (142 and 164 members in R. viridisi and S. longicornis, respectively) and the Environmental Information Processing category (57 and 75 members, respectively) (Supplementary Fig. S10B and Supplementary Table S4). In R. viridisi, the predominant terms included protein kinases [BR:ko01001], membrane trafficking [BR:ko04131], chromosome, and associated proteins [BR:ko03036], pathways in cancer [PATH:ko05200], PI3K-Akt signaling pathway [PATH:ko04151]. Meanwhile, S. longicornis exhibited a similar profile with protein kinases [BR:ko01001], MAPK signaling pathway [PATH:ko04010], pathways in cancer [PATH:ko05200], and membrane trafficking [BR:ko04131] as the most abundant terms.

Using the bioinformatics pipeline, 152 kinases of R. viridisi were identified as orthologous with the C. elegans lethal RNAi phenotype. Of these, 128 mapped to the ChEMBL database of target associations; 10 had equivalents in S. longicornis but were absent in the fish hosts C. viridis and O. niloticus. In addition, we found one kinase specific to R. viridisi. In S. longicornis, 182 kinases mapped to the C. elegans lethal RNAi phenotype. Of these, 157 mapped to the ChEMBL database and 11 had equivalents in R. viridisi but not with C. viridis and O. niloticus (Supplementary Table S5; one kinase in R. viridisi, TRINITY_DN4663_c0_g2_i1.p1, had two orthologous paralogues in S. longicornis).

These 11 kinases of R. viridisi and the 11 of S. longicornis, which we assigned as MDTs (Table 2 and Supplementary Table S5), belong to the major groups CK1, AGC, TK, CAMK, STE, and CMGC, and are listed in the KEGG BRITE categories metabolism (09181), signal transduction (09132), genetic information processing (09182), cellular community—eukaryotes (09144), endocrine system (09152), and immune system (09151). According to the analysis carried out with Orthovenn2, all the kinases classified as MDTs turned out to be homologous to proteins with known 3D structures. Of the 22 MDTs mentioned above, two were unique kinases involved in two separate KEGG pathways (Supplementary Table S5), and three had high values of betweenness centrality (betweenness ≥ 168; Supplementary Table S6). Therefore, these five MDTs were considered priorities for the molecular docking analysis (Table 2 and Supplementary Fig. S11).

The molecular docking was performed using FDA-approved drugs from the ZINC database and five priority MDTs. The MDTs used for the molecular docking analysis are shown in Table 2. Dihydroergotamine was found to have the highest affinity (ΔG = −9.7 kcal/mol) interactions with kinases of both S. longicornis and R. viridisi (Table 3). The scores obtained from the molecular docking are shown in Supplementary Table S7.

The bioinformatic pipeline, the inspection of each protein, and the phylogenetic analysis let us identify and classify the kinomes of S. longicornis and R. viridisi with a high degree of confidence. The number of kinases in R. viridisi and S. longicornis was lower compared to other metazoan organisms (Stroehlein et al., 2018), including parasitic helminths such as cestodes, trematodes, and other flatworms. Given that our identification relied on transcriptomic data from adult organisms, unlike the genomic data used in previous studies, the observed decrease in the number of kinases in monogeneans may be attributed to the discrepancy in data sources. All the major kinase groups were represented by fewer members in monogeneans than in other platyhelminths, except for CK1, which was represented similarly in monogeneans, trematodes, cestodes, and other metazoans (Stroehlein et al., 2018).

In both monogenean species studied here, members of CMGC and AGC were the most abundant kinases. Similarly, CMGC is the principal group in most cestodes and trematodes, followed by AGC in the case of cestodes. In the trematode S. mansoni, the abundance of CMGC kinases has been linked to essential processes such as control cell proliferation, organelle replication, and segregation (Andrade et al., 2011), whereas in the cestode T. solium, these have been associated with hermaphroditism (Arora et al., 2020). This dominance of CMGC and AGC among platyhelminth species is distinct from other metazoa (Stroehlein et al., 2018). Therefore, higher representation of these major groups of kinases might tie them to crucial roles in platyhelminths.

The workflow facilitated the identification of MDTs. Although the transcriptomes of R. viridisi and S. longicornis were obtained from adult organisms, it is essential to acknowledge the lack stage-specific transcriptomes, which prevents definitive conclusions about the expression of the predicted MDT kinases at all stages of development. However, it is notable that most kinases play fundamental roles in crucial signaling processes throughout the parasite life cycle, such as in S. haematobium, where 214 of 258 transcribed kinases show expression in all sexes and life stages studied (Stroehlein et al., 2015), underscoring the fundamental functions these kinases likely play throughout various stages of parasite development.

CAMK, AGC, and TK were the most represented major groups in the 22 MDTs. This is consistent with the results of Giuliani et al. (2018), whose bioinformatics pipeline also found CAMK and AGC kinases as the most represented drug targets. Coincidentally, members of the CAMKL family were the most abundant in both studies—five in ours, four in the aforementioned study. The kinase GSK (group CMGC) of R. viridisi and S. longicornis was proposed here as an MDT, which is consistent with the result of a comparative chemogenomics study in S. mansoni (Caffrey et al., 2009). Members of GSK are involved with cell proliferation, regulation of glycogen synthase, and participation in the Wnt signaling pathway, which is an integral component of embryonic development (Grimes & Jope, 2001). Inhibition of GSK3 stops S. mansoni laying eggs (Morel et al., 2014). Therefore, the GSK family members represent potential drug targets in monogeneans.

To find candidate drugs to be repurposed against these parasites, we performed molecular docking analysis with 2000 FDA-approved drugs and the five priority MDTs. The drug–kinase pairs with the highest affinities were dihydroergotamine–BRSK and dihydroergotamine–MEKK1. Dihydroergotamine is used to treat migraine through vasoconstriction in the cranial vascular bed by its agonist activity at 5-HT, α-adrenergic, and dopaminergic receptors (Ramírez-Rosas et al., 2013). It actively inhibits the growth of the parasite Plasmodium falciparum (Weisman et al., 2006). Although experimental studies in mammals indicate that dihydroergotamine acts by activating α-adrenergic and 5-HT receptors, our cheminformatic study showed that in monogeneans, the kinases BRSK and MEKK1 might be targets of this drug. Therefore, the efficacy of dihydroergotamine against these parasites should be tested in future experiments.

The interaction between ergotamine and BRSK kinase was one of the strongest in our study. Ergotamine is an α₁-selective adrenergic drug used to treat migraine and has affinity for a wide variety of receptors, such as 5-HT, dopamine, and noradrenaline receptors, and induces contraction of the peripheral, pulmonary, cerebral, and coronary arteries (Dahlöf & Maassen Van Den Brink, 2012). In a study performed by Chan et al. (2018), ergotamine demonstrated remarkable efficacy in combating S. mansoni infection. In that study, schistosome-infected mice exhibited reduced parasite load and egg production when treated with ergotamine. Moreover, in the host, the drug mitigated the organ damage induced by the infection and lowered mortality rates. Although the repeated use of ergotamine in humans is not feasible due to secondary effects (Chan et al., 2018), it is unknown whether its mechanism of action would lead to similar side effects in fish. Therefore, future research should be directed toward validating the effect of ergotamine as a possible anti-monogenean agent. The activity of ergotamine on BRSK kinases, in addition to its known activity at G-protein-coupled receptors, and its selectivity, remains to be biochemically determined.

Another of the strongest interactions was between lomitapide and MEKK1 kinase. Lomitapide inhibits lipid transfer by directly binding to the microsomal triglyceride transfer protein in the liver and intestines, and is used to treat homozygous familial hypercholesterolemia in adults (Rader & Kastelein, 2014). Experimental evaluation of lomitapide against malaria parasites demonstrated its inhibition of β-haematin formation and parasite growth (de Sousa et al., 2020). In addition, this drug inhibits hemozoin formation in the parasite, the same mechanism by which traditional antimalarial drugs work (de Sousa et al., 2020).

This is the first study aimed at predicting and classifying kinases in monogeneans. From our bioinformatics pipeline, 22 kinases were identified as MDTs, of which five were proposed as leads. We suggest that drugs with conceivable anti-monogenean activity—as hinted at by a molecular docking analysis—could be experimentally evaluated for repurposing, particularly dihydroergotamine, ergotamine, and lomitapide given their high affinity with the monogenean kinases BRSK and MEKK1. Further research in this field is essential to enhance the control of parasitic diseases in finfish aquaculture.