Intestinal Transcriptomic and Histologic Profiling Reveals Tissue Repair Mechanisms Underlying Resistance to the Parasite Ceratonova shasta

Background: Myxozoan parasites infect fish worldwide causing significant disease or death in many economically important fish species, including rainbow trout and steelhead trout (Oncorhynchus mykiss). The myxozoan Ceratonova shasta is a parasite of salmon and trout that causes ceratomyxosis, a disease characterized by severe inflammation in the intestine resulting in hemorrhaging and necrosis. Populations of O. mykiss that are genetically fixed for resistance or susceptibility to ceratomyxosis exist naturally, offering a tractable system for studying the immune response to myxozoans. The aim of this study was to understand how steelhead trout that are resistant to the disease respond to C. shasta once it has become established in the intestine and identify potential mechanisms of resistance. Results: Sequencing of intestinal mRNA from resistant steelhead trout with severe C. shasta infections identified 417 genes differentially expressed during the initial stage of the infection compared to uninfected control fish. A strong induction of interferon-gamma and interferon-stimulated genes was evident, along with genes involved in cell adhesion and migration. A total of 11,984 genes were differentially expressed during the late stage of the infection, most notably interferon-gamma, interleukin-6, and immunoglobulin transcripts. A distinct hardening of the intestinal tissue and a strong inflammatory reaction in the intestinal submucosa including severe hyperplasia and inflammatory cell infiltrates were observed in response to the infection. The massive upregulation of caspase-14 early in the infection, a protein involved in keratinocyte differentiation might reflect the rapid onset of epithelial repair mechanisms, and the collagenous stratum compactum seemed to limit the spread of C. shasta within the intestinal layers. These observations could explain the ability of resistant fish to eventually recover from the infection. Conclusions: Our results suggest that resistance to ceratomyxosis involves both a rapid induction of key immune factors and a tissue response that limits the spread of the parasite and the subsequent tissue damage. These results improve our understanding of the myxozoan–host dialogue and provide a framework for future studies investigating the infection dynamics of C. shasta and other myxozoans.


Introduction
Myxozoans are parasitic cnidarians characterized by a complex two-host, two-spore lifecycle that alternates between a myxospore stage that infects the invertebrate host (annelids or bryozoans) and an actinospore stage that infects the vertebrate host (generally fish) [1]. Myxozoans are widely distributed with over 2400 known species and have adapted to both freshwater and marine hosts [2]. Infections in the fish host are often asymptomatic, however, certain myxozoan species are known to cause severe pathology or death of the fish host. Mortality and morbidity as a result of myxozoan infections is becoming more common as aquaculture continues to expand worldwide, and more fish species become intensively managed by humans [3]. For example, throughout Asia and the Mediterranean, cultivation of turbot (Scophthalmus maximus) and several perciform fish is severely limited by enteromyxosis, caused by the myxozoan parasites Enteromyxum scophthalmi and E. leei, respectively [4,5]. Similarly, pharyngeal myxobolosis caused by Myxobolus hunghuensis is a limiting factor in the development of gibel carp (Carassius auratus gibelio) aquaculture in China [6]. Globally, over 35 species of marine fish are affected by Kudoa thyrsites, which encysts in muscle and causes postmortem myoliquefaction of the surrounding tissue, resulting in large economic losses due to poor fillet quality [7].
Rainbow trout and its anadromous form, steelhead trout (Oncorhynchus mykiss), are among the most widely cultivated species of fish and are negatively affected by several myxozoan pathogens, including K. thyrsites, Myxobolus cerebralis (the causative agent of whirling disease) [8], and Tetracapsuloides bryosalmonae, which is a related malacosporean and causes proliferative kidney disease [9]. Within the Pacific Northwest of the United States and Canada, O. mykiss and related salmonids are severely impacted by the intestinal myxozoan Ceratonova shasta (syn. Ceratomyxa shasta) which causes ceratomyxosis, a disease characterized by hemorrhaging and necrosis of the intestine and potentially death of the fish host. C. shasta has been linked to population level declines of Pacific Northwest wild fish stocks [10,11], and exerts such a selective pressure on the fish host that endemic populations became genetically fixed for resistance to the disease [12][13][14][15][16][17]. However, the parasite is not established in all river systems where salmonids are native and allopatric salmonids are highly susceptible to the parasite, and exposure to a single spore is capable of causing ceratomyxosis and mortality [18,19]. This leads to an almost binary resistance phenotype, where mortality either occurs after exposure to one spore (susceptible hosts), or thousands (resistant hosts) [11,20].
The combination of a highly virulent myxozoan, along with naturally occurring resistant and susceptible strains of the fish host, offers an attractive system for studying the immune response to myxozoan parasites. Our research group previously conducted a comparative transcriptomic analysis of resistant and susceptible phenotypes of steelhead trout exposed to a low dose of C. shasta, sufficient to cause mortality in the susceptible phenotype [21]. RNA-seq analysis revealed a downregulation of genes involved in the interferon-gamma (IFN-gamma) signaling pathway in the gills, the site of initial parasite invasion, of both phenotypes at 1 day post exposure (dpe) to C. shasta. By 7 dpe, resistant fish had effectively contained the infection, having either low or undetectable levels of C. shasta in their intestine, with several immune genes upregulated. In contrast to the resistant fish, susceptible fish had a significant parasite burden in their intestine, but no genes with known immune functions were upregulated, highlighting a lack of early recognition in these fish. The parasite continued to replicate exponentially in susceptible fish, and sequencing of intestinal mRNA from these fish at 14 and 21 dpe revealed an interferon-gamma driven T H 1 response and a suppression of the T H 17 response. This intense T H 1 response failed to ameliorate the progress of the disease and the intestines of these fish progressively deteriorated. A comparable dataset from resistant fish with a similar parasite burden in their intestine could help us discern why the susceptible fishes' immune response failed to offer any protection.
To understand how fish with a resistant phenotype respond to the parasite once it becomes established in the intestine, we challenged resistant steelhead trout with a dose of C. shasta sufficient to cause pathology. Accomplishing this task in the laboratory presented a significant hurdle given the difficulty of generating large quantities of actinospores (the infective parasite stage for the fish) and the extremely high resistance threshold of these fish. For this reason, we chose to focus on two key timepoints in the infection, representing the early and late immune response. We anticipated that these fish would respond with a stronger transcriptomic signal due to the higher parasite burden and that the early timepoint would confirm the results of our previous study: that resistant fish mount an earlier immune response to the parasite. Additionally, we believed it would provide a clear view of the early innate immune response, undistorted by the changes in gene expression that arises later in the infection due to the breakdown of the intestinal structure. The late timepoint would help characterize the adaptive immune response, as well as provide a dataset for comparison with susceptible fish from our previous study. We established severe infections by continuously exposing resistant steelhead trout to a high concentration of C. shasta for 5 days to achieve a high parasite burden in the intestine. Intestines were collected at 7 and 21 dpe and changes in gene expression were analyzed by RNA-seq, and the histopathological response was also investigated at both timepoints.

Infection of Resistant Steelhead Trout
The first clinical sign of C. shasta infection occurred at 19 dpe, when all the exposed fish stopped feeding, and thus anorexia became evident. At the 21 dpe sampling timepoint, the intestines of all the exposed fish were grossly swollen, with an apparent increase in capillarization and a hemorrhagic appearance. The tissue was rigid and difficult to cut and later homogenize in liquid nitrogen. Of the three exposed fish that remained after the 21 dpe sampling timepoint, one succumbed to the infection at 32 dpe. Its intestine appeared bloody and swollen and a vent swab revealed abundant mature myxospores. The remaining two fish resumed feeding again at 37 dpe and continued to do so until they were euthanized at 60 dpe. The parasite burden of the exposed fish was quantified by qPCR. The prevalence of infection was 100% among the exposed fish and at 7 dpe they had an average Cq of 27.0 ± 2.0, which increased to 15.2 ± 1.6 at 21 dpe ( Figure 1). For reference, the 1 and 1000 actinospore standards have Cq values of 34 and 24, respectively. All of the control fish were negative by qPCR. mount an earlier immune response to the parasite. Additionally, we believed it would provide a clear view of the early innate immune response, undistorted by the changes in gene expression that arises later in the infection due to the breakdown of the intestina structure. The late timepoint would help characterize the adaptive immune response, as well as provide a dataset for comparison with susceptible fish from our previous study We established severe infections by continuously exposing resistant steelhead trout to a high concentration of C. shasta for 5 days to achieve a high parasite burden in the intestine Intestines were collected at 7 and 21 dpe and changes in gene expression were analyzed by RNA-seq, and the histopathological response was also investigated at both timepoints

Infection of Resistant Steelhead Trout
The first clinical sign of C. shasta infection occurred at 19 dpe, when all the exposed fish stopped feeding, and thus anorexia became evident. At the 21 dpe sampling timepoint, the intestines of all the exposed fish were grossly swollen, with an apparen increase in capillarization and a hemorrhagic appearance. The tissue was rigid and diffi cult to cut and later homogenize in liquid nitrogen. Of the three exposed fish that re mained after the 21 dpe sampling timepoint, one succumbed to the infection at 32 dpe. Its intestine appeared bloody and swollen and a vent swab revealed abundant mature myx ospores. The remaining two fish resumed feeding again at 37 dpe and continued to do so until they were euthanized at 60 dpe. The parasite burden of the exposed fish was quan tified by qPCR. The prevalence of infection was 100% among the exposed fish and at 7 dpe they had an average Cq of 27.0 ± 2.0, which increased to 15.2 ± 1.6 at 21 dpe ( Figure  1). For reference, the 1 and 1000 actinospore standards have Cq values of 34 and 24, re spectively. All of the control fish were negative by qPCR.

Histopathology
Histological observation evidenced a substantial hyperplasia in the lamina propriasubmucosa of the infected fish at 21 dpe, which was not observed in control fish (Figure 2). At 7 dpe, intestines of infected fish maintained their basic tissue structure, though cell proliferation and infiltration in the submucosa were patent and localized detachment of the epithelium from the lamina propria occurred. These signs were significantly magnified in the intestines of infected fish at 21 dpe, which presented large areas of detached epithelia and a severe submucosal hyperplasia. At this timepoint, trophozoites as well as disporoblasts with mature spores were observed invading the submucosa, but they did not reach beyond the stratum compactum ( Figure 3). This collagen sheath, which lies beneath the lamina propria, seemed to act as a barrier, preventing parasite from spreading into deeper layers. Inflammatory cell foci in infected fish mainly consisted of lymphocyte-like cells in the submucosa and infiltrating the epithelium, and granulocytes. The latter are expected to be eosinophils, though eosinophilia was very weakly distinguished in tissue sections due to fixation with Dietrich's solution.
Interestingly, the fish that survived until 60 dpe were able to largely regenerate their intestinal structure. There was evidence of some localized epithelial detachment ( Figure 2d) and large clusters of granulocytes were present along the submucosa (Figure 3d). However, the parasite was not completely cleared by 60 dpe, and C. shasta stages were still found, mainly in the intestinal lumen ( Figure 2e).
A pattern of increasing goblet cell abundance with increasing time post-exposure was observed in both experimental groups. These goblet cells contained neutral mucins (PAS+) (Figure 3).

Histopathology
Histological observation evidenced a substantial hyperplasia in the lamina propriasubmucosa of the infected fish at 21 dpe, which was not observed in control fish ( Figure  2). At 7 dpe, intestines of infected fish maintained their basic tissue structure, though cell proliferation and infiltration in the submucosa were patent and localized detachment of the epithelium from the lamina propria occurred. These signs were significantly magnified in the intestines of infected fish at 21 dpe, which presented large areas of detached epithelia and a severe submucosal hyperplasia. At this timepoint, trophozoites as well as disporoblasts with mature spores were observed invading the submucosa, but they did not reach beyond the stratum compactum ( Figure 3). This collagen sheath, which lies beneath the lamina propria, seemed to act as a barrier, preventing parasite from spreading into deeper layers. Inflammatory cell foci in infected fish mainly consisted of lymphocyte-like cells in the submucosa and infiltrating the epithelium, and granulocytes. The latter are expected to be eosinophils, though eosinophilia was very weakly distinguished in tissue sections due to fixation with Dietrich's solution. Intestinal structure of unexposed control (a,b) and infected (c-e) resistant steelhead trout at 7 dpe to Ceratonova shasta (a,c), at 21 days (b,d) and 60 dpe (e). Note the mild (c) and severe (d) submucosal hyperplasia in recipient fish at the two early time points (double arrowheads). At 60 dpe, the intestinal tissue structure has recovered and resembles that of controls (e). Note the presence of parasite stages in the intestinal lumen (insert). Giemsa staining. Scale bars = 50 µm, except in insert = 5 µm. E; epithelium, SM; intestinal submucosa; SC; stratum compactum, SG; stratum granulosum, M; muscularis.
Interestingly, the fish that survived until 60 dpe were able to largely regenerate their intestinal structure. There was evidence of some localized epithelial detachment ( Figure  2d) and large clusters of granulocytes were present along the submucosa (Figure 3d). However, the parasite was not completely cleared by 60 dpe, and C. shasta stages were still found, mainly in the intestinal lumen (Figure 2e).
A pattern of increasing goblet cell abundance with increasing time post-exposure was observed in both experimental groups. These goblet cells contained neutral mucins   Intestinal histopathology of resistant steelhead trout exposed to Ceratonova shasta. At 21 days post exposure (dpe) (a,b,e) parasite stages (asterisks) present in the intestinal submucosa (SM) reach the stratum compactum (SC) but do not invade the underlying stratum granulosum (SG) and muscularis (M). Note that the lamina propria-submucosa is invaded by the parasite and has lost its tissue structure (e). Inflammatory cell proliferation and infiltration at 7 dpe consisting of lymphocyte-like cells (c). Cluster of granulocytes (black arrows) including one evident eosinophil (white arrow) in the submucosa at a 60 dpe (d). Mucins stained with periodic acid Schiff (PAS) (magenta) in 7 dpe (f), 21 dpe (g) and 60 dpe (h) intestines. Note the high abundance of PAS+ goblet cells at the two latter time points. Giemsa staining (a-e). Scale bars = 20 µm. E; epithelium. Intestinal histopathology of resistant steelhead trout exposed to Ceratonova shasta. At 21 days post exposure (dpe) (a,b,e) parasite stages (asterisks) present in the intestinal submucosa (SM) reach the stratum compactum (SC) but do not invade the underlying stratum granulosum (SG) and muscularis (M). Note that the lamina propria-submucosa is invaded by the parasite and has lost its tissue structure (e). Inflammatory cell proliferation and infiltration at 7 dpe consisting of lymphocyte-like cells (c). Cluster of granulocytes (black arrows) including one evident eosinophil (white arrow) in the submucosa at a 60 dpe (d). Mucins stained with periodic acid Schiff (PAS) (magenta) in 7 dpe (f), 21 dpe (g) and 60 dpe (h) intestines. Note the high abundance of PAS+ goblet cells at the two latter time points. Giemsa staining (a-e). Scale bars = 20 µm. E; epithelium.

Sequencing of Intestinal mRNA
The sequencing of intestinal mRNA from resistant fish at 7 and 21 dpe generated 519 million single-end reads. The percentage of reads that could be uniquely mapped to a single locus in the rainbow trout genome varied between groups, ranging from 40.4% to 75.4% (Table 1). Principal component analysis showed a high degree of uniformity between samples within a group and clearly separated infected fish from uninfected controls with the first principal component explaining 63% of the total variation ( Figure 4). Hierarchical clustering also clearly separated the infected samples from 21 dpe ( Figure 5) and all the samples showed a high degree of correlation (r > 0.96), suggesting no outliers were present.

Differential Gene Expression in the Intestine at 7 dpe and Analysis of Key Gene Expression
The sequencing reads from 7 dpe were mapped to 32,959 genes in the rainbow trout genome, and 417 differentially expressed genes (DEGs) were identified in infected fish relative to their time-matched controls, including 338 (81%) upregulated genes and 78 (19%) downregulated genes ( Figure 6 and Supplementary S1). Among the DEGs, 30 genes were strongly upregulated (Log2-FC > 4), including C-C motif chemokine 4-like (285.6 fold), two paralogs of C-C motif chemokine 13-like (46.7, 79.9 fold), three copies of C-C motif

Sequencing of Intestinal mRNA
The sequencing of intestinal mRNA from resistant fish at 7 and 21 dpe generated 519 million single-end reads. The percentage of reads that could be uniquely mapped to a single locus in the rainbow trout genome varied between groups, ranging from 40.4% to 75.4% (Table 1). Principal component analysis showed a high degree of uniformity between samples within a group and clearly separated infected fish from uninfected controls with the first principal component explaining 63% of the total variation ( Figure 4). Hierarchical clustering also clearly separated the infected samples from 21 dpe ( Figure 5) and all the samples showed a high degree of correlation (r > 0.96), suggesting no outliers were present. Table 1. Average number of sequencing reads generated per sample for each group and the corresponding number of reads that could be mapped to one locus in the rainbow trout genome.

Group
Average

Differential Gene Expression in the Intestine at 7 dpe and Analysis of Key Gene Expression
The sequencing reads from 7 dpe were mapped to 32,959 genes in the rainbow trout genome, and 417 differentially expressed genes (DEGs) were identified in infected fish relative to their time-matched controls, including 338 (81%) upregulated genes and 78 (19%) downregulated genes ( Figure 6 and Supplementary S1). Among the DEGs, 30 genes were strongly upregulated (Log2-FC > 4), including C-C motif chemokine 4-like (285.6 fold), two paralogs of C-C motif chemokine 13-like (46.7, 79.9 fold), three copies of C-C motif chemokine 19-like (  Based on the DEGs functional description and their associated GO terms, 114 DEGs known to have functions involved in immune responses were grouped into ten categories and presented in Table 2. Based on the DEGs functional description and their associated GO terms, 114 DEGs known to have functions involved in immune responses were grouped into ten categories and presented in Table 2.

GO Enrichment 7 dpe
GO enrichment analysis was conducted on the 338 upregulated genes to gain insight into the biological functions of these DEGs, as well as the involved molecular processes. This identified 115 significantly enriched GO terms among the upregulated genes. ClueGo analysis clustered these terms into networks revolving around GO terms for "defense response to other organism", "defense response", "innate immune response", and "interferon gamma-mediated signaling pathway" (Figure 7). No specific GO enrichment was found among the 78 downregulated genes.

Differential Gene Expression in the Intestine at 21 dpe and Analysis of Key Genes
The sequencing reads from 21 dpe mapped to 34,286 genes and 11,984 DEGs were identified in infected fish relative to their time-matched controls. Of these, 5801 (48%) were upregulated and 6183 (52%) were downregulated. In addition to the increased number of DEGs at 21 dpe, the magnitude of the differential expression also greatly increased, with 1,149 DEGs being strongly upregulated (Log2-FC > 4) (Figure 8 and Supplementary S1).
The most highly upregulated genes were those related to the immune system, cell communication, and tissue remodeling: interleukin-6-like (7118.6 fold), two paralogs of gap junction C × 32.  Enriched GO terms were grouped into functionally related nodes using ClueGO, a Cytoscope plugin. Nodes are colored and grouped according to related functions and labelled by the most significant term of the group. Node size corresponds to the FDR-adjusted p value of each GO term and is specific to each graph.

Differential Gene Expression in the Intestine at 21 dpe and Analysis of Key Genes
The sequencing reads from 21 dpe mapped to 34,286 genes and 11,984 DEGs were identified in infected fish relative to their time-matched controls. Of these, 5801 (48%) were upregulated and 6183 (52%) were downregulated. In addition to the increased number of DEGs at 21 dpe, the magnitude of the differential expression also greatly increased, with 1,149 DEGs being strongly upregulated (Log2-FC > 4) (Figure 8)

GO Enrichment 21 dpe
GO analysis of the 5801 upregulated genes revealed that they were enriched for 1135 GO terms, predominantly immune-related. ClueGO analysis clustered the biological process GO terms into five networks revolving around the terms "transmembrane receptor protein tyrosine kinase signaling pathway", "immune response-activating signal transduction", "cytokine secretion", "immune response-activating cell surface receptor signaling pathway", "T cell differentiation", regulation of leukocyte migration", and "positive regulation of actin filament polymerization" (Figure 9). 21, 10, x FOR PEER REVIEW 13 of 28 . Figure 9. Gene ontology (GO) enrichment of the genes upregulated in the intestine of resistant steelhead trout at 21 days post exposure to Ceratonova shasta. Enriched GO terms were grouped into functionally related nodes using ClueGO, a Cytoscope plugin. Nodes are colored and grouped according to related functions and labelled by the most significant term of the group. Node size corresponds to the FDR-adjusted p-value of each GO term and is specific to each graph.
The same analysis of the 6189 downregulated genes found 491 enriched GO terms which ClueGO analysis clustered into networks based on metabolic and energy producing GO terms, including "generation of precursor metabolites and energy", and "carboxylic acid catabolic process" (Figure 10). Enriched GO terms were grouped into functionally related nodes using ClueGO, a Cytoscope plugin. Nodes are colored and grouped according to related functions and labelled by the most significant term of the group. Node size corresponds to the FDR-adjusted p-value of each GO term and is specific to each graph.
The same analysis of the 6189 downregulated genes found 491 enriched GO terms which ClueGO analysis clustered into networks based on metabolic and energy producing GO terms, including "generation of precursor metabolites and energy", and "carboxylic acid catabolic process" (Figure 10). Enriched GO terms were grouped into functionally related nodes using ClueGO, a Cytoscope plugin. Nodes are colored and grouped according to related functions and labelled by the most significant term of the group. Node size corresponds to the FDR-adjusted p-value of each GO term and is specific to each graph.

Analysis of Key Immune Genes Differentially Expressed at 21 dpe
Based on the GO enrichment profiles and the most highly upregulated genes, a set of 89 key genes related to the innate and adaptive immune response that were differentially expressed in infected steelhead trout at 21 dpe were identified and are presented in Table  3.  Enriched GO terms were grouped into functionally related nodes using ClueGO, a Cytoscope plugin. Nodes are colored and grouped according to related functions and labelled by the most significant term of the group. Node size corresponds to the FDR-adjusted p-value of each GO term and is specific to each graph.

Analysis of Key Immune Genes Differentially Expressed at 21 dpe
Based on the GO enrichment profiles and the most highly upregulated genes, a set of 89 key genes related to the innate and adaptive immune response that were differentially expressed in infected steelhead trout at 21 dpe were identified and are presented in Table 3.

Comparison to the Transcriptome of Susceptible Fish at 21 dpe
To understand the difference between the response of resistant fish and their susceptible counterparts, we compared the differential gene expression results from 21 dpe in this study to those previously generated by our research group for susceptible fish with a similar parasite burden at 21 dpe [21]. The experimental conditions of this study were designed to align with those of our previous one, and while other factors that influence gene expression cannot be ruled out, we believe this comparison will provide insight into the broader details of how resistant fish respond compared to their susceptible counterparts. Comparison of the 11,984 DEGs identified in resistant fish at 21 dpe with the DEGs found in susceptible fish at 21 dpe identified 7820 genes differentially expressed in both phenotypes ( Figure 11). 28 of the 7820 were expressed in opposite directions in resistant fish compared to susceptible fish (supplemental). Among the 28 genes, 2 paralogs of "60 kDa heat shock protein, mitochondrial" were the only genes with known immune functions and were downregulated in resistant fish (−2.3 and −2.5-fold). The overall expression of the 7820 shared DEGs was highly positively correlated (r = 0.94, p < 2.2 × 10 −16 ), and the average difference in fold change was 1.76 ± 1.73. GO enrichment analysis of the shared DEGs revealed 453 enriched GO terms, primarily related to metabolic processes, innate immune response, and cell adhesion. Analysis of the genes that were differentially expressed only in resistant fish revealed a similar profile with 302 enriched GO terms. However, adaptive immune system processes were much more prominent, including many terms that were not present among the shared DEGs or those unique to susceptible fish. These included "adaptive immune response", "adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains", "lymphocyte mediated immunity", and "T cell mediated immunity." However, adaptive immune system processes were much more prominent, including many terms that were not present among the shared DEGs or those unique to susceptible fish. These included "adaptive immune response", "adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains", "lymphocyte mediated immunity", and "T cell mediated immunity." The number of secreted IgT and IgM transcripts was quantified for both resistant and susceptible fish at 21 dpe to compare the antibody response between the two phenotypes. This revealed a much stronger induction of both immunoglobulins in resistant fish, particularly IgT, which was upregulated over 300-fold (table 4).  The number of secreted IgT and IgM transcripts was quantified for both resistant and susceptible fish at 21 dpe to compare the antibody response between the two phenotypes. This revealed a much stronger induction of both immunoglobulins in resistant fish, particularly IgT, which was upregulated over 300-fold (Table 4).

Discussion
The dialogue between host and parasite is complex, involving not only the immune system, but also changes in the host's metabolism and the cellular structure of the infected tissues. RNA-seq offers a broad, non-targeted approach to understand this interaction, and can identify putative resistance genes for future functional studies. Ceratonova shasta offers an attractive biological model for studying the host-parasite dynamic, given that it infects well-studied and economically important salmonids that naturally occur on the phenotypic extremes of resistance. In this study, we sequenced intestinal RNA from resistant steelhead trout that were exposed to a high dose of C. shasta to elicit a vigorous immune response. Intestinal samples from 7 and 21 dpe were sequenced to further our understanding of the innate and adaptive immune response to this parasite. We found that by 7 dpe resistant fish were already mounting a vigorous IFN-gamma driven T H 1 response that was accompanied by remodeling of the intestinal tissues. This continued at 21 dpe, where signs of a strong antibody response were evident along with a possible transition to a T H 2 response. Despite the overwhelming parasite dose (>30,000 spores) these fish received and the damage to the intestinal structure that resulted, some fish were able to regenerate their intestinal tissues and resume feeding as normal.
In our previous study comparing the transcriptomic response of resistant and susceptible steelhead trout exposed to a moderate dose of C. shasta, we observed that despite having a significant parasite burden in their intestine by 7 dpe, susceptible fish did not respond to the infection [21]. Although the exposure was not sufficient to establish an active infection in the intestines of the resistant fish in that study, we observed an upregulation of several immune genes in the intestine at 7 dpe, suggesting that a more rapid immune response was occurring in resistant fish. The results presented here confirm this, with resistant fish mounting a vigorous immune response at 7 dpe, with 338 genes upregulated. This response was predominantly an IFN-gamma driven T H 1 response, with upregulation of numerous interferon-stimulated genes, as well as MHC class I genes and those involved in antigen processing and presentation. This would suggest that an adaptive immune response is already being mounted against the parasite. We also observed upregulation of NLRC5 and GTPase IMAP family member 4 (GIMAP 4). NLRC5 is a cytosolic pathogen recognition receptor involved in MHC class I-dependent immune responses [22] and was found to be upregulated in the gills and intestine of resistant fish in our previous study. In that same study, two paralogs of GIMAP 4, a protein involved in T-lymphocyte development [23], were the most highly upregulated genes in the gills of resistant fish at 1 dpe. Intriguingly, the exact same paralogs were the most highly upregulated immune genes in the intestine of susceptible fish at 14 dpe. The early induction of these genes in resistant fish compared to susceptible fish strongly suggests a role for these genes in C. shasta resistance. Although different paralogs of NLRC5 and GIMAP 4 were identified in this study, our findings further support their involvement in C. shasta resistance.
It is also important to highlight that despite being exposed to much higher concentration of C. shasta actinospores, and over a longer period of time, 2-fold less parasite DNA could be detected in the intestine of resistant fish at 7 dpe (Cq of 27.0 ± 2.0 for resistant fish vs. a Cq of 24.8 ± 0.8 for susceptible fish). This further supports the hypothesis that resistant fish can recognize the parasite and elicit an immune response early in the infection that is capable of either directly killing the parasite or inhibiting its replication. It remains to be determined precisely when and where this occurs, but it may be occurring in the blood vessels during C. shasta migration from the gills to the intestine. When resistant and susceptible fish have been exposed in parallel to C. shasta, no difference in parasite burden has been detected in the gills, and both phenotypes had a suppressed IFN-gamma response in this tissue [21,24]. Elimination of the parasite in the intestine does not seem to occur until much later in the infection (25+ dpe), and between 7 and 21 dpe the amount of parasite DNA in the intestine seems to follow an exponential growth curve in both resistant and susceptible fish. Our observations suggest that initial resistance either occurs in the blood or when the parasite first reaches the intestine.
While it was not surprising to find a strong induction of IFN-gamma, something that is commonly seen in C. shasta infected salmonids and other myxozoan infected fish [24][25][26][27], it was perhaps surprising that the early immune response was so strongly focused on this one pathway. IFN-gamma is the signature T H 1 cytokine and mediates the response to viruses and other intracellular pathogens. C. shasta is considered an extracellular pathogen and certainly occupies that role in the intestine. This raises the question as to how or why a pathway geared towards intracellular pathogens is invoked and why a cytosolic pathogen recognition receptor (NLRC5) would be upregulated. This may be explained by C. shasta having an intracellular phase early in the infection. The only evidence for this comes from a study by , which reported the invasion and migration of fluorescently stained C. shasta [28]. The authors noted a potential intracellular phase in the endothelium of blood vessels starting at 3 dpe and suggested the use of electron microscopy to resolve the question of an intracellular developmental stage of the parasite. This would not be without precedent, as other myxozoans are known to have an early intracellular developmental stage or even having an intracellular final location. M. cerebralis initially replicates intracellularly in the epithelia of the epidermis of trout [29], and Sphaerospora molnari is intracellular in the gills of carp, prior to extracellular proliferation in the blood [30], and Kudoa species develop within muscle fibers. In the latter case, it has been suggested that the parasite remains undetected by the host immune system due to their intracellular location, firmly enveloped by the remnants of the muscle fiber and the surrounding connective network of the endomysium [31].
Whether an early intracellular phase is required for parasite development or represents a form of immune evasion remains to be determined. An intriguing possibility, at least for C. shasta, is that going intracellular causes the host to initiate a cytotoxic T H 1 response that is ineffective against the extracellular stages in the intestine. This is evidenced by the upregulation of antigen presentation machinery and MHC class I molecules, which are responsible for the presentation of cytosol-derived peptides to CD8 + cytotoxic T lymphocytes (CTLs). Activated CTLs recognize and directly kill infected host cells, which is vital for controlling infections caused by intracellular pathogens [32]. A mistargeted immune response would explain why the rapid induction of IFN-gamma in the intestine at 7 dpe failed to slow parasite proliferation, with significantly more parasite DNA being detected at 21 dpe.
If IFN-gamma is playing a role in resistance to C. shasta in the intestine, the actual mechanism remains unclear. IFN-gamma mediates parasite clearance by inducing the production of nitric oxide (NO) in classically activated macrophages [33,34]. This antiparasitic effect has been demonstrated for numerous parasites, including the well-studied intestinal parasites Entamoeba histolytica, Cryptosporidium parvum, and Giardia spp. [34][35][36]. In the present study, we did not observe upregulation of inducible nitric oxide synthase, which is responsible for NO production in macrophages, and instead observed upregulation of arginase at both 7 and 21 dpe. Arginase, a marker for alternatively activated macrophages, utilizes the same substrate as nitric oxide synthase (L-arginine) and redirects it to the formation of polyamines and proline, which are necessary for collagen synthesis and wound healing [37]. In addition to depleting the necessary substrate for NO production, upregulation of arginase also inhibits the expression of inducible nitric oxide synthase. It is also curious that we observed upregulation of arginase at 7 dpe in the absence of any T H 2 cytokines, which typically drive its expression. The protozoan parasites Trypanosoma brucei and T. cruzi have both been shown to induce host expression of arginase in a T H 2-cytokine independent manner as a means of suppressing NO production [38,39]. Our observations may indicate either a protective host response geared towards controlling inflammation and maintaining tissue integrity, or a pathogenic strategy to escape IFN-γ driven NO production. Given how important the interplay between classically activated NO producing macrophages and alternatively activated arginase expressing macrophages are to the resolution of parasitic infections, an in-depth analysis of the macrophage populations during C. shasta infection would greatly benefit our understanding of the host immune response.
In addition to developing an earlier immune response, resistant fish have a much different reaction to the parasite at the tissue level. At 7 dpe, evidence of tissue remodeling is present with the upregulation of junctional proteins, arginase, and the massive upregulation of caspase-14, which plays a terminal role in keratinocyte differentiation [40]. Keratins are cytoskeletal proteins classically used as epithelial cell markers during injury and disease in vertebrates [41,42] and which have also been described in non-cornified mucosal epithelia of teleosts [43]. In mice, keratinocyte growth factor has been shown to ameliorate drug-induced intestinal mucosal disruption through induction of epithelial repair and tight junction protein expression, together with the inhibition of increased epithelial permeability [44]. The early changes in gene expression observed in this study, including the upregulation of caspase-14, point towards an active repair process of the intestinal barrier disrupted by the parasite. The gross pathology observed at 21 dpe, where the intestine took on a rigid, leathery appearance and was mechanically resistant to homogenization in liquid nitrogen might be related to the thickening of the lamina propria-submucosa due to the observed hyperplasia. The increased vascularization of the organ that we observed would contribute to the tissue repair process. By 21 dpe, no parasite stages were observed beyond the stratum compactum, which appeared to act as a parasite barrier. Thus, parasite proliferation is limited in resistant steelhead trout to the mucosal layers of the intestine, where cell regeneration occurs on a daily basis favoring lesion recovery. This is likely aided by the massive upregulation of IL-10 at 21 dpe, which would facilitate the resolution of inflammation and subsequent tissue repair. This is similar to what is observed in gilthead sea bream that become resistant to reinfection with E. leei, where IL-10 is associated with the resolution of infection [45]. The tissue response of resistant steelhead trout in this study strongly contrasts with what is observed in susceptible fish, where all layers of the intestine become infected and it becomes a soft, spongy mass that loses its overall structure [28,46,47]. How the stratum compactum becomes a physical barrier and whether it contributes to the intestinal hardening, deserves further studies. It does appear that the ability of resistant fish to maintain their intestinal structure in the face of parasite replication is likely a critical factor in resisting C. shasta induced mortality and would explain why previous studies have observed organized, tissue-level responses to C. shasta in resistant fish, but not in susceptible fish [48][49][50][51].
Differences in the intestinal epithelial integrity between hosts have also been observed in studies of the intestinal myxozoans Enteromyxum scophthalmi and E. leei [52]. Turbot is highly susceptible to E. scophthalmi, suffering serious intestinal lesions and barrier dysfunction, leading to high morbidity and mortality rates. Gilthead sea bream (Sparus aurata L.), on the other hand, experience low mortality rates and are able to better maintain their intestinal epithelial integrity even when heavily infected by E. leei [53]. Similar to C. shasta, E. leei is able to infect a wide range of fish species with varying degrees of host susceptibility [4]. For sharpsnout seabream (Diplodus puntazzo), differences in humoral immune factors have been suggested to play a role in their high susceptibility to E. leei [54], however no transcriptomic comparison with gilthead sea bream has been conducted. It should also be noted that gilthead sea bream that survive primary exposure to E. leei acquire a protective immunity to reinfection, that is associated with increased IgM and IgT expression relative to naïve fish [45]. Similar to gilthead sea bream, the resistant fish in this study were better able to maintain their intestinal integrity and had higher expression of IgM and IgT compared to susceptible fish at this same point in the infection.
Given how different the intestinal response of resistant fish is at 7 dpe compared to susceptible fish, it is surprising how similar the response is at 21 dpe. 7820 genes were differentially expressed in both phenotypes and their expression levels were highly correlated. In addition to numerous metabolic and cell junction genes, this shared response includes several immune factors (IL-6, IL-8, IL-10, IFN-γ) that have been found to be upregulated in other studies of C. shasta infected salmonids [24,25,55]. While we cannot rule out temporal differences in the expression of these genes being a critical factor in resistance, something that almost certainly occurs given the delayed parasite recognition of susceptible fish, their differential expression alone does not explain resistance vs. susceptibility to C. shasta. Examination of the genes that are only differentially expressed in resistant fish at 21 dpe revealed a much larger role for the adaptive immune response in these fish, particularly the B cell response. Resistant fish have significantly more heavy and light chain transcripts upregulated, as well as secreted IgM and IgT transcripts at this time. Again, this is likely influenced by the delayed parasite recognition of susceptible fish, but it suggests an earlier and stronger antibody response at play in resistant fish. It has been demonstrated that salmonids are capable of generating IgM and IgT that is specific to C. shasta [56]. However, the effectiveness of this antibody response in reducing mortality from C. shasta remains unclear. A study of susceptible rainbow trout exposed to C. shasta showed upregulation of IgM and IgT in these fish, but it failed to improve their condition and 100% mortality occurred [55]. The authors noted that this may be due to the response coming too late, after the intestine has been severely damaged. Our observations would support this, as resistant fish were better able to maintain their intestinal structure and mounted an earlier adaptive immune response, including the observed foci of lymphocytelike cells.
Differences in the T cell response are also evident at 21 dpe, most notably among T H 17 cytokines which are an important aspect of the gut mucosal barrier that help prevent the dissemination of bacteria [57]. In our previous study, we observed a strong downregulation of IL-17 family cytokines in susceptible fish at 21 dpe, whereas in the resistant fish in this study, no IL-17 cytokines were downregulated and IL-17A was upregulated. Whether this is directly related to the C. shasta is unclear and it may be more related to overall gut health and the invasion of opportunistic bacteria. More interesting is the possible transition to a T H 2 response that may be occurring late in the infection. This adaptive immune response is driven by IL4/13, which we observed massive upregulation of at 21 dpe, along with upregulation of the transcriptional regulators GATA3 and STAT5. The T H 2 response is associated with wound healing, the suppression of T H 1-driven inflammation, and the elimination of helminth parasites [58,59]. Although much smaller than helminth worms, C. shasta is also an extracellular intestinal parasite and the mechanisms used to clear helminths (B cell activation, eosinophil recruitment, mucus hypersecretion, increased cellular turnover) may help eliminate or otherwise create an unfavorable environment for C. shasta.
Comparison of our finding with other studies on myxozoan infections suggest that an effective T cell response is critical for reducing pathology. Similar to what we observed in C. shasta infected rainbow trout, when strains of rainbow trout that are resistant and susceptible to whirling disease are compared, the resistant strain has an earlier and more effective immune response, with upregulation of IFN-gamma sooner in the infection and a stronger T cell response during the initial stages [26,60]. The failure of this early immune response in susceptible fish leads to a sustained inflammatory response that fails to control the infection and likely contributes to tissue damage. The pathology of proliferative kidney disease in rainbow trout is associated with an early imbalance of T H cytokines and a dysregulated B cell response [9]. An in-depth analysis of T H cytokines in gilthead sea bream during E. leei infection revealed a strong induction of T H 1 and T H 17 cytokines in the intestine [61]. The ability of Atlantic salmon (Salmo salar) to resolve infection with K. thyrsites, and develop resistance to reinfection, is linked to a cytotoxic T cell response [7]. Further examination of the T cell response during C. shasta infection, including determination of the cellular kinetics and the specific T cell subsets responding to the infection, would help elucidate the cellular basis to resistance.

Conclusions
In this study we employed RNA-seq to explore the complex dynamic between host and parasite and to understand how resistant steelhead trout are able to overcome an actively progressing C. shasta infection in their intestine. Our results confirmed that early parasite recognition is critical for resistance and that initial invasion of the intestine by the parasite elicits a strong IFN-gamma driven adaptive immune response in resistant fish. This is accompanied by remodeling of the host intestinal tissue, which is likely vital for maintaining the overall intestinal structure in the face of extensive parasite replication and the tissue inflammation that results. This early adaptive immune response leads to a vigorous B cell response later in the infection, characterized by strong antibody production, particularly of IgT. Based on these results, we propose that a core immune response to C. shasta exists among resistant and susceptible fish and that temporal differences in the expression of key immune factors (IFN-γ, IL-6, GIMAP4, IgT), resulting from earlier parasite recognition by resistant fish, largely explains the different infection outcomes for these fish. In conjunction with this, resistant fish have a different response to the parasite at the tissue level with the stratum compactum playing an important role in limiting parasite spreading, which might be related to the induction of keratinization. Given that most RNA-seq studies have been conducted on fish that are susceptible to myxozoans, we believe the present study offers a valuable framework for putting those, and future studies, in perspective.

Experimental Fish
Resistant steelhead trout were collected from the Round Butte Hatchery (OR, USA) and bred as previously described [62]. The fish were reared on 13.5 • C specific-pathogen free (SPF) well water at the Oregon State University (OSU) John L. Fryer Aquatic Animal Health Laboratory (AAHL) in Corvallis, OR, USA, and fed a commercial diet daily (BioClark's Starter, Bio-Oregon, Longview, WA, USA).

Parasite Challenge
Laboratory cultures of the invertebrate host Manayunkia occidentalis are maintained at the AAHL and serve as source of C. shasta actinospores for laboratory challenges. The annelids are housed in indoor mesocosms receiving flow-through UV-treated treated river water and spore production is routinely monitored [63]. The parasite challenge was initiated in March 2019, when spore production from the genotype IIC mesocosm reached 30,000 spores per day, a sufficient dose to meet or exceed their threshold of resistance [11,20]. 21 resistant steelhead trout (average 103.1 g ± 7.1 g) were placed in a 100-liter tank that received effluent from the mesocosm. Water temperature was increased from 14.5°C to 17.5°C ± 0.4°C over 2 days by addition of 18°C SPF well-water. This temperature was chosen as it is representative of the river water temperatures that out-migrating salmonids experience when they are exposed to C. shasta, and it aligns with our previous study of the transcriptomic response of steelhead trout. An additional 21 Round Butte steelhead trout were transferred into an identical tank setup but received effluent from the control mesocosm which contains uninfected M. occidentalis. 5 days after the exposure began, the fish (treatment and control) were transferred into six 25-liter tanks (7 fish per tank) that were randomly assigned and supplied with 18°C SPF well-water.

Tissue Sampling
Fish were sampled at 7-and 21-days post exposure (dpe), with exposure being defined as their initial placement in the exposure tank. At each timepoint, 3 fish from each of the 6 tanks were euthanized with an overdose of MS-222 (tricaine methanosulfonate, Argent Laboratories, Redmond, WA, USA) for a total of 18 fish per timepoint (9 exposed, 9 controls). The entire intestine was removed from each fish and placed in either RNAlater (2 out 3 fish per tank) or Dietrich's fixative (1 out of 3 fish per tank). The samples collected in RNAlater were immediately stored at 4°C and then transferred to -80°C after 24 hours. At each timepoint, the fish were sampled at the same time of day to avoid changes in gene expression due to circadian rhythms [64]. After the fish were sampled at the 21 dpe timepoint, the remaining fish in each tank were consolidated into two tanks (exposed and control) to eliminate distress due to isolation. The fish were monitored until 60 dpe, at which time they were euthanized with an overdose of MS-222.

Sample Processing
The intestine samples collected in RNAlater were homogenized in liquid nitrogen using a porcelain mortar and pestle. Following this, 25 mg of homogenized tissue from each sample underwent RNA extraction using the RNeasy Mini Kit (Qiagen, catalog number 74104) according to the manufacture's protocol. An additional 25 mg of tissue underwent DNA extraction using the DNeasy Blood and Tissue Kit (Qiagen, catalog number 69506). The extracted DNA was eluted in 30 µl of Buffer AE that was applied to the spin column twice to achieve a higher concentration. The concentration and purity and of the extracted RNA and DNA was assessed with a NanoDrop ND-1000 UV-Vis Spectrophotometer.
The amount of parasite DNA present in the intestine was determined using C. shastaspecific qPCR assay [63] and 100 ng of extracted DNA from each sample was assayed in triplicate wells through 40 cycles using an Applied Biosystems StepOnePlus Real-Time PCR System. A sample was considered positive for C. shasta if all three wells fluoresced, and the sample was rerun if the Cq standard deviation between wells was greater than 1. On each qPCR plate, a positive control, a negative control (molecular grade water), and a standard dilution curve equivalent to 1, 10, 100, and 1000 actinospores was included.
Intestinal samples in Dietrich's fixative were routinely processed for histology, embedded in Technovit 7100 resin (Heraeus Kulzer, Wehrheim, Germany) and 3 µm sections were stained with Giemsa or with periodic acid Schiff (PAS). Observations and microphotographs were made with a Leitz Dialux 22 light microscope connected to an Olympus DP70 camera.

Library Prep and Sequencing
Intestinal mRNA from 8 samples at each timepoint (4 treatment, 4 control) were submitted to the Center for Genome Research and Biocomputing at OSU for library preparation and sequencing. The integrity of the RNA was confirmed by running each sample on an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). A totol of 1 µg of RNA was used for library preparation using the Illumina TruSeq™ Stranded mRNA LT Sample PrepKit according to the manufacturer's instructions (Cat. No. RS-122-2101, Illumina Inc. San Diego, CA, USA). Library quality was checked with a 4200 TapeStation System (Agilent Technologies, USA) and quantified via qPCR. The libraries were sequenced on two lanes of an Illumina HiSeq 3000 as 100-bp single-end runs.
Differentially expressed genes (DEGs) between treatment and control fish were identified using the negative binomial Wald test in DESeq2 and were considered significant if they had a Benjamini-Hochberg false discovery rate (FDR) adjusted p value < 0.05 and an absolute log 2 (fold change) > 1. Annotation of the DEGs and gene ontology (GO) enrichment was conducted with OmicsBox (v 1.1.135, https://www.biobam.com/omicsbox/, accessed on 12 August 2020) [71]. To obtain high quality annotations, the blast e-value cutoff was set at 1 × 10 −5 and genes were preferentially annotated with the SWISS-PROT database [72] followed by the NCBI nonredundant database and the 'Vertebrata' taxonomy filter was applied. All genes detected in this study were used as the background for GO enrichment and up-and downregulated genes were analyzed separately. Enriched GO terms and their FDR-adjusted p-values were imported into Cytoscope (v 3.7.2, https://cytoscape.org/, accessed on 20 August 2020) [73] for visualization with the ClueGo (v 2.5.6, http://www.ici.upmc.fr/cluego/, accessed on 20 August 2020) [74] plugin, which clusters the GO terms into functionally related networks. O. mykiss was chosen as the organism for Ontologies/Pathways and the GO Term Fusion option was used to merge GO terms based on similar associated genes. Volcano plots were constructed with the R package EnhancedVolcano (v 1.0.1, https://bioconductor.org/packages/release/bioc/ html/EnhancedVolcano.html, accessed on 12 September 2020) [75].
DEGs found in this study were compared to DEGs previously identified in susceptible steelhead trout exposed to C. shasta [62]. Genes that were differentially expressed in both phenotypes at 21 dpe were identified and the Pearson correlation coefficient between their expression levels was calculated in R. To analyze the antibody response against C. shasta, sequenced reads from both resistant and susceptible at 21 dpe were mapped against the coding sequence for the secreted forms of IgM (GenBank: S63348.1) and IgT (GenBank: AY870263.1) using Salmon (v 0.10.0, https://salmon.readthedocs.io/en/latest/, accessed on 21 September 2020) [76]. The output was imported into DESeq2 to estimate fold-change between exposed and control fish.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/pathogens10091179/s1, Table S1: All differentially expressed genes (DEGs) identified in this study and the associated gene ontology (GO) enrichment terms.  Institutional Review Board Statement: The Institutional Animal Care and Use Committee of Oregon State University reviewed and approved all husbandry practices used in this study (IACUC protocol #5010).

Informed Consent Statement: Not applicable.
Data Availability Statement: All sequencing data generated in this study will be made publicly available on the NCBI's short read archive (SRA) by the time of publication. All other relevant data are contained within the article.