Pan-phylum In Silico Analyses of Nematode Endocannabinoid Signalling Systems Highlight Novel Opportunities for Parasite Drug Target Discovery

Hidden Markov Model and BLAST Analysis

A Hidden Markov Model (HMM)-based approach has previously been reported for the identification of flatworm GPCRs (49). Predicted protein datasets were downloaded from WBP v14 for all nematodes with publicly available genome data (134 genomes; see File SI 2 ). Predicted protein datasets were concatenated for use as a predicted protein database for HMMERv3.3 HMM-searches. Profile HMMs for nematodes were constructed using predicted protein alignments of all candidate EC-effector protein homologs in C. elegans, C. briggsae, C. brenneri, C. japonica and C. remanei (see File SI 1 ). Multiple Sequence Alignments (MSAs) were generated using EMBL-EBI Clustal-Omega [https://www.ebi.ac.uk/Tools/msa/clustalo; (50)]. The hmmsearch function was employed to identify putative EC proteins within the nematode predicted protein datasets (see File SI 2 ) using default settings. Due to the volume of genomic data generated, the highest confidence hits were selected for each protein based on an inclusion threshold of E-value ≤0.01 and/or a score of ≥150.

Putative EC protein sequences identified via hmmsearch were then used as queries in reciprocal BLASTp searches of WBP (https://parasite.wormbase.org/Multi/Tools/Blast; default settings) and NCBI non-redundant (https://blast.ncbi.nlm.nih.gov; default settings) databases. Queries that failed to return a putative ECR or ECE hit within the top 4 BLAST results were excluded from further downstream analyses. All BLASTp searches returning negative hits (hits outside of the outlined inclusion criteria and those that returned no hits) were confirmed negative via tBLASTn searches; this approach also mitigated the impact of poor genome quality (false negatives) on our analyses where relevant.

Post-BLAST Sequence Analysis

Key EC ligand binding residues and functional motifs for mammalian and invertebrate ECRs and ECEs were identified from the published literature (see File SI 3 ) and all positive hits were examined visually for the presence of any key residues/motifs via multiple sequence alignments using EMBL-EBI Clustal-Omega [https://www.ebi.ac.uk/Tools/msa/clustalo (50)]. The presence of known protein family or structural domains in ECE BLAST hits was analysed using InterProScan [https:www.ebi.ac.uk/interpro/search/sequence-search, (51)]. Putative ECR hits were analysed for the presence of GPCR transmembrane (TM) domains using EMB TMpred server [https://embnet.vital-it.ch/software/TMPRED_form.html (52)]. Any putative ECE hits which lacked the required family/protein domains for designation as an ECE, or putative ECR sequences which possessed ≤ 3 TM domains, were excluded from further analysis. ECRs and ECEs were analysed for the presence of conservative substitutions of key residues or within binding motifs using WebLogo3 (53, 54) (see Supplementary Table 2 ).

Phylogenetic Analysis

MUSCLE was used to create multiple sequence alignments (MSAs) of protein sequences for all positive EC protein hits using MEGA X (55). For ECRs, alignments were manually edited to include only TM domains, for ECEs only functional domains were included in analysis. Functional domains for ECEs and TM domains for ECRs were identified via the NCBI Conserved Domains Database [CDD; https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi (56)]. Maximum likelihood (ML) phylogenetic trees were constructed using PhyML [http://www.phylogeny.fr (57)] from the domain only MUSCLE MSAs with default parameters and branch support assessment using the approximate likelihood ratio test (aLRT) with “SH-like” parameters. Trees were exported from PhyML in Newick format and were drawn and annotated using the Interactive Tree of Life [iTOL; https://itol.embl.de (58).

Endocannabinoid Signalling System Effectors Are Broadly Conserved Across Phylum Nematoda

Nematode EC-effectors (ECEs responsible for synthesis and degradation of EC ligands, and putative ECRs) appear to be broadly conserved across phylum Nematoda (Figure 2A). ECEs display greater conservation than putative ECRs across all nematodes, with the exception of the ECE nape-1/2 which appears to be absent from 37% of nematodes, many of which are representatives of clades 8 and 9 (Figures 2A, B). The more conserved profile of ECEs versus ECRs is consistent with the requirement for specific ECEs in the biosynthesis and degradation of EC-ligands and the potential for redundancy among ECRs, which has been documented in other systems (74, 75).

Genes encoding the key ECEs responsible for 2-AG synthesis and metabolism, DAGL-2 and ABHD-12 respectively, are co-conserved in 87% of nematodes (see Figure 2), indicating that a significant proportion of nematode species are likely to possess a functional, canonical, 2-AG synthesis and degradation pathway. Most (95%) nematode species examined possess an gene ortholog of the 2-AG biosynthesis enzyme DAGL-2, while the gene encoding ABHD-12, responsible for catalysing 2-AG metabolism is present in 88% of species. Parascaris equorum and Globodera pallida do not appear to possess either a dagl-2 or abhd-12 homolog, however their genome quality is lower as indicated by CEGMA/BUSCO scores [see Figure 2; (76, 77)].

Genes encoding NAPE and FAAH, the primary enzymes responsible for the synthesis and metabolism of AEA, are co-conserved in 55% of nematode species, suggesting that a significant proportion nematodes have the ability to metabolise AEA (Figure 2). This is corroborated by studies that have identified AEA in several nematodes, including C. elegans and N. brasiliensis, via mass-spectrometry (18, 39). While 88% of species encode at least one FAAH homolog, the NAPE-encoding gene is conserved in only 63% of nematodes (either nape-1 or -2). Notably, many of the species that lack a nape homolog are filarial nematodes, Strongyloides species, or members of the Diplogasteroidea superfamily (Figure 2). Whilst the absence of nape in some species may be explained by genome quality, in other species with robust genome data an alternative AEA synthesis pathway may exist (see below). Previous studies have demonstrated the presence of two, functionally divergent, nape orthologs in C. elegans that occupy adjacent genomic positions (nape-1, IV:3739520.3740880; nape-2, IV:3735470.3738925; [WormBase; (78)], share 73% sequence identity, and display complete conservation of the NAPE-PLD signature sequence (43). Our pan-phylum analysis confirms that C. elegans is the only nematode species that possesses two distinct NAPE-encoding genes. In all other species that encode NAPE the same gene ID was returned for both the nape-1 and nape-2 BLASTp searches (one positive nape hit was considered a positive return for both nape-1 and nape-2 and was designated nape-1/2; see Figure 2 and File SI 3 ). The presence of a single nape-1/2 in all of the nematodes examined, including other Caenorhabditis species, suggests that nape-1 and -2 may have arisen as a result of a relatively recent gene duplication event in C. elegans (43).

Of the seven putative ECRs included in this study NPR-19 and NPR-32 have been most closely linked to ECS in nematodes (16, 17, 27). Our in silico analysis reveals that 79% of the nematode species investigated in this study possess NPR-19 and 89% possess NPR-32, underscoring their putative importance to nematode biology. In addition to NPR-19 and -32, OCR-2 an ortholog of the human transient receptor potential vanilloid channel (TRPV) (79), has also been closely linked to ECS (25). OCR-2 regulates signal transmission and thus modulates several C. elegans behaviours (19, 80). It is interesting to note that 100% of nematode species examined in this study possess a gene encoding OCR-2, also suggesting a significant role in nematode biology. NPR-9 has been implicated in locomotion, regulation of innate immune responses, roaming and foraging behaviours in C. elegans (81–83), while GPR-55, the human ortholog of nematode NPR-9, is known to interact with human CB1 and CB2 to form functionally important heteromers (84, 85). The lower level of npr-9 conservation across phylum Nematoda revealed here (npr-9 conserved in 60% of species examined; Figure 2) may indicate a less conserved functional role for this putative ECR in nematodes.

Nematode EC-Effector Conservation Profiles Display Clade Specific Trends

Our data demonstrate distinct conservation patterns of EC-effectors across nematode clades. Whilst clade 9 and 10 nematodes exhibit the highest degree of EC-effector conservation, clade 8 species display the most reduced complements (Figure 2). This could be, in part, explained by lower CEGMA/BUSCO scores for clade 8 genomes, such that the profiles presented here may not be a true representation of EC-effector complements in this clade. Variable genome quality is an inherent caveat to in silico approaches such as those employed in this study however, the use of tBLASTn for all negative BLASTp returns can help to mitigate this in addition to continued improvements in genome quality (86).

For other clades there appear to be ECR specific trends. For example, npr-32 is broadly conserved across all clades, pointing towards a more conserved function for this receptor pan-phylum. Other ECRs display more clade-specific, restricted profiles, including nhr-49 which appears to be entirely absent from clade 2 (Figure 2A). Clade 2 nematode genome assemblies are high quality [as indicated by CEGMA/BUSCO scores; (48)] and likely provide a true reflection of the nhr-49 profile. Notably, npr-19 is also entirely absent from clade 12 species, indicating that alternative ECRs may contribute to the ECS pathway in these nematodes. Indeed, clade 12 nematodes exhibit broad conservation of other putative ECRs (e.g. npr-9, -32, nhr-49, octr-1, ser-4, ocr-2; see Figure 2A).

Filarids Possess Distinctive EC-Effector Profiles

Whilst there appear to be limited differences between free-living and parasitic nematode EC-effector profiles in general [96% ECR conservation in free-living nematodes (FL) vs 92% in parasitic species; 91% ECE conservation in free-living nematodes vs 83% in parasitic species], filarial nematodes display distinct EC-effector gaps. For example, filarial species completely lack genes encoding NPR-9 and NAPE-1/2, and have a significantly reduced FAAH-2 encoding gene profile (Figure 2A). The absence of NAPE-1/2 points towards the presence of an alternative AEA synthesis pathway in filarids, which may also be the case for other species e.g. Pristionchus spp. and Strongyloides spp. that likewise lack NAPE-1/2 encoding genes (see below; Figure 2A).

Nematodes That Lack NAPE-1/2 Possess Putative Alternative AEA Synthesis Enzymes

Our data indicate that 41 nematode species lack genes encoding the nematode NAPE-PLD ortholog NAPE-1/2, the enzyme primarily responsible for AEA synthesis (Figure 2A) (34). In vertebrates, two additional pathways have been implicated in the synthesis of AEA: (i) hydrolysis of NAPE by ABHD-4 forming the intermediate glycerophosphoanandamide (Glycero-p-AEA) which, following further hydrolysis by glycerophosphodiester phosphodiesterase 1 (GDE-1), results in the formation of AEA (see Figure 3A ii); (ii) hydrolysis of NAPE by phospholipase A2 (PLA-2) to form the intermediate N-acyl-1-acyl-lyso-PE (lyso-NAPE) and subsequently, following hydrolysis by lysophospholipase D (Lyso-PLD), results in the formation of AEA (see Figure 3A iii) (87–90). To determine if alternative AEA synthesis pathways exist in nematodes that appear to lack the classical NAPE-1/2 AEA biosynthesis pathway, we mined all available genome data for abhd-4, gde-1, pla-2 and lyso-PLD.

Our data reveal that 39 of the 41 species that lack NAPE-1/2 possess at least one putative alternative AEA synthesis pathway (see Figure 3B, File SI 3 ). 95% of these species possess both abhd-4 and gde-1 (alternative pathway shown in Figure 3A ii), while 56% encode both PLA-2 and Lyso-PLD (alternative pathway shown in Figure 3A iii); 56% of nematodes encode the enzymes for both alternative AEA synthesis pathways. Therefore these data suggest that nematodes which lack the classical NAPE-1/2 biosynthesis pathway predominantly synthesise AEA via ABHD-4 and GDE-1. However, within clades 8, 9 and 10 there are examples of species which may have the ability to synthesise AEA via either alternative pathway e.g. Strongyloides species and several filarial nematodes (Figure 3B).

Interestingly in mammals, in addition to PLA-2, other enzymes have been linked to the synthesis of lyso-NAPE (38), for example, ABHD-4 can remove an acyl group from NAPE resulting in the creation of lyso-NAPE (89). Therefore the presence of abhd-4 in 100% of the nematodes investigated in this study suggests that some nematodes that possess lyso-PLD but lack pla-2 could compensate by employing ABHD-4 in the synthesis of lyso-NAPE.

Interestingly, it appears that a subset of the species that do encode NAPE-1/2 (20 key species selected for their relevance to human, veterinary or plant health and broad clade representation), also encode all of the alternative AEA synthesis enzymes (abhd-4, gde-1, pla-2 and lyso-PLD), with the exception of Trichuris muris which appears to lack pla-2 (see File SI 4 ). It is unclear whether NAPE-1/2 is indeed the predominant AEA synthesis pathway in these species or, if the alternative pathways also contribute to AEA production.

In silico evidence for the presence of putative, alternative, AEA biosynthesis enzymes in a range of therapeutically relevant nematodes strongly suggests that NAPE-1/2 independent pathways may contribute to AEA synthesis in these species. Further analysis, including mass spectrometry to isolate AEA, will begin to unravel the importance of alternative AEA biosynthesis pathways in these nematodes.

Nematode FAAH Homologs Display Conservative Substitutions in a Key AEA Binding Site

The mammalian AEA hydrolysis enzymes FAAH-1 and FAAH-2 possess four key residues required for catabolic activity (FAAH-1: K142, M191, S217, S241; FAAH-2: K131, C180, S206, S230) (42, 91, 92). Analysis of nematode FAAH-1 homologs identified in this study (see Figure 2) revealed that >90% possess the key mammalian binding site residues K142, S217 and S241, while 83% of identified nematode FAAH-2 homologs possess K131, S206 and S230 (see Figures 4A–C). However, 87% of nematode FAAH-1 homologs display a conservative substitution (methionine for leucine) at position 191 while 80% of nematode FAAH-2 homologs substitute cysteine for leucine at position 180 (Figures 4A–C). M191 has been implicated in mammalian FAAH-1 EC-derivative binding, with studies demonstrating lipophilic interactions between this residue and partial cannabinoid receptor agonists, linking it to AEA binding (93), however the significance of C180 (FAAH-2) is less clear (94–96).

Figure 4

Nematode FAAH homologs conserve key functional domains and motifs, but display conservative substitutions in a key AEA binding site. (A) Protein sequence alignment of Homo sapiens fatty acid amide hydrolase (FAAH-1), H. sapiens FAAH-2, Caenorhabditis elegans FAAH-1 and C. elegans FAAH-2. Amino acids are highlighted in the same colour if > 60% of residues are conserved. Legend denotes conserved amino residues, key ligand binding residues in FAAH-1 and FAAH-2, and the amidase signature domain. AA denotes amino acid; H. sapiens denotes Homo sapiens, C. elegans denotes Caenorhabditis elegans. (B) Amino acid sequence-logo representing sequence diversity between key residues in nematode FAAH-1 homologs (consensus) vs H. sapiens FAAH-1 [O00519, FAAH1_HUMAN] and, (C) Amino acid sequence-logo representing sequence diversity between key residues in nematode FAAH-2 homologs vs H. sapiens FAAH-2. Key EC binding site residues were derived from H. sapiens FAAH-1 and -2 and are detailed in the top column of each table. Colours indicate hydrophobicity of amino acid residues (hydrophilic residues are blue, neutral residues are green and hydrophobic residues are black).

Mammalian FAAH-1 and -2 also possess a highly conserved 130 bp amidase signature domain that enables enzyme characterisation (97, 98); this domain is conserved in 94% of the nematode species examined in this study (Figure 4A).

While these data demonstrate that nematodes possess homologs for FAAH enzymes that display broad conservation with vertebrates, the presence of a distinct substitution in a key AEA binding site across many nematodes may highlight the potential for drug target selectivity towards parasitic nematode species.

Nematodes Possess Homologs for the Mammalian 2-AG Degradation Enzyme ABHD-12

In mammals, MAGL is primarily responsible for 2-AG degradation and thus the termination of EC signalling (31, 99). Additional 2-AG degradation enzymes, ABHD-6 and -12, have also been reported (33, 100) however, ABHD-6 is thought to play a less significant role in 2-AG degradation (31, 99).

The nematode literature presents conflicting data on the identity of the 2-AG degradation enzyme; indeed prior to this study it was unclear whether nematodes possess a true ortholog of MAGL (19, 27) or, if the nematode 2-AG degradation enzyme is actually an ABHD-12 ortholog (11, 18). It is interesting to note that previous work indicates that magl homologs are absent in some nematode species (18), which has been confirmed here using our pan-phylum in silico approach (Figure 2A). Indeed, in our analyses only two nematode species (of the 109 in this study) returned an magl-like sequence within the top 5 BLASTp hits (Strongylus vulgaris, SVUK_0001964001; Steinernema scapterisci, L892_g30127.t1), both of which failed to meet E-value inclusion criteria. The absence of MAGL in nematodes indicates the presence of an alternative 2-AG degradation enzyme.

In light of the limited role for ABHD-6 in mammalian 2-AG degradation and the absence of MAGL across nematodes (as reported here and in previous studies), we focused our attention on ABHD-12 as a putative alternative to MAGL in nematodes. Pan-phylum analysis of nematode genomes identified orthologs for ABHD-12 in 88% of nematode species examined. Significantly, phylogenetic analyses of these putative abhd-12 homologs demonstrated that 99% of nematode BLAST returns for abhd-12 cluster more strongly with human abdh-12 than human magl (Figure 5 and Supplementary Figure 1 ) suggesting that the nematode 2-AG hydroxylase enzymes identified here, and originally designated as magl in C. elegans [Y97E10AL.2; (19)], are orthologous with abhd-12.

Nematode NPR-19 and -32 Orthologs Possess Key Functional Motifs and EC Binding Residues

NPR-19 has been identified as a putative EC-GPCR in C. elegans and appears to modulate key aspects of nematode biology including axon regeneration, locomotion, modulation of noiceception and feeding, and therefore may represent a promising anthelmintic target (27).

78% of the nematodes examined in this study possess an npr-19 ortholog (see Figure 2 and Figure 6A). Phylogenetic analysis demonstrates that nematode npr-19 orthologs failed to cluster with human CB1 and CB2, confirming that NPR-19 is not the direct ortholog of the human EC-GPCRs receptors (CB1 and CB2) (Figure 6A and Figure SI 1 ). Further analysis of the nematode npr-19 orthologs revealed that whilst the nematode NPR-19 consensus sequence has only 23% similarity to human CB1, several of the known human CB1 EC binding residues (N46, D88, F189, L193, F379, S383) (101, 102) are conserved (Figure 6B). This indicates that nematode NPR-19 orthologs are likely to possess EC ligand binding capacity and aligns with previous work in C. elegans (27). When examined at the clade level, sequence analysis showed that two conservative substitutions are present at position K192 (82% of clade 2 species substitute K192 for N192), while species in clades 8, 9, 10 and 12 substitute K192 for D192. K192 forms a hydrogen bond with the amide oxygen of AEA, implicating this residue in ligand binding (101) (see Figure 6B).

Figure 6

Maximum likelihood phylogeny of nematode npr-19 homologs. (A) 85 nematode NPR–19 homologs are shown in addition to Homo sapiens CB1 and CB2 [P21554 (CNR1_HUMAN), P34972 (CNR2_HUMAN)], H. sapiens GPR-55 [Q9Y2T6 (GPR55_HUMAN)] and several Caenorhabditis elegans biogenic amine receptors (serotonin [SER-1,-3,-6; WBGene00004776, WBGene00004778, WBGene00021897] dopamine [DOP-1-3; WBGene00001052, WBGene00001053, WBGene00020506] and tyramine [TYRA-2 and -3; WBGene00017157, WBGene00006475]). Non-NPR–19 homologs are marked with a red asterisk (*). Outer ring denotes nematode clade and coloured circles represent branch support values. Tree was generated from an alignment trimmed to include functional domains. Branch supports indicate statistical support from approximate likelihood ratio test (aLRT). (B) Amino acid sequence-logo demonstrating sequence diversity between nematode NPR-19 and NPR-32 orthologs (consensus) and H. sapiens CB1 [P21554 (CNR1_HUMAN)]. Known vertebrate endocannabinoid binding and GPCR motifs/residues are indicated in the top row of the sequence logo table, transmembrane regions 1-7 are indicated by orange boxes and numbers, amino acid colours indicate hydrophobicity of amino acid residues (hydrophilic residues are blue, neutral residues are green and hydrophobic residues are black).

NPR-32 is implicated in C. elegans axon regeneration (17) and, in addition to NPR-19, is believed to be a putative EC-GPCR (27). In this study we identified 97 NPR-32 homologs (Figure 2 and Supplementary Figure 1 ) which share (consensus sequence) only 21% identity with the human EC-GPCR CB1, but conserve several key residues (N46, D88, K115, F189; Figure 6B) that are believed to be important for ECS function (17).

In addition, CB1 possesses a “toggle switch” (residue W356), a putative molecular hinge that interacts with F200 to change the form and state of the receptor which in turn aids EC ligand binding (103, 104). In nematodes the “toggle switch” W356 is conserved in >85% of NPR-19 and -32 orthologs, whereas NPR-19 F200 is substituted for W200 in 59% of species and NPR-32 F200 is substituted for L200 in 90% of species, (see Figure 6B). The significance of these observations will be revealed through molecular docking studies, crystal structure analysis, and functional genomics in relevant parasite species, and will inform the role and importance of these receptors in nematode ECS biology.

The authors wish to thank WormBase ParaSite for helpful assistance with transcriptome resources.