The transcriptome of the invasive eel swimbladder nematode parasite Anguillicola crassus

Nematode samples, RNA extraction, cDNA synthesis and Sequencing

Raw sequencing reads are archived under study-accession number SRP010313 in the NCBI Sequence Read Archive (SRA; http://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP010313) [24]. All samples were sequenced using the FLX Titanium chemistry, except for the Taiwanese female sample T1, which was sequenced using FLX standard chemistry, to generate between 99,000 and 209,000 raw reads per sample. For the L2 library, which had a larger number of non-A. crassus, non-Anguilla reads, we confirmed that these data were not laboratory contaminants by screening Roche 454 data produced on the same run in independent sequencing lanes.

Trimming, quality control and assembly

Raw sequences were extracted in FASTA format (with the corresponding qualities files) using sffinfo (Roche/454) and screened for MINT adapter sequences using cross-match [25] (with parameters -minscore 20 -minmatch 10). Seqclean [26] was used to identify and remove poly-A-tails, low quality, low complexity and short (<100 base) sequences. All reads were compared to a set of screening databases using BLAST [27] (expect value cutoff E <1e-5, low complexity filtering turned off: -F F). The databases used were (a) a host sequence database comprising an assembly of the An. japonica Roche 454 data, a unpublished assembly of An. anguilla Sanger dideoxy sequenced expressed sequence tags (made available to us by Gordon Cramb, University of St Andrews) and transcripts from EeelBase [28], a publicly available transcriptome database for the European eel; (b) a database of ribosomal RNA (rRNA) sequences from eel species derived from our Roche 454 data and EMBL-Bank; and (c) a database of rRNA sequences identified in our A. crassus data by comparing the reads to known nematode rRNAs from EMBL-Bank. This last database notably also contained cobiont rRNA sequences. Reads with matches to one of these databases over more than 80% of their length and with greater than 95% identity were removed from the dataset. Screening and trimming information was written back into sff-format using sfffile (Roche 454). The filtered and trimmed data were assembled using the combined assembly approach [23]: Two assemblies were generated, one using Newbler v2.6 [22] (with parameters -cdna -urt), the other using Mira v3.2.1 [29] (with parameters–job=denovo,est,accurate,454). The resulting two assemblies were combined into one using Cap3 [30] at default settings and contigs were labeled by whether they derived from both assemblies (high confidence assembly; highCA), or one assembly only (lowCA; for a detailed analysis of the assembly categories see the supporting Methods file). The superset of highCA contigs, lowCA contigs and the remaining unassembled reads defines the set of tentatively unique genes (TUGs).

Post-assembly classification and taxonomic assignment of contigs

We rescreened the assembly for host and other contamination by comparing it (using BLAST) to the three databases defined above, and also to NEMBASE4, a nematode transcriptome database derived from whole genome sequencing and EST assemblies [31, 32]. For each contig, the highest-scoring match was recorded, if it spanned more than 50% of the contig. We also compared the contigs to the NCBI non-redundant nucleotide (NCBI-nt) and protein (NCBI-nr) databases, recording the taxonomy of all best matches with expect values better than 1e-05. Sequences with a best hit to non-Metazoans or to Chordata within Metazoa were excluded from further analysis.

Protein prediction and annotation

Protein translations were predicted from the contigs using prot4EST (version 3.0b) [33]. Proteins were predicted either by joining single high scoring segment pairs (HSPs) from a BLAST search of uniref100 [34], or by ESTscan [35], using as training data the Brugia malayi complete proteome [36] back-translated using a codon usage table derived from the BLAST HSPs, or, if the first two methods failed, simply the longest ORF in the contig. For contigs where the protein prediction required insertion or deletion of bases in the original sequence, we also imputed an edited sequence for each affected contig. Annotations with Gene Ontology (GO), Enzyme Commission (EC) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were inferred for these proteins using annot8r (version 1.1.1) [37], using the annotated sequences available in uniref100 [34]. Up to 10 annotations based on a BLAST similarity bitscore cut-off of 55 were obtained for each annotation set. The complete B. malayi proteome (as present in uniref100) and the complete C. elegans proteome (as present in WormBase v.220) were also annotated in the same way. SignalP V4.0 [38] was used to predict signal peptide cleavage sites and signal anchor signatures for the A. crassus transcriptome and for the proteomes of the two model nematodes. InterProScan [39] (command line utility iprscan version 4.6 with options -cli -format raw -iprlookup -seqtype p -goterms) was used to obtain domain annotations for the highCA contigs. We recorded the presence of a lethal RNAi phenotype in the C. elegans ortholog of each TUG using the biomart-interface [40] to WormBase v. 220 using the R package biomaRt [41].

Single nucleotide polymorphism analysis

We mapped the raw reads to the complete set of contigs, replacing imputed sequences for originals where relevant, using ssaha2 (with parameters -kmer 13 -skip 3 -seeds 6 -score 100 -cmatch 10 -ckmer 6 -output sam -best 1) [42]. From the ssaha2 output, pileup files were produced using samtools [43], discarding reads mapping to multiple regions. VarScan [44] (pileup2snp) was used with default parameters on pileup files to output lists of single nucleotide polymorphisms (SNPs) and their locations.

In the 10,496 SNPs thus defined, the ratio of transitions (ti; 6,908) to transversion (tv; 3,588) was 1.93. From the prot4EST predictions, 7,189 of the SNPs were predicted to be inside an ORF, with 2,322 at codon first positions, 1,832 at second positions and 3,035 at third positions. As expected, ti/tv inside ORFs (2.39) was higher than outside ORFs (1.25). The ratio of synonymous polymorphisms per synonymous site to non-synonymous polymorphisms per non-synonymous site in this unfiltered SNP set (dn/ds) was 0.45, rather high compared to other analyses. Roche 454 sequences have well-known systematic errors associated with homopolymeric nucleotide sequences [45], and the effect of exclusion of SNPs in, or close to, homopolymer regions was explored. When SNPs were discarded using different size thresholds for homopolymer runs and proximity thresholds, the ti/tv and in dn/ds ratios changed (Additional file 1: Figure S1). Based on this SNPs associated with a homopolymer run longer than 3 bases within a window of 11 bases (5 bases to the right, 5 to the left) around the SNP were discarded. There was a relationship between TUG dn/ds and TUG coverage, associated with the presence of sites with low abundance minority alleles (less than 7% of the allele calls), suggesting that some of these may be errors. Removing low abundance minority allele SNPs from the set removed this effect (Additional file 1: Figure S2). For enrichment analysis of GO terms associated with positively selected TUGs we used the R package GOstats [46].

Using Samtools [43] (mpileup -u) and Vcftools [47] (view -gcv) we genotyped individual libraries for each of the master list of SNPs. Genotype- calls were accepted at a phred-scaled genotype quality threshold of 10. In addition to the relative heterozygosity (number of homozygous sites/number of heterozygous sites) we used the R package Rhh [48] to calculate internal relatedness [49], homozygosity by locus [50] and standardised heterozygosity [51] from these data. We confirmed the significance of heterozygote-heterozygote correlation by analysing the mean and 95% confidence intervals from 1000 bootstrap replicates estimated for all measurements.

Gene expression analysis

Read-counts were obtained from the bam files generated for genotyping using the R package Rsamtoools [52]. LowCA contigs and contigs with fewer than 32 reads over all libraries were excluded from analysis. Libraries E1 and L2 had very low overall counts and thus we excluded these libraries from analysis. The statistic of Audic and Claverie [53] as implemented in ideg6 [54] was used to contrast single libraries. Differential expression between libraries from male versus female nematodes was accepted for genes that differed in expression values between all the female libraries (E2, T1 and T2; see Table 1) versus the male (M) library (p <0.01), but had no differential expression within any of the female libraries at the same threshold. Differential expression between libraries from nematodes of European An. anguilla and Taiwanese An. japonica origin was accepted for genes that differed in expression values between library E2 and both libraries T1 and T2 (p <0.01), but showed no differences between T1 and T2.

Overrepresentation analyses

The R package annotationDbi [55] was used to obtain a full list of associations (along with higher-level terms) from annot8r annotations prior to analysis of GO term overrepresentation in gene sets selected on the basis of dn/ds or expression values. The R package topGO [56] was used to traverse the annotation graph and analyse each node term for overrepresentation in the focal gene set compared to an appropriate universal gene set (all contigs with dn/ds values or all contigs analysed for gene expression) with the “classic” method and Fisher’s exact test. Terms for which an offspring term was already in the table and no additional counts supported overrepresentation were removed. Mann-Whitney u-tests were used to test the influence of factors on dn/ds values. To investigate multiple contrasts between groups (factors) Nemenyi-Damico-Wolfe-Dunn tests were used, and for overrepresentation of one group (factor) in other groups (factors) Fisher’s exact test was used.

General coding methods

The bulk of analysis (unless otherwise described) presented was carried out in R [57] using custom scripts. For visualisation we used the R packages ggplot2 [58] and VennDiagram [45].

Sampling A. crassus

One female A. crassus and one male A. crassus were sampled from an An. japonica aquaculture with high infection loads in Taiwan, and an additional female was sampled from an An. japonica caught in a stream with low infection pressure adjacent to the aquaculture. A female nematode and pool of L2 were sampled from An. anguilla in the river Rhine, and one female from A. anguilla from a lake in Poland. All adult nematodes were replete with host blood. To assist in downstream filtering of host from nematode reads, we also sampled RNA from the liver of an uninfected Taiwanese An. japonica.

Assembly and post-assembly screening

We screened the complete assembly for remaining host contamination, and identified 3,441 TUGs that had significant, higher similarity to eel (and/or chordate; EMBLBank Chordata proteins) than to nematode sequences [32]. Given the identification of cercozoan ribosomal RNAs in the L2 library, we also screened the complete assembly for contamination with transcripts from other taxa.

1,153 TUGs were found with highest significant similarity to Eukaryota outside of the kingdoms Metazoa, Fungi and Viridiplantae. These contigs matched genes from a wide range of protists from Apicomplexa (mainly Sarcocystidae, 28 hits and Cryptosporidiidae 10 hits), Bacillariophyta (diatoms, mainly Phaeodactylaceae, 41 hits), Phaeophyceae (brown algae, mainly Ectocarpaceae, 180 hits), Stramenopiles (Albuginaceae, 63 hits), Kinetoplasitda (Trypanosomatidae, 26 hits) and Heterolobosea (Vahlkampfidae, 38 hits). Additionally 298 TUGs had best, significant matches to genes from fungi (e.g Ajellomycetaceae, 53 hits) and 585 TUGs had best, significant matches to genes from plants. Outside the Eukaryota there were significant best matches to Bacteria (825 TUGs; mostly to members of the Proteobacteria), Archaea (8 TUGs) and viruses (9 TUGs). No TUGs had significant, best matches to Wolbachia or related Bacteria known as symbionts of nematodes and arthropods. All TUGs with highest similarity to sequences deriving from taxa outside Metazoa were excluded. The final, screened A. crassus assembly has 32,525 TUGs, spanning 12,733,095 bases (of which 11,372 are highCA-derived, and span 6,575,121 bases). All analyses reported below are based on this filtered dataset.

Annotation

For 32,418 of the A. crassus TUGs a protein was predicted using prot4EST [33] (Table 2). An apparently full-length open reading frame (ORF) was identified in 353 TUGs, while for 29,877 the 5’ ends and for 24,277 the 3’ ends were missing. In 13,383 TUGs the corrected sequence with the imputed ORF was slightly changed compared to the raw sequence by insertions or deletions necessary to obtain a continuous reading frame. Using BLAST we determined that 9,556 had significant similarity to C. elegans proteins, 9,664 TUGs matched B. malayi proteins, and 11,620 TUGs had matches in NEMPEP4 [31, 32]. Comparison to the UniProt reference identified 11,115 TUGs with significant similarities. We used annot8r [37] to assign GO terms to 6,511 TUGs, EC numbers for 2,460 TUGs and KEGG pathway annotations for 3,846 TUGs (Table 2). Additionally 5,125 highCA contigs were annotated with GO terms through InterProScan [39]. Nearly one third (6,989) of the A. crassus TUGs were annotated with at least one identifier, and 1,831 had GO, EC and KEGG annotations (Figure 1).

We compared our A. crassus GO annotations for high-level GO-slim terms to the annotations (obtained in the same way) for the complete proteome of the spirurid filarial nematode B. malayi and the complete proteome of C. elegans (Figure 2). The occurrence of GO terms in the annotation of the partial transcriptome of A. crassus was more similar to that of the proteome of B. malayi (0.95; Spearman correlation coefficient) than to the that of the proteome of C. elegans (0.9).

Despite the lack of completeness at the 5’ end suggested by peptide prediction, just over 3% of the TUGs were predicted to be secreted (920 with signal peptide cleavage sites and 65 signal peptides with a transmembrane signature). Again these predictions are more similar to predictions using the same methods for the proteome of B. malayi (742 signal peptide cleavage sites and 41 with transmembrane anchor) than for the proteome of C. elegans (4,273 signal peptide cleavage sites and 154 with transmembrane anchor).

By comparison to RNAi phenotypes for C. elegans genes [59, 60] likely to be orthologous to A. crassus TUGs, 6,029 TUGs were inferred to be essential (RNAi lethal phenotype in C. elegans).

Similar patterns were observed for conservation assessed at different stringency, and when assessed across all TUGs, except that a higher proportion of all TUGs were apparently unique to A. crassus.

Proteins predicted to be restricted to Nematoda and novel in A. crassus were significantly enriched in signal peptide annotation compared to conserved proteins, proteins novel in Metazoa and novel in Spirurina (Fisher’s exact test p <0.001 ; Figure 3). The proportion of lethal RNAi phenotypes was significantly higher for C. elegans presumed orthologs of TUGs conserved across kingdoms (97.23%) than for orthologs of TUGs not conserved across kingdoms (94.59%; p <0.001, Fisher’s exact test).

Identification and analysis of single nucleotide polymorphisms

Single nucleotide polymorphisms (SNPs) were called using VARScan [44] on the 1,100,522 bases of TUGs that had coverage of more than 8-fold available. SNPs predicted to have more than 2 alleles, or that mapped to an undetermined (N) base were excluded, as were SNP likely to be due to base calling errors close to homopolymer tracts and SNP calls resulting from apparent rare variants.

Signal peptide containing proteins have been shown to have higher rates of evolution than cytosolic proteins in a number of nematode species. A. crassus TUGs predicted to contain signal peptide cleavage sites showed a non-significant trend towards higher dn/ds values than TUGs without signal peptide cleavage sites (p=0.22; two sided Mann-Whitney-test). Orthologs of C. elegans transcripts with lethal RNAi phenotype are expected to evolve under stronger selective constraints and the values of dn/ds showed a non-significant trend towards lower values in TUGs with orthologs with a lethal phenotype compared to a non-lethal phenotypes (p=0.815, two-sided U-test).

The genome-wide representativeness of these 199 SNP markers for the whole genome in population genetic studies was confirmed using heterozygosity-heterozygosity correlation [48]: mean internal relatedness = 0.78, lower bound of 95% confidence intervals from 1000 bootstrap replicates (cil) = 0.444; mean homozygosity by loci = 0.86, cil = 0.596; standardised heterozygosity = 0.87, cil= 0.632.

Differential gene expression

Gene expression was inferred by the unique mapping of 252,388 (71.49%) of the raw reads to the fullest assembly (including all assembled contigs as a “filter”; total contigs/all TUGs in Table 2). Non-A. crassus contigs, and all contigs with fewer than 32 reads overall were excluded. Thus 658 TUGs were analysed for differential expression using ideg6 for normalisation and the statistic of Audic and Claverie [53] for detection of differences. Of these TUGs, 54 showed expression predominantly in the male library, 56 TUGs were more highly represented in the female library (Additional file 5), 56 TUGs were primarily expressed in the libraries from Taiwan, and 22 TUGs were overrepresented in the European library (Additional file 6).

TUGs annotated as acyltransferase were upregulated in the European libraries. However, the expression patterns for other TUGs with overrepresented terms connected to metabolism did not show concerted up or down-regulation. Thus for the term “steroid biosynthetic process”, 2 TUGs were downregulated and 3 contigs upregulated in European nematodes. No enrichment of signal peptide positive TUGs, of TUG conservation categories, or TUGs with C. elegans orthologs with lethal or non-lethal RNAi-phenotypes was identified. Significantly elevated dn/ds was found for TUGs differentially expressed in European versus Asian nematodes (Fisher’s exact test p=0.007; also both up- or down-regulated were significant). TUGs overexpressed in the female libraries showed elevated levels of dn/ds (Fisher’s exact test p=0.041), but contrast male overexpressed genes showed decreased levels of dn/ds (Fisher’s exact test p=0.014).