Functional insights into the infective larval stage of Anisakis simplex s.s., Anisakis pegreffii and their hybrids based on gene expression patterns

Taxonomic identification and selection of specimens for RNA sequencing

Larvae specimens in the same development stage (third stage larvae L3) of Anisakis were obtained from the FAO 27 area of fishing distribution into a general survey carried out by our laboratory. L3 were extracted from hosts following published procedures [34, 37, 38]. L3 were rinsed in 0.9% saline solution, placed in an antibiotic solution for 30 min and then rinsed in bi-distilled water for 1 h before molecular identification and RNA extraction. The larvae were individualized (from each host) and for each specimen the caudal part was used for DNA extraction and PCR amplification for species identification; this caudal part and the rest of the body were separately stored at -80 °C until required. Species identification in that survey was performed following the taxonomic criteria of [23, 39] by using the ITS1 region of the nuclear ribosomal DNA (rDNA). Breifly: individual Anisakis spp. L3 stage juveniles was placed in an Eppendorf tube after previously having a small part of the caudal region removed to allow molecular identification of each individual. DNA was extracted and purified using the Speedtools Tissue DNA Extraction Kit (Biotools) following manufactures instructions. Molecular identification was carried out for each individual using PCR-RFLP. The forward primer A 5’-GTCGAATTCGTAGGTGAACCTGCGGAAGGATCA-3′ and reverse primer B 5’-GCCGGATCCGAATCCTGGTTAGTTTCTTTTCCT -3′ [39] were used in reactions containing, 10 mm Tris-HCl (pH 8.3), 1.5 mm MgCl2, 50 mm KCl, 200 mm each of dATP, dCTP, dGTP and dTTP and 1 unit of DNA polymerase (Biotools B and M labs, S.A. Madrid, Spain). Initial denaturalization was carried out for a period of 2 min at 94 °C followed by 35 cycles of 94 °C for 1 min, 58 °C for 1 min, 72 °C for 1 min followed by a final 7 min extension at 72 °C. Amplified DNA fragments were digested with the restriction enzymes HhaI and HinfI (New England Biolabs, Massachusetts, MA, USA) following manufacturer’s instructions. Restriction fragments were separated by electrophoresis in Tris- Borate-EDTA (TBE) buffered 2.5% agarose gel, stained with SYBERsafe and visualized with UV illumination.

Sequence annotations

The transcriptomes of A. simplex s.s., A. pegreffii and the hybrids were annotated using the NCBI tool blastx [55] with an evalue of 10e-⁵ as cutoff threshold and the NCBI non-redundant protein (NR) and the Eukaryotic Orthologous Groups (KOG) databases [56] as reference subjects. Gene Ontology (GO) terms as well as their evidence code [57] and Enzyme Commission (EC) numbers [58] were assigned to annotations using correlation tables among GenBank accessions and GO and EC information. InterPro domain [59] annotations were also annotated using identical methodology. Transcripts that had no informative blastx matches were further used to search the NCBI nucleotide database (NT) using NCBI BLAST tool blastn and identical evalue cut-off. The annotation protocol is automatized as a pipeline cited in the section below “pipelines”. Proteins (open reading frames) longer than 75 residues were predicted from transcripts using Transeq from the EMBOSS package [60] and OrfPredictor [61]. Signal 4.1 [62] and TMHMM v.2.0 [63] were used to predict secretory signal peptides and transmembrane regions. Metabolic pathways were retrieved from the KEGG web site [64] using EC annotations as queries.

Reference transcriptome

Minimus2 from the Amos package [65] was used to merge the 3 assemblies into a merged transcriptome. Remaining redundancy was eliminated using CD-HIT [66] at a similarity threshold of 0.95. The merged transcriptome after redundancy filter resulted in 75,380 consensus transcripts that were annotated via BLAST as described above for transcriptomes using the NR database.

Differential expression and enrichment analyses

Bowtie2 [67] was used to map the two RNAseq pair-end libraries sequenced per each Anisakis sample on the merged transcriptome. In all cases over 97% of the reads were successfully mapped on the reference (data not shown). Next, Corset [68] was used to parse the bam files resulting from the mapping step and extract a cluster file grouping the 75,380 consensus transcripts into 74,751 transcript clusters and a count file summarizing the read counts obtained per each cluster from all pair-end libraries mapped to the reference merged transcriptome. The count file was used as input to EdgeR [69] to perform 3 differential expression tests at the whole-transcriptome level: “hybrids vs. A. pegreffii”, “hybrids vs. A. simplex s.s.” and “A. pegreffii vs. A. simplex s.s.”. The SI and LI pair-end libraries prepared from each sample were considered as technical replicates. Transcripts were considered as differentially expressed when they had a p-value < 0.05 after False Discovery Rate (FDR) correction applied by EdgeR using the Benjamini-Hochberg method [70] (referred as just FDR through the rest of the manuscript). Transcripts clusters were assigned the annotations of their consensus transcripts using the GPRO worksheet [71]. Clusters from transcripts with no NR/NT annotations were considered in differential expression analyses only if they had non-zero read counts in at least 2 of the 3 samples.

GO and metabolic pathway enrichment analyses were performed using GOseq [72] and the results from differential expression as input, following the GOseq manual indications for non-native transcriptomes. Concretely, we used 4 input files: 1) a file summarizing the consensus transcripts found as differentially expressed as a FDR < 0.05 as a cutoff in differential expression analyses; 2) a file with all assayed consensus transcripts; 3) a file with the GOs or with the metabolic maps associated to the consensus transcripts (depending on the analysis); and 4) a file with the sequence size of each assayed cosensus transcript. A priori, enriched GO terms and metabolic pathways supported by P-values < 0.05 in the resulting Wallenius distribution, were considered as significant. A 5% FDR correction was performed on the P-values obtained from the enrichment analyses using the Benjamini-Hochberg method to outline the most enriched GOs and pathways.

With the aim to also asses differential expression of genes encoding for excreted or secreted (ES) compounds read counts of transcript clusters with annotation of secretory leader peptide from Signal-P domain or with annotation of GO terms for extracellular compartments were exported to a new count file. Counts in this file from transcripts sharing the same protein descriptions (from NR annotations) were combined using a custom php script called combine_counts.php that is available in AnisakisDB (see the section below AnisakisDB). We combined information from the Uniprot browser [73] and the literature to filter the transcript set to include only those verified as coding for ES products or with confirmed secretory signal peptides but unknown sub-localization. The final count file had 356 non-redundant protein descriptions associated with genes encoding putative ES products. Three differential expression analyses (“hybrids vs. A. pegreffii”, “hybrids vs. A. simplex s.s.” and “A. pegreffii vs. A. simplex s.s.”) were performed using EdgeR, with the small and large insert size libraries of each transcriptome considered as technical replicates (as above).

A third differential expression study focusing on the distinct families of allergens detected in the transcriptomes of A. simplex s.s., A. pegreffii and the hybrids, was also performed. We retrived 509 allergen sequences belonging to 150 allergen families (AllFams) in fungi, plants and animals from the Uniprot protein knowledgebase [73]. We have also used additional information provided by the AllFam database of allergen families [14], and the WHO/IUIS Allergen Nomenclature Database [13, 14]. These 509 allergen sequences were used to identify (via blastx search) homologs in the merged transcriptome. Nine hundred and thirty seven consensus transcript were identified as potential allergen candidates based on a blast e-value threshold of 2e^− 06. Using identical methods to those described above for ES protein genes, read counts for these 937 candidates were collated for each AllFam family of allergens in a count file that was then used as input to perform 3 differential expression analyses using EdgeR as above.

Pipelines

The above described protocols for pre-processing, de novo assembly, annotation and differential expression were executed using the DeNovoSeq and the RNAseq pipelines provided by the GPRO suite [70].

AnisakisDB

We have reconstructed the transcriptomes of the L3 species A. simplex s.s. and A. pegreffii and their hybrids producing a collection of 290,296 transcripts spanning more than 200 Mb. We have also reconstructed a consensus transcriptome of 75,380 consensus sequences. Sequences have been annotated using different databases and classificatory systems (NR/NT, KOG, InterPro, GO, EC) thus creating abundant material (annotations and other material such as scripts, statistics graphical representations and the differential expression and enrichment files). To house the aforesaid material and make it freely available for readers, we constructed a web-available database, AnisakisDB. The database was programmed using Laravel 5 framework for public client platform and a CodeIgniter 2 framework for server platform, both based on the PHP programming language [74], a MySQL database [75] and an Apache server [76] hosted in a Linux environment. The database portal includes a blast search powered by the NCBI BLAST package [77]. Within AnisakisDB, sequences, annotations and backups of all performed analyses of differential expression and enrichment are distributed in seven sections (“Transcriptomes”, “Blast”, “Download”, “COII markers”, “Secretome” and Allergome) managed by an intuitive menu and additional utilities to retrieve or compare sequences and annotations. In addition, within the section “Transcriptomes” we have implemented different web site sections providing access to 5 independent Venn Diagrams created based on the GI, KOG, InterPro, GO, EC annotations used to create the Venn diagram of Fig. 1. In each Venn-based web site, clicking on any intersected number the Venn opens a dialog summarizing the annotations that correspond to the intersection. The Section “COII markers” is another web site where the COII-based phylogeny inferred for selection of specimens has been graphically implemented, also as a dynamic tree representation. AnisakisDB is freely and online accessible at www.anisakis.mncn.csic.es.

Transcriptome de novo assembly and annotation

We sampled Anisakis L3 from a single ocean region (FAO 27) and fish species (M. merluccius), to avoid issues that might arise from response to different paratenic hosts and minimize other environmental and phylogeographic confounding factors. L3 were identified according the taxonomic approaches of [23, 39] by using the ITS1 region from rDNA. The phylogenetic tree of 141 COII sequences from 100 populations plus another 41 sequences obtained from GenBank showed all A. simplex s.s and all A. pegreffii samples clustering in separate sister clades (Additional file 1). The tree also showed that all COII sequences from hybrids fell into one or another of these two clades, allowing immediate classification of hybrids as to their maternal species. RNA was extracted from pools of 10 A. simplex s.s. 10 A. pegreffi and 10 hybrids (five mothered by A. simplex s.s. and five mothered by A. pegreffii). Two paired-end RNA-Seq datasets were generated from each pooled sample, forming technical replicates of the transcriptome. Four transcriptome assemblies were generated, one for each sample (A. simplex s.s., A. pegreffii and hybrids) and one that merged the 3 sample-derived ones (Table 1). The 3 sample-derived transcriptomes are similar in span and complexity. The assembly of A. simplex s.s. was longer and had an elevated number of predicted transcripts, apparently because of a greater number of isoforms (30,366 more than A. pegreffii assembly and 45,059 more than the hybrid assembly). We think this might be due to sequencing technical issues, elevated levels of alternative splicing (i.e. isoforms) or elevated heterozygosity in the transcriptome of A. simplex s.s. as the number of unigenes is in contrast quite equilibrated (4657 unigenes more than A. pegreffii assembly and 6989 more than the hybrid assembly).

Expression patterns and enrichment of metabolic pathways and gene ontologies

The merged (consensus) transcriptome reconstructed in this study based on the three Anisakis transcriptome consist of 75,380 consensus sequences, which were used as a mapping reference to compare the expression patterns at the whole-transcriptome level of the 3 taxa. Of the 75,380 consensus transcripts, 32,999 have no annotation and therefore correspond to unknown sequences expressed by the at least two of three taxa (they might be unknown proteins, potential non-coding RNAs and/or even mobile elements). The remaining 42,381 consensus transcripts correspond to coding sequences annotated based on 12,511 different proteins (according to the number of non redundant GI annotations) thus meaning that any biological interpretation of differential expression and enrichment based on the merged transcriptome is based on the expression patterns of 12,511 coding genes (around 68% of the pan-transcriptome here characterized based on the three transcriptomes). After mapping the reference transcriptome, the 75,380 consensus sequences of the merged transcriptome were grouped into 74,751 transcript clusters (see Methods) correcting any potential bias in isoforms as detected in the assembly of A. simplex s.s. The Biological Coefficient of Variation (BCV) and the average dispersion of the assayed samples were inferred and evaluated to asses to suitability of the fastq libraries, resulting in quite appropriate values (BCV = 0.4359 and dispersion = 0.19005) for a study of differential expression. A similar conclusion resulted from plotting BCV against the average log of counts per million of mapped reads. Additionally, a multidimensional scaling (MDS) plot based on logFC of the differences among taxa sources and the sequenced replicates per taxa shows that while the SI and LI libraries used as taxon replicates fall close to one another, the replicates of different taxa fall widely separated (Fig. 2b). This result anticipates important differences among species and hybrids. The first plot dimension correspond to the differences among taxa sources (A. simplex s.s., Hybrids, A. pegreffii) that provides the major source of variability. The second plot dimension corresponds to the sequenced replicates per taxa (small and medium insert size fastqs) that present practically no differences with the exception of the hybrids where the 2 replicates are slightly more heterogeneous because they were created as a pool of 5 hybrids close to A. simplex s.s., and 5 other close to A. pegreffii. Consistent with this perspective, of these 74,751 transcript clusters, 8239 resulted as differentially expressed in the analysis performed between the hybrids and A. pegreffii using FDR < 0.05 as cutoff, 23,549 were differentially expressed in the analysis performed between the hybrids and A. simplex s.s. using same FDR, and 24,813 were differentially expressed between A. pegreffii and A. simplex s.s. with same FDR (Fig. 2a and Additional file 3). In general, these analyses of differential expression based on the merged transcriptome suggest that the A. pegreffii parental expression pattern predominates over A. simplex s.s. in the hybrids, in agreement with population genetic studies of Moroccan Anisakis [36].

The 272 ECs identified as common in the 3 transcriptomes link to 53 metabolic pathways [64] related to nutritional requirements and biosynthesis, energy metabolism, metabolism of xenobiotics, environmental information processing, signal transduction and genetic information processing. Ten of these 53 pathways were shown as differentially enriched in at least one of the 3 enrichment comparisons after correction of the pvalue at FDR < 0.05 (Additional file 4). As shown in Fig. 3, two of these pathways (oxidative phosphorylation and methane metabolism) were differentially enriched in all comparisons thus indicating that the main metabolic differences among A. simplex s.s., A. pegreffii and their hybrids concern to the energy metabolism. Both A. pegreffii and A. simplex s.s. present overexpression with respect to the hybrids as for transcripts associated to methane metabolism while differentially expressed transcripts involved in oxidative phosphorylation are overexpressed in A. pegreffii with respect to the 2 other taxa. Figure 3 also reveals that the enrichment of the glycosaminoglycan degradation pathway was significant between the performed comparisons of the hybrids and A. pegreffii against A. simplex s.s. but not when the hybrids and A. pegreffii are compared. The degradation of glycosaminoglycans pathway mainly concerns to transcripts coding for hyalurononglucosaminidases, which are known to play role in the extracellular matrix as soluble components and polyelectrolytes interacting with growth factors and other transient components of the extracellular matrix [78]. In this case, the enrichment derives from a strong overexpression of the transcripts coding for these enzymes in the hybrids and A. pegreffii with respect to A. simplex s.s. Reversely, the pathway of alanine aspartate and glutamate metabolism is shown as significantly enriched when the hybrids and A. simplex s.s. are compared with A. pegreffii but not when they were compared each to other thus suggesting that the A. pegreffii parental does not predominate in all expression patterns of hybrids. The remaining pathways (pyruvate metabolism, aminoacyl-tRNA biosynthesis, streptomycin biosynthesis, starch and sucrose metabolism, pentose phosphate pathway, and cysteine and methionine metabolism) were found as differentially enriched only in 1 of the 3 performed comparisons, predominantly when A. pegreffii was compared against any of the 2 other taxa (Fig. 3). Results sugest that A. pegreffii and the hybrids exhibit for these pathways more differentially expressed transcripts than those that would be expected by chance. In fact, the ratio between differentially expressed and assayed transcripts associated to the enrichment of each pathway is predominantly lower when the hybrids are compared with A. pegreffii than when they are compared with A. simplex s.s. or when A. simplex s.s. is compared with A. pegreffii. Finally, it is worth to also discuss the detection of streptomycin biosynthesis pathway, which is associated with the over expression of 7 transcripts encoding for 2 type of enzymes, inositol-3-phosphate synthases (lyases) and dTDP-4-dehydrorhamnose 3,5-epimerases (isomerases) in A. pegreffii and A. simplex s.s. when compared to the hybrids. The expression of the genes coding for these enzymes might have been induced by the antibiotic treatment applied on L3 specimens before the extraction (see methods) thus apparently revealing some deficiencies in the hybrids with respect to their parents of for this metabolic pathway. However, we should to also take into account that Inositol-3-phosphate synthase also participates in the metabolism of inositol phosphates, mainly to produce membrane receptor playing role in signalling, cell growth, apoptosis, cell migration, endocytosis, and cell differentiation [79], while dTDP-4-dehydrorhamnose 3,5-epimerases are also associated to the nucleotide sugar metabolism and with the biosynthesis of dTDP-rhamnose to produce rhamnose. The biosynthesis of dTDP-rhamnose is a new pathway involved in hypodermal development or in the production of the cuticle or surface coat in larval stages of C. elegans and other nematodes [80]. We thus think that inositol-3-phosphate synthases and dTDP-4-dehydrorhamnose 3,5-epimerases of anisakids might link not only to the metabolism of secondary metabolites but also to signalling, nucleotide sugar metabolism and larval development.

Analyses of GO enrichment were also performed based on the 1211 GOs shared by the 3 transcriptomes and revealed differential enrichment of 91 GOs at FDR < 0.05 in at least on of the 3 performed comparisons (Additional file 4). Forty of these 91 correspond to biological processes terms, 37 are molecular functions and 14 cellular components (Fig. 4). Thirteen of these 91 were significant in all comparisons and correspond with functions, processes and sub-localizations associated with the energy metabolism and this mainly includes oxidases, oxidoreductases, dehydrogenases and other enzymes and binding proteins with assigned roles in redox processes, and electron carriers probably related to the transmembrane electrochemical gradient as well as roles for heme-biding proteins (Cytochrome P450s). We think that the latter might also be related with sterol and heme recruitment. As other nematodes, anisakids probably recruit exogenous sterol and heme recruited from host sources to implement sterol and heme dietary needs. Consistent with this argument, the 3 transcriptomes addressed in this study present annotations of prenyltransferases, squalene oxidases and sterol regulatory element-binding proteins (SREBP); prenyltransferases play role in the synthesis of pyrophosphate and geranylpyrophosphate, a major metabolic chokepoint in de novo sterol biosynthesis pathway [81] while squalene and derivatives are the main source of host sterol and SREBP is known to stimulate sterol biosynthesis. The enrichment of the remaining GOs shown in Fig. 4 was significant in 1 or 2 of the 3 performed comparisons and refers to molecular functions biological processes and cellular components related to the enriched pathways shown and discussed above but also with transposition (mainly integrases and proteases encoded by Bel/Pao LTR retrotransposons), transport, growth and development, locomotion, and with the extracellular region (ES proteins, allergens and cuticle collagens). In broad terms and for almost all evaluated GOs (with several exceptions) the ratio between differentially expressed and assayed transcripts is again lower if the hybrids are compared with A. pegreffii than when the hybrids are compared with A. simplex s.s. or when the latter is compared with A. pegreffii. In turn, the comparison between the hybrids and A. pegreffii presents largely more enrichment than the others; of the 91 GOs shown as differentially enriched, 67 terms are significant in the aforesaid comparison between the hybrids and A. pegreffii while 27 GOs are significant when the hybrids and A. simplex s.s compared and 34 when the latter is compared with A. pegreffii. Interestingly, the most prominent differences in GO enrichment observed between the hybrids and the 2 parents derive from the differential expression of transcripts related with the inositol polyphosphate 5-phosphatase activity involved in inositol phosphate-mediated signaling. While these 2 GO terms are significantly enriched when they are compared with their parental species but no when the parents are compared each to other. Indeed, we think that the deregulation of the inositol polyphosphate 5-phosphatase enzyme and how it affects the whole inositol phosphate-mediated signalling pathway it is worth to be investigated, moreover if we take into account that hybrids are not able to reach adult stages and that inositol phosphate-mediated signalling has been associated to longevity, growth and reproduction in nematodes and other invertebrates [82]. The enrichment of GOs is also consistent with the perspectives provided by the differential expression analyses and the enrichment of metabolic pathways that together supports the idea of that although the alleles of the A. pegreffii parental apparently predominate over those of A. simplex s.s. in the expression patterns of L3 larval stage hybrids the latter are intermediate biological entities with metabolic and functional requirements, quantitatively different to those of the two parental species.

The set of metabolic pathways and GOs (i.e. the functional profiles) of A. simplex s.s., A. pegreffii and their hybrids are consistent with the aerobic intermediary metabolism of nematodes thus supporting the idea that almost all pathways typically observed in C. elegans and other nematodes are present also in these anisakids. A graphical idealisation of the intermediary metabolism of nematodes is represented in Fig. 5 according to [81–85] and the most representative pathways and GOs detected in the 3 Anisakis transcriptomes have been highlighted. Our conclusion is that the A. pegreffii parental predominates over A. simplex s.s. in the expression of hybrids albeit the latter, and A. pegreffii also present significant differences between them indicating that hybrids are intermediate biological entities with metabolic and functional requirements quantitatively different to those of their parental. The hybrids thus offer an excellent opportunity to investigate the mechanisms of speciation in Anisakis including the novel biological and ecological processes where hybrid genotypes might have a relevant role or significant lacks. Hybridization among A. simplex s.s. and A. pegreffii has indeed an important role in the evolutionary biology of these 2 species through gene introgression [86]. While hybrids do not appear to reach the definitive host, allowing development to reproductive adults, they reflect the dynamics of moving areas of sympatry. Hybrids may thus be crucial in understanding the microevolutionary processes active in A. simplex s.l., in host-parasite ecology (including the range of paratenic hosts exploited), and in whale migrations (as interbreeding must take place in whales as final host).

Genes encoding for excretory/secretory proteins

To specifically focus on the potential sub-localizations of the products encoded by differentially expressed transcripts associated to enrichment we created a box plot based on the logFCs of transcripts associated to enriched cellular component GOs (Fig. 6). The analysis shows the relationships of over/under expression observed among the 3 taxa with respect to products sub-localized in mitochondrion, respiratory chain, ATP transporting factors and myosin complexes with ATP-dependent roles in muscle contraction, motility and the extracellular region. A significant number of outliers (due to their extremely high logFCs) correspond to transcripts whose product is sub-localized in the extracellular region thus indicating relevant differences with respect to genes encoding likely ES proteins. All ascaridids, Anisakis included, release ES products to the host environment. Infective L3 penetrate the intestinal wall of the fish host, and may encyst on the surface of the internal organs or migrate to the muscular tissues. The migration may occur not only when the host is alive [87–91] but also after death [92–94]. These processes are likely to be mediated by secreted enzymes. ES proteins have additional functional roles connected to nutrition, infectivity, allergy, immune evasion or pathogenicity [34, 95, 96] playing important roles in the host-parasite interaction. In the aforeshown analyses performed at the whole-transcriptome level we showed that part of the differential enrichment between taxa corresponds to transcripts encoding likely ES proteins (above). According to our annotations these transcripts correspond to the expression of at least 356 non-redundant gene descriptions of which we create a count file of mapped reads to explore how different are the expression patterns of the investigated taxa with respect to ES genes (see Methods). Nearly three-quarters (266/356) were found to be significantly differentially expressed in one or more comparisons under a FDR < 0.05 as threshold cutoff (Additional file 5). More specifically, 151 genes were found as differentially expressed when comparing hybrids and A. pegreffii. Of these, 79 genes were overexpressed in hybrids and 72 overexpressed in A. pegreffii. Comparing hybrids with A. simplex s.s., identified 191 differentially expressed genes, of which 107 were overexpressed in hybrids and 84 in A. simplex s.s.. Comparing A. simplex s.s. and A. pegreffii, 100 transcripts were overexpressed in A. pegreffii and 82 in A. simplex s.s. A heatmap with clustering of the logFC per comparison of the 70 most differentially expressed ES genes (Fig. 7) permits to visualize four major patterns of expression. Rows (i.e. the gene expression patterns) split in 4 Clusters of expression designated as clusters 1,2,3 and 4. Cluster 1 groups genes over-expressed in A. pegreffii compared to both A. simplex s.s. and hybrids, cluster 2 is a single gene, overexpressed in hybrids. Cluster 3 is defined by genes over-expressed in A. simplex s.s. compared to both A. pegreffii and hybrids, and genes belonging to Cluster 4 were predominantly over-expressed in both A. pegreffii and hybrids compared to A. simplex s.s. ES genes with allergenic properties are highlighted with red asterisks. In particular, the analysis shows that in mixed infections different species and hybrids of Anisakis may produce different kinds and levels of ES proteins, or in other words that the modulation of the host environment may depend on the mix of larvae from species and hybrids. Considering that hybrids are not able to complete their life cycle and that no colony of L3 hybrids has been yet described in nature parasiting alone their hosts (as far as we know hybrids have been only reported accompanied by their parental species) it is tempting to speculate based on our results that hybrids might be able to survive in host environments only when they are accompanied by the parental species. Facilitation may also derive from genes overexpressed or downexpressed in the hybrids.

Detection and classification of allergens

Anisakiasis is a consequence of allergic reaction to several allergenic molecules expressed and/or secreted by Anisakis in their intermediate hosts. Briefly, 15 allergens (Ani s 1 to Ani s 14 and Ani s troponin C) are formally considered from A. simplex s.s. according to the AllFam database of allergen families and the WHO/IUIS Allergen Nomenclature Database [13, 14]. Five of the 15 Anisakis formal allergens - Ani s1, and Ani s 2 and Ani s 3, Ani s 7, and Ani s 12 - are considered to be major allergens (one to which > 50% of allergic patients immunologically react and present unusual levels of IgE and IgG in serum [97, 98]); Ani s 1 belongs to the animal Kunitz serine protease inhibitor family of protease inhibitors (AF003); Ani s 2 is a myosin heavy chain family (AF100) [99] that includes the highly cross-reactive mite paramyosin major allergens [100]; Ani s 3 is part of the tropomyosin family (AF054), which have been identified as minor inhalative allergens in arthropods and as major food allergens in crustaceans, mollusks and Anisakis [101]; and Ani s 7 and Ani s 12 are 2 proteins with unknown biological functions. The remaining Anisakis allergens are minor allergens; Ani s 4 belongs to the cystatjn family (AF005) [102] of reversibly binding cysteine protease inhibitors; Ani s 6 belongs to the cysteine-rich trypsin inhibitor-like domain family [103] (AF027); Ani s 5, Ani s 8 and Ani s 9 are members of the nematode-specific SXP/RAL-2 family (AF137), which have unknown function [104, 105]; Anis s 10, Ani s 11 and Ani s 14 similarly remain functionally unclassified [106–108]; Anis s troponin C is a member of the EF hand family (AF007), which has roles in calcium signaling and binding [109]; and Ani s 13 belongs to the globin family (AF009) [110]. Yet more interesting, recent proteomic analyses have identified up to 28 potential new allergens in Anisakis that can be classified into distinct allergen families [15] and a more recent transcriptomics study [16] has identified several new allergens from A. simplex s.s. and A. pegreffii. Additional potential allergens (Ani s 24 kDa, Ani s CCOS3, Ani s cytochrome B, Ani s FBPP, Ani s NADHDS4L, Ani s NARaS and Ani s PEPB) have been also considered according to previous studies [111–113]. All these molecules are classified or can be classified as food allergens because they are present in marine fish. The characterization of allergens in anisakids is still open and that will aid diagnosis and possible alleviation of the serious symptoms of anisakiasis. A complete allergome (all allergen proteins) catalogue for anisakids is required.

Special importance should be paid to proteins related to mechanisms of parasitism and pathogeny of these nematodes. According to a recent transcriptomic study [114] in the pharingeal region of A. simplex s.s. there are at least 226 transcripts which are potentially involved in invasion and host pathogen interplay. The same study found at least 339 transcripts with a similar role in A. pegreffii. Regarding this, 13 immunoreactive proteins from the esophageal glands of A. simplex s.s. have been described as potential allergens (data not published; PARASITE project [115].