Molecular and Cellular Response to Experimental Anisakis pegreffii (Nematoda, Anisakidae) Third-Stage Larval Infection in Rats

Background: Anisakiasis is a zoonotic disease caused by accidental ingestion of live Anisakis spp. third-stage larvae present in raw or undercooked seafood. Symptoms of this emerging infectious disease include mild-to-severe abdominal pain, nausea, and diarrhea. Some patients experience significant allergic reactions. Aims: In order to better understand the onset of anisakiasis, we aimed to: (i) histopathologically describe severe inflammatory/hemorrhagic infection site lesions in Sprague-Dawley rats experimentally infected with Anisakis pegreffii larvae; and (ii) qualitatively and quantitatively characterize the transcriptomes of affected tissues using RNA-Seq. Methodology: The experiment was performed on 35 male rats, sacrificed at 5 time points (6, 10, 18, 24, and 32 h post-infection). Gastric intubation was performed with 10 A. pegreffii larvae (N = 5 infected rats per time point) or 1.5 ml of saline (external control N = 2 rats). 16 pools, seven for muscle tissues and nine for stomach tissues, were created to obtain robust samples for estimation of gene expression changes depicting common signatures of affected versus unaffected tissues. Illumina NextSeq 500 was used for paired-end sequencing, while edgeR was used for count data and differential expression analyses. Results: In total, there were 1372 (855 up and 517 down) differentially expressed (DE) genes in the Anisakis-infected rat stomach tissues, and 1633 (1230 up and 403 down) DE genes in the muscle tissues. Elicited strong local proinflammatory reaction seems to favor the activation of the interleukin 17 signaling pathway and the development of the T helper 17-type response. The number of DE ribosomal genes in the Anisakis-infected stomach tissue suggests that A. pegreffii larvae might induce ribosomal stress in the early infection stage. However, the downstream pathways and post-infection responses require further study. Histopathology revealed severe inflammatory/hemorrhagic lesions caused by Anisakis infection in the rat stomach and muscle tissues in the first 32 h. The lesion sites showed infiltration by polymorphonuclear leukocytes (predominantly neutrophils and occasional eosinophils), and to a lesser extent, macrophages. Conclusion: Understanding the cellular and molecular mechanisms underlying host responses to Anisakis infection is important to elucidate many aspects of the onset of anisakiasis, a disease of growing public health concern.


INTRODUCTION
Tissue-dwelling intestinal nematodes encompass many diverse species, such as members of genera Trichuris (whipworms), Strongyloides (threadworms), Necator, Ancylostoma (hookworms), Ascaris, Anisakis, and Trichinella. Each species occupies a distinct host microenvironment, affecting, in turn, the type and effectiveness of the host immune response, the latter also being greatly influenced by the parasites' life cycle and behavior (1).
Typically, the host immune response to helminth infection is characterized by the T helper 2 (Th2) response, which involves the production of cytokines interleukin 3 (IL3), IL4, IL5, IL9, and IL13, eosinophilia, goblet, and mast cell hyperplasia, as well as alternatively activated macrophages (2). Despite the dominant Th2 phenotype, the central players in the course of helminth infection are naïve CD4 + T cells that can differentiate into several different regulatory and effector lineages (3). In addition to the pool of naïve T-cell precursors, the relatively plastic nature of Th subsets enables them to switch phenotype, mainly under the influence of varied cytokines (4). This ebb and flow of Thcell populations in response to helminths has a major role in determining if the elicited immune kinetics will contribute to host protection or pathological onset of inflammation (3).
It is assumed that the immune response to tissue-dwelling intestinal nematodes fails to be modulated or down-regulated in an atypical or accidental host (5). An interesting experimental model of accidental infection in humans is rat anisakiasis caused by third-stage (L3) A. pegreffii larvae (6). In humans, L3 Anisakis larvae are either unable to molt to the subsequent stage (L4) and reach reproductive maturity or molting to L4 occurs extremely rarely (7). However, in the course of infection, patients can suffer significant clinical symptoms (8). The estimated frequency of anisakiasis in Japan is approximately 2000 cases/year and in South Korea 200 cases/year, while for some European countries estimates range from 20 (9) to 500 cases/year (10). However, a recently published quantitative risk assessment analysis estimated the risk of anisakiasis from the consumption of anchovies to be between 7700 and 8320 cases annually in the Spanish population (11). Thus, suggesting a considerable underestimation of the current number of cases in Europe. In addition, a multi-criteria decision analysis of foodborne parasites by European region ranked Anisakidae within the top 10 and top 5 in Northern and South-Western Europe, respectively (12).
Anisakiasis is caused by ingestion of live Anisakis L3 larvae in raw or undercooked seafood and it can develop according to the following scenarios: (i) expulsion of larvae with stool and/or vomit with no associated lesions of the gastrointestinal mucosa; (ii) attachment/embedding/penetration of larvae, most commonly affecting gastrointestinal tissues, or less commonly tissues encountered during larval migration through the peritoneal/pleural cavities, mesentery, liver, pancreas, lymph nodes, ovary, and subcutaneous tissues (13)(14)(15). Therefore, the disease is commonly classified as gastric, intestinal, and ectopic (or extra-gastrointestinal) (13,14). In addition, ingestion of live Anisakis L3 can lead to sensitization and subsequent allergic reactions in humans. In a small number of cases, sensitized individuals may react even after ingestion of dead larvae (16). This clinical form is termed allergic/gastroallergic anisakiasis.
In a well-developed case of anisakiasis, fibrosis and eosinophilic infiltration into tissues surrounding the parasite are the most distinctive features of the local inflammatory response. Eosinophilic granulomas are observed in humans and infection models (17)(18)(19)(20). Other cell lineages usually reported in histopathological examinations include mast cells, lymphocytes, and plasma cells (18). Additional component of the innate response includes Toll-like receptors (TLRs) on epithelial cells and activated dendritic cells (17).
To better understand accidental Anisakis infection in humans and subsequently develop effective intervention strategies, it is important to understand molecular and cellular mechanisms underlying host responses to this uncommon infection. Rodent models of intestinal nematode infection have proved relevant in mapping the cellular and molecular basis of mammalian protective mechanisms (21). In addition, previous research confirmed that Anisakis infection of rats leads to similar pathologies as humans (16). Therefore, the aim of this study was to: (i) describe histopathological changes at the infection site of A. pegreffii larvae in Sprague-Dawley rats; and (ii) qualitatively/quantitatively characterize the transcriptomes of affected tissues using RNA-Seq.

Animal Ethics
All animal experiments and protocols were approved by the Ethical Committee of the School of Medicine at the University of Split (registry number 2181-198-03-04-18-004), as well as the Veterinary and Food Safety Office of the Ministry of Agriculture (registry number 525-10/0255-16-7). Rat experiments were performed at the University of Split Animal Facility (permit number HR-POK-019) where they were raised and housed in pairs, in plastic cages with sawdust and corn bedding. The animals were kept in a controlled environment: food and water ad libitum, temperature 22 ± 1 • C, with a 12 h light/dark cycle. The animals were separated in individual cages 24 h prior to the experiments and deprived of food 24 h before experimental infection.

Larvae Collection and Rat Infections
A. pegreffii larvae were collected from blue whiting Micromesistius poutassou, freshly caught in the C1 fishing zone of the Adriatic Sea (FAO 37.2.1), provided by a trusted dealer on the morning of the preliminary and experimental infection and delivered on ice. Gastric intubation and probe preparation protocols were adopted from the literature (20). Briefly, actively moving larvae were washed several times in physiological saline solution and checked under an Olympus BX 40 light microscope (Olympus Corp., Shinjuku, Tokyo, Japan) to confirm no cuticle damage had occurred during collection. Selected larvae were placed in previously prepared gastric probes. Each probe contained 10 larvae.
A preliminary experiment was performed to assess the temporal dynamics of Anisakis larvae infection. Samples were collected from 15 female Sprague-Dawley rats (average weight 197 ± 13.6 g, 3 animals per time point) at 6, 10, 24, 48, and 72 h post-infection.
Based on the results of the preliminary experiment, the duration of the in vivo experiment was set at 32 h post-infection. The experiment was performed on 35 male Sprague-Dawley rats in total (average weight 207 ± 20.1 g), with seven rats sacrificed per time point (6,10,18,24, and 32 h post-infection). The experimental procedure remained the same.
Seven rats per group were administered a mixture of anesthetic and analgesic; Ketaminol (Richter Pharma AG, Wels, Austria), 50-100 mg/kg and Xylapan (Vetoquinol UK Ltd, Buckingham, UK), 5-10 mg/kg, by intraperitoneal injection, with additional administration of Ketaminol (Richter Pharma AG, Wels, Austria), 50-100 mg/kg, if toe pinch reflex was detected. Next, gastric intubation was performed with 10 A. pegreffii larvae (infected N = 5 rats per time point) or 1.5 ml of physiological saline solution (external control N = 2 rats). All animals were euthanized by decapitation, following an overdose of anesthesia by > 150 mg/kg Ketaminol (Richter Pharma AG, Wels, Austria) at designated time points. Tissue samples and recovered A. pegreffii larvae were collected, washed in physiological saline solution and immediately stored in Tri Reagent (Ambion Inc., Invitrogen, Carlsbad, CA, USA) at −80 • C.

Tissue Preparation for Semi-Thin Sections
For histopathological analysis, small fragments of the Anisakisinfected stomach (N = 5) and abdominal muscle tissues (N = 2) were collected and fixed in 4% paraformaldehyde in 0.1 M phosphate buffered saline (PBS) on ice. Tissue samples were postfixed in 1% aqueous osmium tetroxide for 1 h, dehydrated in ascending series of acetone, and embedded in Durcupan resin (Honeywell-Fluka, Morris Plains, NJ, USA). Semi-thin 0.5 µm sections were cut, stained with 1% toluidine blue, and examined under an Olympus BX 40 light microscope (Olympus Corp., Shinjuku, Tokyo, Japan). Images were captured with an Olympus Camedia camera (Olympus Corp., Shinjuku, Tokyo, Japan) and assembled with Photoshop CS 5 software (Adobe Systems, San Jose, CA, USA).

RNA and DNA Extraction
Total RNA was extracted using Tri Reagent (Ambion Inc., Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. Stomach and muscle tissues were selected when severe inflammatory/hemorrhagic lesions with or without migrating Anisakis larvae were observed. In addition, adjoining unaffected stomach and muscle tissues from the same rats were collected as the internal controls (referred to as unaffected tissues or control) and used in all statistical analyses for comparison with affected tissues. Details on stomach and muscle tissue samples are provided in Supplementary Table 1. Muscle tissues include abdominal, dorsal, intercostal muscle, and thoracic diaphragm. Stomach encompassed greater curvature, fundus, pyloric antrum, and pyloric canal. Hereafter, we will refer to our sequenced pools as "stomach" and "muscle" (tissues) when referring to RNA-Seq. In addition, RNA extraction was performed from the same tissues of uninfected rats (the external controls, referred to as uninfected, to avoid confusion with internal controls).
Additionally, DNA was extracted from 14 larvae found in the process of stomach/muscle penetration to confirm Anisakis species. DNA was extracted from Tri Reagent (Ambion Inc., Invitrogen, Carlsbad, CA, USA) after RNA extraction (data not shown) following the manufacturer's protocol and used for PCRbased restriction fragment length polymorphism (PCR-RFLP) analysis of the ribosomal DNA (rDNA) internal transcribed spacers (ITS-1 and ITS-2). RFLP pattern characteristic of Anisakis simplex (sensu stricto) × A. pegreffii putative hybrid (620-370-300-250 bp) (22) were observed for a single larva. All other larvae were confirmed as A. pegreffii according to an established RFLP pattern (370-300-250 bp) (22).

cDNA Preparation and Illumina Next-Generation Sequencing
Total RNA was dissolved in 20-40 µl of RNase/DNase free water (Merck Millipore, Billerica, MA, USA) and shipped on dry ice to the laboratory for advanced genomics at the Ruder Bošković Institute, Croatia, which provided the sequencing service. RNA concentration, purity, and integrity were determined using a 2100 BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit 3.0 (Thermo Fisher Scientific, Waltham, MA, USA). Based on sample quality and specific lesion sites, 16 pools of at least three biological replicates, seven for muscle tissues and nine for stomach tissues, were created to obtain robust samples for estimation of gene expression changes depicting common signatures of affected vs. unaffected tissues (Supplementary Table 1). Following the manufacturer's protocol, the cDNA library was prepared using TruSeq Stranded mRNA kit (Illumina, San Diego, CA, USA) and subsequently sequenced using Illumina NextSeq 500 platform (Illumina, San Diego, CA, USA) over four lanes.

RNA-Seq Raw Reads Pre-processing and Mapping
The mean number of paired-end reads generated per sample was 31.2 million (range from 16.4 to 37.6 million reads). The quality assessment of the reads derived from different lanes using FASTQC (Babraham Bioinformatics, Babraham Institute, Cambridge, UK) indicated good data quality and no lane effects. Thus, the reads were joined into two paired FASTQ files per sample. Higher quality was observed for forward reads, as generally reported for Illumina sequencing (23). Trimmomatic (24) was used to trim Illumina adapter sequences, perform the sliding window clipping of reads (quality threshold 20, window size 4), and to remove reads shorter than 30 bases. Subsequently, PRINSEQ (25) was used to remove low-complexity reads (entropy threshold 70) and reads with more than 10% ambiguous (N) nucleotides. On average, 87% of paired-end reads survived this procedure and were mapped to Rattus norvegicus (v6) genome, Ensembl release 91 (26) using STAR (27). STAR (27) was also used to perform read counting. Approximately 90% of the reads were uniquely mapped. All data were submitted to NCBI Sequence Read Archive (SRA) under the accession number of SRP150499 (project ID PRJNA475982).

Differential Expression Analyses
EdgeR (28,29) package for R software ver. 3.4.2 (30), Bioconductor release 3.6 (31) was used for the count data exploratory analyses, log2-counts-per-million (logCPM) normalization, and inference of differential gene expression. Based on different sample variances and discrete clustering after principal component analyses (PCA), samples from stomach and muscle tissues were analyzed separately. Samples from uninfected rats were not included in statistical analyses, as they were represented by a single pool per dataset and were used only for comparative purposes in illustrations.
Low expression features (with less than one count per million per sample and at least three in any experimental condition) were filtered out prior to analysis. The dataset was mapped to Entrez Gene identifiers using rat annotation package org.Rn.eg.db (32). Statistical testing for differential expression was initially performed using a paired design with FDR cutoff of 0.05 and subsequently, relative to a minimum foldchange threshold of 2, using the edgeR implementation of t-tests relative to a threshold (TREAT) (33) method under the generalized linear framework. TREAT is an extension of the empirical Bayes moderated t-statistic presented by Smyth (34) and it achieves reliable p-values and FDRs for finding genes with differential expression that is biologically meaningful by including the fold-change threshold of interest in a formal hypothesis test. In this manner, we focused on genes with stronger regulation. Heatmap representations of the expression patterns of differentially regulated genes were plotted using ComplexHeatmap package (35).

Enrichment Analysis
The test for over-representation of Gene Ontology (GO) terms within differentially expressed genes (FDR < 0.05, |LogFC| >> 1) was conducted using edgeR's in-built goana function with FDR cut-off set at 0.05, taking into account gene length bias. The top 10 terms for stomach and top six for muscle associated with Biological Processes were selected and plotted using GOplot package (36) based on ggplot2 (37). The functional analysis of KEGG signaling, metabolic, and disease pathways (38) was conducted using GAGE (Generally Applicable Geneset Enrichment), as implemented in the gage package (39). Log (base 2)-fold changes resulting after the linear model fit with edgeR were used as per gene score, and a two-sample t-test was run to compare mean log-fold changes per set relative to mean background log-fold change calculated using all filtered genes. Tests were run for up-, down-, and bi-directionally perturbed pathways using the default q-value cut-off of 0.1. Briefly, in results, we focused only on the metabolic and signaling sets with up-or downregulation that satisfied more stringent statistical significance criteria (q-value < 0.01). Pathview package (40) was used to visualize log-fold changes recorded for differentially expressed genes (FDR < 0.05) on their respective KEGG pathway maps, for pathways deemed significantly perturbed in GAGE analyses.

Time Course of Infection
The preliminary experiment showed that the A. pegreffii infection duration (period from application of larvae until expulsion from the host) in rats was relatively short ( Table 1). Six h post-infection, a fraction of the larvae had already left the host (between 80 and 90% were recovered from within the host in preliminary and experimental infections, respectively) and 24 h post-infection approximately 30% of applied L3 were found within the host. Larvae showed no synchronized migratory behavior correlated with time post-infection. Two larval clearance routes were observed. Route one was passage through the digestive tract facilitated by peristalsis with no apparent tissue damage. Route two was by penetration through the stomach, small intestine, caecum, or large intestine, with or without subsequent migration through the abdominal wall muscle. In a few cases, L3 larvae were found inside the pelvic cavity, both in preliminary and experimental infections. Epidermal penetration was also observed, but only in the preliminary experiment. In the majority of cases, L3 larvae were observed penetrating different parts of the stomach mucosa and abdominal muscles. Therefore, these tissues were the main subject of our downstream molecular and histopathological analyses.

Histopathology of A. pegreffii Experimental Infection
In the stomach, A. pegreffii larvae were found in the early phase of migration through the gastric wall, causing compression and necrosis of surrounding parietal and zymogenic cells in the mucosa ( Figure 1A). Several strata of polymorphonuclear leukocytes (predominantly neutrophils and occasional eosinophils), and to a lesser extent, macrophages, were observed infiltrating penetration site. However, these cells appeared not to be in the direct contact with larval cuticle, probably due to tissue shrinkage during sample processing. Moreover, the rupture of both basement membrane and muscularis mucosae occurred at the site of larval migration. Close to the two adjoining arteries in the submucosa, a large necrotic area admixed with scant inflammatory infiltrate and mild hemorrhage was observed (Figures 1A,B, asterisk). However, extensive hemorrhages, visible at gross pathological examination, were also noticeable in the submucosal layer. The submucosa was predominantly infiltrated with polymorphonuclear leukocytes (predominantly neutrophils), and macrophages, interspersed between connective tissue fibers, fibroblasts, and fibrocytes. Extravasated polymorphonuclear leukocytes concentrated toward the necrotic area adjoining blood vessels and the migrating larva ( Figure 1C). Macrophages were seen in deeper submucosa strata, closer to the large submucosal hemorrhage ( Figure 1D).
In the abdominal musculature, penetrating A. pegreffii larva caused tissue fragmentation and necrosis along its migration path (Figure 2A). The perimysial connective tissue exhibited extracellular edema and moderate hemorrhage with inflammatory infiltration, the latter composed of polymorphonuclear leukocytes (predominantly neutrophils), and macrophages ( Figure 2B). A few polymorphonuclear leukocytes were also seen in the endomysium or between myocytes at the site of larval penetration. Higher magnification revealed structural disintegration and necrosis of muscle fibers surrounding the larva, with an accumulation of abundant cellular debris mixed with numerous polymorphonuclear leukocytes and macrophages ( Figure 2C).

Differentially Expressed Genes
Consistent with the non-synchronized migrating behavior of L3 larvae during the experiment and no preference for penetration site within inspected tissues; somewhat higher variability was observed between samples within experimental groups using PCA (Supplementary Figure 1). Nevertheless, statistical analysis detected common gene regulation patterns associated with Anisakis infection in the stomach and muscle tissues.  (Figure 3). The expression profiles of the uninfected samples for both datasets were consistent with control samples (Figure 3). These include S100 proteins, namely S100a8 (4.68 logFC in stomach and 4.62 in muscle) and S100a9 (3.53 logFC in stomach and 4.41 in muscle). S100A8 is a calcium-and zinc-binding protein with an important role in the regulation of inflammatory processes and the immune response and is also involved in orchestrating chemotaxis and adhesion of neutrophils. S100A8/A9 is predominantly present as calprotectin, displaying a plethora of intra-and extracellular functions. S100A8/A9 proinflammatory activity includes promotion of cytokine and chemokine production, and evidence for this activity can be supported by a number of expressed Ccl (Ccl2, Ccl3, and Ccl7) and Cxcl (Cxcl1, Cxcl2, and Cxcl6) chemokines (Figure 3) in Anisakis-infected tissues. Another highly expressed gene in both stomach and muscle tissues was matrix metallopeptidase 3 (Mmp3, 4.43 logFC in stomach and 2.89in muscle), which is involved in degradation of fibronectin, laminin, collagens and cartilage proteoglycans, and activation of procollagenase.

Functional Interpretation of Transcriptomic Changes Induced by Anisakis Infection
From the 259 overrepresented GO terms (FDR < 0.05) associated with the 41 top DE genes for Anisakis-infected rat stomach, the top five ranked were response to stress, defense response, inflammatory response, response to cytokine, and neutrophil chemotaxis (Figure 4, Supplementary Table 3). The presence of hemoglobin-related genes [hemoglobin alpha 1 (Hba1) and 2 (Hba2) and hemoglobin subunit beta-1 (Hbb)] was directly related to response to stress. Several regenerative family member genes (Reg3a, Reg3g, Reg3b, and Reg1a), with a particularly pronounced upregulation in stomach (Figure 4 and Supplementary Table 2), were also associated with other inflammation, immune, and defense-related sets. The remainder of the top five GO terms had many genes in common. In addition to S100a8/S100a9 and Ccl/Cxcl chemokines, interesting transcripts encompassed lipocalin 2 (Lcn2), involved in multiple processes such as apoptosis and innate immunity, and chitinase 3-like 1 (Chi3l1), which plays a role in Th2 inflammatory response and IL13-induced inflammation.
The top five ranked overrepresented GO terms for the Anisakis-infected rat muscle, out of 494 selected (FDR < 0.05), were very similar to the stomach, but ranked differently: defense response, inflammatory response, response to external stimulus, granulocyte migration, and neutrophil migration (Figure 5 and Supplementary Table 3). The response to stress was also found in the overrepresented GO terms list for the muscle, although ranked 21st according to FDR (Supplementary Table 3). The largest set of DE genes was found in the response to external stimulus, followed by the defense response ( Figure 5). Some of the genes were associated only with the response to external stimulus term, such as activated leukocyte cell adhesion molecule (Alcam), which promotes T-cell activation and proliferation via its interactions with Cd6. However, most of the genes were common to the top five GO terms. In addition to the omnipresent S100a8/S100a9 and Ccl/Cxcl chemokines, two other potent proinflammatory cytokines are worth mentioning; interleukin 1β (Il1β) and Il6, the latter known as the inducer of the acute phase response.
The complete KEGG pathway analysis results are shown in Supplementary Table 4, including metabolic, signaling, and disease sets for the infected stomach and muscle tissues undergoing Anisakis larval migration. Metabolic and signaling sets with up-or downregulation that satisfied more stringent statistical significance criteria (q-value < 0.01) are shown in Figure 6. The top five overrepresented KEGG metabolic and signaling pathways upregulated in the rat stomach included 03010 Ribosome, 04610 Complement and coagulation cascades, 04060 Cytokine-cytokine receptor interaction, 04657 IL17 signaling pathway and 04640 Hematopoietic cell lineage (Figure 6 and Supplementary Table 4). Except for 03010 Ribosome (Figure 7), the rest of the top five KEGGs were also found in the infected muscle, however, not all appeared in the top five list (Supplementary Table 4). In addition to, 04060 Cytokine-cytokine receptor interaction and 04640 Hematopoietic cell lineage, the top five upregulated KEGG pathways for rat muscle included 04062 Chemokine signaling pathway, 04514 Cell adhesion molecules (CAMs), and 04145 Phagosome (Figure 6 and Supplementary Table 4). The only downregulated KEGG pathway in the rat muscle was 00190 Oxidative phosphorylation. According to KEGG functional categories, the greatest number of differentially perturbed pathways belongs to the immune system, followed by signal transduction and signaling molecules and interaction (Figure 6). The pathways with the greatest number of DE genes (N = 12) included in both the stomach and muscle were 04657 IL17 signaling pathway (Figures 6, 8) and 04060 Cytokine-cytokine receptor interaction (N = 10). 04657 IL17 signaling pathway is shown only for muscle data (Figure 8), as there were only slight differences between the two datasets.

DISCUSSION
This is the first RNA-Seq transcriptomic analysis of vertebrate host tissues affected by Anisakis larval migration. We chose Sprague-Dawley rats as the host model to simulate the course and pathogenesis of accidental infection in humans. In addition, molecular identification of the anisakid confirmed that A. pegreffii was the parasite in our study, which is important to consider, as it is the most frequent causative agent of anisakiasis in Mediterranean countries (18,19,41,42). An in vivo and in vitro study in the same animal model confirmed the pathogenic potential of A. pegreffii The study also found that A. pegreffii larvae tend to remain in the stomach of rats much longer than Anisakis simplex (sensu stricto) (6). Consistent with our study, the authors concluded that the migratory behavior of L3 larvae was non-synchronous and L3 larvae showed no preference for penetration site within inspected tissues. However, the authors did not report muscle penetration, probably due to the shorter time-span of the experiment (24 h) (6). While RNA-Seq has previously been applied to identify and characterize putative novel parasite allergens in Anisakis (43) and to reveal the pathogenic mechanisms associated with infection in two genus members, A. simplex (s. s.) and A. pegreffii (44), no description of the host transcriptomic response has yet been given.
Even though anisakids are not natural human parasites and they are not evolutionary adapted to each other (5), both human and murine models of Anisakis infection have been predominantly characterized in the context of Th1, Th2, or Th1/Th2 immune responses inferred through different serological studies aimed at defining the kinetics of allergic reactions to the nematode after different sensitization regimes (16,45,46). By assessing immune response traits in L3-infected Anisakis-naïve rats, we have identified an acute proinflammatory host response, which seems to be predominantly characterized by the activation of a Th17 lineage of effector T helper lymphocytes. This lineage is known for the production of IL17 (47), a cytokine that induces secretion of chemokines and antimicrobial peptides in diverse tissues and cell types (fibroblasts, endothelial cells, epithelial cells, keratinocytes, and macrophages), consequently leading to the recruitment of neutrophils and an enhanced proinflammatory reaction (48). Indeed, we found both neutrophils and macrophages mixed with necrotic cells and cellular debris in our study. Necrosis can stimulate inflammation due to leakage of intracellular components into surrounding tissue (49). Such a scenario can be more pronounced in tissues with high metabolic activity, such as the stomach epithelium, which is rich with parietal and zymogenic cells producing hydrochloric acid and various digestive enzymes, respectively (50). Necrosis and eosinophils are usually present in the center of granulomas caused by Anisakis larval remains, together with foreign-body giant cells and connective tissue. Subsequently, the cellular infiltration encompasses primarily lymphocytes, as memory T and B cells are generated during primary infection of the host. Thus sensitized, hosts can experience allergic reactions upon re-exposure to Anisakis (16).
In previous allergy models of Anisakis infection in C57BL/6 mice, significantly increased Th17-related cytokines IL6 and IL17A were observed after subchronic treatment (7 times over 15 days) with recombinant Ani s 1 (rAni s 1) allergen (51). In vivo results were further supported by proliferation and activation of IL6 and IL17A secretion in splenocytes treated with rAni s 1 or rAni s 9, suggesting that IL17 might play a critical role in the Anisakis-associated allergic reaction (51). However, in another study, Wistar rats exposed orally to fresh and frozen Anisakis larvae (treated twice over 7 days and sampled after 21 days) had a significantly higher level of plasma IL17. Notably, intraperitoneal injection of crude larval extract did not induce significant IL17 production (52). This might suggest that live larvae trigger different immune recognition and activation pathways than crude extracts or even secretory/excretory products (53). This hypothesis is supported by a recent study where dendritic cell (DC) cytokine and chemokine secretion was conditioned by live A. pegreffii larvae or crude extract (54). Interestingly, the autologous CD4 + T cells that were stimulated by DCs differentiated in the presence of live larvae or crude extract but failed to produce IL4, IL17, and IL10. Observing no IL17 production in both scenarios, the authors concluded that DCs may contribute to establishing localized inflammation at the earliest stage of infection, but their role in Th2/Th17 or T regulatory (Treg) polarization remains inconclusive. Similarly, Th17 and IL17 signaling pathways seem to be strongly induced in our study, without regulation of the Il17 gene.
Conversely, Th17 cells are known to produce cytokines that are specialized in promoting responses against extracellular bacteria and fungi (55). However, Trichinella spiralis, another zoonotic nematode suggested as the best terrestrial counterpart of the tissue-dwelling marine Anisakis, can elicit both Th2and Th17-related immune responses during its intestinal phase, which occurs after the first week of infection (56). The life cycle of T. spiralis in humans begins with the consumption of meat containing the first-stage muscle larvae. After larval release from the cyst and penetration into epithelial cells of the small intestine, the intestinal phase begins, where larvae molt, grow and reproduce (57). However, at the onset of nurse-cell formation (induced by the parasite approximately 1-week post-infection), Th2 and Th17 responses are inhibited by Treg cell recruitment to the peripheral lymph nodes near the nurse cells 2-weeks post-infection (56). This shows that Th2 and Th17 play crucial roles in the initiation of the inflammatory response to nematode infection, as observed in our study, which can subsequently develop in another direction as infection becomes chronic. Interestingly, we found upregulated diseases-related KEGG pathways in the rat stomach and muscle that are principally associated with strong inflammatory responses, such as 05323 Rheumatoid arthritis and 05321 Inflammatory bowel disease (Supplementary Table 4), corroborating the observed Th17 response in experimental Anisakis infection. As Th17 cells are potent inducers of tissue inflammation, they are also associated with the pathogenesis of many autoimmune diseases and human inflammatory conditions. In experimental animal models, Th17 cells are responsible for autoimmune encephalitis, collageninduced arthritis, colitis, and psoriasis (58)(59)(60). Furthermore, the numerous bacteria-related diseases significantly upregulated in our stomach and muscle dataset, such as rno05150 Staphylococcus aureus infection and rno05132 Salmonella infection (Supplementary Table 4), are most likely due to the extensive study of the Th17 response in bacterial infections.
The frequent appearance of the bacteria-related host immune responses in our results posits another possible explanation. The host intestinal mucosal immune system is a unique environment that must restrain immune responses directed against commensal microbes and dietary antigens while also protecting mucosal tissue during nematode parasite infection (1). It appears that commensal microbiota are physically separated from intestinal surfaces and this barrier is critical for limiting immune activation and maintaining homeostasis and mutualistic host-microbe associations (61). Furthermore, perturbed spatial FIGURE 5 | GOChord plot of top six ranked overrepresented GO terms belonging to the Biological Process subontology for Anisakis-infected rat muscle. The genes are linked to their assigned terms via colored ribbons. Genes are ordered according to the observed log-fold change (logFC), which is displayed in descending intensity of red squares displayed next to the selected genes.
relationships between microbiota and host correlate with disease states such as inflammatory bowel disease (62). RegIIIγ, a secreted antibacterial lectin that we identified as the top DE gene in the rat stomach tissues (logFC 6.25), is essential for maintaining a ∼50 µm zone that physically separates the microbiota from the small intestine epithelium (63). In addition, loss of the host-bacterial segregation in RegIIIγ −/− mice results in increased bacterial colonization of the intestinal epithelium and activation of intestinal adaptive immune responses by the microbiota (63). It is plausible that the tissue damage induced by Anisakis penetration into the stomach mucosa disrupted the delicate host-microbe microenvironment, resulting in the host immune response being targeted toward bacteria crossing the mucosal barriers, whether of autochthonous origin or carried over by contact with Anisakis larvae.
Another important inflammation-related pathway upregulated in our study is the KEGG 04657 IL17 signaling pathway, implicated by 12 DE genes in both the stomach and muscle tissues of rats. IL17A and IL17F have protective roles in the host mucosal barrier defense mechanism against certain pathogens. Although mostly related to bacterial infections (64), they have also been associated with filarial nematodes (65) as well as other nematode infections (66). Both stomach and muscle datasets had S100 proteins (S100a8 and S100a9), matrix metallopeptidase (Mmp3 and Mmp13), and Cxcl and Ccl chemokines (Cxcl1, Cxcl2, Ccl2, Ccl7, and Ccl20) from the IL17 signaling pathway. Il6 was present only in the muscle data and lipocalin 2 (Lcn2) was present only in the stomach dataset. S100A8/A9 proinflammatory functions were manifested through leukocyte recruitment, increased cytokine and chemokine production, and regulation of leukocyte adhesion and migration in our study (67). Nematode infections have previously been reported as the reason behind neutrophil recruitment in response to high S100A8 and S100A9 protein expression levels. For example, neutrophil recruitment was significant in a transient lung inflammatory response caused by the migratory phase of Litomosoides sigmodontis filarial infective larvae in BALB/c mice (68). The intracellular and extracellular functions of S100a8/a9 in the rat stomach and muscle were further supported by observation of neutrophil/leukocyte chemotaxis and migration, cytokine and chemokine production, upregulation of TLR and NF-κB signaling pathways (67), supporting an intensive inflammatory onset timed by 32 h post-infection in our study.
Interestingly, no evidence to support active wound healing was found within the transcriptomes of stomach and muscle tissues. However, it is possible that among the myriad MMPs that appeared in the DE stomach and muscle tissue gene lists (Mmp3, Mmp10, Mmp13, and Mmp19, to name some of the most significant), some were involved in tissue remodeling and repair in response to the injury (69). We abstain from inferring a generalized assumption on the subject, as a specific MMP secreted by one cell type (for example, a macrophage) probably performs a different function than the same enzyme produced by another cell type (for example, an epithelial cell) (70). However, MMP functions encompass inflammation and immunity, acting on pro-inflammatory cytokines, chemokines, and other proteins (70). In addition, it is believed that MMP13 plays a role in wound healing via a mechanism that probably involves activation of transforming growth factor beta 1 (TGFβ1) and degradation of connective tissue growth factor (CTGF), as well as in keratinocyte migration during wound healing (67). In our study, however, we believe that hemorrhagic lesions and IL1β and IL6 could have delayed the wound healing process (71) or the sampling occurred at too early in the process to fully capture it. Similarly, the appearance of fewer eosinophils than neutrophils is most likely the result of sampling in the peracute phase, which occurs prior to the typical eosinophilic infiltration of the tissue surrounding the parasite and before the Th2 response. According to the literature, eosinophilia does not typically develop until several days after the onset of clinical symptoms in anisakiasis (72).
Furthermore, we observed an intensive neutrophil and macrophage infiltration into the damaged tissues, which may have potentially biased consequences. Initially, their beneficial action results in removing the infectious agent, clearing debris, and expressing factors that promote wound healing. However, during the healing process, the control of inflammation may be crucial for effective tissue repair (73), with macrophages having an important role in the resolution of tissue damage, as observed in the course of viral infection (74). In addition, IL17 may also indirectly contribute to the early stages of the tissue repair process, as it serves as a mediator of successive neutrophil recruitment and activation (75). However, prolonged production of the IL17 also potentiates inflammation and tissue damage (76).
IL6 might be an important factor in Anisakis epidemiology (51) and it appeared in several differentially perturbed pathways in our study. In addition to its role in activation and differentiation of macrophages, lymphocytes, and terminal differentiation of B cells, IL6 also regulates acute and chronic inflammation (77). Murine studies have identified TGFβ and IL6 as key cytokines driving Th17 differentiation (78,79). Additionally, pro-inflammatory cytokines such as IL1β and tumor necrosis factor α (TNFα) can increase the efficiency of Th17 differentiation (80). Both Il6 and Il1β were DE in the Anisakis-infected rat muscle. However, the role of Il6 in our study should be interpreted with caution, as IL6 can also act as an anti-inflammatory myokine, which is elevated and discharged into the bloodstream in response to the muscle contraction (81). Hyperactivation of IL6/STAT3 signaling in the muscle may stimulate the synthesis and systemic release of the acute-phase response proteins such as serum amyloid A (SAA) and fibrinogen, which may amplify catabolic signals in the muscle, turning muscle into a key player in innate immunity in experimental cancer-induced muscle wasting in mice (82). Whether Anisakis infection triggers pro-or anti-inflammatory IL6 responses remains to be elucidated.
In addition to the remarkable proinflammatory response based on Th17 and IL17 signaling pathways and DE genes included in these pathways, we also observed the stress response. The highest-ranking GO term in rat stomach was the response to stress, but the strongest q value was unexpectedly designated to the KEGG 03010 Ribosome pathway. This feature was specific to Anisakis-infected stomach. Although none of the ribosomal genes in our study showed large expression differences, they were found consistently upregulated (FDR < 0.05) across the stomach dataset (Supplementary Table 2). A visualization of the KEGG Ribosome pathway differentially perturbed in the Anisakisinfected rat stomach is shown in Figure 7. This perturbation suggests that Anisakis infection could have triggered the onset of ribosomal stress during penetration of the stomach mucosa. Interestingly, most ribosomal proteins have been implicated in host immune responses by boosting immune signaling or facilitating pathogen proliferation under various circumstances (83).
In the host-parasite context, ribosomal stress has only been addressed in a single study of pig infection with Cryptosporidium parvum merozoites (early asexual stages), where both host and parasite responses were assessed by RNA-Seq. Interestingly, transcripts of the C. parvum-infected pig intestinal cells did not reveal stress-and apoptosis-related genes. However, genes encoding ribosomal functions were highly enriched in both host and C. parvum cells, although to a greater extent in the parasite than the host (84). The authors pointed out that transcriptomic changes observed might be influenced by rapid parasite proliferation forcing a high energetic and FIGURE 8 | Visualization of KEGG IL17 signaling pathway genes differentially perturbed (FDR < 0.05) in rat muscle infected with A. pegreffii. Plotted using Pathview package for R/Bioconductor. metabolic load on host cells. However, they did not associate it directly with ribosomal genes. Supported by the number of DE ribosomal genes in the Anisakis-infected rat stomach tissues and in light of their recently discovered activities in host immune responses, we argue that A. pegreffii larvae could induce such a response at the early infection stage in specific tissues. The downstream pathways and whether it is eventually balanced in the post-infection stages remains to be elucidated.

CONCLUSION
Our research represents an enhancement in understanding the early infection mechanisms of Anisakis larvae in an accidentalinfection model. The molecular results showed a strong localized proinflammatory response, favoring the activation of the IL17 signaling pathway and the development of the Th17type response, with a possible co-occurrence of ribosomal stress. This is further supported by the cellular findings showing severe inflammatory/hemorrhagic lesions. Gaining an improved understanding of the effector mechanisms at the host-parasite interface lays the foundations for more focused mechanistic studies of non-adapted host-parasite interactions.

AUTHOR CONTRIBUTIONS
IM and ŽT conceived and designed the study. Acquisition, analysis, and interpretation of data were performed by all authors (IBu, JH, AV, and IM: experimental infection and tissues sampling; IBu and JH: RNA isolation; AV: DNA isolation and RFLP; JH and IBo semi-thin sections; ŽT and IBu: RNA-Seq data analyses). The manuscript was drafted by IBu and all authors contributed to drafting sections of their area of expertise and revised it critically. All authors have read and approved the final content of the version to be published.