De novo transcriptome sequencing and analysis of Anisakis pegreffii (Nematoda: Anisakidae) third-stage and fourth stage larvae

Abstract Anisakis pegreffii is known as one of the causes of a fish-borne zoonosis, anisakidosis. Despite its significant public health and food hygiene impacts, little is known of the pathogenesis, genetic background of this parasite, at least partly due to the lack of genome and transcriptome information. In this study, RNA-seq and de novo assembly were conducted to obtain transcriptome profiles of the A. pegreffii third and fourth larvae. The third stage larvae (APL3) were collected from chub mackerel and the fourth stage larvae (APL4) were obtained by in vitro culture. In total, 47,243 and 43,660 unigenes were expressed in APL3 and APL4 transcriptomes. Of them, 18,753 were known and 28,490 were novel for APL3, while 18,996 were known and 24,664 were novel for APL4. The most abundantly expressed genes in APL3 were mitochondrial enzymes (COI, COII, COIII) and polyubiquitins (UBB, UBIQP_XENLA). Collagen-related genes (col-145, col-34, col-138, Bm1_54705, col-40) were the most abundantly expressed in APL4. Mitochondrial enzyme genes (COIII, COI) were also highly expressed in APL4. Among the transcripts, 614 were up-regulated in APL3, while 1,309 were up-regulated in APL4. Several protease and protein biosynthesis-related genes were highly expressed in APL3, all of which are thought to be crucial for invading host tissues. Collagen synthesis-related genes were highly expressed in APL4, reflecting active biosynthesis of collagens occurs during moulting process of APL4. Of these differentially expressed genes, several genes (SI, nas-13, EF-TSMT, SFXN2, dhs-27) were validated to highly transcribed in APL3, while other genes (col-40, F09E10.7, pept-1, col-34, VIT) in APL4. The biological roles of these genes in vivo will be deciphered when the reference genome sequences are available, together with in vitro experiments.

The family Anisakidae comprises the nematode species whose adult stages can be found in aquatic animals, while the third stage larvae (L3) generally exist in the body cavity, visceral organs and muscles of various fish and squid species. Although humans are not the final hosts of these nematodes, the larval forms of anisakid nematodes, particularly those of the genera Anisakis and Pseudoterranova are known to be associated with human infection by the ingestion of raw or undercooked fish or cephalopods harboring these larvae (Mattiucci and D'Amelio, 2014). Infection with these nematodes is therefore considered to be a threat to public health due to their zoonotic potential and the presence of larvae in fish products cause aesthetic problems, reducing commercial value. Many studies have also revealed that alive or dead anisakid larvae can cause allergic reactions to humans, which are often associated with elevated level of immunoglobulins E (Audicana and Kennedy, 2008;Nieuwenhuizen and Lopata, 2013).
Of genus Anisakis, Anisakis simplex sensu lato (s.l.) has been traditionally considered as the main causative agent of anisakidosis (Audicana and Kennedy, 2008). But since molecular approaches have been introduced for identification of anisakid nematodes, it was proved that there are three species in A. simplex complex and of them, the two sibling species A. simplex sensu stricto (s.s.) and A. pegreffii were shown to be the causative agent of human infection (Mattiucci and D'Amelio, 2014 and the references therein). Several in vitro and in vivo studies have also demonstrated that A. pegreffii has pathogenic potentials to humans, as well as A. simplex (s.s.) (Suzuki et al., 2010;Jeon and Kim, 2015). Therefore, more emphasis should be placed on studying epidemiology and pathogenicity of A. pegreffii as well.
Anisakis pegreffii is known to be widely distributed in the Austral Region between 30°N and 55°S, particularly in the Mediterranean Sea and the East Asia (Umehara et al., 2006;Mattiucci and Nascettii, 2008;Suzuki et al., 2010;Bak et al., 2014;Jeon et al., 2016). Many interemediate/paratenic hosts and definitive hosts of A. pegreffii are known to be shared by A. simplex (s.s.), and recombinant genotypes of these two sibling Anisakis species are reported to sympatrically exist in the waters mentioned above (Abollo et al., 2003;Umehara et al., 2006;Farjallah et al., 2008;Cavallero et al., 2014). In several human clinical cases, Anisakis larvae removed from patients were identified as A. pegreffii by molecular methods, reflecting their pathogenic potentials to humans (Mattiucci et al., 2013).
Recent impressive progress in genome-wide analyses of socio-economically important nematode parasites helped us to better understand the genetic information of them covering the general biology, hostparasite interaction, pathogenicity and development of drug or vaccine candidates. Moreover, the increase of sequencing data have opened a new era in comparative studies in parasitic nematodes, thereby offering many useful genetic information from different species or closely related species with which to compare pathogenesis, conduct differential diagnosis and a large scale epidemiology (Cantacessi et al., 2015, Lv et al., 2015. Currently, 134 parasitic or free-living nematodes genomes or transcriptomes are available in databases (http://www.wormbase.com, accessed on November 25, 2019) and the numbers will be increased. For anisakid nematodes, several transcriptomic analyses have been focused on A. simplex, a well-known cause of human anisakidosis or on comparative studies with a sibling species, A. pegreffii (Baird et al., 2016;Cavallero et al., 2018;Llorens et al., 2018), but not on A. pegreffii itself. In this study, we obtained transcriptomes of A. pegreffii L3 (APL3) and in vitro-induced L4 by RNA-seq and de novo assembly. Then, we compared transcriptomes of these two different stages of A. pegreffii to provide detailed information on differences regarding their biology and pathogenicity against their hosts.

Sample preparation
Alive APL3 were obtained from chub mackerel (Scomber japonicus), which is known to be the one of the main intermediate/paratenic hosts for A. pegreffii in Korea (Bak et al., 2014). Freshly caught chub mackerel were purchased from a local fish market located in the east coast of Korea and immediately transported to the laboratory. The body cavity was longitudinally cut and the viscera were carefully examined for collecting the nematodes.
The collected nematodes were placed in a Petri dish filled with sterile PBS and washed several times to remove any tissue debris. Then, they were observed under a stereomicroscope and actively moving nematodes without any injury were selected for in vitro culture to obtain A. pegreffii fourth stage larvae (APL4). Each larva was treated with antibioticantifungal solution as described elsewhere before starting in vitro culture (Iglesias et al., 2001). Then, individual larvae were placed in a sterile 24 well tissue culture plate (SPL, Korea) with 1 ml of culture medium (pH 4.0) in each well, at 37˚C with 5% CO 2 atmosphere. The composition of culture medium and culture conditions were followed as described by Iglesias et al. (2001). The culture medium was changed twice a week and the larvae were checked every day. Each individual L3 were considered to be molted to L4 when the sheath were found by stereoscope, and APL4 is known to be induced within approximately five days in these conditions (Jeon and Kim, 2015).
Each individual APL3 and APL4 were washed with sterile PBS, then individually placed in a 1.5 ml Eppendorf tube and stored at −80˚C until use. Each larva was identified by PCR-RFLP and subsequent sequencing for the mitochondrial cox2 gene according to the method described elsewhere (D'Amelio et al., 2000;Nadler and Hudspeth, 2000). The caudal end of each larva was cut and used for molecular identification, and the rest of them were used for next generation sequencing. The larvae identified as A. pegreffii were used for further analysis.

Library construction and next generation sequencing
Total RNA (1 μ g per each sample set) was prepared by homogenizing approximately 50 larvae per each developmental stage with Trizol, following the manufacturer's protocol (Invitrogen, USA). Prior to mRNA isolation, total RNA samples were treated with DNAse. And the RNA yield and integrity was measured by NanoDrop 1000 (Thermo Scientific, Wilmington, USA) and BioAnalyzer 2100 (Agilent Technology, Santa Clara, USA), respectively.
The TruSeq stranded mRNA sample preparation kit (RS-122-2101, Illumina, San Diego, USA) was used for preparing mRNA sequencing libraries. Poly-Acontaining mRNA was purified from 1 μ g of total RNA, by using Oligo dT attached magnetic beads. Then the purified mRNA was disrupted into short fragments, and first-stranded cDNAs were synthesized using SuperScript∏ -reverse transcriptase (Invitrogen, USA) and random hexamers. cDNA with adaptors ligated to both ends were enriched by PCR. The cDNA library size and quality were electrophoretically evaluated by using the Agilent DNA 1000 kit on a BioAnalyzer 2100. The libraries were subsequently sequenced on an Illumina HiSeq 2500. Image analysis was performed using the HiSeq control software version 2.2.58. Raw data were processed and base calling was performed using the standard Illumina pipeline (CASAVA version 1.8.2 and RTA version 1.18.64).

De novo assembly and functional annotation
The complete sequences of two sample sets were subjected to the adaptor and quality trimming by FastQC method. The contigs were obtained by the assembly of preprocessed sequences assembly using Trinity (Haas et al., 2013). Then, they were clustered using the CD-HIT and TGIGL, to obtain the unique genes as a reference transcriptome (Pertea et al., 2003;Li and Godzik, 2006). The preprocessed reads were further mapped to the reference transcriptome using BWA and the expression patterns (i.e. FPKM) of each genes/transcripts were obtained using RSEM method (Li and Durbin, 2010;Li and Dewey, 2011). Normalization of the FPKM values was conducted by edgeR using TCC R package (Sun et al., 2013). The reference transcriptomes were annotated by mapping against to NR database using BLASTx and SwissProt/UniProt database using InterProScan. The protein coding regions were also annotated with TransDecoder (Haas et al., 2013). The complete annotations were merged together to improve the annotations of each genes/transcripts using in-house scripting.

Orthologous cluster analysis
In total, six genomes of nematodes were selected for whole-genome orthologous cluster analysis, i.e. Ancylostoma ceylanicum, Brugia malayi, Trichinella spiralis, Caenorhabditis elegans, A. simplex along with A. pegreffii in our study. Complete proteins sequences of these five nematodes were downloaded from WormBase (http://www.wormbase.org) for those selected nematode species and the orthologous clusters were analyzed using OrthoVenn2 with default options (Xu et al., 2019).

Differentially expressed genes (DEGs) analyses and GO enrichment analysis of DEGs
For DEG analyses, complete FPKM and edgeR scores were taken, and the statistical significance (P value and Q value) and fold changes (FC) were calculated, based on the normalized values. The filter cut-offs (i.e., P ≤ 0.05, Q ≤ 0.05 and FC ≥ 2) were used to select the DEGs/transcripts from the given pairs. Gene Ontology (GO) terms for the identified genes from DEG analysis were determined using the GOA database (Barrell et al., 2009).

Real-time PCR validation
DEGs in either APL3 or APL4 were validated by realtime PCR analysis. cDNA targeting highly ranked DEGs in each ASL3 and ASL4 transcriptomes were synthesized by using PrimeScript™ RT reagent Kit with gDNA Eraser(Perfect Real Time; Takara, Japan), according to the manufacturer's instructions. As a reference gene, GAPDH (Glyceraldehyde-3phosphate dehydrogenase) gene was selected. Total RNA was extracted from each APL3 and APL4 samples as mentioned above; 1 μ l of total RNA extracted using Trizol (Sigma. USA) was mixed with 2 μ l of 5 X gDNA Eraser Buffer, 1 μ l of gDNA Eraser, 6 μ l of RNase-free dH 2 O, and incubated at 42˚C for 2 min. Then, cDNA was synthesized by adding 4 μ l of  X PrimeScript™ Buffer 2, 1 μ l of PrimeScript™ RT Enzyme Mix I, 1 μ l of RT Primer Mix, 4 μ l of RNasefree dH 2 O, and incubated at 37˚C for 15 min.
The selected A. pegreffii gene sequences and the internal control primers were retrieved from Genbank and gene-specific primers were designed using Primer Expression 3.0 software (Applied Biosystems, USA). The sequences of the primers used in this study are listed in Table 1. The real-time PCRs were set up using SYBR Premix Ex TaqII (Tli RNase H Plus, Takara, Japan) on 96 well plates in 20 μ l reaction volumes consisting 2 ul of 10-fold diluted cDNA as templates, 1.6 μ l of primers (10 μ M), 0.4μ of 50 X ROX reference dye II and 6.0 μ l of dH 2 O. The reaction was conducted with an Applied Biosystems 7500 Real Time PCR System (Applied Biosystems, USA) and cycling conditions involved 15 min of initial pre-denaturation step at 95˚C, then 45 cycles as follows: 15 sec at 95˚C, 15 sec 60˚C, and 30 sec at 72˚C. GAPDH gene was used as internal control and each gene was run in triplicates. Relative expression was calculated using the Δ Δ CT method (Livak and Schmittgen, 2001) and significant differences were identified with the Student's t-test at P < 0.01.

Identification of putative allergens
To make a list of putative allergens from A. pegreffii, the contigs of APL3 and APL4 transcriptomes were compared with the sequence data retrieved from the AllergenOnline Database (Goodman et al., 2016; http:// www.allergenonline.com.; version 19), using BLASTp (e-value cut off: <1e-05, identity matching cut off: >70%). In addition, the gene expression level of putative allergens from APL3 and APL4 were compared by DEG analysis with the same conditions as above.

Transcriptome profiles of A. pegreffii L3 and L4
In total, 57,831,158 (4,395,168,008 bases) and 61,963,202 (4,709,203,352 bases) raw reads were obtained for APL3 and APL4 by Illumina sequencing. The generated raw reads were deposited in the Sequence Read Archive database of NCBI (accession number: PRJNA602791, PRJNA602795).
After cleansing and removing low quality reads (Q < 30), 53,790,942 (4,079,160,122 bases, APL3) and 56,726,910 clean reads (4,301,151,997 bases, APL4) were obtained (Table 2). In total, 53,383 unigenes were obtained by de novo assembly, with 22,387 known genes and 30,996 novel genes. The average length of unigenes was 1,091 base pairs (data not shown).
In total, 47, 243 and 43,660 genes were expressed in APL3 and APL4, respectively (> fpkm 1.0). Of these genes, 28,490 and 24,664 genes were novel for APL3 and APL4 (Table 3). The predicted gene number of A. pegreffii and several other representative nematodes were compared.

Orthologous cluster analysis
The number of orthologous transcripts of A. pegreffii and those of other representative nematodes were compared. In total, 8,344 A. pegreffii transcripts were orthologous to any of those of the reference nematode groups and 2,549 were common to all the five references (data not shown). Most of these 2,549 transcripts are thought to encode proteins involved   in developmental and common biological process of nematodes. In total, 1,022 transcripts did not match with any sequence of the five reference nematode groups, suggesting they may include those putatively encoding novel gene sequences of A. pegreffii.
Highly expressed genes and DEG analyses of A. pegreffii L3 and L4 transcriptomes The abundance of each gene in APL3 and APL4 transcriptome was calculated; ubiquitin genes (UBB, UBIQP_XENLA) were highly expressed in APL3, while collagen genes (col-34, col-138, col-40) were highly expressed in APL4. Several mitochondrial enzyme genes (COI, COII, COIII) were highly expressed both in APL3 and APL4 (Tables 4 and 5). Different gene expression profiles between APL3 and APL4 transcriptomes were investigated by DEG analyses. In total, 1,923 genes were up-regulated either in APL3 or APL4. Among them, 614 were upregulated in APL3 while 1,309 were up-regulated in APL4 (data not shown). The top 20 up-regulated genes either in APL3 or APL4 were listed in Tables 6 and 7.

GO enrichment analyses of DEGs in A. pegreffii L3 and L4 transcriptomes
DEGs were subjected to GO analyses. In total 1,184 DEGs were mapped to either of three main categories (cellular component, biological process and molecular functions). GO analysis of the DEG revealed that 14 cellular component, 41 biological process, and 10 molecular function categories (P < 0.001) (data not shown). Cell, cellular process, and binding were the subcategories to which the largest number of genes belongs to each three categories (Fig. 1).

Discussion
A. simplex complex consist of three morphospecies: A. simplex (s.s.), A. pegreffii and A. berlandi (formerly named as A. simplex C) (Mattiucci et al., 2018), with different genetic structure, biology and geographical distribution. Of these three closely related species, A. simplex (s.s.) and A. pegreffii are known to sympatrically occur in several areas such as the Iberian Atlantic coast and Japanese Sea waters, where the putative hybrids have been reported as both of larvae and adults (Abollo et al., 2003;Umehara et al., 2006;Farjallah et al., 2008;Cavallero et al., 2014). Given that these two species frequently occur in commercially important marine fish species, and humans can become infected by eating raw or undercooked fish and cephalopod harboring these larvae, it is necessary to investigate their pathological potentials of both species. Several in vitro and in vivo studies demonstrated that both of these two species LGALS8  are invasive to cause harmful effects to humans (Suzuki et al., 2010;Jeon and Kim, 2015;Lee et al., 2017). Their pathologic potentials at transcriptomic levels were also recently compared (Cavallero et al., 2018) and partially differentially upregulated repertoires of transcripts were found in each species, but their detailed pathological meanings still await further study.
A. simplex has been known as the major cause of human infection, and therefore many publications exist regarding their pathological process. For example, the existence of several proteases was confirmed by either their biochemical activities or the sequences (Sakanari et al., 1989;Sakanari and McKerrow, 1990;Hotez et al., 1994;Quiazon et al., 2013). Currently several A. simplex transcriptome studies are available (Cavallero et al., 2018;Kim et al., 2018;Llorens et al., 2018), which made the systematic approach possible. However, A. pegreffii, a congeneric species of genus Anisakis is less investigated in terms of their genetic and biological characteristics; there is no reference genome available, thereby we conducted de novo transcriptome assembly with L3 and in vitroinduced of L4 of A. pegreffii, and compared their characteristics in this study.
Mitochondrial enzymes-related (COI, COII, COIII) and polyubiquitin-related genes (UNN, UBIQP_ XENLA) were highly expressed in L3, while collagenrelated genes 138,Bm1_54705) were highly expressed in L4. Mitochondrial cytochrome c oxidase is the key enzyme of aerobic cell respiration in all eukaryotes and many prokaryotes, and different level of this enzyme at different developmental stages of a parasitic nematode Ascaris suum, having the highest activity in L2 (the second stage larvae) and the lowest activity in adult stage, was reported (Takamiya et al., 1996). Adult anisakid nematodes are known to inhabit gastrointestinal tracts of marine vertebrates (Mattiucci et al., 2018), which probably have a low or fluctuating oxygen tension. Therefore, relatively lower activity of cytochrome oxidase would be expected with the development of anisakid nematodes larvae. But, the atmospheric oxygen concentrations of the hosts' microhabitats for anisakid nematodes are currently unknown and should be investigated to make clear discussion.

JOURNAL OF NEMATOLOGY
Collagens are ubiquitous structural proteins constituting the cuticle, a multi-layered flexible exoskeleton for protecting nematodes from adverse environmental conditions (Page, 2013). A new cuticle is synthesized for each developmental stage of nematodes, with many enzymes being involved in this process. Large families of the cuticle collagen genes are known in C. elegans, and also in many parasitic nematodes (Page, 2013 and the references therein). In our study, collagen-related genes were abundantly transcribed in L4, suggesting new cuticles are actively synthesized.
Differential metabolic profiles among different developmental stages of several parasitic nematodes were shown by transcriptome analysis. A filarial nematode Brugia malayi showed different transcriptome repertoires of the overrepresented protease activity in L3 and the overrepresented structural component in L4 (Choi et al., 2011). A strongylid nematode Haemonchus contortus also showed the genes associated with oxygen transport and heme binding were up-regulated in L3, and the genes associated with body morphogenesis and various metabolic process were up-regulated in L4 (Laing et al., 2013). Parasitic nematodes should adapt themselves to different environments with differing food sources  (Cavallero et al., 2018;Kim et al., 2018;Llorens et al., 2018), but information on transcriptomic changes during the life cycle are still scarce. In this study, we selected the highly ranked DEG either in L3 or in L4 and validated them by real-time PCR. Of 13 highly ranked differentially transcribed genes in L3, five genes (SI, nas-13, EF-TSMT, SFXN2, dhs-27) were validated by real-time PCR. For L4, five genes F09E10.7,VIT) were validated by real-time PCR from 11 highly ranked differentially transcribed genes. Some metabolic enzymes-related genes (e.g. SI, dhs-27) were more highly transcribed in L3 than L4, while some (e.g. PEPT-1) were more highly transcribed in L4 in this study. Many of these transcripts have orthologues in other species, which are so far uncharacterized and therefore await further study.
Allergic (Nieuwenhuizen and Lopata, Baird et al., 2016). Currently, 15 allergens in total have been described from A. simplex (s.s.) according to the WHO/IUIS Allergen Nomenclature Database (Pomés et al., 2018), but the exact number and characteristics of potential allergen molecules in anisakid nematodes are yet to be defined. Recently Baird et al. (2016) and Llorens et al. (2018) analyzed allergens from larval Anisakis transcriptomes and described a number of putative allergens including previously reported Anisakis allergens. In this study, we listed 56 transcripts matching to consensus sequences from AllergenOnline database (e-value cut off: <1e-05, identity matching cut off: >70%). Of them, 42 transcripts were matched with those of allergens from other organisms. This might be related with cross-reactivity of the sera from Anisakis-induced allergic patients to allergens from other sources, as described by Llorens et al. (2018). 16 transcripts from A. pegreffii transcriptomes in this study were matched with novel A. simplex (s.s.) allergens, suggesting some of them commonly exist in Anisakis species. DEG analysis for these putative allergen sequences revealed that there were no considerable differences in log FC value between APL3 and APL4, except for one transcript corresponding to Ani s 11-like protein 2 precursor gene, which needs further validation. Different expression of allergen profiles depending on the life cycle stages of anisakid nematodes are yet to be reported.
While 15 allergens (Ani s 1 to Ani s 14 and Ani s troponin C) have been described from anisakid nematodes (Pomés et al., 2018), 8 Anisakis allergens including 2 unassigned ones were found from the AllergenOnline Database in the A. pegreffii transcriptomes in this study. Several explanations can be provided for this result; some of these Anisakis allergens might be unregistered in the database when browsed. Otherwise, some of the transcriptome profiles in our study might have been inevitably missed out because we cut a small portion of caudal end of each individual A. pegreffii larvae for molecular identification and conducted transcriptome analysis with the rest of them. Cavallero et al. (2018) recently conducted transcriptome analysis with whole A. simplex (s.s.) and APL3, and specifically focused their pharyngeal tissue transcriptome. Our methodological approach might have affected our transcriptome dataset profiles, consequently underestimated the total number of allergens in our data set.
There are several classes of proteases in nematodes (e.g. cysteine, serine, aspartic, and metalloproteases), with a variety of functions, and some of them are thought to be well conserved across the nematode species such as molting, embryogenesis or cellular survival. On the other hand, proteases of parasitic nematodes are also involved in diverse adaptive functions including host tissue penetration, larval migration, immune evasion, host tissue digestion, and so forth  and the references therein). A: the real-time PCR results of five genes with significantly higher expression level in APL3 compared to APL4; B: the real-time PCR results of five genes with significantly higher expression level in APL4 compared to APL3; C and D are 14 additional genes selected to compare the expression levels of ASL3 and APL4 (C: eight gene for APL3 compared to APL4, D: six gene for APL4 compared to APL3). GAPDH genes were used as the reference genes. The results are given as mean of samples (n = 3) and expressed as fold change in mRNA expression. *P < 0.05. Different parasitic nematode species and different developmental stages in a certain nematode species may have different repertoire of proteases for adapting themselves to a specific niche, and comparative analyses of these differences would be essential to understand their pathogenesis and develop potential control strategies. For example, transcriptome analysis of a parasitic nematode Angiostrongylus cantonensis revealed that the infective third stage larvae highly expressed metallo-, aspartic-, and cysteine protease, while the fifth stage larvae highly expressed cysteine-, aspartic-, and serine proteases (Chang et al., 2014), suggesting these proteases have different roles in different developmental stages of parasites. Moreover, transcript profiles of aspartic protease genes (e.g., 113,155) in an entomopathogenic nematode Steinernema carpocapsae were characterized and different expression profiles of each transcript were found in the different developmental stages (Balasubramanian, Nascimento, Ferreira, Martinez, and Simões, 2012;Balasubramanian, Toubarro, Nascimento, Ferreira, and Simões, 2012;Balasubramanian and Simões, 2013). Of genes investigated in this study, many cysteine protease (Atg4c, atg4a, clp-1 G1), metalloprotease (TM104_CAEBR, ccpp-6, dpy-31), aspartic protease (asp-6 G1, asp-6 G2), hyaluronidase (HYAL1), enolase (enol-1)-related genes were found to be highly transcribed in APL3 by real-time PCR, while clp-1 G2 was validated to be highly transcribed in APL4 (Fig. 2). The third stage larvae of parasitic nematodes should successfully establish infection in the final hosts, and have many strategies to evade host defense mechanisms. The proteases secreted in third stage larvae of parasitic nematodes are therefore thought to work as one of the strategies to evade host defense mechanisms. Several proteases of A. simplex have been reported (Sakanari et al., 1989, Sakanari andMcKerrow, 1990;Hotez et al., 1994). Lee et al. (2017) reported strong activities of matrix metalloproteinases and serine proteases in excretory/secretory product of A. pegreffii infective third stage larvae, which is thought to support our findings. In particular, Cavallero et al. (2018) recently analyzed pharyngeal tissue transcriptomes of 2 sibling Anisakis species (A. simplex (s.s.) and A. pegreffii) and found upregulated transcripts involving pathogenicity (e.g. proteolysis, anesthesia, hemostasis inhibitors, anticoagulant, virulence factors, and immunomodulatory peptides) in the pharyngeal transcriptomes, when compared with the corresponding whole larvae transcriptomes. They suggested that these DEGs are potentially involved in host tissue invasion and pathogenesis. In this study, we also found different DEG repertoires between APL3 and APL4, mainly proteolytic enzymes and cuticle synthesis as mentioned above.
While we report comparative transcriptomic analyses of different developmental stages of A. pegreffii in this study, there are still much gaps in our knowledge of many DEGs genes either in L3 or in L4. Currently, the draft genome of A. pegreffii is not available. The draft whole genome of a sibling species A. simplex was recently uploaded (https://www.ncbi. nlm.nih.gov/assembly/GCA_900617985.1), yet to be completely annotated. Thus, the number of total orthologue genes in A. pegreffii transcriptome and DEGs either in L3 or in L4 might be underestimated. More precise information on the predicted gene functions of A. pegreffii transcriptomes will be obtained when the reference genome sequences are available and in vitro experiments are conducted.