Genomics of Loa loa, a Wolbachia-free filarial parasite of humans.

Loa loa, the African eyeworm, is a major filarial pathogen of humans. Unlike most filariae, L. loa does not contain the obligate intracellular Wolbachia endosymbiont. We describe the 91.4-Mb genome of L. loa and that of the related filarial parasite Wuchereria bancrofti and predict 14,907 L. loa genes on the basis of microfilarial RNA sequencing. By comparing these genomes to that of another filarial parasite, Brugia malayi, and to those of several other nematodes, we demonstrate synteny among filariae but not with nonparasitic nematodes. The L. loa genome encodes many immunologically relevant genes, as well as protein kinases targeted by drugs currently approved for use in humans. Despite lacking Wolbachia, L. loa shows no new metabolic synthesis or transport capabilities compared to other filariae. These results suggest that the role of Wolbachia in filarial biology is more subtle than previously thought and reveal marked differences between parasitic and nonparasitic nematodes.

Introduction

Filarial nematodes dwell within the lymphatics and the subcutaneous tissues of up to 170 million people worldwide and are responsible for notable morbidity, disability and socioeconomic loss[1]. Although eight filarial species infect humans, only five cause significant pathology -- Wuchereria bancrofti, Brugia malayi, and Brugia timori, the causative agents of lymphatic filariasis; Onchocerca volvulus, the causative agent of ‘river blindness’ or onchocerciasis; and Loa loa, the African eyeworm; L. loa affects an estimated 13 million people and causes chronic infection most often characterized by localized angioedema (Calabar swelling) and/or subconjunctival migration of adult worms across the eye (“African eyeworm”). Complications of infection include encephalopathy, entrapment neuropathy, glomerulonephritis and endomyocardial fibrosis[2]; L. loa is restricted geographically to equatorial West and Central Africa, where its deerfly vector (Chrysops spp.) breeds; L. loa microfilariae (L1) are acquired by flies from human blood and subsequently develop into infective larvae (L3) before being reintroduced into a human host during a second blood meal (Supplementary Fig. 1). While L. loa is the least well-studied of the pathogenic filariae, it has been gaining prominence of late because of the severe adverse events (encephalopathy and death) associated with ivermectin treatment[3] in mass drug administration campaigns in West and Central Africa.

L. loa was targeted for genomic sequencing for two reasons. First, in contrast to other pathogenic filariae, L. loa lacks the α-proteobacterial endosymbiont Wolbachia. The obligate nature of Wolbachia symbiosis in W. bancrofti, B. malayi, and O. volvulus has been inferred by studies in which antibiotics (e.g. doxycycline) that target Wolbachia (but not the worm itself) have shown efficacy in treating humans with these infections[4, 5]. Through genomic analysis, Wolbachia have been hypothesized to provide essential metabolic supplementation to their filarial hosts[6, 7]. The absence of the Wolbachia endosymbiont in L. loa suggests that either there has been lateral transfer of important bacterially-encoded genes or that the obligate relationship between the endosymbiont and its filarial host is dispensable, at least under certain circumstances. Understanding the comparable adaptations of L. loa was considered essential to gain insight into the potential impact of the endosymbiont[8]. Second, as the most neglected of the pathogenic filariae but one gaining increasingly more clinical prominence, understanding the host parasite relationship as it relates to the severe post-treatment reactions typical of both Wolbachia-containing and Wolbachia-free filarial parasites became of paramount importance.

Thus, we generated a draft genome sequence of L. loa generated and produced a refined gene annotation aided by transcriptional data from L. loa microfilariae. We also generated draft genome sequences of two of the most pathogenic (and Wolbachia-containing) filarial species, W. bancrofti and O. volvulus. This approach enabled us to define more comprehensively the genomic differences between L. loa and other filarial parasites.

Genome assemblies and repeat content

The nuclear genome of Loa loa consists of five autosomes plus a sex chromosome. Using 454 whole genome shotgun sequencing, L. loa was sequenced to 20x coverage and assembled into 5774 scaffolds with an N50 of 172 Kb and total size of 91.4 Mb (Table 1). The W. bancrofti and O. volvulus genomes derived from single adult worms (an unsexed juvenile adult worm for W. bancrofti and an adult male worm for O. volvulus) were sequenced to 12x and 5x coverage, respectively (Table 1). Because of the low coverage of the O. volvulus genome, it was not included in further analysis. While the assembly sizes of the L. loa and B. malayi genomes are comparable (91.4 Mb versus 93.7 Mb), the scaffold N50 of the L. loa genome is almost twice that of B. malayi, making the L. loa genome assembly the most contiguous assembly of any filarial nematode to date. The filarial genomes differ widely in repeat content (Table 1, Supplementary Tables 1–14, Supplementary Note) with the L. loa genome being more repetitive than W. bancrofti but less than B. malayi.

Table 1

Genome features of filarial worms and their Wolbachia endosymbionts; Wolbachia genome abbreviations stand for Wolbachia of B. malayi (wBm), Wolbachia of W. bancrofti (wWb), and Wolbachia of O. volvulus (wOv).

Organism	Coverage	Sequence(Mb)	Scaffolds	Scaffold N50(Kb)	% GC	% Repetitive	% LowComplexity	# of Genes
L. loa	20x	91.4	5774	172	31.0	9.3	1.7	14907a
W. bancrofti	12x	81.5	25884	5.16	29.7	6.2	3.9	19327a
O. volvulus	5x	26.0	22675	1.27	32.5	--	--	--
B. malayi	9x	93.7	8180	94	30.2	12.1	1.1	18348
wBm	11x	1.08	1	--	34.2	--	--	805
wWb	2x	1.05	763	1.62	34.0	--	--	--
wOv	2x	0.44	341	1.51	32.8	--	--	--

Due to fragmentation of the genome assemblies, the true W. bancrofti gene count is estimated to be 14,496–15,075 genes, while the true L. loa gene count is estimated to be 14,261 genes (Supplementary Note).

As nuclear Wolbachia transfers (nuwts) have been identified in all Wolbachia-colonized and Wolbachia-free filarial nematodes examined[9], we expected to find similar transfers in the L. loa genome. However, a BLAST-based search of the assembled L. loa genome did not reveal any large transfers of Wolbachia DNA. A more sensitive read-based analysis determined that the L. loa genome does not have any large, recent transfers (> 500 bp, Supplementary Note). It does however have small, presumably “older”, transfers supporting the hypothesis that L. loa was once colonized by Wolbachia but subsequently lost its endosymbiont (Supplementary Table 15 and Supplementary Fig. 2). Of the transfers definitively of Wolbachia ancestry, and not of possible mitochondrial ancestry, there is no evidence that these are functional in L. loa (Supplementary Note).

Gene content and synteny

Initial gene sets were produced for both L. loa and W. bancrofti based on a combination of gene predictors with refinements to the L. loa annotation based on RNA sequence (RNA-Seq) data (see Methods). The final L. loa gene set contains 14,907 genes, 70% of which are supported by RNA-Seq (Table 1, Supplementary Tables 16, 17). The W. bancrofti genome is predicted to encode 19,327 genes (Table 1, Supplementary Note). The filarial genomes show a high degree of synteny (Figure 1), with 40% and 13% of L. loa genes being syntenic relative to B. malayi and W. bancrofti, respectively. Nearly all syntenic breaks between filarial genomes occur at scaffold ends (Supplementary Fig. 3B), suggesting the percentage is limited by assembly contiguity and the true level of synteny is much higher. When the L. loa genome is compared to that of C. elegans, orthologs from a single L. loa scaffold map predominantly to a single C. elegans chromosome (Figure 1). However, only 2% of all L. loa genes were syntenic relative to C. elegans (Supplementary Fig. 3A), supporting the hypothesis that genome rearrangements during filarial evolution were mostly intra-chromosomal[7]. An intermediate level of synteny (12%) is seen between L. loa and related non-filarial parasite Ascaris suum (Supplementary Fig. 3A).

Figure 1

Synteny between filarial worms and C. elegans. (A) Gene distribution on the C. elegans genome. Black-edged vertical bars represent each one of the six C. elegans chromosomes, labeled accordingly in the bottommost part of the panel. Horizontal colored boxes within each bar indicate the location and strand of C. elegans genes (leftmost column = plus strand; rightmost column = minus strand). The color of the lines designates each chromosome and serves as a color-based legend for panel B; (B) Gene distribution on the twelve longest L. loa scaffolds. Scaffolds are represented by black-edged horizontal bars, and identified by labels on the left. Vertical colored boxes indicate the position and strand of each gene (uppermost column = plus strand; bottommost column = minus strand). The color-coding indicates the chromosome where each respective ortholog in C. elegans is located. Grey colored boxes represent either genes without orthologs in C. elegans or genes with two or more homologs in distinct C. elegans chromosomes; (C) Distribution of L. loa scaffold 4 orthologs on the C. elegans W. bancrofti and B. malayi genomes. The scaffolds and chromosomes with best matches to L. loa scaffold 4 based on whole genome alignment are depicted here. Each row contains one or more black-edged horizontal bars representing either chromosomes (C. elegans) or scaffolds (L. loa B. malayi and W. bancrofti) from each sequenced genome. Purple boxes indicate position and strand of genes. Grey projections connect orthologous genes across organisms.

More than half of the genes encoded by the L. loa and W. bancrofti genomes could be assigned functional categories, PFAM domains, Gene Ontology (GO) terms, and/or Enzyme Commission (EC) numbers (Supplementary Fig. 4, Supplementary Tables 16, 18). Relative to other filarial genomes, the L. loa genome is enriched (p < 0.05, Fisher’s exact test) for numerous domains, including pyridoxamine 5’-phosphate oxidases that synthesize vitamin B (Figure 2; see below). The L. loa genome is also enriched for numerous chemoreceptors suggesting that L. loa may be capable of more complex interactions with its host environment than other filarial worms (see below; Supplementary Note). An RNA helicase domain involved in viral DNA replication is enriched in the L. loa genome; this domain was likely horizontally transferred to the L. loa genome from cyclovirus infection (Supplementary Note). While not statistically significant, the L. loa genome encodes many more hyaluronidases (6) than the B. malayi or W. bancrofti genomes (2 each). Hyaluronidases are often involved in tissue penetration and could allow L. loa to move more readily through human host tissue, as L. loa adults are highly mobile whereas B. malayi and W. bancrofti adults are commonly tethered to the lymphatic endothelium.

Figure 2

Enriched and depleted PFAM and TIGRfam domains in each filarial genome relative to the other two. All domains significantly enriched (p < 0.05, Fisher’s exact test) are shown; red indicates enriched while blue indicates depleted. Numbers of identified domains are given in each box. Broad functional categories representing each domain are shown to the left.

The genome of W. bancrofti is enriched (p < 0.05, Fisher’s exact test) for genes with domains related to cellular adhesion and the extracellular matrix (e.g., cadherins, laminins, and fibronectins). Whether these are important in mediating the fibrosis associated with lymphatic filarial disease (e.g. elephantiasis, lymphedema) in W. bancrofti infection[10] or with establishing an anatomical niche within the afferent lymphatics where the adults reside awaits clarification.

Gene products associated with immunologic responses

Each of the filarial parasites interacts with both its definitive mammalian and intermediate arthropod hosts (Chrysops spp. in the case of L. loa) during its life cycle (Supplementary Fig. 1). The parasite is thought not only to have its own innate immune system to protect itself from microbial pathogens, but also to have evolved mechanisms to exploit and/or subvert host and vector defense mechanisms. Although adaptive immune molecules such as immunoglobulins or Toll-like receptors (TLRs) are absent in L. loa and other filarial nematodes, L. loa, like other filariae, appears to possess a primordial Toll-related pathway (Supplementary Table 19, Supplementary Note). The innate immune system encoded by the L. loa genome also includes C-type lectins, galectins, jacalins, and scavenger receptors; L. loa contains a number of lipopolysaccharide binding proteins, proteins implicated in modulating the effects of host bacteria or microbial translocation products. Like B. malayi, the L. loa (and W. bancrofti) genomes do not encode antibacterial peptides described in C. elegans and Ascaris suum[7], suggesting that these molecules are either dispensable in filariae or too divergent to detect.

Analysis of L. loa genes identified a number of human cytokine and chemokine mimics/antagonists including: genes encoding macrophage migration inhibition (MIF) family signaling molecules, transforming growth factor-β and their receptors, members of the interleukin (IL)-16 family, an IL-5 receptor antagonist, an interferon regulatory factor, a homolog of suppressor of cytokine signaling (SOCS)-7, and two members of the chemokine-like family (Supplementary Table 19). In addition, the L. loa genome encodes 17 serpins and 7 cystatins shown to interfere with antigen processing and presentation to T cells[11], 2 indoleamine 2,3 dioxygenase (IDO) genes that encode immunomodulatory proteins implicated in strategies of immune subversion, and a number of the Wnt family of developmental regulators that typically modulate immune activation. The L. loa genome encodes proteins that have sequences similar to human autoantigens (Supplementary Note). Although some of these can also be found in the other filariae, the slight expansion of these in L. loa suggests that antibodies induced by L. loa infection may be more autoreactive than those induced by other parasites.

Protein kinases

In addition to elucidating host-pathogen interactions, pathogen genomes can be evaluated for potential drug targets such as protein kinases. We therefore annotated protein kinases in the L. loa genome and compared them to those in other nematode genomes (Supplementary Tables 20–23, Supplementary Fig. 5). Numerous differences were observed between filarial and non-parasitic nematode kinases, particularly regarding meiosis. The widely-conserved TTK kinase (Mps1), which plays a key role in eukaryotic meiosis[12], is present in L. loa and absent in C. elegans. By contrast, filarial nematodes lack the nearly universally conserved RAD53-family kinase chk-2, which is present in C. elegans. In most eukaryotes RAD53 is involved in initiating cell-cycle arrest when DNA damage is detected but in C. elegans is essential for chromosome synapsis and nuclear rearrangement during meiosis[13]. This reciprocal difference suggests that meiosis in filarial parasites may be regulated in a manner more like typical eukaryotes than C. elegans (Supplementary Note). Six L. loa protein kinases are orthologous to targets of drugs currently approved for humans (Supplementary Table 23), including the tyrosine kinase inhibitor imatinib that has been shown to kill schistosomes [14] and Brugia parasites of all stages at concentrations ranging from 5 to 50 uM (unpublished data). Therefore repurposing already approved drugs that target these kinases may be promising in treating filarial (and other helminth) infections[15].

Nematode phylogenomics

To examine the evolution of filarial parasites in the context of other nematodes, we estimated a phylogeny from 921 single copy core orthologs across nine nematode genomes using maximum likelihood, parsimony, and Bayesian methods. All methods converged on a single topology with 100% support (either bootstrap values or posterior probabilities) at all nodes (Figure 3). This phylogeny indicates that Meloidogyne hapla occupies a position basal to a clade of Rhabdtina (i.e; C. elegans, C. briggsae, and Pristionchus pacificus) and to the Spirurina (i.e. filarial worms and Ascaris suum). While these results contrast with previous studies based on ribosomal subunits that placed M. hapla closer to Rhabdtina than the filarial worms[16, 17], our analysis utilized a larger gene set and has higher nodal support values.

Figure 3

Phylogenomic analysis of nematodes. Maximum likelihood, parsimony, and Bayesian methods all estimated an identical phylogeny using the concatenated protein sequences of 921 single copy orthologs. Adjacent to each node are likelihood bootstrap support values/parsimony bootstrap support values/Bayesian posterior probabilities. The distribution of genes in ortholog clusters is shown to the right of the phylogeny. Core genes are encoded by all genomes, shared genes are encoded by at least two but fewer than all genomes, and unique genes are found only in one genome. Orthologs specific to the non-parasitic nematodes (C. elegans C. briggsae, and P. pacificus) and filarial nematodes are also highlighted. Of the 6,280 L. loa genes with no functional assignment, 3,665 are unique to L. loa and 1,158 are filaria-specific.

Relative to the genomes of non-parasitic nematodes, numerous orthologs were identified as unique to the filarial parasites (Figure 3). Proteins encoded by the filarial genomes show enrichment of immunogenic domains such as extracellular and cell adhesion domains, and in a metabolic context are enriched for trehalase domains involved in trehalose degradation (q < 0.05, Fisher’s exact test; Supplementary Fig. 6). Trehalose is known to be involved in the protection of nematodes from environmental stress[18] and could potentially play a key role in filarial survival. Trehalose levels have also been associated with increased lifespan in C. elegans[19] and might support the idea that an increased use of trehalose by filarial nematodes could be related to their relatively long lifespan.

The filarial genomes lack a wide array of seven transmembrane G-protein-coupled chemoreceptors (7 TM GPCRs; Supplementary Fig. 6). Profiling of 7 TM GPCRs revealed a pattern of progressive loss of many families in the transition from non-parasitic to parasitic lifestyles (Figure 4). For example, filarial nematodes and T. spiralis completely lack the Str superfamily, including odr-10, known to be involved in detection of volatiles[20], and Kin-29, a protein kinase that regulates Str expression in C. elegans[21]. If the Str superfamily is more broadly involved in odorant detection, this could explain why these molecules are lacking in filarial nematodes and T. spiralis parasites that only live in aqueous environments, while they are retained in A. suum and M. hapla, which are exposed to volatiles in part of their life cycle. Only the Srab, Srx, Srsx, and Srw families were conserved across all nematodes, suggesting that these GPCRs mediate vital nematode functions.

Figure 4

Phylogenetic profile of chemoreceptors in nematode genomes. Both the 7 transmembrane domain G-protein-coupled chemoreceptors (GPCRs) and guanylate and adenylate cyclases are shown.

Filarial genomes are also depleted in both soluble and receptor guanylate cyclases – these genes are involved in the regulation of environmental sensing and complex sensory integration functions (Figure 4). However, gcy-35 and gcy-36, which are involved in the detection of molecular oxygen in solution[22], are encoded in the filarial genomes. Protein kinase profiling revealed 18 receptor guanylate cyclases which are present in C. elegans but not in filarial worms, including environmental sensors gcy-14 and gcy-22 (Supplementary Table 23). Depletion of these and other kinases involved in olfactory and gustatory sensing, including kin-29, suggests that the environments of filarial nematodes are less complex in terms of chemosensory inputs than are those inhabited by non-parasitic nematodes (Supplementary Note). The L. loa genome does, however, encode significantly more chemoreceptors than do other filarial nematodes (p < 0.05, Fisher’s exact test), which may be related to the increased mobility of L. loa adult worms.

Phylogenetic profiling of metabolism

Previous genomic analysis identified five biosynthetic pathways (heme, riboflavin, FAD, glutathione and nucleotides) present in Wolbachia but missing from its relatives, e.g; Rickettsia. These Wolbachia-encoded pathways were hypothesized to provide metabolites needed by their filarial hosts[6]. As L. loa lacks Wolbachia, it was theorized that the L. loa genome must encode genes to replace these pathways, potentially laterally transferred from Wolbachia to an ancestor of L. loa. However, no transfers relating to these metabolic functions were apparent (see above). Thus, we generated complete metabolic pathway reconstructions for nine nematode and four Wolbachia genomes (Table 2, Supplementary Tables 24, 25) to determine how L. loa acquires these metabolites and placed the results in an evolutionary context. Unexpectedly, none of the five “complementary” pathways differ between L. loa and the other filarial nematodes, calling into question the role of these pathways in filarial-Wolbachia symbiosis.

Table 2

Phylogenetic profiles of biosynthesis pathways hypothesized to be involved in the Filaria-Wolbachia symbiosis. Conservation each pathway across nematodes and Wolbachia are shown in Supplementary Tables 24 and 25, respectively. Pathways are labeled as complete (+), partial (+/–), or absent (–); Wolbachia genome abbreviations stand for Wolbachia of B. malayi (wBm), Wolbachia of D. melanogaster (wMel), Wolbachia of C. pipiens (wPip) and Wolbachia of W. bancrofti (wWb).

BiosynthesisPathway	C. elegans	C. briggsae	P. pacificus	M. hapla	T. spiralis	A. suum	B. malayi	W. bancrofti	L. loa	wBm	wMel	wPip	wWb
Heme	–	–	–	–	–	–	–a	–a	–a	+	+	+	+
Riboflavin	–	–	–	–	–	–	–	–	–	+	+	+	+
FAD	+	+	+/–	+/–	+	+	+	+	+	+	+	+	+/–
Glutathione	+	+	+	+	+	+	+	+	+	+	+	+	+
Purines	+	+	–	–	+/–	+	–	–	–	+	+	+	+
Pyrimidines	+	+	+	+	+/–	+	+/–	+/–	+/–	+	+	+	+

all filarial worms encode a ferrochelatase, the last enzyme in heme synthesis (Supplementary Note).

Furthermore, in only two pathways (heme and nucleotide synthesis) did the filarial genomes differ from the other nematodes. The FAD and glutathione pathways are complete in all nematode genomes, while the riboflavin pathway is missing from all nematode genomes. The heme biosynthesis pathway, previously reported to be absent in B. malayi[7], is missing from not only the filarial worms but from all nematode genomes characterized to date. Experimental work on C. elegans (also Wolbachia-free) has shown that it cannot synthesize heme de novo[23]; B. malayi has been previously noted as having a single member of the heme synthesis pathway, ferrochelatase (an enzyme that catalyzes the last step in heme synthesis[7] [Supplementary Note]). The gene encoding ferrochelatase is also present in the L. loa and W. bancrofti genomes, but absent in all other nematode genomes, including A. suum. It is possible that this gene in filarial nematodes is not involved in heme synthesis, but rather in an alternate, unknown pathway.

Like B. malayi, both L. loa and W. bancrofti lack the ability to synthesize nucleotides de novo. All three filarial genomes lack the majority of the proteins involved in the purine synthesis pathway as well as the first enzyme involved in the pyrimidine synthesis pathway (Table 2, Supplementary Table 24). Other nematodes have also lost portions of these pathways; the purine synthesis pathway has been largely lost in P. pacificus and M. hapla, while the first two enzymes in the pyrimidine synthesis pathway have been lost in T. spiralis. These multiple and likely independent losses could underscore a general flexibility in the need for de novo nucleotide synthesis in nematodes. All nematodes, including the filariae, have complete sets of purine and pyrimidine interconversion pathways (Supplementary Table 24), implying that they could generate all necessary nucleotides from a single purine and pyrimidine source, a concept supported by experimental data in B. malayi[24]. Filarial genomes do encode two purine-specific 5’ nucleotidases for salvage whereas all other nematodes encode only one; the extra copy in the filariae appears to arise from a single gene duplication event and diverged significantly from the ancestral gene (Supplementary Fig. 7). Additionally, we profiled known nematode and Wolbachia transporters linked to to these pathways, and found no evidence of differences between filarial and non-filarial nematodes nor among Wolbachia endosymbionts (Supplementary Note and Supplementary Fig. 8). Given the uniformity of these pathways across nematodes and the apparent lack of any related transfers of Wolbachia DNA to the L. loa genome, it is likely that the symbiotic role played by Wolbachia in filarial nematodes either lies outside these pathways or involves more subtle metabolic supplementation rather than the wholesale provision of unproduced metabolites.

The only metabolic pathway found to differ in gene content between L. loa and other nematodes with sequenced genomes is vitamin B6 synthesis and salvage. Most nematode genomes encode single copies of the two enzymes involved in vitamin B6 salvage, while the L. loa genome encodes five copies of the second enzyme, pyridoxal 5’-phosphate synthase (Supplementary Note). This pathway also differed among Wolbachia genomes. While both of the insect Wolbachia genomes also encoded two genes involved in the synthesis of vitamin B6 (pdxJ and pdxK), neither of the filarial Wolbachia genomes did (the difference between wBm and wMel having been noted previously[6]). If the filarial Wolbachia endosymbionts need to acquire vitamin B6 exogenously, this could explain a metabolic need of Wolbachia fulfilled by the nematode. However, with that hypothesis in mind, it is unclear why L. loa, the one pathogenic filarial nematode without Wolbachia, would encode a greater number of vitamin B6 salvage genes than either B. malayi or W. bancrofti. We could not exclude differences in pyroxidine transporters as no orthologs of known transporters could be identified in either nematode or Wolbachia genomes (Supplementary Note).

Conclusion

The study of some of the nematode genomes has already provided great insight into the genomic structure, biology and evolution of this major division of nematode parasites. With the release of the genome of L. loa, a human pathogen and parasitic nematode that does not contain Wolbachia, we have been able provide surprising insights into the dispensability of this endosymbiont that deepen the mystery surrounding the “essential nature” of Wolbachia for many filarial worms.

Through large-scale genomic comparisons within the phylum Nematoda, we have not only been able to define molecules and pathways that are either Loa loa- or filaria-specific but also, by comparison with the non-parasitic nematodes (e.g; C. elegans), we have gained a first glimpse into the nature of parasitism itself. Moreover, this effort has identified new targets for intervention that should aid programs aimed at control and elimination of these important but neglected parasites.

Online Methods

Sequencing and Assembly

For L. loa, 5×105 microfilariae (mf) were purified during a therapeutic apheresis from a patient with loiasis infected in Cameroon seen at the NIH under protocol 88-I-83 (NCT00001230). A single unfertilized adult W. bancrofti worm was obtained under ultrasonic guidance (as part of protocol NCT00339417) in Tienegabougou, Mali. A single adult O. volvulus male was isolated from a surgically-removed subcutaneous nodule in Ecuador after collagenase digestion. Genomic DNA for all samples was prepared using the Qiagen genomic DNA kit (Qiagen, Gaithersburg, MD). DNA obtained from W. bancrofti and O. volvulus was amplified using the Qiagen Repli-g midi kit. For L. loa, W. bancrofti, and O. volvulus, approximately 50, 10, and 5 ug of DNA was used for genomic sequencing, respectively. For L. loa, 454 shotgun fragment and 3Kb jumping sequencing libraries were prepared and sequenced as previously described[25]. Only fragment libraries were constructed for W. bancrofti and O. volvulus. Assemblies were then generated using Newbler version 2.1 (Roche 454 Life Sciences). Given the overall low coverage of the W. bancrofti and O. volvulus assemblies (5–12x), no bias normalization was done for the whole-genome amplified sequence data. Also for the W. bancrofti and O. volvulus assemblies, contigs were screened by BLASTing against Genbank’s non-redundant nucleotide database (NT) using a cutoff of 1e-25 and minimum match length of 100 bp, and all contigs where the top match was to Wolbachia were removed. Any contigs remaining in the nematode assembly that had secondary matches to Wolbachia were screened manually to ensure that no large chimeric contigs had been generated and retained. Unassembled reads were also screened for Wolbachia sequence using the same BLAST parameters and database. Unassembled reads identified as Wolbachia, along with reads underlying the contigs identified as Wolbachia, were assembled together using Newbler version 2.1 to generate the Wolbachia genome assemblies.

Repeat Content Analysis

Repeat content was identified using RepeatScout[26] followed by RepeatMasker using both nematode repeats from RepBase v17.06[27] and the output from RepeatScout. Only hits with Smith-Waterman score above 250 were maintained. Additional repeats were then identified based on abnormally high read coverage in the genome assemblies, utilizing genome sequence scanning with hysteresis triggering. Positions with read depth 20 times the mode of the read depth distribution switched the “collapsed reads” state to ON during the scanning process, while positions with read depth lower than 10 times the mode switched the “collapsed reads” state to OFF. Only regions longer than 100 nucleotides were reported. Read mapping was performed by runMapping application of the Newbler suite[28]. The output was converted to SAM file format by seq.Newbler2SAM option of GLU package. Only the best alignment of each read was kept. Read depth was calculated by the genomeCoverageBed program of BEDTools suite[29].

RNA Sequencing

RNA was prepared from one million L. loa mf purified from the blood of a patient. Under liquid nitrogen, mf were disrupted by a stainless steel piston apparatus. Total RNA was extracted using the RNeasy Kit (Qiagen, Valencia, CA, USA). A non-strand-specific cDNA library for Illumina paired-end sequencing was prepared from ~37 ng of total RNA as previously described[30] with the following modifications. RNA was treated with Turbo DNase (Ambion, TX) and fragmented by heating at 80 °C for 3 min in 1x fragmentation buffer (Affymetrix, CA) before cDNA synthesis. Sequencing adapter ligation was performed using 4,000 units T4 DNA ligase (New England Biolabs, MA) at 16°C overnight. Following adapter ligation the resulting library was cleaned, size selected twice using 0.7x volumes of Ampure beads (Beckman Coulter Genomics, MA), enriched using 18 cycles of PCR, and cleaned using 0.7x volumes of Ampure beads (Beckman Coulter Genomics, MA). The resulting Illumina sequencing library was sequenced with 76 base paired-end reads on an Illumina GAII instrument (v1.8 analysis pipeline) following the manufacturer’s recommendations (Illumina, CA).

Identification of transfers (nuwts)

An initial search of the Wolbachia of B. malayi genome against the L. loa genome was done using BLASTN with a cutoff of 1e−5. After this assembly-based search, nuclear Wolbachia transfers (nuwts) were identified through a screen of the L. loa sequencing reads as being >80% identical to Wolbachia sequences over 50% of the read. Searches were refined to examine reads with >50 bp match to Wolbachia and were manually curated to remove spurious matches that had a nematode ancestry. Reads matching the bacterial rRNA were removed as they could arise from any bacterial genome that might be contaminating the sample. Regions of homology <50 bp were included if they were detected through analysis of an adjacent region with homology over >50 bp. All of the reads containing nuwts were mapped back to the L. loa genome to identify the consensus sequence and the relationship confirmed using BLASTN to NT. Phylogenetic analysis was conducted on nucleotide sequences of predicted nuwts using RAxML[30].

Annotation

Genes for both L. loa and W. bancrofti were predicted using a combination of ab initio gene prediction tools as previously described[31]. We also used TBLASTN to search the genome assembly against protein sequences of the following species: C. elegans, C. briggsae, Schistosoma mansoni, S. japonicum, and B. malayi (downloaded from Genbank on February 16, 2010). The top BLAST hits are used to construct GeneWise[32] gene models. In addition, we generated gene models using available EST data from L. Loa, W. bancrofti, O. volvulus, and B. malayi (downloaded from GenBank on December 2, 2009). All of these models were used as input into EVM[33] to generate combined gene predictions. To incorporate the L. loa RNA-Seq data, we aligned all RNA-Seq reads to the L. loa genome using BLAT[34]. Next, we use the Inchworm module of the Trinity package[35] with default settings in genome-guided mode to assemble the reads into EST-like transcripts. These transcripts along with the models from EVM into PASA[33] were used for gene model improvement. Gene sets were subsequently filtered to remove repeats, including genes overlapping rRNA, tRNA, or output from RepeatScout[26] or TransposonPSI. Every annotated gene was given a locus ID of the form LOAG_##### (L. loa) or WUBG_##### (W. bancrofti). Pfam domains within each gene were identified using Hmmer3[36], while gene ontology terms were assigned using BLAST2GO[37]. Secretion signals and transmembrane domains were identified using SignalP 4.0[38] and TmHmm[39], respectively. Core eukaryotic genes were identified using CEGMA[40].

Identification of Fragmented Genes

Fragmented W. bancrofti genes were associated to their putative intact orthologs in L. loa or B. malayi by unidirectional BLAST of W. bancrofti peptides against peptides from the reference genome (L. loa or B. malayi); W. bancrofti proteins with less than 80% similarity to the reference, based on query length, and an e-value higher than 1×10−10 were disregarded. A gene was considered fragmented if its length in W. bancrofti was at least 50% shorter than its respective ortholog. The number of reference genome orthologs with multiple assigned fragments in W. bancrofti was then used to extrapolate a corrected gene count for W. bancrofti. An identical analysis was done for L. loa genes by comparison to B. malayi.

Synteny Analysis

Whole genome alignments of C. elegans, B. malayi and W. bancrofti against L. loa were performed by progressive Mauve[41] with default parameters. The extent of the alignment between a pair of sequences was defined as the length spanning all their respective colinear blocks. For each comparison, chromosomes or scaffolds having the longest alignment against L. loa scaffold number 4 (100 scaffolds from W. bancrofti, 30 scaffolds from B. malayi and C. elegans chromosome III) were selected for visualization. For the systematic evaluation of synteny, pairwise syntenic blocks between the genomes of L. loa, C. elegans, B. malayi and W. bancrofti were defined by DAGchainer[42] with minimum numbers of colinear genes set to 3.

Gene Clustering and Phylogenetic Analysis

We built a comparative set of genomes including those sequenced in this study and Pristionchus pacificus (from www.pristionchus.org), Caenorhabditis elegans (release 224 from WormBase), C. briggsae (CAAC00000000), Meloidgyne hapla (from www.pngg.org), Brugia malayi (release 230 from WormBase), Ascaris suum (published release from WormBase), and Trichinella spiralis (ABIR00000000). Genes were clustered using OrthoMCL with a Markov inflation index of 1.5 and a maximum e-value of 1×10−5 [43]. Amino acid sequences of orthologs present as single copies in all genomes were aligned using MUSCLE[44] and concatenated. We then estimated phylogenies from this dataset using three methods. Parsimony bootstrapping analysis was conducted using PAUP[45] using unweighted characters and 1,000 bootstrap replicates. For maximum likelihood analysis, we first selected the JG model[46] using ModelGenerator[47], and then used the PROTCATJG model in RAxML[30] with 1,000 bootstrap replicates. For Bayesian analysis, we used MrBayes[48] with a mixed amino acid model and gamma distributed rates. We ran the analysis with 1 chain for 1 million generations, sampling every 500 generations, and discarding the first 25% of samples as burn-in. Enrichment analyses were conducted using Fisher’s exact test and multiple comparisons were corrected using false discovery rate[49].

Kinase Classification

Initial sets of protein kinases were identified by orthology with annotated C. elegans kinases. Kinases without orthologs were identified in a search of the proteome against a protein kinase HMM derived from an alignment of Dictyostelium protein kinases[50] using a cutoff score of −66. Low-scoring sequences were additionally screened for conservation of known protein kinase sequence motifs. All protein kinases were classified using a controlled vocabulary[51, 52], and classifications of filarial kinases with C. elegans orthologs were mapped from the curated set from the KinBase database. Kinases without orthologs in C. elegans were searched against the curated set using BLAST and classified if the top three hits agreed. Orthology across all nematodes was then used to identify potentially missed kinases and ensure consistent classification.

Metabolic Reconstruction

In addition to the nine nematode genomes listed above, we utilized three additional Wolbachia genomes from B. malayi (AE017321), Drosophila melanogaster (AE017196) and Culex pipiens (AM999887). Metabolic pathways were characterized using Pathway Tools[53]. Metabolic reconstruction was performed using EFICAz2[54] to assign Enzyme Commission (EC) numbers for each enzyme. EC numbers and gene names were used as input to the Pathologic software[55] with transport-identification-parser and pathway-hole-filler options set to assign MetaCyc[56] pathways for each organism. The full set of metabolic pathways for each genome are available at the WormCyc database.