Polymorphism in ABC transporter genes of Dirofilaria immitis.

Dirofilaria immitis, a filarial nematode, causes dirofilariasis in dogs, cats and occasionally in humans. Prevention of the disease has been mainly by monthly use of the macrocyclic lactone (ML) endectocides during the mosquito transmission season. Recently, ML resistance has been confirmed in D. immitis and therefore, there is a need to find new classes of anthelmintics. One of the mechanisms associated with ML resistance in nematodes has been the possible role of ATP binding cassette (ABC) transporters in reducing drug concentrations at receptor sites. ABC transporters, mainly from sub-families B, C and G, may contribute to multidrug resistance (MDR) by active efflux of drugs out of the cell. Gene products of ABC transporters may thus serve as the targets for agents that may modulate susceptibility to drugs, by inhibiting drug transport. ABC transporters are believed to be involved in a variety of physiological functions critical to the parasite, such as sterol transport, and therefore may also serve as the target for drugs that can act as anthelmintics on their own. Knowledge of polymorphism in these ABC transporter genes in nematode parasites could provide useful information for the process of drug design. We have identified 15 ABC transporter genes from sub-families A, B, C and G, in D. immitis, by comparative genomic approaches and analyzed them for polymorphism. Whole genome sequencing data from four ML susceptible (SUS) and four loss of efficacy (LOE) pooled populations were used for single nucleotide polymorphism (SNP) genotyping. Out of 231 SNPs identified in those 15 ABC transporter genes, 89 and 75 of them were specific to the SUS or LOE populations, respectively. A few of the SNPs identified may affect gene expression, protein function, substrate specificity or resistance development and may be useful for transporter inhibitor/anthelmintic drug design, or in order to anticipate resistance development.

Introduction

Dirofilaria immitis, a filaroid nematode that causes dirofilariasis or heartworm disease, is arguably the most important parasite that infects dogs in North America (Wolstenholme et al., 2015). This mosquito-borne nematode affects 30 other mammal species including cats, wild canids, felids and occasionally humans (Lee et al., 2010, Genchi et al., 2011). The disease has worldwide distribution but is more prevalent in temperate, tropical and sub-tropical regions of the world (Simón et al., 2009). For over 2 decades, prevention of the disease has relied upon a single class of drug, the macrocyclic lactones (MLs). MLs were first used as heartworm preventives in 1987. Concern about possible loss of efficacy (LOE) to this class of drug was reported in 2005 (Hampshire, 2005). Recent studies have unequivocally confirmed ML resistance in D. immitis (Bourguinat et al., 2011a, Bourguinat et al., 2015, Pulaski et al., 2014). The ATP binding cassette (ABC) transporters have been implicated in ML resistance (see review Lespine et al., 2012). Specific studies on the role of ABC transporters and ML resistance include studies in Caenorhabditis elegans (James and Davey, 2009, Ardelli and Prichard, 2013, Bygarski et al., 2014), Haemonchus contortus (Blackhall et al., 1998; Xu et al., 1998; Sangster et al., 1999, Prichard and Roulet, 2007), Onchocerca volvulus (Eng and Prichard, 2005, Ardelli and Prichard, 2007, Bourguinat et al., 2008, Osei-Atweneboana et al., 2011, Nana-Djeunga et al., 2014), Teladorsagia circumcincta (Dicker et al., 2011), and for benzimidazole resistance in H. contortus (Blackhall et al., 2008), triclabendazole resistance in Fasciola hepatica (Mottier et al., 2006) and praziquantel resistance in Schistosoma mansoni (Messerli et al., 2009). In D. immitis, a diplotypic genotype in Dim-pgp-11 (ABC-B class of the transporter superfamily) gene was found to correlate strongly with a loss of response to ivermectin (Bourguinat et al., 2011a, Bourguinat et al., 2011b).

Members of the ABC transporter superfamily, particularly P-glycoproteins (PGPs) and multidrug resistance proteins (MRPs) have broad substrate specificity (e.g., antibiotics, antimalarials, herbicides, antifungals and other chemotherapeutic agents) (Higgins, 2007). ABC transporters, as efflux pumps, influence the bioavailability and disposition of drugs through active efflux of compound out of cells, thereby contributing to a phenomenon called multidrug resistance (MDR) (Lespine et al., 2012, Ardelli, 2013). Overexpression of these transporters has been suggested as a MDR conferring mechanism (Ardelli, 2013). In nematodes, multidrug resistance is suggested to be conferred by members from subfamily B and C, though many reports specifically implicate the role of PGPs in the resistance mechanism (Prichard et al., 2012 Lespine et al., 2012). A study in B. malayi has shown that transporters from subfamily A and G may also have a role in resistance, as they were also overexpressed along with PGPs and MRPs following treatment with ivermectin and moxidectin (Stitt et al., 2011, Tompkins et al., 2011).

Although ABC transporters may play a role in developing and maintaining anthelmintic resistance, it is believed that their main physiological function is to protect neurons and other tissues in which they are expressed, from a broad spectrum of toxins (Prichard and Roulet, 2007, Ardelli, 2013). Unlike vertebrates, which have a vascular system contained within epithelial cells which line the blood capillaries and other blood vessels where PGPs and MRPs are commonly expressed, nematodes do not have a similar barrier cell layer expressing ABC transporters. The number of efflux transporters is greater in nematodes, compared to mammals, as they may have evolved this diversity as part of their protective mechanism for a variety of different tissues and specific cells. ABC transporters can regulate anthelmintic efficacy by modulating the intracellular or intra cell membrane concentration of lipophilic drugs. Therefore, a deletion or inhibition of such transporters may partially overcome resistance (Lespine et al., 2012, Ardelli, 2013, Greenberg, 2013a, Greenberg, 2014a). For example, deletion or disruption of mammalian and nematode Pgp genes leads to increased sensitivity to ivermectin (Schinkel et al., 1994, Janssen et al., 2013). Also, substrates of PGP that are competitive inhibitors (e.g., verapamil), or that block the transport function directly, have restored anthelmintic sensitivity (Ardelli and Prichard, 2013). In B. malayi, ABC transporter inhibitors have potentiated sensitivity to ivermectin (Tompkins et al., 2011). Similar reversal mechanisms have been shown for H. contortus resistant to ivermectin (Molento and Prichard, 1999, Lifschitz et al., 2010, Bartley et al., 2012). Inhibitors of MRP, such as MK571 and buthionine sulfoxamine, have also reversed resistance (Prichard et al., 2012). In nematodes, ABC transporters are believed to be essential for the survival of nematodes as they are known to be involved in a wide range of processes. For example, in C. elegans, they play critical roles in apoptotic cell corpse removal (Wu and Horvitz, 1998), dauer formation (Yabe et al., 2005), RNAi (Sundaram et al., 2006), directed sperm motility (Kubagawa et al., 2006), apart from resistance to toxins (Broeks et al., 1995) and heavy metals (Broeks et al., 1996, Vatamaniuk et al., 2005), and resistance to pathogens (Mahajan-Miklos et al., 1999). Taking all of this into consideration, ABC transporter genes may prove to be attractive therapeutic targets for new or repurposed MDR reversing agents, and/or for new anthelmintics (Lespine et al., 2012, Ardelli, 2013, Greenberg, 2014b).

Knowledge of allelic variants in ABC transporter genes may be useful information to have during drug design. Knowledge about polymorphism in these genes may help ensure that a new drug will be active against all the allelic forms of these targets. Some polymorphism in MDR transporter genes may influence their level of expression and/or the functional characteristics (Gerloff, 2004), and such knowledge may be useful to understand the mechanism of resistance. Also, heterogeneity of a drug target may determine the drug selection process (Prichard, 2001) and knowledge of polymorphism in drug transporters may help to anticipate resistance development. In order to address these possibilities, the objective of this study was to identify putative ABC transporter genes/subunits from subfamilies A, B, C and G in D. immitis, by comparative genomic approaches, and to analyze the genes for single nucleotide polymorphisms (SNPs).

Materials and methods

Identification of all of the ABC-A, -B, -C and -G transporter genes of D. immitis

The genome of D. immitis is not fully annotated (Bourguinat et al., 2015) and therefore a complete inventory of ABC transporter genes is not yet available. In a recent study (Bourguinat et al., 2016), eight ABC transporter genes including a pseudogene from D. immitis were described. To identify the remaining, un-annotated homologs of ABC transporter genes in D. immitis, we followed a similar approach as previously described (Mani et al., 2016). Briefly, all the nucleotide sequences that encode genes or subunits of ABC transporter subfamilies (A, B, C and G) 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 used in the nucleotide BLASTt server v2.2 (http://nematodes.org/genomes/dirofilaria_immitis/) to indicate putative protein type. These BLAST searches were helpful to locate each of the putative ABC transporter 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 from whole genome sequencing data

Whole genome data from diverse pooled samples of 122 adult D. immitis, from 17 dogs from the USA, Grand Canary, Grenada and Italy, and approximately 32,000 D. immitis microfilariae from 4 dogs from USA, described previously (Bourguinat et al., 2015, Mani et al., 2016) were used. More information on the D. immitis samples and their classification into ML susceptible (SUS) and loss of efficacy (LOE) populations is given in the Supplementary table ST1. The adult worms were mixed sex populations characterized as susceptible to macrocyclic lactone heartworm preventives, while the microfilariae were termed loss of efficacy samples as it was suspected that the dogs may have become infected despite being on macrocyclic lactone heartworm preventives, as previously described (Bourguinat et al., 2015). Some of the adult female worms were gravid. This would only enhance the possibility for detecting polymorphism in the samples. The method followed for synchronized file generation has been described elsewhere (Mani et al., 2016). Paired-end reads were assembled, filtered and trimmed. D. immitis nuclear genome v2.2 (http://nematodes.org/genomes/dirofilaria_immitis/) was used as the reference genome. BAM files that consisted of aligned sequences from each population against the reference genome were generated. The program Popoolation2 was used to generate a synchronized file that contained nucleotide frequencies for each population at every base in the reference genome. This information was obtained for each gene in a concise format after filtering for base quality.

Nomenclature and classification of ABC transporter genes identified

Previously 8 ABC-B transporter genes in D. immitis had been identified, named and classified based on phylogenetic analysis (Bourguinat et al., 2016). For the rest, the naming was either based on the closest homolog in C. elegans, identified during a BLASTP search, or as given during annotation of the D. immitis nDi.2.2.2 nuclear genome. Fifteen ABC transporter genes/subunits of D. immitis were classified, based on the classifications of ABC transporter homologs in C. elegans (Sheps et al., 2004) and B. malayi (Ardelli et al., 2010).

Identification of SNP within ABC transporter genes

A similar methodology, as described by Mani et al. (2016) was followed to identify the location of a SNP within a GFF (Generic Feature Format) file that contained annotations of the nuclear genome of D. immitis nDi.2.2.2 was used to call whether a locus was in an exonic or intronic region of the gene. The type of SNP, whether synonymous (sSNPs) or non-synonymous (nsSNPs) (missense or nonsense), and the respective amino acid changes were assessed using Sequencher software 4.10.1 (Gene Codes Corporation, Ann Arbor, MI, USA). Transmembrane helices of these membrane proteins were predicted using online web tools such as Phobius (http://phobius.sbc.su.se/) and HMMTOP (http://www.enzim.hu/hmmtop/index.php). For those of the ABC transporter genes whose transmembrane domain and ATP-binding cassette domain were previously predicted (Bourguinat et al., 2016), the topological structures were displayed using the TOPO2 transmembrane protein display software (http://www.sacs.ucsf.edu/cgi-bin/open-topo2.py). Such topological models incorporated every SNP identified in this study.

Focus on polymorphism in potential drug targets

ABC transporters with potential drug target capabilities were highlighted, based on similar filtering methodology used for B. malayi (Kumar et al., 2007) and D. immitis (Godel et al., 2012, Mani et al., 2016). One such filter was evidence of detrimental effects of its homolog during RNAi in C. elegans. Such detrimental effects could be, for example, embryonic lethal/larval arrest, shortened life span, locomotion variant, organism development variant, sluggish/fainter, slow growth, egg size defective, pharyngeal pumping variant, as known from the Wormbase website (http://legacy.wormbase.org/). The other type of filter used was the absence of a BLASTP hit with an E-value below 10−5, in the predicted proteomes of Homo sapiens and Canis lupus familiaris. To study the impact of identified SNPs, especially non-synonymous SNP (nsSNPs), on the secondary structure of an ABC transporter, a secondary structure prediction algorithm called PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred/) (McGuffin et al., 2000, Mani et al., 2016) was used.

Results

ABC transporter genes/subunits of all nematodes extracted from databases

A total of 169 nucleotide sequences that could encode ABC transporter genes or subunits from subfamilies A, B, C and G were extracted from available databases such as NCBI, NEMBASE4, Wormbase and Broad Institute. The number of nucleotide sequences, of related nematodes, extracted for analysis is given in Table 1 (complete information of sequences is available in Supplementary data file S1). Of the total nucleotide sequences extracted, 13 of them belonged to ABC-A, 84 to ABC-B, 29 to ABC-C and 43 to ABC-G encoding genes/subunits.

Table 1

Number of nucleotide sequences of ABC transporter genes/subunits collected per nematode from NCBI, Wormbase, NEMBASE4 and Broad Institute.

Nematodes	Number of gene/subunit sequences extracted for polymorphism analysis
Brugia malayi	25
Loa loa	20
Onchocerca volvulus	6
Wuchereria bancrofti	10
Caenorhabditis brenneri	5
Caenorhabditis briggsae	23
Caenorhabditis elegans	40
Caenorhabditis remanei	23
Haemonchus contortus	9
Others (Ascaris suum, Onchocerca ochengi, Toxocara canis, Strongyloites ratti, Dirofilaria immitis, etc.)	8
Total	169

Number of ABC transporter genes identified in the D. immitis genome

By the complementary genomic approaches, we have identified 17 genes or subunits that encode 15 unique ABC transporter genes (see Table 2). Out of the 15, eight ABC transporter genes have already been described (Bourguinat et al., 2016). These genes, identified on the basis of PCR amplification, sequencing and bioinformatics tools, included three PGP encoding genes (Dim-pgp-3, Dim-pgp-10, Dim-pgp-11), a pseudogene and two half-transporter genes (Dim-haf-1, Dim-haf-4) from the ABC-B subfamily, and two half transporter genes (Dim-haf-5.1, Dim-haf-5.2) from the ABC-C subfamily (Bourguinat et al., 2016). We identified 7 more ABC transporter genes, two (Dim-abt-2, Dim-abt-4) from sub-family A, three (Dim-mrp-1, Dim-mrp-5, Dim-mrp-7) from sub-family C and two (Dim-wht-4, Dim-wht-8) ABC transporter genes from sub-family G.

Table 2

ABC transporter genes identified in the genome of Dirofilaria immitis.

ABC transporter class	Putative gene name	Location on nDi.2.2 genome
A	Dim-abt-2	nDi.2.2.scaf04723:166..1096nDi.2.2.scaf05671:1..844
	Dim-abt-4	nDi.2.2.scaf01473:143..1713nDi.2.2.scaf01473:3942..7078
B Full transporters	Dim-pgp-3*	nDi.2.2.scaf00046:226523..237100
	Dim-pgp-10*	nDi.2.2.scaf00048:63901..79265
	Dim-pgp-11*	nDi.2.2.scaf00004:79181..88673
	Pseudogene*	nDi.2.2.scaf00041: 17947..27954
Half transporters	Dim-haf-1*	nDi.2.2.scaf00023:388438..395064
	Dim-haf-4*	nDi.2.2.scaf00101:45240..51224
C Full transporters	Dim-mrp-1	nDi.2.2.scaf00010:187828..203859
	Dim-mrp-5	nDi.2.2.scaf00052:175763..187462
	Dim-mrp-7	nDi.2.2.scaf00004:629508..641785
Half transporters	Dim-haf-5.1*	nDi.2.2.scaf00496.1:20854..25667
	Dim-haf-5.2*	nDi.2.2.scaf00496.2:20854..25667
G	Dim-wht-4	nDi.2.2.scaf00589:2368..14928
	Dim-wht-8	nDi.2.2.scaf00252:36959..42309

The table shows newly identified ABC transporter genes from subfamilies ABC-A, ABC-B, ABC-C and ABC-G, updated from previously reported numbers (indicated by *, Bourguinat et al., 2016).

Some of the ABC transporter genes identified by Bourguinat et al. (2016) were different from those of their respective annotations in the D. immitis genome (http://nematodes.org/genomes/dirofilaria_immitis/). For example, Dim-pgp-3, Dim-pgp-10, Dim-haf-1, Dim-haf-4 and Dim-haf-5 were differently annotated in the genome. In such cases, every polymorphic position within those genes has been called, based on the annotations by Bourguinat et al. (2016), as the authors had used bioinformatic tools and PCR amplification to confirm each gene. For the rest of the ABC transporter genes, the position of every SNP was called as per gene annotations in the D. immitis genome. The BLASTP search-filter, although applied to help identify nematode-specific ABC transporters (and their associated polymorphisms), did not affect any of the results on polymorphism.

Polymorphic patterns in ABC transporter genes of D. immitis

Out of 17 ABC transporter genes/subunits that encoded 127,275 base pairs (bp), sequence coverage was only available for 119,033 bp. This number of base pairs, that covered 0.14% of ∼84.2 Mb sized nuclear assembly of D. immitis, was analyzed for polymorphism. As previously described (Mani et al., 2016), a position was called polymorphic if there was a 20% difference between the frequency of 2 or more nucleotides. Using this 20% threshold level, we identified 231 SNPs in 119,033 bp with a nucleotide diversity rate of 0.19% (see Table 3). Of the total SNPs, 67 were present across all of the populations, while 89 and 75 were specific to SUS and LOE populations, respectively. We identified 193 SNPs in intron regions, with an intronic SNP rate of 1/412 bp, and the remaining 38 SNPs were present in exons, with SNP rate of 1/1041 bp. Twenty-one of the SNPs in exons caused changes in amino acids (See Table 3). Polymorphic rates calculated for the pooled SUS populations, from different locations, are given in Table 4.

Table 3

Summary of Dirofilaria immitis ABC transporter gene/subunit sequences investigated for SNP analysis.

Parameter	Results
Number of ABC transporter genes/subunits studied	17
Total bases (bp) encoded	127,275
Total bases (bp) covered for SNP analysis	119,033
Intronic bases	79,466
Exonic bases	39,567
Number of SNPs in introns	193
Number of SNPs in exons	38
Number of missense causing SNPs in exons	21

Table 4

Genetic variability of Dirofilaria immitis ABC transporter genes based on geographical locations.

Geographical locations	No. of SNPs found	Polymorphic rate (%)a
USA	102	65.8
Grand Canary (Spain)	98	63.2
Grenada	99	63.8
Italy	38	24.5

Number of SNPs identified (as percentage) among susceptible populations in each country, based on total SNPs identified.

Polymorphism in ABC-B genes

Three full-length transporter Pgps, a pseudogene and 2 half transporters from ABC-B subfamily of D. immitis were investigated for polymorphism. Together, they contained 59 SNPs within 55,588 bps covered with an exonic SNP rate of 1/1457 and an intronic SNP rate of 1/797. Dim-pgp-3, Dim-pgp-10, Dim-pgp-11 and a pseudogene were identified in scaffolds nDi.2.2.scaf00046, -48, −4 and −41 respectively, as previously described (Bourguinat et al., 2016). These authors identified half-transporter genes Dim-haf-1 and Dim-haf-4 in scaffolds nDi.2.2.scaf00023 and nDi.2.2.scaf00101, respectively. In Dim-pgp-3, we located 10 SNPs, three of which were in exon regions. A polymorphism at the 645th codon of the gene changed aliphatic amino acid isoleucine to nucleophilic amino acid threonine (Table 5). This SNP was predicted to be located at the linker region of the transporter (see Supplementary data file S3, Fig. S3A). Two synonymous types of polymorphism occurred at amino acid positions 322 and 834 of the gene. Dim-pgp-10 was the least polymorphic (0.05%) among the Pgp genes with 8 SNPs. Of the two SNPs localized in exon regions, one at amino acid position 176 caused an amino acid variant with either serine if the corresponding codon was TCT or phenylalanine if the corresponding codon was TTT (Table 5). This locus was predicted to be located at the cytoplasmic region between TM2 and TM3 (see Supplementary data file S3, Fig. S3B). A second polymorphism (GCT→GGT) identified at 317th codon of the protein was identified with A317G, predicted to be present at the fifth transmembrane region of the protein (see Supplementary data file S3, Fig. S3B). Dim-pgp-11 (Bourguinat et al., 2016) previously implicated in loss of response to ivermectin (Bourguinat et al., 2011a, Bourguinat et al., 2011b) was found with 18 SNPs. Of those, 8 were only detected in SUS and 5 were only observed in LOE populations. Two (R230P, K1203R) of five exonic SNPs in Dim-pgp-11 were non-synonymous types of polymorphism (Table 5). The SNP K1203R, found previously as an exonic part of the diplotypic “GG-GG” genotype of D. immitis microfilariae that correlated with LOE response to ivermectin (Bourguinat et al., 2011b), was present predominantly in LOE pooled populations as expected (see Supplementary data file S2). Two adjacently placed polymorphic sites at the 230th position of the amino acid were known. A second position change at the 230th codon caused a change from arginine to proline whereas the third position nucleotide change did not affect the amino acid coded. An ABC transporter-encoding gene identified in scaffold nDi.2.2.scaf00041 was called as a pseudogene as it contained multiple genetic lesions that could prevent translation to a functional protein (Bourguinat et al., 2016). Though the gene appears to be non-functional, polymorphism in pseudogenes can be important (Macphee et al., 2002, Bourguinat et al., 2016). Polymorphic analysis of the pseudogene allowed us to identify 20 SNPs, out of which, three were in exonic regions. A deletion type of polymorphism at codon 135 caused a stop codon at the downstream flanking region across the populations. At position 293 of the pseudogene, a polymorphic variant with either an alanine or valine was present (see Supplementary data file S2). Among the half ABC-B transporter genes, Dim-haf-4 had two SNPs in the intron regions of the gene, one exclusively in the LOE populations and the other (deletion polymorphism) across all populations. A single intronic polymorphic locus in the Dim-haf-3 gene was found to be present only in the SUS populations (see Supplementary data file S2).

Table 5

SNP analysis of ABC-B transporter genes/subunits in Dirofilaria immitis.

Putative gene name	SNP position in the scaffolda	Nucleotide change	Effect of polymorphic codon on protein
Dim-pgp-3	232,474	ATT↔ACT	I645T
Dim-pgp-10	74,93776,446	GCT↔GGTTCT↔TTT	A317GS176F
Dim-pgp-11	79,76687,087	AAG↔AGGCGA↔CCA	K1203RR230P

Nucleotides that change are shown in bold.

Scaffold number for each gene given in Table 2.

Polymorphism in ABC-C genes

Three MRP genes (Dim-pgp-2, Dim-pgp-5, Dim-pgp-7) and two half-transporter genes (Dim-haf-5.1, Dim-haf-5.2) (Bourguinat et al., 2016) from the ABC-C subfamily, that together covered 44,689 bp of the nDi2.2 D. immitis genome, were used for SNP investigation. With a total of 75 and 20 SNPs in the intron and exon regions of these genes, respectively, an overall polymorphic rate of 0.2% was determined. The calculated intronic SNP rate was 1/378 and the exonic SNP rate was 1/818. A homolog of the mrp-1 gene, identified in scaffold nDi.2.2.scaf00010, had 38 SNPs, out of which 22 were specific to SUS populations and 8 were specific to LOE populations. Two of the five exonic SNPs caused amino acid change, namely T24P and L1185F, whereas the remaining three SNPs were synonymous types of polymorphism. The SNP L1185F caused a change in predicted secondary structure of the protein (Table 6, Fig. 1). A gene identified in scaffold nDi.2.2.scaf00052 was determined to be a homolog of mrp-5. Dim-mrp-5 had 30 SNPs with estimated polymorphic rate of 0.25%; the most polymorphic amongst the mrp genes studied. Eleven polymorphic sites were specific to SUS populations and another 14 were specific to LOE populations. Seven out of the eight exonic SNPs at amino acids positions 411, 412, 634, 638, 640, 1230 and 1295 of mrp-5, were identified to cause amino acid changes (Table 6), out of which, the SNP A411G was predicted to cause a change in secondary structure (Fig. 1). Dim-mrp-7 in scaffold nDi.2.2.scaf00004 was identified with 16 SNPs, four of them were in exonic regions. All four exonic SNPs were non-synonymous, namely, W258G, K263N, L280S and G1100C (Table 6). As shown in Fig. 1, two of the SNPs, K263N and G1100C would cause a change in secondary structure, as predicted by PSIPRED. The half transporter gene, Dim-haf-5.1 was identified with 11 SNPs. One out of three in the exon regions was non-synonymous, causing a change in amino acid from serine (AGC) to arginine (AAC) at position 61 of the predicted amino acid sequence. This SNP is anticipated to be in the cytoplasmic region, before the TM1 (see Supplementary data file S3, Fig. S3D). Two other SNPs, at positions 84 and 126 of the amino acid sequence, were synonymous types of polymorphism.

Table 6

SNP analysis of ABC-C transporter genes/subunits in Dirofilaria immitis.

Putative gene name	SNP position in the scaffolda	Nucleotide change	Effect of polymorphic codon on protein
Dim-mrp-1	193,274202,979	TTG↔TTTACC↔CCC	L1185FcT24P
Dim- mrp-5b	176,681176,986181,925181,943183,983183,986	CTA↔GTAGTG↔GGGTGT↔TTTAGG↔AAGGAT↔GGTGCT↔GGT	L1295VV1230GC640FR634KD412GA411Gc
Dim-mrp-7	631,503631,520631,867638,339	TGG↔GGGAAA↔AACTTA↔TCAGGC↔TGC	W258GK263NcL280SG1100Cc
Dim-haf-7	21,168	AGC↔AAC	S56N

Nucleotides that change are shown in bold.

Scaffold number for each gene given in Table 2.

Homolog of the gene in C. elegans has a detrimental RNAi phenotype.

PSIPRED predicted change in secondary structure of the protein due to SNP identified.

Fig. 1

PSIPRED predicted changes in secondary structures of ABC transporters due to SNPs identified by this study. The PSIPRED server annotates the query sequence with secondary structure cartoons ( for helix (H), for coil (C)). Confidence level of prediction for each position in a sequence is given by histograms above the cartoon schematics. Only the flanking region of every polymorphic locus (highlighted in yellow) that caused a change in predicted secondary structure of a gene is shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Polymorphism in ABC-G genes

Two ABC-G transporter genes (Dim-wht-4, Dim-wht-8) together covered 17,912 bp of the genome. With 70 SNPs in introns and 5 in exons, they had an overall polymorphic rate of 0.4%. The wht-4 gene, which spanned 12,560 bp in the scaffold nDi.2.2.scaf00589, had at least 69 SNPs, with a calculated polymorphic rate of 0.54%, which was the highest among all the ABC transporter genes studied. Among the SNPs identified, 33 of them were only in LOE populations and 14 of them were specific to SUS populations. All of the 5 SNPs in exons were synonymous, of which 4 occurred only in SUS populations. In another wht gene (Dim-wht-8), identified in scaffold nDi.2.2.scaf00252, all 6 SNPs discovered were in intron regions (Supplementary data file S2) and 5 were specific to SUS populations.

Discussion

Multidrug transporters such as P-glycoproteins are members of the ABC transporter superfamily of efflux transporters and they have been associated with drug resistance in helminths, including nematodes (Lespine et al., 2012, Ardelli, 2013) and trematodes (Mottier et al., 2006, Greenberg, 2013b). These transporters may influence anthelmintic drug concentrations at the target sites of the parasite and the co-administration of MDR inhibitors (e.g., verapamil, valspodar) with anthelmintics has enhanced efficacy in vivo (Lespine et al., 2012). In addition, these transporters are likely to play key roles in worm physiological functions such as excretion, reproduction and modulation of host responses (Greenberg, 2014a, Greenberg, 2014b). Therefore, ABC transporters may prove to be attractive therapeutic targets for agents that can act as inhibitors, or as anti-nematode agents in their own right. Recently, there has been an interest to exploit these transporters for resistance reversal. For resistance mediated by MDR transporters, different strategies for reversal may be considered. New compounds may either act as MDR inhibitors, directly inhibiting efflux, or as inhibitors of transcriptional regulation that may modulate gene expression at the transcription level or at the level of mRNA stability, at the post-transcriptional level (Lespine et al., 2012).

Knowledge of the polymorphism of ABC transporter genes of D. immitis, especially from sub-families A, B, C and G, implicated in drug resistance (Ardelli, 2013), may be relevant and provide useful information during transport inhibitor or reversing agent drug design. Such knowledge might also be helpful to assess potential resistance development. To study polymorphism, we wished to identify all putative ABC transporter genes/subunits of D. immitis by comparative genomic approaches. We followed this approach as the genome of D. immitis was not fully annotated (Bourguinat et al., 2015). We found that the genome of D. immitis contains at least 15 ABC transporter genes; 2 abt genes from sub-family A, 3 pgp, 1 pseudogene and 2 half transporters from sub-family B, 3 mrp, 2 half transporters from subfamily C, and 2 white-like genes from subfamily G. Ardelli (2013) predicted 33 ABC transporter genes (including 8 pgp, 8 haf and 5 mrp genes) in B. malayi, a closely related filarial nematode; whereas another filarial nematode, O. volvulus was reported to have 7 ABC transporter genes (including 1 pgp, 3 haf and 1 mrp genes) (Ardelli, 2013).

The polymorphic rate was calculated for each SUS pooled population and it was found that the rate was similar between samples from the USA, Grand Canary and Grenada. However, a low polymorphic rate was determined in samples from Italy and this can be explained by the fact that those worms were from a single dog. This hypo-variability of ABC transporter genes is in agreement with the finding of low genetic variability among ion channel genes (Mani et al., 2016). Moreover, polymorphism studied at the microsatellite level (Belanger et al., 2011) and at the whole genome level (Godel et al., 2012) also suggested low genetic variability among D. immitis populations. The latter study, however, used two single D. immitis worms, one each from Italy and the USA to draw this conclusion. This trend of hypovariability would be consistent with recent common ancestry of European and American heartworm populations due to migration.

D. immitis Pgps could be of interest as targets for multidrug resistance reversal, based on evidence that overexpression of Pgps plays a role in resistance in other nematode species and that resistance mediated by Pgp genes in a related nematode can be reversed after the introduction of inhibitors. A significant increase in the expression profiles of Cel-pgp-11, Cel-haf-1, Cel-mrp-1, Cel-mrp-5, Cel-mrp-7 after IVM treatment (Yan et al., 2012, Ardelli and Prichard, 2013, Bygarski et al., 2014) was reported in an IVM-resistant strain, compared to the wild-type Bristol N2 strain. Sensitivity to ivermectin was increased by MDR inhibitors in B. malayi (Tompkins et al., 2011). P-glycoproteins of D. immitis could also be interesting targets for new anthelmintic drugs, considering the physiological roles their homologs play in other nematodes, including C. elegans. Deletion of some of the Pgp genes had a measurable effect on phenotypes of C. elegans (Janssen et al., 2013). For example, C. elegans deletion strains pgp-3 and pgp-10 were found to have significantly higher levels of pharyngeal pump rate than the wild-type strain, suggesting that these genes affected feeding in these strains. Also, the pgp-11 gene is known to be essential for larval development of C. elegans (Janssen et al., 2013). These findings suggest that PGP-34, PGP-10, and PGP-11 of C. elegans have phenotypic effects on resistance, feeding and development of the parasite. If we suppose similar gene product activities for Dim-pgp-3, Dim-pgp-10, Dim-pgp-11, in D. immitis, then these genes could be of potential interest as targets for candidate drugs.

Six ABC-B transporter genes (Dim-pgp-3, Dim-pgp-10, Dim-pgp-11, 1 pseudogene, Dim-haf-1, Dim-haf-4), previously identified and named (Bourguinat et al., 2016) were found to have at least 59 SNPs, an overall polymorphic rate of 0.1%. In Dim-pgp-3, a missense variant (I645T) with an aliphatic amino acid, isoleucine, or a nucleophilic amino acid, threonine, predicted to be in the linker region of the protein (Supplementary data file S3, Fig. S3A). It has been shown that the linker region of ABC transporters mediates ubiquitination and control of protein turnover (Kölling and Losko, 1997, Sato et al., 2009), two processes implicated in mechanisms of drug resistance in Plasmodium (Hunt et al., 2007, Deplaine et al., 2011). The pgp-3 homolog in C. elegans, Cel-pgp-3, found to be expressed in the intestine and excretory cell (Zhao et al., 2004), has a role in feeding of the parasite (Bygarski et al., 2014). Considering all this, the pgp-3 of D. immitis could be a potential target for a reversal mechanism. In the Dim-pgp-10 gene, 8 SNPs, including 2 exonic SNPs, were identified. In C. elegans, pgp-10 is expressed in the intestine (Zhao et al., 2004) but its role is not yet known. Further studies need to be done to understand the implications of the SNPs identified in Dim-pgp-10.

Another Pgp gene of D. immitis, Dim-pgp-11, was found to be very polymorphic as it contained at least 18 SNPs. Of them, five were exonic including one of interest, K1203R, found mainly in LOE populations. The same SNP (K1203R (AAG→AGG)) in the Dim-pgp-11 gene has been captured in different forms by two separate studies, validating our approach to identify SNPs in ABC transporter genes. It was identified as position 11 that corresponded to the 1st GG of the GG-GG genotype of a partial Dim-pgp-11 sequence (GenBank accession number HM596853) (Bourguinat et al., 2011b). In another study, it was known as LOCUS_48,992_B, identified by comparing read frequencies of four SUS and four LOE D. immitis populations (Bourguinat et al., 2015). The NBDs and Walker A and Walker B sites are very conserved across ABC-B proteins, and recently it has been shown (Verhalen et al., 2017) that changing E1197 to 1197Q; which makes the NBD region just adjacent to the Walker B site more polar and more basic, impairs ATP turnover in mammalian P-glycoprotein, and thus likely transport function. A similar change occurring with the K1203R SNP, found in Dim-PGP-11, similarly introduces a more polar, very basic group adjacent to the Walker B site which may modulate ATP turnover and transport activity. Alignment of this region in a number of organisms indicates that K1203 is highly conserved. Thus the mutation to 1203R may represent an important change affecting the functioning of the Dim-PGP-11 transporter.

In the pgp-11 of D. immitis, two additional SNPs at the second and third position of codon 230 were identified, with the former causing an amino acid change from arginine (basic amino acid) to proline (cyclic amino acid). This non-synonymous SNP was predicted to be present at the cytoplasmic region between TM4 and TM5, immediately after the TM4 region (see Supplementary data file S3, Fig. S3C). The crystal structure of a eukaryotic homolog, CmABCB1 (from Cyanidioschyzon merolae) revealed that the TM4 helix acts as a “gatekeeper” of the substrate entrance pore and is therefore critical for efficient drug binding and transport (Kodan et al., 2014). In the mouse Pgp, it was suggested that TM4 forms one of the portals open to the inner leaflet for drug entry and so is believed to play an important role in selecting drugs from the lipid bilayer (Li et al., 2014). If we assume a similar scenario, any SNP present in and around the TM4 helix of PGP-11 of D. immitis, may play a role in substrate selection, binding and transport; therefore the SNP R230P in Dim-pgp-11 maybe of interest during drug design. Another nsSNP, K1203R, predicted to be present immediately before the Walker B motif (ILLDE) in the second nucleotide binding domain of the transporter (see Supplementary data file S3, Fig. S3C) may be of possible interest for further study as this polymorphic locus was implicated in the LOE phenotype to IVM treatment (Bourguinat et al., 2011b). D. immitis PGP-11 is only 36% identical to PGP-1 of C. elegans for which a higher resolution (3.4 ºA) crystal structure is available (Jin et al., 2012) and has been shown to bind ML anthelmintics (David et al., 2016) and therefore caution is required in using the structural knowledge of Cel-PGP-1 to compute functional implications of identified SNPs in Dim-PGP-11. In C. elegans, the pgp-11 gene is known to be expressed in the intestine and excretory cell (Zhao et al., 2004). Moreover, the homolog of pgp-11 is essential for survival in C. elegans, as the gene is required for larval development (Janssen et al., 2013). In D. immitis, a diplotypic GG-GG genotype of the gene strongly correlated with loss of ivermectin response phenotype of mfs (Bourguinat et al., 2011a, Bourguinat et al., 2011b). Taken together, the pgp-11 of D. immitis may be a gene of potential interest for drug discovery.

The MRPs may also be of interest as candidate drug targets. Their protective role in worms is well known as they were found to reduce the toxicity of heavy metals (Broeks et al., 1996) and MLs (Lespine et al., 2012). There is an increased expression of MRPs in wild-type (Lespine et al., 2012) and ivermectin-resistant strains (James and Davey, 2009) of C. elegans after exposure to MLs. IVM is a known substrate of mammalian MRPs (Lespine et al., 2006). Considering these factors, MRP homolog genes in D. immitis, namely Dim-mrp-1, Dim-mrp-5 and Dim-mrp-7 may be interesting drug targets. In the mrp-1 gene, a polymorphism at codon 24 changed the amino acid from threonine (ACC) to proline (CCC). Another SNP with either a leucine or phenylalanine at position 1185 of the amino acid sequence was identified. The homolog of mrp-1, in C. elegans, regulates dauer larva formation (Yabe et al., 2005) in addition to its protective role against heavy metals (Broeks et al., 1996). Therefore Dim-mrp-1 (and its associated SNP) may be an interesting target for drug design.

In C. elegans, the likely role of the MRP-5 protein, as a heme exporter, was shown and targeted deletion of mrp-5 in the intestine caused embryonic lethality in the mutant worms (Korolnek et al., 2014). MRP-5 of C. elegans is thus essential for viability and therefore could be an interesting drug target. The gene product of Dim-mrp-5 may also be a target of interest. However, Dim-mrp-5 was found to be highly polymorphic, with at least 30 SNPs (eight of them in the exon regions) and so there is a higher chance of resistance development against a drug targeting this gene. In Dim-mrp-7 of D. immitis, 4 SNPs were identified in exonic regions, all of which resulted in non-synonymous types of polymorphism. Some of the SNPs identified within mrp genes, namely L1185F (mrp-1), A411G (mrp-5), K263N and G1100C (mrp-7) are predicted to cause a change in secondary structure (see Fig. 1). Those SNPs might have interesting effects on their respective protein function and therefore further studies are needed to understand their implications. In the Dim-mrp-7 (ABCC10) homolog of humans, two SNPs (one in intron 4 and another in exon 12) were associated with renal tubular dysfunction (Kathawala et al., 2014).

We believe that some of the SNPs identified, including SNPs in non-coding sequence, may have effects on gene expression, mRNA stability, function, substrate specificity, or drug resistance. We are aware that some polymorphic sites may have been missed in this study, due to lack of coverage during sequencing, small sample size and limited sample diversity, or due to the use of the 20% threshold value for SNP calling. Further studies need to be done to confirm each SNP in the field and to understand its implications.

Conclusion

Information on the degree of heterogeneity of D. immitis has been limited and this study adds to that information. The study lists all genes detected from A, B, C and G sub-families of ABC transporters. The merit of the study lies in the diverse large number of samples (122 adult worms and ∼32,000 mfs) used for SNP genotyping, with 231 SNPs identified that may be of relevance for MDR inhibitor or antinematodal drug design. This study also provides an updated inventory of ABC transporter genes in D. immitis.