Polymorphism in ion channel genes of Dirofilaria immitis: Relevant knowledge for future anthelmintic drug design.

Dirofilaria immitis, a filarial parasite, causes cardiopulmonary dirofilariasis in dogs, cats and wild canids. The macrocyclic lactone (ML) class of drugs has been used to prevent heartworm infection. There is confirmed ML resistance in D. immitis and thus there is an urgent need to find new anthelmintics that could prevent and/or control the disease. Targeting ion channels of D. immitis for drug design has obvious advantages. These channels, present in the nematode nervous system, control movement, feeding, mating and respond to environmental cues which are necessary for survival of the parasite. Any new drug that targets these ion channels is likely to have a motility phenotype and should act to clear the worms from the host. Many of the successful anthelmintics in the past have targeted these ion channels and receptors. Knowledge about genetic variability of the ion channel and receptor genes should be useful information for drug design as receptor polymorphism may affect responses to a drug. Such information may also be useful for anticipation of possible resistance development. A total of 224 ion channel genes/subunits have been identified in the genome of D. immitis. Whole genome sequencing data of parasites from eight different geographical locations, four from ML-susceptible populations and the other four from ML-loss of efficacy (LOE) populations, were used for polymorphism analysis. We identified 1762 single nucleotide polymorphic (SNP) sites (1508 intronic and 126 exonic) in these 224 ion channel genes/subunits with an overall polymorphic rate of 0.18%. Of the SNPs found in the exon regions, 129 of them caused a non-synonymous type of polymorphism. Fourteen of the exonic SNPs caused a change in predicted secondary structure. A few of the SNPs identified may have an effect on gene expression, function of the protein and resistance selection processes.

Introduction

Dirofilaria immitis is a mosquito-borne-filarial nematode that causes dirofilariasis or heartworm disease in dogs and cats, and occasionally infects humans. For approximately the last 25 years, prevention of heartworm infection has been solely dependent on a single drug class, the macrocyclic lactones (MLs). Two sub-groups of MLs, namely the avermectins - ivermectin (IVM) and selamectin, and the milbemycins – moxidectin and milbemycin oxime - are used as heartworm preventatives (Lespine et al., 2012). However, there have been reports of loss of efficacy (LOE) of these MLs against heartworm in the USA, especially in those areas where heartworm challenge is high (Hampshire, 2005). Recent studies have confirmed, using in vivo efficacy studies and genetic analysis, actual ML resistance in D. immitis (Bourguinat et al., 2011a, Bourguinat et al., 2015, Pulaski et al., 2014). There also exists some evidence of ML resistance in the human parasite, Onchocerca volvulus (Osei-Atweneboana et al., 2007, Osei-Atweneboana et al., 2011, Nana-Djeunga et al., 2012, Nana-Djeunga et al., 2014, Pion et al., 2013), a filarial nematode which is closely related to D. immitis. Considering the range of nematode parasites in animals and humans, anthelmintic resistance occurs against virtually all of the major families of broad spectrum anthelmintics, including the MLs, the benzimidazoles, levamisole and monepantel (Kotze et al., 2014). In the absence of effective vaccines or means to prevent infection, both treatment and prophylaxis for most helminth parasites relies mainly on the use of these drugs. Since resistance is an increasing, serious concern there is an urgent need to develop new classes of anthelmintic drugs.

Current programs of anthelmintic drug discovery, especially those that operate in industrial settings, focus primarily on the discovery of new drugs for veterinary indications: mainly gastrointestinal nematodes of livestock and companion animals, and canine heartworm (Woods et al., 2007, Woods et al., 2011, Woods and Knauer, 2010, Geary et al., 2015). Furthermore, it is worthy of note that many of the anthelmintics used in human medicine were originally developed for the veterinary sector (Geary et al., 2015). In terms of economic factors, in the USA alone, with over 80 million dogs, and with heartworm drugs costing US $75–100/dog per year (Godel et al., 2012), there exists an attractive market for heartworm drug development. Any new pharmaceutical for heartworm control may also benefit efforts to control filarial and potentially other parasitic infections in humans (Wolstenholme et al., 2015).

The neuromuscular system of D. immitis shows great potential as a drug target; nematodes have well-developed neuromuscular systems that control motility, navigation, feeding, mating and responses to environmental cues, which are essential for their survival, development and reproduction (White et al., 1986, Perry et al., 2004, Greenberg, 2014). Interference with any of these activities can paralyse, kill or stop reproduction in the parasite (Greenberg, 2014). Ion channels and associated receptors that underlie neuromuscular systems are targets for a wide variety of naturally occurring toxins and synthetic compounds (Camerino et al., 2007). Ion channels are targets for many nematocidal drugs currently in the market (Wolstenholme, 2011). Furthermore, any new drug that acts on these receptors can be easily screened, in a semi-automated manner, for anthelmintic activity based on its effect on worm motility (Wolstenholme, 2011). Ion channels of nematodes may be classified by the type of ions that are allowed to pass through them (e.g. cations (Na+, K+ and Ca++) or anions (Cl−)), and by the type of gating: voltage gated (sodium, calcium or potassium channels) or ligand gated (in which the ligand may be an amino-acid such as glutamate, ɤ-amino butyric acid (GABA), or a biogenic amine such as serotonin, tyramine or dopamine) (Wolstenholme, 2011, Greenberg, 2014). Among the ligand-gated ion channels (LGICs), the cysteine-loop (cys-loop) superfamily includes cation-permeable acetylcholine receptors as well as anion-selective channels gated by GABA, glutamate, 5-hydroxytryptamine, dopamine or tyramine (Raymond and Sattelle, 2002, Hobert, 2013, Ringstad et al., 2009, Lees et al., 2012). A complete inventory of the homologs of these channel genes, in D. immitis, is so far unknown.

Due to current screening activities based on worm motility, there are good chances that ion channels and receptors in D. immitis will be targets for new anthelmintics. Genetic variability, in the form of single nucleotide polymorphism (SNPs), insertions or deletions (indels), in ion channel and receptor genes may create structural changes in the protein products. Such changes may alter the physiochemical or structural properties, disrupting folding, affecting stability or function of the protein, or making it totally a null receptor for a drug effect (Lahti et al., 2012). Genetic changes within ion channel genes may also modulate their expression levels and such changes, though not examined here in this study, could be a source of variability in channel properties (Mulley et al., 2005). Knowledge of possible polymorphism in ion channels and associated receptors in D. immitis may be relevant for drug design processes. In addition, such knowledge can be used to make sure that any new drug is active against all of the allelic forms of the target, including allelic variants found in LOE populations. Furthermore, heterogeneity in an ion channel drug target may facilitate resistance development (Prichard, 2001, Bourguinat et al., 2011c) and knowledge of that heterogeneity may be informative in anticipating possible resistance development. The objectives of this study were to identify all the putative ion channel genes/subunits in D. immitis by comparative genomic approaches and to analyze them for polymorphism.

Materials and methods

Identification of all the ion channel genes of D. immitis

The nDi.2.2 genome of D. immitis is in draft stage (Godel, 2012) and the gene annotation is not complete (Bourguinat et al., 2015). A complete inventory of ion channel genes of D. immitis is not available. To identify all the putative ion channel genes/subunits in the genome of D. immitis, a complementary approach was followed. Protein-encoding genes from the assembled D. immitis contigs were predicted and validated by Godel et al. (2012) using three parallel strategies; i) prediction with the ab initio gene finders SNAP (Korf, 2004) and Augustus (Stanke and Morgenstern, 2005) using the training set of Brugia malayi, ii) direct alignment to B. malayi proteome and iii) alignment to the RNA-seq assembly (Godel, 2012, Kumar, 2013). Thus, all the ion channel genes so far annotated and available in the GFF (Generic Feature Format) file format under Nuclear Annotation Freeze nDi.2.2.2 in the D. immitis website (http://nematodes.org/genomes/dirofilaria_immitis/) were used. To study the un-annotated homologs of D. immitis ion channel genes, nucleotide sequences of all the genes or subunits that belong to the family of cys-loop LGIC, voltage-gated (e.g., potassium, calcium) and other channel types of all nematodes were extracted from available databases such as NCBI (http://www.ncbi.nlm.nih.gov/), Wormbase (http://www.wormbase.org/#01-23-6), Broad Institute (https://www.broadinstitute.org/) and NEMBASE4 (http://www.nematodes.org/nembase4/). These sequences were then blasted (BLASTN 2.2.25) (Altschul et al., 1997) in the nucleotide blast server v2.2 (http://nematodes.org/genomes/dirofilaria_immitis/) to locate each of the putative ion channel genes in the scaffolds of the D. immitis nuclear genome (version nDi.2.2.2) (http://salmo.bio.ed.ac.uk/cgi-bin/gbrowse/gbrowse/nDi.2.2.2/).

Synchronized file generation to assess polymorphism in ion channels

Pooled worm samples and the method followed for synchronized file generation for genetic variability analysis were as previously described (Bourguinat et al., 2015) except that the reference genome used was different in this study. Briefly, our study samples included a total of 122 worms from 17 ML susceptible dogs from the USA (Missouri isolate), Gran Canary (Spain), Grenada (West Indies) and Italy (Po Basin, Northern Italy). All worms from each country population were pooled, resulting in 4 pools for susceptible worms. Other phenotypic populations included ∼8000 microfilariae (mfs) from each of 4 ML-LOE dogs, originally from four different locations in the USA (Mechanicsville, VA; New Orleans, LA; Haywood County, TN; Monroe, LA). Loss of efficacy to MLs in such dogs was assessed, as described in Bourguinat et al. (2015). Mfs from each individual LOE dog were analyzed as a pool, resulting in 4 pools for LOE samples. All 8 pooled samples were subjected to whole genome sequencing using the HiSeq2000 platform from Illumina® by Illumina Next Generation Sequencing (The McGill University and Genome Quebec Innovation Centre) and Bam files, corresponding to the alignment of the reads from each population against the reference genome, were generated. D. immitis nuclear genome v2.2 (http://nematodes.org/genomes/dirofilaria_immitis/) was used as the reference genome for the alignment. The program, PoPoolation2, which allows comparison of nucleotide frequencies of two or more populations and identifies any significant differences in nucleotides at a position among populations (Kofler et al., 2011), was used to generate a synchronized file. After filtering for base quality, the synchronized file contained nucleotide (read) frequencies for every position for every population along the length of the reference genome. Based on the population read frequencies obtained for every gene/subunit in the synchronized file, a locus was considered to be polymorphic and different if a change in nucleotide frequency at that locus was >15% different between populations, a threshold limit set on the basis of the number of reads and the base pairs covered for SNP genotyping.

Nomenclature, classification of identified ion channels

Naming of the identified ion channel genes/subunits was based on their respective orthologs in Caenorhabditis elegans during BLAST searches. These names may change in the future based on phylogenetic analysis. The gene description for the ion channels was the same as given during annotation; available as nuclear annotation freeze nDi.2.2.2 in the website http://nematodes.org/genomes/dirofilaria_immitis/. Classification of each ion channel of D. immitis into sub-groups/sub-classes was done as previously described for the neuronal genome of C. elegans (Hobert, 2013).

Assessment of position, type and impact of identified polymorphism

A GFF file containing the nuclear annotations of the nDi.2.2.2 version of the D. immitis genome (http://nematodes.org/genomes/dirofilaria_immitis/) was used to identify whether a SNP was located in an intron or an exon region of a transcript. Protein sequence predicted for every gene annotation, available as file “Protein predictions from nDi.2.2.2” in the same website, was used to locate the position of amino acid changes within a protein. For the un-annotated genes, all the exon regions of a gene were predicted, combined and translated to confirm the full length of a gene. The type of single nucleotide polymorphism, synonymous (nSNP) or non-synonymous (nsSNP), was determined using Sequencher software 4.10.1 (Gene Codes Corporation, 2010). PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred/), a general purpose secondary structure predictor (SSP) and a highly accurate SSP tool for ionic channels (Konopka et al., 2009) was used to predict the effect of non-synonymous mutations in the secondary structure of a gene/subunit.

Focus on potential drug targets and their polymorphism

To highlight any potential D. immitis drug target (and its polymorphic forms), a similar filtering methodology was used as was followed for B. malayi (Kumar et al., 2007) and D. immitis (Godel et al., 2012). One of the filters used was that their respective orthologous gene should have had deleterious effects following gene knockout studies in C. elegans. Such deleterious effects could be embryonic lethal/larval arrest, shortened life span, locomotion variant, organism development variant, sluggish/fainter, slow growth, egg size defective, or pharyngeal pumping variant, as shown in the Wormbase website (http://legacy.wormbase.org/). A second type of filter used was the absence of a BLASTP (http://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome) hit with an E-value below 10−5, in the predicted proteomes of Homo sapiens and Canis lupus familiaris (Godel et al., 2012, Kumar, 2013).

Results

Number of ion channel genes/subunits studied for SNP genotyping

A total of 1249 nucleotide sequences for all the known ion channel genes/subunits of both free living and parasitic nematodes were extracted from NCBI and nematode databases (Supplementary file 1, Table S1). Complete information of these extracted nucleotide sequences is available in Supplementary file 2, Table S2 (worksheet “Details of extracted sequences”). We identified 224 genes/subunits in total (Supplementary file 2, Table S2, worksheet “Ion channel genes in D. immitis”) that covered 965,735 base pairs (bps) (Table 1); about 1.1% of the ∼84.2 Mb sized nuclear assembly of D. immitis (Godel, 2012). All of the 224 genes/subunits included for SNP genotyping were either identified from BLAST hits in the D. immitis nucleotide blast server (nDi.2.2) with extracted nucleotide sequences of ion channels of related nematodes or those ion channels were already annotated in the D. immitis genome (http://salmo.bio.ed.ac.uk/cgi-bin/gbrowse/gbrowse/nDi.2.2.2/).

Table 1

A summary of ion channel gene/subunit sequences investigated for the Dirofilaria immitis SNP analysis.

Parameter	Results
Number of ion channel genes/subunits studied	224
Total bases (bp) covered	965,735
Intronic bases	698,720
Exonic bases	267,015
Number of SNPs in introns	1508
Number of SNPs in exons	254
Number of missense causing SNPs in exons	129

Polymorphic pattern in ion channel genes of D. immitis

Among the eight pooled samples studied for polymorphism within the 224 genes/subunits of ion channels in D. immitis, we identified 1762 SNPs (Supplementary file 3, Table S3) with a nucleotide diversity rate of 0.18% (Table 1). This number of SNPs was obtained after using 15% as the threshold level for SNP calling at a position. If the alternative nucleotide frequency at a locus was set at ≥20% in order for a position to be considered polymorphic, then 1504 of the 1762 loci were retained. About 85.5% of the SNPs identified were in intron regions with an intronic SNP rate of 1/463 bp, whereas the remaining SNPs were in exons with an exonic SNP rate of 1/1051 bp. One hundred and twenty nine of 254 SNPs in exons caused non-synonymous polymorphism, 14 of which changed secondary structure of the protein as predicted by PSIPRED. The nsSNPs identified in this analysis are in most cases described in detail. The polymorphic rates, calculated for each population as a percentage of the total SNPs identified, are given in Table 2. A similar polymorphic rate (around 70%) was seen in pooled populations from the USA, Grand Canary and Grenada. However, a lower polymorphic rate (43%) was identified in the pooled samples from Italy.

Table 2

Genetic variability within ion channels and receptors genes of Dirofilaria immitis in populations from different locations.

Geographical locations	No. of SNPs found	Polymorphic rate (%)a
USA	1203	75.2
Grand Canary (Spain)	1135	70.9
Grenada	1120	70.0
Italy	688	43.0

Number of SNPs identified (as percentage) in each country based on total SNPs identified.

SNPs in Cys-loop LGIC

The genome of D. immitis was found to contain 81 genes/subunits that could encode 42 unique genes (15 nAChR, 1 ionotropic glutamate receptor (iGluR), 6 ACh-gated, 6 GABA-gated, 5 glutamate-gated, 2 biogenic amine-gated and 7 genes in a diverse group). From all of the genes encoding Cys-loop LGICs studied for SNP analysis, 410 SNP loci were identified. Of the 64 SNPs in exon regions, 39 were nsSNP. In this superfamily of cys-loop LGICs, an overall polymorphic rate of 0.17% was noticed with 0.14 and 0.33% in cation-selective nAChR and ionotropic glutamate receptors, respectively. Anion-selective channels such as acetylcholine-gated, GABA-gated, glutamate gated and the biogenic amine-gated subgroup were found to have polymorphic rates of 0.08, 0.32, 0.18 and 0.15%, respectively (Table 3).

Table 3

Summary of SNPs and SNP rate calculated for each sub-family of cys-loop ligand-gated ion channel genes.

Ion selectivity	Sub-groups	No. of subunits studied	Base pairs covered	No. of SNPs identified			Polymorphic rate (%)
				Introns	Exons	Total
Cation	nAChR-type	34 (15)	88,768	107	19	126	0.14
Cation	iGluR	2 (1)	18,180	53	7	60	0.33
Anion	ACC	6 (6)	37,509	26	3	29	0.08
Anion	Aminergic	3 (2)	33,752	40	9	49	0.15
Anion	GABA-gated	14 (6)	21,198	59	8	67	0.32
Anion	GluCl	8 (5)	18,564	26	8	34	0.18
Unknown	diverse	14 (7)	28,484	35	10	45	0.16

No. of unique genes are indicated in brackets.

SNPs in the exon regions of nAChR-types of the cys-loop LGIC superfamily, such as acr-8, acr-11 and acr-16 were recorded (Table 4). The gene acr-8 identified in the scaffold nDi.2.2.scaf00069 had a deletion mutation at the 65th amino acid position and the percentage of deletions was 33% in susceptible populations compared to 64% in LOE populations. The other subunit of acr-8 in the same scaffold was identified to have three SNPs with R364G in particular causing a change in predicted secondary structure of the protein (Fig. 1). In the homolog of coelomocyte-specific gene cup-4, identified in the scaffold nDi.2.2.scaf01506, a SNP was identified at position 141 between aliphatic isoleucine and acidic amide asparagine. In C. elegans, the cup-4 gene was shown to be required for efficient endocytosis of fluids by coelomocytes (Patton et al., 2005) and the loss-of-function RNAi studies showed increased sensitivity to oxidative stress and reduced lifespan (Park et al., 2009). Among the acetylcholine gated chloride channel genes studied, both acc-3 and lgc-47 were identified with the SNPs Y311D and Y256D, respectively. D. immitis has five GluCl encoding genes namely glc-2, glc-3, glc-4, avr-14 and avr-15; three of these genes had SNPs in exon regions. The homolog of the avr-14 gene in D. immitis had two SNPs, C382F and C385F. The presence of these two adjacently positioned SNPs was predicted to cause a change in secondary structure in the flanking region (Fig. 1). The glc-2 homolog was found with two SNPs, Q90E and V381A, with the latter causing a change in the predicted secondary structure. Another GluCl gene homolog of glc-4 had a SNP at T375K. The glutamate-gated cation channel (ionotropic glutamate receptor) gene glr-1 homolog in D. immitis, was identified with two SNPs, N287S and P809T. The former was specific to LOE, and the latter to susceptible populations. GABA-gated channel gene gab-1 had a SNP, I20T identified only in the susceptible populations, whereas the other gene, unc-49 with N57D polymorphism was in both susceptible and LOE populations. SNP loci were also identified in the homologs of the dopamine gated channel gene lgc-53, tyramine gated channel gene lgc-55 and also in uncharacterized members of the cys-loop LGIC superfamily, such as lgc-39, lgc-41, lgc-44 and lgc-45.

Table 4

SNP analysis of cys-loop LGIC genes/subunits of Dirofilaria immitis.

Putative gene name/Description	Coverage in nDi.2.2 scaffold	SNP position in the scaffold	Nucleotide change	Amino acid polymorphism
(lgc-27) ligand-gated ion channel family member	nDi.2.2.scaf00004 (380863..387209)	384,333	TTG↔GTG	L137V
		384,340	GTA↔GGA	V139G
(lgc-30) neurotransmitter-gated ion-channel ligand binding domain containing protein	nDi.2.2.scaf00014 (297209..303637)	299,340	ACT↔AAT	T257N
(acr-8)b nicotinic acetylcholine receptor alpha subunit 8	nDi.2.2.scaf00069 (86715..87299)	87,289	CTT↔C-T	L65Deletion
(acr-8) nicotinic acetylcholine receptor alpha subunit 8	nDi.2.2.scaf00069 (89908..93744)	92,048	TAT↔TTT	Y264F
		93,364	CGA↔GGA	R364Gc
		92,025	TTA↔TTT	L256F
(acr-11)a cre-acr-11 protein	nDi.2.2.scaf00565 (17269..21582)	19,839	CAA↔CCA	Q46P
(acr-16)a acetylcholine receptor subunit alpha-type	nDi.2.2.scaf07899 (1..688)	305	TCA↔TTA	S53L
		337	TTG↔TTT	L43F
		345	CCG↔TCG	P40S
(unc-63) nicotinic acetylcholine receptor alpha subunit	nDi.2.2.scaf00751 (16790..21721)	16,834	CTT↔TTT	L353F
(cup-4)a acetylcholine receptor-like protein	nDi.2.2.scaf01506 (618..6786)	1869	ATT↔AAT	I141N
(acc-3) cre-acc-3 protein	nDi.2.2.scaf00812 (363..7136)	2540	TAT↔GAT	Y311D
(lgc-47) neurotransmitter-gated ion-channel ligand binding domain containing protein	nDi.2.2.scaf00030 (142465..147144)	144,803	TAC↔GAC	Y256D
(glc-4) glutamate-gated chloride channel subunit beta	nDi.2.2.scaf00002 (266371..271808)	267,383	ACG↔AAG	T375K
(glc-2)a Glutamate-gated chloride channel	nDi.2.2.scaf00035 (124763..129718)	125,276	CAA↔GAA	Q90E
		129,561	GTG↔GCG	V381Ac
(avr-14)a glutamate-gated chloride channel	nDi.2.2.scaf00410 (1..5378)	5254	TGT↔TTT	C382Fc
		5263	TGT↔TTT	C385Fc
(glr-1) Ionotropic glutamate receptor	nDi.2.2.scaf00632 (11556..19637)	15,358	AAT↔AGT	N287S
		19,548	CCA↔ACA	P809T
(gab-1) gamma-aminobutyric acid receptor subunit beta-like	nDi.2.2.scaf01694 (74..792)	132	ATC↔ACC	I20T
(unc-49) gamma-aminobutyric-acid receptor beta subunit	nDi.2.2.scaf01074 (11518..13243)	12,418	AAC↔GAC	N57D
(lgc-53) neurotransmitter-gated ion-channel ligand binding domain containing protein	nDi.2.2.scaf00139 (44015..65641)	49,539	CTA↔CCA	L691P
(lgc-53) neurotransmitter-gated ion-channel ligand binding domain containing protein	nDi.2.2.scaf00238 (35622..42385)	36,421	CGA↔GGA	R375G
		36,422	TCG↔GCG	F374L
		36,430	AGA↔AAA	R362K
		36,442	CAA↔AAA	Q368Kc
		36,459	TCC↔ CCC	S372P
(lgc-39) cre-lgc-39 protein	nDi.2.2.scaf02810 (1..1288)	531	CCA↔GCA	P71A
(lgc-41) ligand-gated ion channel family member	nDi.2.2.scaf00010 (375585..379494)	379,448	CAC↔CCC	H335P
		379,483	GCA↔ACA	A347T
(lgc-44) neurotransmitter-gated ion-channel ligand binding domain containing protein	nDi.2.2.scaf01696 (393..4458)	557	CTT↔ATT	L139I
		596	GCA↔ACA	A126T
		4443	GCA↔ACA	A6T
(lgc-45)b cre-lgc-45 protein	nDi.2.2.scaf02553 (1..1867)	1356	GAG↔AAG	E59K
(lgc-55) cre-lgc-55 protein	nDi.2.2.scaf00086 (64805..70168)	67,173	CCT↔ACT	P177T

Homolog genes with detrimental RNAi phenotype in C. elegans.

Absence of a significant BLAST hit (E-value below 10−5) in the predicted proteomes of H. sapiens and C. lupus familiaris.

SNP causes change in secondary structure of the protein as predicted by PSIPRED.

Fig. 1

PSIPRED predicted secondary structure changes in ion channels due to polymorphism. Secondary structure changes predicted for each ion channel are represented by cartoons: for helix (H), for strand (E) and for coil (C). Confidence value for prediction at each position is given as a series of blue bar graphs - . Any change in amino acid due to polymorphism is highlighted in yellow.

Potassium channels

Seventy three of the gene annotations in the nDi.2.2 genome encode 43 potassium channel genes (plus 1 auxiliary subunit), with representatives for each of the three major classes: 2-pass (3 genes), 4-pass (23 genes) and 6-pass (14 genes) transmembrane proteins. With 284 SNPs in these voltage gated channels, a polymorphism rate of 0.09% was observed; lowest among the genes studied (Table 5). A SNP at amino acid position 99 of the inward rectifier potassium channel gene irk-3, caused an amino acid variant between arginine (basic amino acid) if the corresponding codon was CGA and glycine (aliphatic amino acid) if the corresponding codon was GGA. Homologs of potassium channels with 4-transmembrane topology, namely twk-48 and twk-8, were polymorphic with SNPs I171T and M468V, respectively. The latter caused a change in predicted secondary structure. The SNP in the gene twk-8 was in the susceptible populations only, but was well conserved with methionine in the LOE populations. A subunit of the twk-47 gene homolog in D. immitis was identified with a SNP V9A, with predominantly valine at this position in the susceptible populations and alanine in the LOE populations. Two subunits of the gene twk-18, which existed as gene clusters in the same scaffold of nDi.2.2.scaf01557, were found to be polymorphic, with one of the subunits having 4 SNPs, including a SNP that leads to a stop codon. The other subunit had 2 SNPs (M1V and H5D) towards the 5′ end of the gene. Homologs of the six-pass transmembrane channel genes in D. immitis such as exp-2, kvs-4 (Kv2 subfamily), and shw-1, shw-3 (Kv3 subfamily) were found to have a single SNP in each (see Table 6). In the kvs-4 homolog gene, a locus with the SNP T208I was more common in the LOE populations than in the susceptible populations. A couple of SNPs, identified in the potassium channel gene containing a tetramerisation domain (homolog of uncharacterized C. elegans gene F59F3.6), with SNP I194M in particular, being predicted to cause a change in the secondary structure of the protein (Fig. 1). A secondary structure change caused by a SNP (H346R) in a calcium activated potassium channel gene, kcnl-2, was found to be present only in the LOE populations. An auxiliary subunit and a multipass transmembrane protein, unc-93 had a SNP R131Q that was only observed in the susceptible populations.

Table 5

Summary of genes/subunits of voltage-gated and other channel types studied for polymorphism, SNPs identified and the SNP rate calculated among Dirofilaria immitis populations from different geographical locations.

Channel types	No. of subunits studied	Base pairs covered	No. of SNPs identified			Polymorphic rate (%)
			Intron	Exon	Total
Calcium	26 (9)	183,111	278	46	324	0.18
Potassium	72 + 1* (44)	313,602	245	39	284	0.09
Chloride	13 (9)	55,311	224	41	265	0.48
DEG/ENaC	8 (7)	52,637	91	14	105	0.20
TRP	14 (7)	82,860	280	39	319	0.38
Others	8 (7)	28,970	47	11	58	0.20

No. of unique genes are indicated in brackets. *One auxiliary subunit of a potassium channel.

Table 6

SNP analysis of potassium channel genes/subunits of Dirofilaria immitis.

Putative gene name/description	Coverage in nDi.2.2 scaffold	SNP position in the scaffold	Nucleotide change	Amino acid polymorphism
(irk-3) inward rectifier potassium channel 2	nDi.2.2.scaf0048335 (35621..36678)	35,701	CGA↔GGA	R99G
(twk-48) potassium channel subfamily k member 18-like	nDi.2.2.scaf00031 (282982..285733)	284,309	ATC↔ACC	I171T
(twk-8) t family of potassium channels protein	nDi.2.2.scaf00297 (29943..41201)	31,960	ATG↔GTG	M468Vc
(twk-18) cre-twk-18 protein	nDi.2.2.scaf01557 (1..1936)	1832	AAA↔TAA	K103Stop
		1846	ATG↔ATA	M107I
		1869	CGA↔CAA	R115Q
		1875	TCT↔TTT	S117F
(twk-18)b cre-twk-18 protein	nDi.2.2.scaf01557 (1976..6354)	1976	ATG↔GTG	M1V
		1988	CAT↔GAT	H5D
(twk-47)a,b protein twk-47	nDi.2.2.scaf03721 (1..1197)	1172	GTA↔GCA	V9A
(exp-2) expulsion defective family member	nDi.2.2.scaf00006 (218879..222098)	219,314	GTA↔GGA	V24G
(kvs-4) potassium voltage-gated channel subfamily b member 1	nDi.2.2.scaf00301 (583..5688)	4802	ACA↔ATA	T208I
(shw-3) cre-shw-3 protein	nDi.2.2.scaf05035 (1..436)	428	AGC↔ATC	S3I
(shw-1) voltage-gated potassium channel	nDi.2.2.scaf07713 (1..653)	624	GAC↔AAC	D80N
(F59F3.6) k + channel tetramerisation domain containing protein	nDi.2.2.scaf00117 (106270..113892)	113,664	ATT↔ATG	I194Mc
		113,723	GCT↔GTT	A214V
(unc-93)a potassium channel regulatory protein	nDi.2.2.scaf00192 (61319..67673)	62,629	CGA↔CAA	R131Q
(kcnl-2)a small conductance calcium-activated potassium isoform m	nDi.2.2.scaf01340 (1485..8980)	1522	CAT↔CGT	H346Rc

Homolog genes with detrimental RNAi phenotype in C. elegans.

Absence of a significant BLAST hit (E-value below 10−5) in the predicted proteomes of H. sapiens and C. lupus familiaris.

SNP causes change in secondary structure of the protein as predicted by PSIPRED.

Calcium channels

The nDi.2.2 D. immitis genome contained four α1 subunits, one α2∂ subunit, one β subunit and two auxiliary proteins. Out of 26 gene annotations that could encode 9 unique genes, 321 SNPs were observed with a polymorphic rate of 0.18%. Investigation of these calcium channels for SNP genotype revealed the presence of SNPs in the homologs of phylogenetically defined L-type (‘long-lasting’) gene egl-19, Non – L (P/Q) -type gene unc-2, T-type (‘transient’) gene cca-1 and also in nca-2, an α1U branch of invertebrate cation channel (Hobert, 2013) (Table 7). In gene nca-2, a SNP resulting in either glutamate or lysine at position 1619, could be of importance since its RNAi phenotype is embryonic lethal in C. elegans. A SNP predicted to cause a secondary structure change (Fig. 1) was found to be present at amino acid position 147 of the α2∂ subunit gene unc-36. The SNP V147G was present mainly in the susceptible populations. Multiple SNPs were identified in the Ryanodine receptor (RyR) class of intracellular calcium channels. Among the auxiliary proteins of calcium channels, two SNPs (A141T, F144L) in unc-79, and a T1836M SNP in unc-80 were identified. Orthologs of these genes in C. elegans had RNAi detrimental effects (Jospin et al., 2007).

Table 7

SNP analysis of calcium channel genes/subunits of Dirofilaria immitis.

Putative gene name/Description	Coverage in nDi2.2 scaffold	SNP position in the scaffold	Nucleotide change	Amino acid polymorphism
(egl-19) Ca	nDi.2.2.scaf00156 (87..12852)	11,273	ACC↔CCC	T1447P
(unc-2) P/Q – type calcium channel	nDi.2.2.scaf01292 (5521..9462)	7528	TTG↔GTG	L163V
(cca-1) calcium alpha subunit family member	nDi.2.2.scaf02061 (1..2561)	748	CTG↔CCG	L102P
(cca-1) calcium alpha subunit family member (cca-1)	nDi.2.2.scaf02732 (1..1740)	693	CGA↔GGA	R130G
(nca-2)a four domain-type voltage-gated ion channel alpha-1 subunit	nDi.2.2.scaf00024 (97885..111975)	110,036	GAA↔AAA	E1619K
(unc-36) voltage-dependent calcium channel	nDi.2.2.scaf00220 (32908..46275)	43,927	GTG↔GGG	V147Gb
RyRa ryanodine receptor 44f	nDi.2.2.scaf00024 (338721..388427)	345,969	CAA↔AAA	Q4384K
		371,014	GAA↔GGA	E1711G
		371,999	CGA↔GGA	R1587G
		377,552	ACC↔CCC	T1285P
		385,665	CTT↔CCT	L182P
(unc-80)a uncoordinated family member	nDi.2.2.scaf00207 (43636..72270)	54,434	ACG↔ATG	T1836M
(unc-79)a uncoordinated family member	nDi.2.2.scaf01338 (32..946)	923	GCA↔ACA	A141T
		932	TTC↔CTC	F144L
(unc-79)a protein unc-79 homolog	nDi.2.2.scaf04898 (1..1130)	1125	CAA↔CAT	Q2H

Homolog genes with detrimental RNAi phenotype in C. elegans.

SNP causes change in secondary structure of the protein as predicted by PSIPRED.

Chloride channels

Nine membrane localized chloride channels were identified; five from the chloride channel (CLC) family, two from the chloride intracellular channel (CLIC) family and two from the bestrophin-related channel types. Two hundred and sixty five SNPs were identified among these 9 chloride channel genes with a polymorphic rate of 0.48%; highest among the channel types studied. In the CLC-type clh-1 gene homolog, two closely positioned SNPs A300G and V302F were identified. In the other CLC-type gene clh-3, we found two SNPs at positions 915 (Q915K) and 1095 (S1095N). A single SNP N325H was identified in the clh-5 gene homolog. The clh-6 gene was multi-polymorphic as we identified 9 missense types of polymorphism in this gene. Bestrophin-related calcium-activated chloride channel genes, namely best-13 and best-24 were found to be polymorphic (see Table 8). One of the subunits of the best-24 gene had a premature stop codon causing polymorphism, whereas the SNP L75M, found in the other subunit, caused a change in predicted secondary structure.

Table 8

SNP analysis of chloride channel genes/subunits of Dirofilaria immitis.

Putative gene name/Description	Coverage in nDi.2.2 scaffold	SNP position in the scaffold	Nucleotide change	Amino acid polymorphism
(clh-1) protein clh- isoform b	nDi.2.2.scaf01639 (1..5958)	5052	GCT↔GGT	A300G
		5057	GTT↔TTT	V302F
(clh-3/clc-3) voltage gated chloride channel family protein	nDi.2.2.scaf00125 (31151..45497)	31,476	AGT↔AAT	S1095N
		32,409	CAA↔AAA	Q915K
(clh-5/clc-5) chloride channel protein 3	nDi.2.2.scaf00353 (9838..16713)	14,624	AAT↔CAT	N325H
(clh-6/clc-6) chloride channel protein 7 (ccp-7)	nDi.2.2.scaf00051 (224832..229977)	225,891	GTT↔ATT	V628I
		226,548	GTG↔GCG	V549Aa
		226,820	ATT↔GTT	I517V
		226,841	TTT↔GTT	F510V
		227,468	AAT↔GAT	N413D
		228,310	GCG↔GGG	A278G
		228,313	GCT↔GGT	A277G
		228,316	GTG↔GGG	V276G
		228,981	GGT↔AGT	G147S
(best-13) bestrophin 1	nDi.2.2.scaf00401 (31248..33870)	32,314	GCT↔GGT	A136G
		33,830	CGA↔CAA	R384Q
(best-24) bestrophin family protein	nDi.2.2.scaf00816 (5183..11705)	5254	GAA↔AAA	E626K
		9709	CAA↔TAA	Q347Stop
(best-24) bestrophin family protein	nDi.2.2.scaf01764 (346..4516)	852	TTG↔ATG	L75Mc

SNP causes change in secondary structure of the protein as predicted by PSIPRED.

Other type of channels

From the DEG/ENaC (DEGenerin/Epithelial Sodium Channels) protein family of D. immitis, we identified 8 annotations that could encode 7 genes. The calculated polymorphic rate within this group of sodium channels was 0.2%. One SNP in each of the amiloride-sensitive sodium channel genes, namely unc-105, acd-1, acd-5 and del-10, was identified (see Supplementary file 4, Table S4). Double SNPs, F33L and D148N were seen in del-10. The TRP (Transient Receptor Potential) superfamily of cation channels were 7 in number in the D. immitis genome. These canonical TRP channels carried 319 SNPs in both introns and exons, with a polymorphism rate of 0.38%. Both the first and second SNPs were found in codon 214 of the gene spe-41, changing alanine to valine or serine. In the trp-4 gene, 3 SNPs, T83P, A538T and S894A, were identified; the first predicted to cause a secondary structure change from helix to coil. Another TRP channel gene, ced-11, was found to be highly polymorphic with 15 variable loci identified in the exon regions alone. So far, the genome of D. immitis was found to have three cation-selective cyclic nucleotide gated (CNG) ion channel genes, namely tax-2, tax-4 and che-6, that respectively encode β, α, and α/β subunits (Smith et al., 2013). A single polymorphic site T106I was found in the β-type subunit gene, tax-2, whereas the other subunit gene, che-6, was also highly polymorphic (see Supplementary file 3, Table S3). The gene ncs-4 (“neuronal calcium sensor”) belonging to the calcium binding protein family had a SNP at amino acid position E169K. Also, a single SNP, V124I was found in the cation channel, and 2 SNPs, F30L and S66G, were observed in an unnamed voltage-dependent anion channel. The presence of these amino acid changes were predicted to cause secondary structure changes in the protein.

Discussion

Parasitic nematode associated diseases cause serious health issues in millions of humans (Hotez et al., 2007, Hotez et al., 2008, Brooker et al., 2010, Lustigman et al., 2012) and animals (both domestic and livestock) each year (Bird and Kaloshian, 2003, McKellar and Jackson, 2004, Charlier et al., 2014). MLs have been the drugs of choice to treat and/or to prevent many parasitic nematode infections, including heartworm in dogs (Geary, 2005, Omura, 2008, Prichard et al., 2012). In recent times, many of the existing anthelmintics, including the MLs, have begun to face problems of emerging resistance. Therefore, improved diagnostics, new drugs and effective vaccines are important goals for efficient control and prevention of parasitic diseases. For the animal health sector, there has been an increased focus on the discovery of new drugs for canine heartworm (Geary et al., 2015). In the context of discovering new anthelmintics, it is of interest that ion channels in the neuromuscular systems of nematodes have been successful drug targets for many of the anthelmintics discovered in the past, and remain targets of choice for discovery of new anthelmintics (Wolstenholme, 2011, Greenberg, 2014).

In this study, we employed comparative genomic approaches, by using ion channel gene sequences from related nematodes, to localize their homologs in the nDi.2.2 genome of D. immitis. A total of 224 genes/subunits that are likely to encode 126 unique ion channels and receptors genes were identified. The genome of D. immitis has 42 cys-loop LGIC (15 cationic nAChR, 1 cationic iGluR, 19 anionic chloride channels and 7 diverse channel genes). This is similar to the 44 receptor genes (21 nAChR-like genes and 23 chloride channel subunits) reported in the closely related B. malayi genome (Scott and Ghedin, 2009), but fewer compared to 102 LGIC genes (52 nAChR, and 50 chloride channel genes) in C. elegans (Jones and Sattelle, 2004). We report that the genome of D. immitis may have at least 43 voltage-gated potassium channels; close in number to the 36 in B. malayi (Scott and Ghedin, 2009) but fewer than the 72 reported in C. elegans (Hobert, 2013) and 80 in mammals (Pardo and Stühmer, 2014). Six subunits encoding calcium channel genes (also 2 auxiliary subunits) were identified in D. immitis, with four encoding α1 subunits, one each encoding α2∂, β subunits. In C. elegans, 5 genes encode α1, 2 β, and 2 α2∂ subunits (Hobert, 2013). Among the chloride channels, 5 were identified from the CLC type, 2 CLICs and 2 genes from the bestrophin-related genes, compared to 6 genes of the CLC and 26 bestrophin-related genes found in C. elegans (Schriever et al., 1999, Hobert, 2013). We also identified 7 genes that belong to DEG/ENaC/ASIC channels and another 7 belonging to the TRP-type channels. This compares with 30 and 23 genes, respectively, in C. elegans (Hobert, 2013).

We aimed at investigating SNPs in these ion channel genes that cover 1.1% of ∼84.2 Mb nuclear assembly of D. immitis. Within these ion channel genes, we report an overall polymorphic rate of 0.18%. For each of the polymorphisms, especially for those in the exonic regions, we studied its likely effects, such as a change in amino acid, and change in predicted secondary structure of the protein as it could be relevant for anthelmintic drug development. Furthermore, this study identified all the allelic forms in these ion channels. This information may be used to help ensure that a prospective drug candidate is active against different genetic forms. As we also used genomic data from pooled LOE populations for genotyping, the populations’ specific genetic variants were also identified. This information may be helpful to check the effectiveness of a new drug, even against ML LOE populations. Moreover, the polymorphic knowledge of ion channel targets may help to anticipate resistance development, as the more heterogeneous a drug target is, the greater the potential for resistance to develop (Prichard, 2001).

We analyzed polymorphic rates per geographical location and found rates of 75.2, 70.9, 70.0 and 43.0% in pooled populations from the USA, Grand Canary, Grenada and Italy respectively. A low polymorphic rate in the sample of worms from Italy may be explained by the fact that the worms were from a single dog, in contrast to worms coming from other countries in which several dogs donated worms. Thus, barring the Italian samples, we saw a trend of low genetic variability among ion channel and receptor genes between populations from the USA, Grand Canary and Grenada. This hypo-variability, though not anticipated, is in agreement with the population genetic results studied at the microsatellite level (Belanger et al., 2011) and also at the whole genome level (Godel et al., 2012). However, Godel et al. (2012) used only two different D. immitis isolates, one from Italy and the other from the USA to draw this conclusion. Our findings could support the hypothesis that heartworm disease is a New World disease (Bowman and Atkins, 2009). However, canine heartworm was reported for the first time in Italy in 1626 (Birago, 1626) but is known to have been present in the USA from 1847 (Osborne, 1847). On the other hand, the result of our study is in contrast to reports suggesting that D. immitis has considerable genetic variability (Geary et al., 2011, Bourguinat et al., 2011c). These earlier studies were based on only 5 gene sequences in individual adult worms from the USA (field and laboratory samples) and Japan (field samples). It must also be kept in mind that in the current study we did not analyze polymorphism when it occurred in less than 15% of the reads of a sequence. Thus, additional polymorphism may occur at low levels. Other parasitic nematodes have marked genetic variability such as Ascaris lumbricoides (Criscione et al., 2007), Haemonchus contortus (Cerutti et al., 2010), Telardorsagia circumcincta and Mazamostrongylus odocoilei (Grillo et al., 2007). However, these studies were conducted on whole genomes and some genes may be more conserved than others. Therefore, it is not appropriate to compare genetic variability in a subset of genes in one organism with different genes in other organisms. D. immitis appears to be a nematode with a relatively low level of heterogeneity in ion channel genes. Genome wide studies need to be done on worms from different locations to make definitive conclusions on the overall degree of heterogeneity in D. immitis.

Genetic changes either in the form of insertions/deletions (indels) or mutations in the intron region of a gene can influence splicing (Kubota et al., 2011), structure and function of a protein, and the possibility of resistance selection (Barrère et al., 2014). For example, an indel of 63 bp present in intron 2 of the H. contortus gene, Hco-acr-8 was found to be responsible for the truncated isoform of the protein Hco-ACR-8b and this genetic change correlated significantly with levamisole resistance (Barrère et al., 2014). Also, an intronic SNP out of the diplotypic GG-GG genotype in Dim-pgp-11 of D. immitis (Bourguinat et al., 2016) occurs at the intron region of the 3′ end of the gene. This diplotypic genotype strongly correlated with in vitro (Bourguinat et al., 2011b), and in vivo responses of mf to IVM (Bourguinat et al., 2011a). Accordingly, a few of the 1508 SNPs (82% of total polymorphism) identified within intron regions of the 244 genes/subunits of ion channel genes could have implications for gene expression, response to drugs (Wang et al., 2011), resistance development (Barrère et al., 2014), or other effects, which would need to be confirmed in future investigations. Also, 125 SNPs found within exon regions caused synonymous type variation. There is therefore a possibility that a few of these SNPs might have implications either for aberrant mRNA splicing (Cartegni et al., 2002) or affect mRNA stability that could impact protein expression (Nackley et al., 2006) or change protein conformation. Such effects could lead to changes in protein function (Kimchi-Sarfaty et al., 2007).

This study also identified 129 SNPs that are in exon regions and caused amino acid changes. A few of the nsSNPs may be critical as they were found in ion channels identified as potential drug targets based on two criteria: the detrimental effect of knockout of its homolog gene in C. elegans and absence of a homolog in humans and dogs. The significance of the nsSNPs was analyzed on the basis of changes in predicted secondary structure in association with a change in amino acid.

Glutamate gated chloride channels (GluCls) are proven targets for macrocyclic lactones. GluCls are invertebrate specific and play a key role in locomotion, feeding and sensory input (Greenberg, 2014), and therefore are attractive drug targets. AVR-14, in particular, is an interesting target as it is well conserved in all nematodes studied so far (Laughton et al., 1997, Jagannathan et al., 1999, Yates and Wolstenholme, 2004, Njue et al., 2004, Tandon et al., 2006). GluCls are expressed on motor neurons (Dent et al., 2000, Portillo et al., 2003, Glendinning et al., 2011) that mediate locomotion, feeding, reproduction and secretion/excretion from the excretory pore. Considering this, the AVR-14 of D. immitis (GenBank: CAE46429.1) could act as a drug target; two adjacent SNPs, C382F and C385F, identified in this study could be critically important from a drug design perspective. The homolog of this gene in C. elegans has a locomotion coordination variant phenotype following RNAi (Cook et al., 2006).

Nicotinic acetylcholine receptors (nAChRs), the most common targets of current anthelmintics, are a diverse group of receptors with interesting pharmacology (Wolstenholme, 2011, Greenberg, 2014). nAChRs are pentameric structures with a wide variety of potential subunits that provide them with considerable receptor diversity, and distinct structural properties and pharmacological sensitivities (Greenberg, 2014). Although the possibility of rapid development of resistance, and cross-resistance between drugs that target these receptors, has been anticipated as nematodes may alter receptor sensitivity to different anthelmintics by varying the stoichiometry of subunits (Williamson et al., 2009, Buxton et al., 2014), it was shown that this possibility may not be true at least in the case of neuronal nAChRs targeted by monepantel (Kaminsky et al., 2008). nAChRs get expressed at the neuromuscular junction, nerve ring and in the pharynx of C. elegans (Jones and Sattelle, 2004). Therefore, nematode nAChRs, including those in D. immitis, may be good targets for new anthelmintics. Among the polymorphisms identified in the putative nAChRs genes of D. immitis, SNPs identified in acr-8, acr-11 and acr-16 genes need special mention. Three SNPs (L256F, Y264F, R364G), including a secondary structure change caused by R364G in acr-8, and a deletion-causing SNP in the other subunit of the acr-8 gene, were identified. The acr-11 gene has a detrimental RNAi phenotype in C. elegans and its homolog in D. immitis has been identified with a SNP Q46P. The other nAChR gene, acr-16 was found to have three close SNP loci, at P40S, L43F and S53L.

Inhibitory GABA-gated chloride channels present at the neuromuscular junction of nematodes (Holden-Dye et al., 1989, Richmond and Jorgensen, 1999) mediate the relaxation phase of sinusoidal muscle movement (Accardi et al., 2012) and so any drug-induced activation, for example by the anthelmintic piperazine, a GABA agonist, can cause flaccid paralysis of worms (Martin, 1985). MLs also appear to bind to GABA-gated chloride channels of nematodes (Feng et al., 2002, Brown et al., 2012). Therefore, these channels can be fruitful drug targets. GABA receptors of D. immitis, namely the homologs of gab-1 and unc-49 with identified SNPs I20T and N57D, respectively, may be of interest for drug intervention. The gene unc-49 (and its associated SNP) is worthy of mention as its homolog gene in C. elegans has proven expression in somatic muscle at the neuromuscular junction and therefore plays an important role in locomotion (Bamber et al., 1999). Nematodes possess serotonin, dopamine, tyramine and ACh-gated anion channels not found in mammals (Wolstenholme, 2011, Beech et al., 2013, MacDonald et al., 2014). Therefore, the SNPs in the tyramine-gated receptor coded by lgc-55, and the dopamine-gated receptor coded by lgc-53 could be of interest during drug development. iGluR gene glr-1 of D. immitis, with its SNPs N287S, P809T, could be an interesting drug target, as such receptors have been found to play vital roles in the nematode nervous system (Aronoff et al., 2004, Kano et al., 2008).

Voltage gated ion channels function in response to changes in membrane potential and gate calcium (Cav channels) or potassium (Kv channels) in nematodes. This type of voltage-gated ion channel represents a class of outstanding but underexploited drug targets (Camerino et al., 2007, Davies et al., 2007). To date, only two members of the voltage-gated ion channels, namely Ca2+- activated potassium channel gene slo-1 and a schistosome calcium channel gene have been implicated in the action of anthelmintics: emodepside and praziquantel, respectively (Wolstenholme, 2011). The genome of D. immitis encodes 43 potassium, 6 calcium channel genes but no voltage-gated sodium channels. The potassium channel genes, twk-18 and twk-47 were found to be nematode specific targets and our study has identified SNPs in each of the genes. Further study may be warranted to understand the implications of the identified SNPs, M1V, H5D (twk-18) and V9A (twk-47), on the functionality of the respective proteins. A SNP (H346R), which may be of importance, was found in the small conductance Ca2+ -activated K+ (SK) channel gene kcnl-2, wherein the homolog of this gene is required for regulating the rate of egg-laying in C. elegans (Chotoo et al., 2013). The α1 subunit encoding calcium channel genes, such as nca-2 and egl-19 with SNPs (E1619K and T1447P, respectively) could be good targets for future anthelmintic discovery, considering the fact that their respective homolog genes in C. elegans were found to have RNAi detrimental phenotypes (Kamath et al., 2003, Humphrey et al., 2007). An intracellular Ryanodine receptor (RyR) calcium channel encoded by the unc-68 gene is well conserved in nematodes (Maryon et al., 1996) and is localized in muscles (Hobert, 2013) and in neurons (Liu et al., 2005) of C. elegans. Also, the phthalic acid diamide and anthranilic diamide classes of insecticides act on the RyR receptor (Sattelle et al., 2008). Considering this, the RyR receptor homolog of D. immitis could be an attractive drug target. However, this gene was found to be highly polymorphic (both intronic and exonic SNPs) and so the chances of resistance being selected against drugs that target this gene product could be high. In D. immitis, CLC type chloride channel genes, namely clh-1, -3, -5 and -6, which potentially control the membrane potential of cells, were found to have polymorphic site(s). However, the significance of these SNPs has yet to be studied. There were also two bestrophin-related chloride channel genes identified in D. immitis; the best-24 gene homolog is of interest since the gene is expressed in neurons in C. elegans (Hobert, 2013). For the same reason, SNPs identified in the gene, mainly the one causing a stop codon and the other causing a change in predicted secondary structure, may have interesting effects.

Nematodes possess several other families of ion channels in the nervous system that could also be exploited as drug targets. For examples, the DEG/ENaC (Degenerin/Epithelial Sodium Channel) family of sodium channels have been implicated in mechanotransduction and mechanosensitive behaviour in C. elegans (Syntichaki and Tavernarakis, 2004). Activating this type of sodium channel in the worms was found to cause inappropriate cell death (Driscoll, 1992) due to their constitutive activation (Hong and Driscoll, 1994), and so any agonist drug could be a potential anthelmintic (Wolstenholme, 2011). From this study, SNPs were also found in DEG/ENaC channel genes such as unc-105, acd-1, acd-5 and the worm-specific gene del-10. TRP channels, a superfamily of cationic channels, well represented in both nematodes and schistosomes (Wolstenholme et al., 2011), mediate transduction of sensory stimuli and are important in Ca2+ signaling cascades. Any dysregulation of these channels could interfere with signal transduction and disrupt Ca2+ homeostasis in worms (Greenberg, 2014). SNPs identified in TRP channel genes such as in trp-4, ced-11 and spe-41 may, therefore, have interesting effects.

Finally, the genome of D. immitis codes for three predicted CNG channels compared to six in C. elegans (Kaupp and Seifert, 2002). A SNP, T106I, was identified in the gene tax-2, which could be a potential drug target since its homolog gene regulates thermosensation, chemosensation and neuronal development in C. elegans (Coburn and Bargmann, 1996, Komatsu et al., 1996, Coburn et al., 1998).

SNPs present either in only susceptible or LOE populations were also highlighted in this study. For example, SNPs such as I20T in the gene gab-1, L691P (lgc-53), P809T (glr-1), I717V (mgl-2), A381T (acd-5), V147G (unc-36), L75M (best-24), Q2H (unc-79), to name a few, were specific to susceptible populations. SNPs such as N287S in gene glr-1, S53L in acr-16, T208I in kvs-4, S1095N in clh-3, A538T in trp-4, H346R in kcnl-2 and I36M in che-6 were identified only in the LOE populations. The purpose of highlighting possible phenotypic consequences of SNP variability in populations is to provide a basis to ensure that any new drug is active against all of the allelic variants of the drug target.

Conclusion

Our study describes the extensive profile of ion channel and receptor genes in D. immitis. This study is the first of its kind to determine genetic polymorphism in ion channels in this parasite, to generate information on possible effects on protein structure and function, and to consider these findings in the context of RNAi phenotypes for homologs of the genes in C. elegans. Such information may be useful during drug design and to anticipate polymorphisms which could impact resistance development. Among the total of 1762 SNPs identified, some may affect gene expression, structure and function of the proteins, and resistance selection processes. A merit of this study lies in the large number of samples used for SNP genotyping. Whole genome sequencing data were generated, from 122 ML susceptible adult worms, isolated from 17 dogs from 4 countries, for SNP analysis. In the case of LOE samples, ∼32,000 mf, isolated from 4 ML-LOE dogs from different US states were used (Bourguinat et al., 2015). It is possible that some additional SNPs were not detected in this analysis due to lack of coverage of some sequences in the genome or to the population size and diversity that was analyzed. Moreover, the SNPs that were specific to the LOE populations were called based on a relatively low number of reads, compared to the read frequencies among the susceptible populations, and therefore information specific to the LOE samples should be considered preliminary. Further studies need to be done to confirm the SNPs in field samples and to understand the possible implications of such SNPs in terms of protein structure, function and also any possible interaction of ion channels with potential antiparasitic drugs.