The genome of the heartworm, Dirofilaria immitis, reveals drug and vaccine targets.

The heartworm Dirofilaria immitis is an important parasite of dogs. Transmitted by mosquitoes in warmer climatic zones, it is spreading across southern Europe and the Americas at an alarming pace. There is no vaccine, and chemotherapy is prone to complications. To learn more about this parasite, we have sequenced the genomes of D. immitis and its endosymbiont Wolbachia. We predict 10,179 protein coding genes in the 84.2 Mb of the nuclear genome, and 823 genes in the 0.9-Mb Wolbachia genome. The D. immitis genome harbors neither DNA transposons nor active retrotransposons, and there is very little genetic variation between two sequenced isolates from Europe and the United States. The differential presence of anabolic pathways such as heme and nucleotide biosynthesis hints at the intricate metabolic interrelationship between the heartworm and Wolbachia. Comparing the proteome of D. immitis with other nematodes and with mammalian hosts, we identify families of potential drug targets, immune modulators, and vaccine candidates. This genome sequence will support the development of new tools against dirofilariasis and aid efforts to combat related human pathogens, the causative agents of lymphatic filariasis and river blindness.

The heartworm dirofilaria immitis (Leidy, 1856) is a parasitic nematode of mammals. The definitive host is the dog; however, it also infects cats, foxes, coyotes, and, very rarely, humans (1). Dirofilariasis of dogs is a severe and potentially fatal disease. Adult nematodes of 20 to 30 cm reside in the pulmonary arteries, and the initial damage is to the lung. The spectrum of subsequent pathologies related to chronic heartworm infection is broad, the most serious manifestation being heart failure. Recent rapid spread of D. immitis through the United States and southern Europe (2, 3) is being favored by multiple factors. Global warming is expanding the activity season of vector mosquitoes, increasing their abundance and the likelihood of transmission of the parasite, and there are growing numbers of pets, reservoir animals, and “traveling” dogs (2, 3).

D. immitis is an onchocercid filarial nematode, related to important parasites of humans, such as Onchocerca volvulus, the agent of river blindness. The D. immitis lifecycle is typical for Onchocercidae. Microfilariae, shed into the bloodstream by adult females, are ingested by a mosquito (various species, including Aedes, Anopheles, and Culex spp.) where they develop into third-stage larvae (L3) and migrate to the labium. Feeding by an infected mosquito introduces L3 into the skin. The prepatent period in the newly bitten dog is 6–9 mo, during which the injected larvae undergo two further molts and migrate via muscle fibers to the pulmonary vasculature, where the adult nematodes develop. At present, diagnosis is effective only for patent infections, because it is based on detection of circulating microfilariae or antigens from mature females. Treatment of dirofilariasis is also problematic, because the arsenical melarsomine dihydrochloride, the only adulticide approved by the U.S. Food and Drug Administration, can cause adverse neurological reactions. Treatment carries a significant risk of lethality due to blockage of the pulmonary artery by dead nematodes. No vaccine is available. These issues, together with the alarming increasing spread of D. immitis, prompted the American Heartworm Society to recommend year-round chemoprophylactic treatment of dogs (4) to kill the larval stages before they develop into adults. This requires monthly administration of anthelmintics, predominantly macrocyclic lactones, such as ivermectin, milbemycin, or moxidectin.

Human-infective parasites related to D. immitis cause subcutaneous filariasis and river blindness and are endemic in tropical and subtropical regions around the globe, with an estimated 380 million people affected (5). Improved diagnostics, new drugs, and, ultimately, effective vaccines are sorely needed. The sequencing of the Brugia malayi genome provides a platform for rational drug design, but by itself this single sequence cannot distinguish between idiosyncratic and shared targets that could be exploited for control (6).

Most of the filarial nematodes that cause diseases in humans and animals, including D. immitis, O. volvulus, Wuchereria bancrofti, and B. malayi, have been shown to harbor intracellular symbiotic bacteria of the genus Wolbachia (e.g., refs. 7–10). These bacteria are vertically transmitted to the nematode progeny, via transovarial transmission. In most of the infected nematode species, all individuals are infected (reviewed in ref. 11). Even though the exact role of Wolbachia in filarial biology has not yet been determined, these bacteria are thought to be beneficial to the nematode host. Indeed, antibiotics that target Wolbachia have been shown to have deleterious effects on filarial nematodes, blocking reproduction, inducing developmental arrest, and killing adult nematodes (e.g., refs. 7–9). This has led to development of research projects with the aim of developing anti-Wolbachia chemotherapy as a novel strategy for the control of filarial diseases. Wolbachia has also been implicated in the immunopathogenesis of filarial diseases, with a role in the development of pathological outcomes, such as inflammation and clouding of the cornea that is typical of river blindness (12). The genome of Wolbachia is thus an additional source of potential drug targets (7–10), but a single genome cannot reveal shared vs. unique biochemical weaknesses.

The human pathogenic Onchocercidae do not represent an attractive market for the pharmacological industry, because projected incomes from impoverished communities in developing endemic nations would be unlikely to cover the costs of drug development. The heartworm may hold a possible solution to this problem, because the market potential for novel canine anthelmintics is big, given the costs for heartworm prevention of $75–100/dog/yr and the estimated number of 80 million dogs in the United States (13). Choosing drug targets that are likely to be conserved in related, human pathogenic species may benefit both canine and human medicine. Here we present the draft genome sequences of D. immitis and its Wolbachia endosymbiont (wDi) and use these data to investigate the relationship between nematode and endosymbiont and identify new drug and vaccine targets.

MATERIALS AND METHODS

D. immitis isolates and DNA sequencing

We sequenced two canine D. immitis isolates, one from Pavia, Italy, and the other from Athens, Georgia, USA. The Pavia isolate was established in a laboratory lifecycle after primary isolation from an infected dog. Adult Pavia nematodes used for DNA extraction were recovered after necropsy of dogs infected as a control group in ongoing investigations (permit FR401e/08 from the Veterinary Office Canton de Fribourg, Switzerland). The Athens nematodes used for DNA and RNA extraction were from a naturally infected dog necropsied as part of routine clinical surveillance and were not from an established strain. Genomic DNA was extracted (QIAamp DNA extraction kit; Qiagen, Valencia, CA, USA) from individual adult female nematodes from Pavia and Athens isolates, and RNA was extracted (RNeasy kit; Qiagen) from individual female and male nematodes from Athens. Whole-genome shotgun sequences were generated at The GenePool Genomics Facility (University of Edinburgh, Edinburgh, UK) and at Fasteris SA (Geneva, Switzerland) using Illumina GAIIx and HiSeq2000 instruments (Illumina, Inc., San Diego, CA, USA). Several short insert (100- to 400-bp) paired-end amplicon libraries and long insert (3- to 4-kb) mate-pair amplicon libraries were made, and data from four of these were used in the final assembly (details are given on the Web site http://www.dirofilaria.org). These yielded a raw data total of 28 Gb in 295 million reads [European Bioinformatics Institute (EBI) Short Read Archive, accession number ERA032353; http://www.ebi.ac.uk/ena/]. After trimming low-quality bases (Phred score <20) and filtering out reads with uncalled bases or length <35 b, 271 million reads were used for assembly (Supplemental Table S1).

Nuclear, mitochondrial, and Wolbachia genome assemblies

The short-read data were assembled using ABySS 1.2.3 (14). A number of test assemblies were performed using other assemblers, and a range of parameters was tested within ABySS, and the final, optimal assembly was performed using a k-mer length of 35 and scaffolding with the paired-end data only. Assembly qualities were assessed using summary statistics including maximizing the N50 (the contig length at which 50% of the assembly span was in contigs of that length or greater), maximum contig length, and total number of bases in contigs (see Supplemental Table S1) and using biological optimality assessment, such as maximizing the coverage of published D. immitis expressed sequence tag (EST) sequences and maximizing the number of B. malayi genes matched and the completeness of representation of core eukaryotic genes (using CEGMA; ref. 15). Redundancy due to allelic polymorphism was reduced with CD-HIT-EST (16), merging contigs that were ≥97% identical over the full length of the shorter contig. The mitochondrial genome was assembled by mapping the reads to the published D. immitis mitochondrial genome (17) and predicting a consensus sequence of the mitochondrial genomes of the Athens and Pavia nematodes separately. The wDi genome was assembled by first identifying likely wDi contigs in the whole assembly with BLASTn (18) using all Wolbachia genomes from EMBL-Bank, and then collecting all raw reads (and their pairs; n=6,912,659) that mapped to these putative wDi genome fragments. The reduced set of likely wDi reads was then assembled using an independently optimized ABySS parameter set, using mate-pair information where available. Mitochondrial and wDi contigs were removed from the full assembly to leave the final nuclear assembly.

Transcriptome shotgun sequencing (RNA-Seq) assembly

The preparation of amplicon libraries and RNA-Seq analysis were performed following standard Illumina TruSeq protocols. A total of 11,019,886 (male) and 21,643,293 (female) read pairs of length 54 b were produced on the Illumina GAIIx platform (ArrayExpress accession number E-MTAB-714; ENA study accession number ERP000758). After quality filtering, the remaining 31,396,183 pairs were assembled with Trans-ABySS using k-mer values from 23 to 47 in steps of 4 (Supplemental Table S1).

D. immitis nuclear genome protein-coding gene prediction and analysis

Repeats in the D. immitis genomic assembly were identified and masked using RepeatMasker 3.2.9 (19), including all “Nematoda” repeats in the RepBase libraries (20). The MAKER 2.08 annotation pipeline (21) was used to identify protein-coding genes based on evidence from the RNA-Seq assembly, alignments to the B. malayi proteome (WormBase release WS220; http://wormbase.sanger.ac.uk/), predictions made by the ab initio gene finder SNAP (22), and predictions from the ab initio gene finder Augustus (23) based on the Augustus hidden Markov model (HMM) profiles for B. malayi. MAKER predicted 11,895 gene models, and, with alternative splicing, a total of 12,872 transcripts and peptides. We compared the nuclear proteome of D. immitis with those of four other species for which complete genome data are available and which span the phylogenetic diversity of the phylum Nematoda (B. malayi, Ascaris suum, Caenorhabditis elegans, and Trichinella spiralis). The complete proteomes were compared using all-against-all BLAST, and then clustered using OrthoMCL (24). OrthoMCL clusters were postprocessed to classify clusters by their species content and analyzed with reference to the robust molecular phylogeny of the Nematoda (25). The prediction of D. immitis orthologs from B. malayi, C. elegans, Homo sapiens, and Canis lupus to identify drug targets was performed with InParanoid (26).

Analyses of orthology and divergence in filarial Wolbachia

The wDi genome was annotated with the RAST server (27), an online resource that uses best-practice algorithms to perform both gene finding and gene functional annotation. Selected metabolic pathways were annotated based on enzyme lists from the KEGG Pathway database (28), after an HMM profile was generated for each enzyme (29) from a ClustalW (30) multiple alignment of a redundancy-reduced set of all the manually curated entries in UniProt (31). Analysis of orthology was performed using the BLAST reciprocal best-hits algorithm (32), with the following cutoff values: E value 0.1 and ID percentage 60%. Protein distance for each pair of orthologs was calculated using Protdist in Phylip 3.69 (33) with the Dayhoff PAM matrix option. Proteins were allocated to functional categories using BLAST against the COG database. Protein distances were then analyzed based on COG categories: within each category we calculated the average distances of protein pairs. To evaluate whether some categories were significantly more variable than others, we performed the Kruskal-Wallis test on COG categories containing more than one ortholog pair. The pairwise Mann-Whitney test was then performed to detect pairs of COG categories that displayed significant differences in their average variation.

Identification of Wolbachia insertions in nematode genomes

To identify potential lateral genetic transfers from Wolbachia to the host nuclear genome, the nuclear genome was queried against the 921-kbp wDi genome using the dc-megablast option in BLASTN (NCBI-blast+2.2.25) with default settings. All high-scoring pairs (HSPs) longer than 100 bp with >80% identity were kept. Overlapping HSP coordinates on the nuclear genome were merged, and sequences from these coordinates were extracted to obtain putative nuclear Wolbachia DNA elements. The B. malayi nuclear genome was screened with the B. malayi Wolbachia (wBm) genome in the same way. The small numbers of Wolbachia insertions identified in the nuclear genomes of Acanthocheilonema viteae and Onchocerca flexuosa (34) were surveyed for matches to the wDi and wBm genomes and cross-compared with the insertion sets from the complete D. immitis and B. malayi genomes using reciprocal best BLAST searches and filtering alignments shorter than 100 bp. Reciprocal best BLAST matches were isolated and single-linkage clustered.

RESULTS

Genome assembly of D. immitis and its Wolbachia symbiont

Genome sequence was generated from single individuals of D. immitis isolated from naturally infected dogs, one from Athens, Georgia (USA) and the other from Pavia (Italy). A total of 16 Gb of raw data was retained after rigorous quality checks, corresponding to ∼170-fold coverage of the D. immitis nuclear genome (likely to be ∼95 Mb, similar to related Onchocercidae). The ABySS (35) assembler performed best based on statistical and biological measures (Supplemental Table S1). The mitochondrial and Wolbachia wDi genomes were assembled independently. The final nuclear assembly contained 84.2 Mb of sequence in 31,291 scaffolds with an N50 of 10,584 bases (). The draft genome of wDi consists of 2 scaffolds spanning 0.92 Mb. We identified 99% of previously deposited genome survey sequences putatively from wDi (GenBank accession numbers ET041559 to ET041665) within our wDi assembly. The wDi genome was 16% smaller than that of wBm (1.08 Mb; GenBank accession number NC_006833), and there was significant breakage of synteny between the two genomes, as has been observed between other Wolbachia.

The D. immitis and wDi genomes, the annotations we have made on these, and additional technical details and analyses are available through a dedicated genome browser (http://www.dirofilaria.org).

Lack of genetic diversity between the sequenced D. immitis isolates

Even though the two sequenced D. immitis came from independent isolates from different continents, they showed low genetic differentiation, allowing the raw sequencing data from both nematodes to be coassembled. We mapped the reads from each nematode back to the draft assembly and identified only 32,729 high-quality single-nucleotide variations, a very low per-nucleotide diversity rate of 0.04%. We identified the sequences corresponding to 11 polymorphic microsatellite loci used previously to analyze the D. immitis population structure in North America (36) and genotyped our two isolates in silico by counting the predicted numbers of microsatellite repeats at each locus. Both our nematodes could be classified within the diversity of the eastern United States population. The mitochondrial genomes of the two isolates differed at only 6 sites (and were thus >99.9% identical). Surprisingly, compared with the published, Australian D. immitis mitochondrion (17), both had many shared differences (each was only 99.5% identical to the published D. immitis mitochondrion). Because the ∼70 differences were often clumped and were unique in the published D. immitis mitochondrial genome compared both with our two genomes and with the genomes of five other filarial nematodes, we suggest that many of these are sequencing errors in the published genome.

A metazoan genome without active transposable elements

The D. immitis genome was surveyed for the three main classes of transposable elements [DNA transposons, long terminal repeat (LTR) retrotransposons, and non-LTR retrotransposons] with tBLASTn (18) using the transposon-encoded proteins as queries. No traces of active or pseudogenized DNA transposons or non-LTR retrotransposons were found, but 376 fragments of LTR retrotransposons of the BEL/Pao family (37) were identified. None of these fragments were predicted to be functional, because all contained frame shifts and stop codons in the likely coding sequence. The D. immitis Pao pseudogenes were most similar to Pao family retrotransposons from B. malayi (6). In B. malayi, several of the Pao retrotransposons are likely to be active, because they have complete open reading frames and LTRs. Overall, however, B. malayi has a lower density of Pao elements and fragments (3.4 Pao/Mb, 8.3% of which are predicted to be functionally intact) compared with D. immitis (4.6 Pao/Mb, none of which were intact).

D. immitis nuclear proteome

Protein coding genes were predicted in the nuclear assembly using the MAKER pipeline (21), integrating evidence-based (RNA-Seq and known protein mapping) and ab initio methods. Of the 11,375 gene models, 897 were predicted to generate alternate transcripts (Table 1). The total number of predicted proteins of length ≥100 aa was 10,179, similar to the 9807 predicted in B. malayi. Based on matches to D. immitis ESTs and core eukaryotic genes (15), the D. immitis proteome was likely to be near-complete. Protein-coding exons occupy ∼18% of the genome of D. immitis and 14% of the genome of B. malayi (Table 1), but in C. elegans there are nearly twice as many genes, and exons cover ∼30% of the genome. The median global identity between a D. immitis protein and its best match (as determined by BLASTp) in B. malayi was 75%.

Table 1

Comparison of the genome assemblies of D. immitis, B. malayi, and C. elegans

Characteristic	D. immitis	B. malayi	C. elegans
Assembly size (Mb)	84.2	93.6a	100.3
Protein-coding gene models	11,375	11,434	20,517
Genes per megabase	135	122	205
Predicted proteins	12,344	11,460	31,249
Protein-coding sequence (%)	18.0	13.8	25.4
Median exons per gene	5	5	6
Median exon size (b)	142	139	147
Median intron size (b)	226	213	73
Overall GC content (%)	28.3	30.2	35.4
Exon GC content (%)	37.4	39.4	43.4
Intron GC content (%)	26.6	27.2	32.5

B. malayi data are from the GenBank RefSeq dataset; C. elegans data from the WS230 dataset.

70.8 Mb scaffolds + 17.5 Mb short contigs.

D. immitis proteins were clustered with the complete proteomes of four other nematode species. These clusters were classified and mapped onto the phylogenetic tree of the five species based on the placement of the deepest node that linked the species that contributed members (). The D. immitis proteome included 3199 proteins (31% of the total proteome) that were unique to this species, a proportion similar to that found in B. malayi (27%), but many fewer (and a lower proportion) compared with those for the other species (for example, C. elegans had 63% of its proteome in species-unique clusters). This difference may be partly due to the 850 proteins in clusters uniquely shared by the relatively closely related D. immitis and B. malayi, but these clusters only raise the proportion of proteins in phylogenetically local clusters to 47%.

Conserved and novel genes in D. immitis. The D. immitis proteome was clustered with those of B. malayi, A. suum, C. elegans, and T. spiralis. Clusters were then classified based on the membership from the five species according to the current phylogeny of the phylum Nematoda. A) Pie chart showing the distribution of classification of D. immitis proteins: D. immitis only, singletons and clusters only found in D. immitis; Onchocercidae, clusters with members only from D. immitis and B. malayi; Spiruria, clusters with members only from Onchocercidae and A. suum; Rhabditia clusters with members only from Spiruria and C. elegans; Nematoda, clusters with members from all five species (i.e., Rhabditia and T. spiralis); and other patterns, clusters with members not fitting simply into the phylogenetic schema (probably arising from gene loss, lack of predictions, or failure to cluster in one or more species). B) Cluster numbers and patterns of conservation mapped onto the phylogeny of the five species.

D. immitis genes homologous to known antinematode drug targets

An array of drugs are effective against nematode parasites (). Of these, flubendazole (38), mebendazole (39), levamisole (40), ivermectin, milbemycin, moxidectin, and selamectin (41, 42) have been demonstrated to be active against D. immitis. Many drug targets have been identified, particularly through forward genetics in the model nematode C. elegans (43) (Table 2). Prominent among these targets are neuronal membrane proteins, highlighting the importance of the neuromuscular junction as a hotspot of anthelmintic drug action. D. immitis appears to lack some known targets, notably members of the DEG-3 subfamily of acetylcholine receptors, which contains the presumed targets of monepantel (44). This contrasts with B. malayi, which possesses orthologs of DEG-3 and DES-2 (45). In C. elegans, the target space of levamisole and ivermectin comprises a large number of ligand-gated ion channels. Although these drugs are effective against heartworm, some of these ion channels do not have an ortholog in D. immitis (Table 2), indicating that those present are sufficient to confer drug susceptibility. The identified D. immitis orthologs of the known anthelmintic targets can now be monitored in suspected cases of drug resistance.

Table 2

Candidate drug targets, top-down search: current anthelmintics and their known targets in C. elegans and orthologs in D. immitis

Chemical class	Drug	Target	C. elegans	D. immitis
Benzimidazole	Albendazole	β-Tubulin	BEN-1	DIMM36740
	Flubendazole
	Mebendazole
Imidazothiazole	Levamisole	nACh receptor	LEV-1
			LEV-8
			UNC-29	DIMM30000
			UNC-38	DIMM45965
			UNC-63	DIMM08405
Macrocyclic lactone	Ivermectin	Glutamate receptor	AVR-14	DIMM16610
	Milbemycin		AVR-15
	Moxidectin		GLC-1
	Selamectin		GLC-2	DIMM25280, DIMM21120
			GLC-3
			GLC-4	DIMM22030
		GABA receptor	EXP-1	DIMM57890
			GAB-1
			UNC-49	DIMM33210
Cyclodepsipeptide	Emodepside	K+ channel	SLO-1	DIMM33710
		Latrophilin GPCR	LAT-1	DIMM37270, DIMM37275
			LAT-2	DIMM17690
Aminoacetonitrile derivative	Monepantel	nACh receptor	ACR-23
			DES-2

nAChR, nicotinic acetylcholine; GPCR, G protein-coupled receptor.

New drug target candidates in D. immitis

New potential drug targets were identified in silico through an exclusion-inclusion strategy (46, 47). Starting from the complete set of predicted D. immitis proteins, we excluded proteins that had an ortholog in the dog or human proteome or had multiple paralogs in D. immitis. We included proteins that had a C. elegans ortholog essential for survival or development (based on RNAi phenotypes) and had predicted function as an enzyme or receptor. Among the 20 candidates identified () were several proven drug targets, such as RNA-dependent RNA polymerase (antiviral), apurinic/apyrimidinic endonuclease and hedgehog proteins (anticancer; ref. 48), UDP-galactopyranose mutase (against mycobacteria, ref. 49; and kinetoplastids, ref. 50), sterol-C24-methyltransferase (antifungal; ref. 51), and the insecticide target chitin synthase (52). The D. immitis orthologs of these enzymes may serve as starting points for the development of new anthelmintics.

Immune modulators and vaccine candidates

Filarial nematodes modulate the immune systems of their mammalian hosts to promote their own survival and fecundity, but the exact mechanisms used remain enigmatic. Proteases such as leucyl aminopeptidase and protease inhibitors such as serpins and cystatins have been implicated in disruption of immune signal processing (53), and we identified D. immitis leucyl aminopeptidase, as well as 3 cystatins, and many serpins (). Another route to modulation is through recruitment of nematode homologs of ancient system molecules that have been redeployed in the mammalian immune system, such as TGF-β and macrophage migration inhibition factor (MIF). In D. immitis, we identified 2 MIF genes, orthologs of the MIF-1 and MIF-2 genes of B. malayi and O. volvulus and 4 TGF-β homologs (Table 4). Another proposed route to modulation is by mimicry of immune system signals. We identified a homolog of suppressor of cytokine signaling 5 (SOCS5), a negative regulator of the JAK/STAT pathway and inhibitor of the IL-4 pathway in T-helper cells, promoting TH1 differentiation (54). Several viruses induce host SOCS protein expression for immune evasion and survival (55). Interestingly, SOCS5 homologs were also identified in the animal-parasitic nematodes B. malayi, D. immitis, Loa loa, A. suum, and T. spiralis, but were absent from the free-living C. elegans, the necromenic Pristionchus pacificus, and the plant parasitic Meloidogyne spp. D. immitis and other filarial nematodes (56) may use SOCS5 homologs to mimic host SOCS5. We also identified a homolog of IL-16, a PDZ domain-containing, pleiotropic cytokine (57). In mammals, IL-16 acts via the CD4 receptor to modulate the activity of a wide range of immune effector cells, including T cells and dendritic cells (58). Again, this molecule was only present in parasitic nematodes (including A. suum; ref. 59) and was absent from genomes of free-living and plant parasitic species. We suggest that these molecules and perhaps other mimics of cytokines and modulators belong to the effector toolkit used by filarial nematodes to build an immunologically compromised niche.

Table 4

D. immitis potential immune modulators and orthologs of onchocercid vaccine candidates

D. immitis protein	B. malayi ortholog	Description	Potential
DIMM39040, DIMM39045	BM18548	Pi-class glutathione S-transferase (GSTP)	VC
DIMM29150	BM02625	Tropomyosin (TMY)	VC
DIMM29270	BM00759, BM19824	Fatty acid and retinoic acid binding protein (FAR)	VC
DIMM47055	BM03010	Fructose bisphosphate aldolase (FBA)	VC
DIMM59360		Astacin metalloprotease MP1	VC
DIMM37935, DIMM46475	BM01859, BM09541, BM14520	Chitinase (CHI)	VC
DIMM48695	BM21967, BM08119	Abundant larval transcript 1 (ALT); unknown function (also known as SLAP)	VC
DIMM48700	BM20051	“RAL-2,” unknown function; DUF148 superfamily (also known as SXP-1)	VC
DIMM62215, DIMM45570, DIMM58880	BM03177, BM05783, BM16294	Activation associated proteins [ASP, also known as venom allergen homologs (VAH)]	VC
DIMM58690	BM02480	“OV103” Onchocerca vaccine candidate of unknown function	VC
DIMM12355	BM07484, BM22082	“B8” Onchocerca vaccine candidate of unknown function	VC
DIMM55190, DIMM50565, DIMM48395	BM00175, BM14240, BM04930, BM07956	“B20” Onchocerca vaccine candidate of unknown function	VC
DIMM56580	BM05118	Cysteine proteinase inhibitor 2 (CPI-2)	VC/IM
DIMM18905	BM04900	Cysteine proteinase inhibitor 3 (CPI-3)	VC/IM
DIMM11425	BM21284	Interleukin-16-like (IL16)	IM
DIMM57180	BM00325	Leucyl aminopeptidase (LAP)	IM
DIMM28945	BM06847	Suppressor of cytokine signaling 5 (SOCS5)	IM
DIMM42430	BM07480	Macrophage migration inhibitory factor (MIF-1)	IM
DIMM40455	BM16561	Macrophage migration inhibitory factor 2 (MIF-2)	IM
DIMM23225	BM17713	Transforming growth factor β (TGF) homolog of C. elegans TIG-2	IM
DIMM37585	BM20852	TGF homologue of C. elegans DAF-7	IM
DIMM29335	BM21753	TGF homologue of C. elegans DBL-1/CET-1	IM
DIMM61250	BM18112	TGF homologue of C. elegans UNC-129	IM

B. malayi orthologs are referred to by their designation in WormBase WS230. IM, immune modulator; VC, vaccine candidate.

We surveyed the D. immitis genome for molecules currently proposed as vaccine candidates in other onchocercids (60, 61) and identified homologs for all 14 classes of molecules (Table 4).

Analysis of the wDi genome: the D. immitis-Wolbachia symbiosis

wDi genes were predicted using the RAST online server. We performed an orthology analysis comparing wBm and wDi and found 538 shared proteins. There were 259 (with 8 duplicated) and 329 (with 4 duplicated) unique genes, respectively, for wBm and wDi. COG analysis showed that the total number of genes in each COG category was similar in the two organisms. Analysis of pairwise protein distances between wDi and wBm in different COG categories indicated that there was significant variation (Kruskal-Wallis P=0.00077) and pairwise Mann-Whitney tests identified 2 of the 14 high-level COG categories as having elevated divergence between the two Wolbachia. The COG categories showing elevated divergence were M (cell wall, membrane, and envelope biogenesis) and S (function unknown).

The relationship between filarial nematodes and their Wolbachia endosymbionts is thought to be a mutualistic symbiosis (62), because extended treatment of infected mammals with tetracycline and other antibiotics results in clearance of the nematodes. The bases of this symbiosis remain unclear. It has been proposed that wBm provides B. malayi with additional sources of critical metabolites such as heme and riboflavin (63). We interrogated the wDi genome to examine the symbiont's biochemical capabilities. C. elegans and other nematodes (including B. malayi, and, on the basis of the genome sequence presented here, D. immitis) are deficient in heme synthesis but wBm has an intact heme pathway () and a CcmB heme exporter, suggesting that it may support its host by providing heme. wBm has a complete pathway from succinyl-CoA to heme (one apparently missing component, HemG, may be substituted by a functional HemY). wDi lacks both HemY and HemG (and the recently described HemJ that can perform the same transformation). This step in the heme pathway is apparently absent in other bacteria, and so this may not indicate a nonfunctional heme synthesis pathway. Further anabolic pathways absent in D. immitis but present in wDi are purine and pyrimidine de novo synthesis (Fig. 2).

Figure 2.

Anabolic pathways in Wolbachia and Dirofilaria. Selected pathways were identified by screening the predicted proteomes with HMM profiles representing each enzyme in the pathway using HMMer (29). The proteomes were hierarchically clustered (77) based on city block distance between the vectors consisting of the best scores (represented as a heat plot) obtained against each profile. A complete prediction of D. immitis metabolic pathways is available online at the Draft Genomes page of the Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/catalog/org_list1.html).

wBm is deficient in folate synthesis because it lacks dihydrofolate reductase and dihydroneopterin aldolase. wDi has both these genes, suggesting that it can use dihydroneopterin as an input to folate metabolism. Wolbachia wMel from Drosophila melanogaster has both these enzymes, and they are variably present in other alphaproteobacteria. Whether this pathway contributes to the nematode symbiosis is unclear, but it does highlight another component of Wolbachia metabolism that may be accessible to drug development. Further wDi gene products that might be exploited as drug targets include nucleic acid synthesis and cell division proteins, such as FtsZ and DnaB, the fatty acid synthesis enzymes FabZ and AcpS, components of the Sec protein secretion system, and, possibly, the peptidoglycan synthesis enzymes of the Mur operon. All these are unique proteins in wDi, do not have counterparts in mammals, and are being developed as antibiotic drug targets for bacterial infections (64–68).

Horizontal gene transfer from Wolbachia to host nuclear genomes is common in animals harboring this endosymbiont (63), and it has been proposed that these transfers may confer new functionality to the nuclear genome (34, 69, 70), although this is unlikely (71). We identified 868 elements, spanning 219 kb, of >100 b with ≥80% identity to wDi. The Wolbachia origin of these elements was supported by clustering based on the frequency distribution patterns (Supplemental Fig. S1) of tetramer palindromes (72). We did not identify the putative complex Wolbachia insertion discussed by Dunning-Hotopp et al. (70) involving the antigen Dg2 gene. We found a version of the Dg2 gene in our predicted transcriptome that contained standard nematode introns, but no evidence of the construct previously described that had the introns of Dg2 largely replaced with sequences that match 100% to the wDi genome. It is likely that this sequence is a laboratory or computational artifact, especially because the construct includes a cloning vector sequence in addition to Wolbachia.

Only 9 of our identified elements matched >80% of the length of a wDi open reading frame and were not interrupted by frame-shifting insertions or deletions or stop codons. Only one of these putative lateral gene transfers had a match to a Wolbachia protein of known function (transcription termination factor, NusB). We found no evidence of transcription of this gene in the male and female RNA-Seq data. We applied the same procedure for finding Wolbachia insertions to the B. malayi genome and identified 654 insertions spanning 327 kb. Only 31 pairs of insertions that were probably derived from homologous Wolbachia genes were found (in both of the two genomes). None of these shared insertions had complete open reading frames. Comparison with the Wolbachia insertions in the partial genomes of A. viteae and O. flexuosa, onchocercid nematodes that have lost their symbionts (34), revealed no insertions shared by all four species. Only 48 insertions were shared by 2 species and 5 were shared by 3. The number of shared fragments was as would be expected from homoplasious, random insertion of Wolbachia fragments independently into their host genomes. If ∼25% of the genome was randomly transferred in all species, the number of shared fragments expected by chance would be ∼45 (0.25×0.25×750 fragments). We thus tentatively conclude that, although elements from wDi have transferred to the nuclear genome, there is no evidence of their functional integration into nematode biology.

DISCUSSION

The D. immitis genome sequence described here is only the second to be determined for an onchocercid nematode, despite the social and economic importance of these parasites. Three genomes were cosequenced: the mitochondrial (at ∼4000-fold read coverage of the 13.6-kb genome; this had been determined previously; ref. 17); the genome of the Wolbachia symbiont wDi (at ∼1000-fold coverage of the 0.9-Mb genome); and the nuclear genome (at ∼150-fold coverage of the estimated 95-Mb genome). We used high-throughput, short-read Illumina technology, stringent quality filtering and optimized assembly methods to derive genomes of good draft quality (73). After redundancy reduction, the span of the nuclear assembly was 84.2 Mb, slightly smaller than the 88.3 Mb assembled for B. malayi (6). Overall, although the number of scaffolds was approximately equivalent, the contiguity of the D. immitis genome assembly was lower than that of B. malayi, because of the availability of long-range scaffolding information for the latter species. The predicted nuclear gene set was much smaller than that of C. elegans, but of a size similar to that of B. malayi. The two onchocercid nematodes also have a lower proportion of species-unique proteins. These two differences may be a feature of the Onchocercidae, because the unpublished L. loa genome has only 15,444 predicted proteins (Filarial Worms Sequencing Project, Broad Institute of Harvard and MIT; http://www.broadinstitute.org/). Another possibility is that the richer analytic environment for C. elegans in particular has permitted the identification of many unique genes using biological evidence (such as transcript information). We will continue to develop and improve the assembly and annotation of D. immitis and wDi as additional tools and biological resources become available.

Two peculiarities of the assembled D. immitis genome are striking: the lack of genetic diversity and the lack of active transposable elements. The lack of diversity was convenient, in that it allowed us to pool data obtained from two different D. immitis isolates, one from Pavia, Italy, and the other from Athens, Georgia, USA. Polymorphisms called from the independent sequencing of the two isolates yielded a per-nucleotide diversity of 0.04%. Both sequenced isolates fall within the single eastern United States population defined by microsatellite analyses (36). The hypovariability may be a result of the recent admixture of European and American heartworm populations through movement of domestic animals or arise from the very recent introduction of heartworm into the New World by Europeans (74). The first report of dirofilariasis in the United States dates from only 1847, as opposed to a 1626 observation from Italy. The lack of genetic diversity in the nuclear genome will make identification of mutations conferring drug resistance much easier. The lack of DNA transposons and active retrotransposons in D. immitis is a strong negative result, because active elements are easy to identify (they are present in multiple, highly similar copies). We identified only fragmented and functionally inactivated segments of Pao-type retrotransposons, similar to those found in and probably still active in B. malayi. To our knowledge, this is the first metazoan genome devoid of active transposable elements. The presence of putatively active Pao elements in B. malayi suggests that their loss was an evolutionary recent event in D. immitis.

The Wolbachia wDi genome, with 823 predicted proteins, complements the D. immitis nuclear genome in that it encodes enzymes for anabolic pathways that are missing in the latter, e.g., biosynthesis of heme, purine, or pyrimidines (Fig. 2). In contrast to wBm, wDi also carries the genes for folate synthesis, suggesting that folate too might be supplied by the endosymbiont. However, essential metabolites could also be taken up from the mammalian or insect host, and so it remains to be shown whether such metabolites are actually delivered from wDi to D. immitis. Analysis of orthology between wBm and wDi revealed that both organisms possess many unique genes (approximately one-third of the total gene complement of each genome). The representation of genes in the different COG categories was similar for wBm and wDi, suggesting that most gene losses occurred before the split of the two lineages or that there have been no biases in gene losses/acquisition after the evolutionary separation. Analysis of protein distances revealed that proteins involved in cell wall/membrane biogenesis (COG category M) displayed more variation between the two organisms compared with the other functional categories. It is reasonable to conclude that the interface between the symbiotic bacterium and the host environment is a place where evolutionary rates are elevated, either as part of an arms race underpinning conflict between the two genomes or as a feature of the dynamic exploitation of the interface in adaptation of the symbiosis. In any case, the endosymbiont, being essential for proliferation of D. immitis, represents a target for control of the heartworm. Screening the predicted wDi proteome returned expected antibiotic drug targets such as Fts and Sec proteins, but also the products of the Mur operon required for peptidoglycan synthesis.

Many of the anthelmintics used in human medicine were originally developed for the veterinary sector. We pursued two approaches to identify potential drug targets in D. immitis: top-down, starting from the known anthelmintic targets of C. elegans (Table 2), and bottom-up, narrowing down the predicted D. immitis proteome to a list of essential, unique, and druggable targets (Table 3). Although the majority of the current anthelmintics activate their target (thereby interfering with synaptic signal transduction), the aim of the second approach was to identify inhibitable targets. The criteria applied—presence of an essential ortholog in C. elegans, absence of any significantly similar protein in human or dog, and absence of paralogs in D. immitis—admittedly missed many of the known anthelmintic targets, e.g., proteins that are not conserved in C. elegans or that possess a mammalian ortholog. The aim of the approach was to maximize the specificity of in silico target prediction at the cost of low sensitivity. Our goal was to end up with a manageable, rather than complete, list of unique D. immitis proteins that are likely to be essential and druggable. Some of the candidates identified are worth further investigation, based on their presumed role in signal transduction, e.g., the nematode-specific G protein-coupled receptors or hedgehog proteins (Table 3). Others have already been validated as drug targets in other systems: sterol-C-24-methyltransferase (EC 2.1.1.41) is a target of sinefungin, chitin synthase (EC 2.4.1.16) is the target of the insecticide lufenuron, and the mannosyltransferase bre-3 is required for interaction of Bacillus thuringiensis toxin with intestinal cells (52). The discovery of new D. immitis drug targets would be timely because resistance to macrocyclic lactones has recently been reported from the southern United States (75).

Table 3

Candidate drug targets, bottom-up search

D. immitis protein	Predicted function	B. malayi ortholog	H. sapiens log10 (E)	C. lupus log10 (E)	C. elegans RNAi
Nucleic acid synthesis and repair
DIMM09370	RNA-dependent RNA polymerase	BM06623	0.28	0.11	Lethal
DIMM23395	Apurinic/apyrimidinic endonuclease	BM17151	>1	−0.12	Lethal
Glycosylation and sugar metabolism
DIMM15580	dTDP-4-dehydrorhamnose 3,5-epimerase	BM18305	0.23	0.04	Lethal
DIMM03355	β-1,4-Mannosyltransferase	BM20353	0.95	0.94	Lethal
DIMM44525	UDP-galactopyranose mutase	BM01820	0.36	0.08	Molt defective
DIMM36945	Chitin synthase	BM18745, BM02779	−3.52	−4.00	Lethal
Lipid metabolism
DIMM52545	Lipase	BM01258, BM03783	0.08	−1.60	Lethal
DIMM13730	Sterol-C24-methyltransferase (Erg11)	BM20515	−3.10	−4.00	Lethal
DIMM28375	Methyltransferase	BM18889	−0.03	−0.15	Lethal
Transport
DIMM21065	Aquaporin	BM04673	−2.05	−0.52	Lethal
Signal transduction
DIMM13570	Nuclear hormone receptor		−2.40	−4.52	Lethal
DIMM11130	G protein-coupled receptor	BM19106	−1.06	−0.59	Lethal
DIMM32415	G protein-coupled receptor		−1.26	−1.57	Lethal
DIMM39455	G protein-coupled receptor		−2.70	−1.96	Lethal
DIMM13630	Groundhog protein		>1	0.04	Lethal
DIMM47150	Warthog protein	BM01098	>1	0.78	Lethal
DIMM03220	Warthog protein	BM01043, BM17326, BM08657	>1	0.32	Lethal
DIMM11410	Haloacid dehalogenase-like hydrolase	BM19541	−3.22	−4.00	Lethal
DIMM13420	Apoptosis regulator CED-9	BM01838	0.77	−1.02	Lethal

Potential drug targets were filtered from the predicted D. immitis proteome using the following criteria: 1) presence of an ortholog in C. elegans that has as an RNAi phenotype lethal, L3_arrest, or molt_defective; 2) absence of a significant BLAST match (E>10−5) in the predicted proteomes of H. sapiens and C. lupus familiaris; and 3) predicted function as an enzyme or receptor.

Filarial nematodes modulate the immune systems of their hosts in complex ways that result in an apparently intact immune system that ignores a large parasite residing, sometimes for decades, in tissues or the bloodstream. They may also require intact immune systems to develop properly (76). Often immune responses result in a pathologic condition for the host in addition to parasite clearance, and Wolbachia may exacerbate these responses (12). We identified a wide range of putative immunomodulatory molecules and, in addition, highlight two D. immitis products that may deflect or distract the host immune response: one similar to SOCS5 and the other similar to IL-18. The host-encoded versions of both of these molecules have been implicated in antifilarial immune responses. Strategies for development of a vaccine against filariases depend on delivering the correct antigens to the right arm of the immune system, avoiding induction of dangerous responses, and deflecting or stopping immune suppression by the parasite. We identified homologs of all the current roster of filarial vaccine candidates in our genome, and these can now be moved rapidly into testing in the dog heartworm model. In addition, we defined a large number of potentially secreted D. immitis proteins that may contribute to the host-parasite interaction and also be accessible to the host immune system.

Onchocercid parasites share not only a fascinating biology involving immune evasion, arthropod vectors, and Wolbachia endosymbionts but also a pressing need for new drugs, improved diagnostic methods, and, ideally, vaccines. We hope that the genome sequence of the heartworm presented here will contribute to an increased understanding of its biology and to new leads for control.