Front. Immunol. Frontiers in Immunology Front. Immunol. 1664-3224 Frontiers Media S.A. 10.3389/fimmu.2021.677730 Immunology Original Research Dual-Organ Transcriptomic Analysis of Rainbow Trout Infected With Ichthyophthirius multifiliis Through Co-Expression and Machine Learning Roh HyeongJin 1 † Kim Nameun 1 Lee Yoonhang 1 Park Jiyeon 1 Kim Bo Seong 2 Lee Mu Kun 3 Park Chan-Il 4 Kim Do-Hyung 1 * 1 Department of Aquatic Life Medicine, College of Fisheries Science, Pukyong National University, Busan, South Korea 2 Aquatic Disease Control Division, National Institute of Fisheries Science (NIFS), Busan, South Korea 3 Korean Aquatic Organism Disease Inspector Association, Busan, South Korea 4 Department of Marine Biology & Aquaculture, College of Marine Science, Gyeongsang National University, Tongyeong, South Korea

Edited by: Geert Wiegertjes, Wageningen University and Research, Netherlands

Reviewed by: Jing Xing, Ocean University of China, China; Javier Santander, Memorial University of Newfoundland, Canada; Louise Jørgensen, University of Copenhagen, Denmark

*Correspondence: Do-Hyung Kim, dhkim@pknu.ac.kr

†Present address: HyeongJin Roh, Pathogens and Disease Transfer, Institute of Marine Research (IMR), Bergen, Norway

This article was submitted to Comparative Immunology, a section of the journal Frontiers in Immunology

 08 07 2021 2021 12 677730 08 03 2021 31 05 2021 Copyright © 2021 Roh, Kim, Lee, Park, Kim, Lee, Park and Kim 2021 Roh, Kim, Lee, Park, Kim, Lee, Park and Kim

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Ichthyophthirius multifiliis is a major pathogen that causes a high mortality rate in trout farms. However, systemic responses to the pathogen and its interactions with multiple organs during the course of infection have not been well described. In this study, dual-organ transcriptomic responses in the liver and head kidney and hemato-serological indexes were profiled under I. multifiliis infection and recovery to investigate systemic immuno-physiological characteristics. Several strategies for massive transcriptomic interpretation, such as differentially expressed genes (DEGs), Poisson linear discriminant (PLDA), and weighted gene co-expression network analysis (WGCNA) models were used to investigate the featured genes/pathways while minimizing the disadvantages of individual methods. During the course of infection, 6,097 and 2,931 DEGs were identified in the head kidney and liver, respectively. Markers of protein processing in the endoplasmic reticulum, oxidative phosphorylation, and the proteasome were highly expressed. Likewise, simultaneous ferroptosis and cellular reconstruction was observed, which is strongly linked to multiple organ dysfunction. In contrast, pathways relevant to cellular replication were up-regulated in only the head kidney, while endocytosis- and phagosome-related pathways were notably expressed in the liver. Moreover, interestingly, most immune-relevant pathways (e.g., leukocyte trans-endothelial migration, Fc gamma R-mediated phagocytosis) were highly activated in the liver, but the same pathways in the head kidney were down-regulated. These conflicting results from different organs suggest that interpretation of co-expression among organs is crucial for profiling of systemic responses during infection. The dual-organ transcriptomics approaches presented in this study will greatly contribute to our understanding of multi-organ interactions under I. multifiliis infection from a broader perspective.

 dual-organ RNA-seq Ichthyophthirius multifiliis white spot disease weighted gene co-expression network analysis machine learning Introduction

White spot disease caused by Ichthyophthirius multifiliis is a well-known condition that causes tremendous economic losses of trout farms but also other species of production fish (1–4). I. multifiliis can invade host epithelial layers, settle above the basal lamina, and cause tremendous damage to trout farms (5–7). Many studies have reported significant changes in innate (e.g., C3, macrophage activation with cytokine expression, serum amyloid A; SAA) and adaptive immune responses (e.g., IgT, IgM, and T and B lymphocytes) subsequent to I. multifiliis infection, but these findings are just the tip of the iceberg of immuno-physiological changes (8–15). In recent years, due to the development of advanced sequencing technologies, global transcriptomic analysis has become a robust tool for the study of immuno-physiological characteristics (16, 17). In the context of these trends, Syahputra et al. (4) monitored whole transcriptomic responses in the trout gill after I. multifiliis infection in an attempt to identify immunological mechanisms of response to infection. It is evident that the skin and gills are the primary I. multifiliis infection sites, but internal systemic responses in other immunologically significant organs (e.g., head kidney and liver) are also important (10, 11, 18). For instance, some studies (12, 19) found that several immune-related genes including complement factor, immunoglobulin, cytokine, and acute phase proteins were significantly differed depending on the type of organs (e.g., liver, head-kidney, and spleen). These results emphasize the importance of understanding global multi-organ transcriptomics, but there is currently not enough data on multi-organ transcriptomics responses under I. multifiliis infection in the rainbow trout.

Many studies have focused on methods of interpreting RNA-seq results, such as differentially expressed gene (DEG) analysis, which is able to identify significant alterations in gene expression (20). However, this method is of limited applicability to multi-organ RNA-seq results because DEG analysis is restricted to comparisons between only two groups. Also, an enormous number of DEGs can be retrieved from a dataset based on a given threshold value; accordantly, interpreting the importance of featured genes and their biological relevance can be very difficult (21–23). To overcome such problems, several statistical algorithms based on machine learning and weighted gene co-expression network analysis (WGCNA) are attracting a great deal of attention for their utility in interpreting massive amounts of transcriptomics data. These methods have advantages in efficiently analyzing multi-dimensional and/or multi-organ RNA-seq data and providing novel insights into systemic transcriptomic changes (23–26). However, a multi-organ transcriptomic atlas of any given hematopoietic organ has not yet been established. Here, we use dual-organ RNA sequencing (RNA-seq) approaches to characterize the despotic temporal transcriptome in detail during hematopoietic stem and progenitor cell expansion in the rainbow trout. Tools based on machine learning technology and WGCNA can make high-dimensional data analysis simple and help to identify featured genes, hub genes, and/or key factors (23, 27). If multi-dimensional transcriptomics results are interpreted through several well-known methods, the advantages of all methods can be maximized by compensating for the limitations of each. For instance, Deng et al. (28) combined the DEG and WGCNA methods to make up for the drawbacks of DEGs, which greatly contributed the understanding of the transcriptomic expression in association with the different pathological stages of renal cell carcinoma. Hence, WGCNA would help in understanding the transcriptomic changes in accordance with the I. multifiliis infection stages.

Therefore, the purposes of this study were to: 1) catalogue organ-enriched, organ-specific, and common gene expression profiles in the head kidney and liver, 2) investigate the interactions and relationship between these two organs, and 3) ultimately understand the systemic immuno-physiological changes in rainbow trout under I. multifiliis infection by means of serological and multi-organ transcriptomics data.

 Materials and Methods Fish and Samples

Severe Ichthyophthirius multifiliis infection occurred in a rainbow trout farm in Gangwon province, Korea with tens of thousands of juvenile rainbow trout (average = 61.7 ± 27.4 g). At the time, the water temperature was approximately 15°C. One week after the first death occurred, a 10% cumulative mortality rate was observed in cases of severe I. multifiliis infection in the gills and skin. At this time, 15 trout were randomly anesthetized with an excess of anesthetic agent (MS-222, Sigma) (Ich group). The trophont number was counted in the 1/2 gill arch to determine infection intensity, and all trout in the Ich group were found to be infected by I. multifiliis. Two weeks post-sampling, eight surviving trout with confirmed primary infection and recovery were additionally sampled as the recovery group (R-Ich group). Blood samples were collected from the caudal vein. Likewise, six trout of similar size without a history of I. multifiliis infection were sampled as the control (Con) group. Hematocrit (Ht) and serological factors (GOT, GPT, ALP, GLU, TCHO, TP, BUN, LDH, and Ca) were measured using a dry chemistry analyzer (FUJI DRI-CHEM 3000). Serum electrolyte levels (Na+, K+, and Cl−) were additionally checked using Fuji Dri-chem Na-K-Cl film; the reference fluid contained 100 mmol of sodium chloride (NaCl), 30 mmol of sodium bicarbonate (NaHCO3), 4 mmol of potassium dihydrogenphosphate (KH2PO4), and 0.1% NaN3. Based on the serological results, the osmolality of the serum was calculated using the following formula (29). The head kidney and liver were then collected and stored in RNAlater solution at −80°C.

Calculated osmolality L − 1 = 2 N a + ( mmol   L − 1 ) + ( glucose   ( mg   d l − 1 ) / 18 ) + ( Urea ( mmol   L − 1 ) / 2.8 ) RNA Extraction and RNA-Seq Analysis

Three rainbow trout in each group were randomly selected for extraction of RNAs from the head kidney and liver using TRIzol reagent (Invitrogen, USA) following the manufacturer’s guideline. The concentration of RNA was quantified by Agilent 2100 BioAnalyzer and Nanodrop, and the sample with higher than 6 RNA integrity number RIN was used in this study. Libraries were prepared for 100 bp paired-end sequencing using a TruSeq Stranded mRNA Sample Preparation Kit (Illumina, CA, USA). Briefly, mRNAs were purified and fragmented from 1 μg of total RNA using oligo dT magnetic beads. These fragmented mRNAs were synthesized as single-stranded RNAs. Double-stranded cDNAs were then prepared. After the sequential end repair process, A-tailing, and adaptor ligation, the cDNA libraries were evaluated with an Agilent 2100 BioAnalyzer (Agilent, CA, USA). cDNAs were quantified with a KAPA library quantification kit (Kapa Biosystems, MA, USA) according to the manufacturer’s library quantification protocol. Following cluster amplification of denatured templates, paired-end (2 × 100 bp) sequencing was performed using an Illumina HiSeq 2500 sequencer (Illumina, CA, USA).

 Mapping, Assembly, and Selection of Differentially Expressed Genes

Low-quality reads were filtered according to the following criteria: reads containing more than 10% skipped bases (marked as “N”s), reads containing more than 40% bases with quality scores of less than 20, and reads with average quality scores of less than 20. The whole filtering process was performed using in-house scripts. Filtered reads were mapped onto a trout reference genome (30) using an aligner TopHat (31). Gene expression levels were measured with Cufflinks v2.1.1 (32) using the SwissProt gene annotation database. To improve measurement accuracy, the multi-read-correction and frag-bias-correct options were applied. All other options were set to default values. Differential expression analysis was performed using Cuffdiff (33).

 Pathway Analysis Using DEGs

Selected DEGs were functionally annotated using the KEGG database (34). A Z-score based on the number of increases or decreases in ko_id for each pathway was calculated using the following formula.

Z − score = Number of up regulated   k o id −   Number of down regulated   k o id Number of ko _ id

Based on the expression tendency of Z-score in the head kidney and liver, the pathways that were up- and down-regulated under I. multifiliis infection were classified. The pathways with the same number of up- and down-regulated DEGs as the Con and Ich groups were considered the etc. group. To investigate co-expression between organs, the pathways that showed up-regulation in both organs (PO-PO), up-regulation in the head kidney but down-regulation in the liver (PO-NE), down-regulation in the head kidney but up-regulation in the liver (NE-PO), and down-regulation in both organs (NE-NE) during I. multifiliis infection or recovery were categorized and visualized using Cytoscape (Ver. 3.6.1) (22, 35).

 Profiling of Featured Genes in Both Organs Using a Machine Learning Interface

To excavate the featured genes that differed between the control and Ich groups in the head kidney and liver RNA-seq analysis, Poisson linear discriminant analysis (PLDA) based on a machine learning method was implemented regardless of the organ. Briefly, twelve RNA-seq results from the control and Ich groups from both head kidney and liver were randomly split into eight samples (four samples from the control group and four from the Ich group) for the training set, and the remaining samples were used as the test set. The FPKM count in all transcripts was normalized by the “deseq median ratio normalization method”, calculated by dividing each FPKM sample by the geometric mean of the FPKM (36). The transformed PLDA model was generated using the training set, and featured genes were selected that effectively discriminated between the control and Ich groups. All statistical calculations and machine learning algorithms were performed using the MLSeq and DESeq2 packages in R (Ver. 3.6.2.) (23). Featured genes were annotated by Blast2GO software based on the GO database, KEGG, and SwissProt databases (34, 37). Featured genes were clustered based on GO and KEGG pathway, and visualized by Cytoscape (Ver 3.6.1.) (35).

 Weighted Gene Correlation Network Analysis (WGCNA)

The co-expression analysis of the FPKM counts of all 46,585 genes in 18 samples and the degree of infection (number of I. multifiliis microbes in the gill) were investigated using the WGCNA package in R (38). The similarity of the signed network was calculated using the formula below, as suggested by Langfelder and Horvath (38) and Kim et al. (26). In brief, blockwise modules were created under the following conditions: topological overlap matrix type = “Signed”, merge cut height = 0.25, verbose = 5, and soft thresholding power (β) = 16; the minimum and maximum number of genes in each module was 200 and 2,000, respectively. Also, gene significance (GSi) was calculated using the formula │cor (Xi, degree of I. multifiliis infection)│, and modules with a P-value lower than 0.05 were selected as featured modules in accordance with I. multifiliis infection regardless of target organ. All genes in the featured modules were classified by gene ontology (GO) pathway, and functionally grouped GO assignments were made using ClueGO (39). The functional networks for pathways with P-value < 0.05 were drawn with a medium level of network specificity (0.55) using ClueGO and Cytoscape (35, 39).

Similarity of signed coexpression = 1 + Cor ( Xi , Xj ) 2 Real-Time PCR

To validate sequence quantification, real-time PCR was performed targeting some selected DEGs (SAA, Hap, HSP70, Hmx, CD99, IgT, Wap65-1, and Col). cDNA synthesis and qPCR were performed following a previous study (40). Briefly, cDNAs were synthesized from 1 µg total RNAs. The qPCR mixture contained qPCR Green 2X Mastermix (Mbiotech, Korea) or qPCR probe 2X Mastermix (Mbiotech, Korea) and 400–500 nM of each primer (with 200–250 nM of supplementary probe added if necessary) for a total volume of 25 µl. Detail information on the primer, probe, and incubation conditions is shown in Table S1 .

 Statistical Analysis

All hemato-serological results and osmolality are presented as average ± standard deviation. The differences between the control and I. multifiliis infection and/or recovery groups were tested for statistical significance using Student’s t-test or one way-ANOVA based on Duncan’s analysis. A P-value of less than 0.05 was taken to indicate a significant difference.

 Results Clinical Signs and Hemato-Serology

On average, 61.7 ± 27.4 ciliates per fish were found in the gills of trout in the I. multifiliis group (Ich). Gill aneurysms and hyperplasia were observed in the gills of all fish in the Ich group (data not shown). However, fish in the control (Con) and I. multifiliis recovery group (R-Ich) showed no clinical signs of I. multifiliis infection. The measurements of hematocrit (Ht), serological factors (Glutamic oxaloacetic transaminase, GOT; Glutamic pyruvic transaminase, GPT; Alkaline phosphate, ALP; Blood urea nitrogen, BUN; Glucose, GLU; Total cholesterol, TCHO; Total protein, TP; Lactate dehydrogenase, LDH; and Calcium, Ca) and electrolytes (Na+, Cl−, K+) in the Con, Ich, and R-Ich groups are shown in Table 1 . The Ht readings of fish in the Ich and R-Ich groups were 21 ± 7 and 25 ± 5%, respectively, which are significantly lower than that of the control fish (39 ± 4%). Likewise, some serological indices, such as ALP, TCHO, and TP, were significantly down-regulated in both the Ich and R-Ich groups, but BUN and LDH levels were elevated only in the Ich group. GOT and GPT levels were significantly up-regulated only in the R-Ich group. Although Na+ concentration was constant in all groups in this study, K+ and Cl− levels significantly differed in the Ich group ( Table 1 ). In particular, the K+ concentration in the serum of fish in the Ich group was significantly higher than that of fish in the control group.

Hematological and serological changes among groups (Control, Con; I. multifillis infected group, Ich; and I. multifiliis recovery group, R-Ich).

Group
Con Ich R-Ich
Hematocrit (%) 39 ± 4a 21 ± 7b 25 ± 5b
GOT (U L−1) 421.8 ± 167.8a 601.4 ± 201.1a 1623.9 ± 1227.8b
GPT (U L−1) 24.2 ± 11.3a 27.1 ± 7.9a 63.8 ± 65.1b
ALP (U L−1) 430 ± 197a 121.4 ± 157.4b 152.8 ± 30.7b
BUN (mg dl−1) 2.3 ± 0.5a 3.1 ± 0.7b 2.1 ± 0.3a
GLU (mg dl−1) 101.7 ± 6.6 83.1 ± 32.2 77.1 ± 17.2
TCHO (mg dl−1) 318.7 ± 91.1a 157.3 ± 78.4b 192.6 ± 54.2b
TP (g dl−1) 4.3 ± 0.6a 3.8 ± 0.5b 3.8 ± 0.5b
LDH (U L−1) 1426.7 ± 354.1a 3116.3 ± 1790.1b 1328.5 ± 607.8a
Ca (mg dl −1) 12.5 ± 0.3 11.7 ± 1.2 12.0 ± 0.7
Na+ (mEq L −1) 158.8 ± 8.0 152.5 ± 11.8 152.0 ± 3.0
K+ (mEq L −1) <1.0a 9.1 ± 2.3c 3.8 ± 1.8b
Cl- (mEq L −1) 103.5 ± 13.6a 88.8 ± 9.8b 95.5 ± 2.8ab

Different letters indicated statistically significant differences by Duncan’s multiple range test (P < 0.05).

 Sequencing and Genome-Guided Assembly

Eighteen cDNA libraries (nine head kidney and liver samples) were individually constructed and sequenced using the Illumina HiSeq platform. A total of 889,869,076 raw reads (453,045,574 and 436,823,502 raw reads for the head kidney and liver, respectively) were obtained. On average, 49,437,171 reads were produced for each sample. Of these reads, the majority (97.1%) passed quality control. Finally, 439,458,201 and 424,700,679 clean reads of the head kidney and liver, respectively, were mapped onto the reference genome. Of those, 356,208,296 and 337,667,476 reads were successfully mapped, respectively (average mapping rate = 80.31%) ( Table S2 ). Among 41,978 predicted genes in the rainbow trout, the average numbers of expressed genes in the Con, Ich, and R-Ich groups were 34,061, 31,854, and 34,660, respectively, in the head kidney and 29,875, 28,899, and 30,467 in the liver. Based on the fragments per kilobase per million mapped reads (FPKM) of expressed genes, 3D-PCA plots were prepared for the head kidney and liver, as shown in Figures 1A, B , respectively.

Principal component analysis (PCA) in the head kidney (A) and liver (B). The red and blue numbers indicate up- and down-regulated DEG, respectively, for different comparison sets.

 Selection of DEGs

Genes with a q-value of less than 0.05 were selected as DEGs (41). Six DEG sets, including comparisons of the Con vs. Ich, Con vs. R-Ich, and Ich vs. R-Ich groups in the head kidney and liver, were calculated. For the head kidney, 6,097 DEGs (up-regulation: 3,529; down-regulation: 2,568), 395 DEGs (up-regulation: 208; down-regulation: 187), and 4,111 DEGs (up-regulation: 1,810; down-regulation: 2,301) were found in the comparison of Con vs. Ich, Con vs. R-Ich, and Ich vs. R-Ich, respectively. Likewise, 2,931, 777, and 3,253 genes in Con vs. Ich, Con vs. R-Ich, and Ich vs. R-Ich were identified as DEGs in the liver ( Figure 1 ).

 Pathway Analysis During <italic>I. multifiliis</italic> Infection and Recovery in a Single Organ

All KEGG-annotated DEGs were categorized in hundreds of KEGG pathways. At least one DEG was involved in 367 and 356 different KEGG pathways in the head kidney and liver, respectively. Genes with Z-scores in the positive direction (Z-score > 0) and negative direction (Z-score < 0) in the comparison of the Con vs. Ich groups were considered to be up- and down-regulated KEGG pathways during I. multifiliis infection. Overall, the number of down-regulated pathways (51%; 182/357) was greater than that of up-regulated pathways (46%; 163/357) among pathways with at least one DEG in the head kidney during I. multifiliis infection. However, the liver showed opposite results: 51% (183/356) and 40% (144/356) up- and down-regulated pathways, respectively, among expressed KEGG pathways ( Table S3 ). The top 10% most profoundly altered pathways were selected by calculating the sum of |Z-score| in Con vs. Ich after excluding human disease-relevant pathways. This produced a total of 16 and 18 featured up- and down-regulated pathways, respectively, in the head kidney ( Figure 2 and Table S4 ). Metabolism (ko01100), RNA transport (ko03013), and the spliceosome (ko03040) were the first, second, and third most highly up-regulated pathways, respectively, while focal adhesion (ko04510), Rap1 signaling (ko04015), and ECM-receptor interactions (ko04512) were the most down-regulated pathways in the head kidney ( Table S4 ). Likewise, 18 and 12 pathways were selected as the top 10% up- and down-regulated pathways, respectively, in the liver during I. multifiliis infection ( Figure 3 and Table S5 ). The proteasome (ko03050), protein processing in the endoplasmic reticulum (ko04141), and the lysozyme (ko04142) were the top three up-regulated pathways, while endocrine resistance (ko01522), MAPK signaling (ko04013), and circadian rhythm (ko04710) were the three most down-regulated pathways in the liver ( Table S5 ).

Top 10% enriched KEGG pathways comprised of 16 up- and 19 down-regulated pathways in the head-kidney of rainbow trout during I. multifiliis infection. The red and blue letters indicate up- and down-regulated KEGG pathways, respectively, during I. multifiliis infection. More detailed information is available in Table S4 .

Top 10% enriched KEGG pathways comprised of 18 up- and 12 down-regulated pathways in the liver of rainbow trout during I. multifiliis infection. The red and blue letters indicate up- and down-regulated KEGG pathways, respectively, during I. multifiliis infection. More detailed information is available in Table S5 .

 Co-Expression Analysis With Transcriptomic Responses in the Head Kidney and Liver Under <italic>I. multifiliis</italic> Infection and Recovery

In order to investigate the shared and independent DEGs and pathways between organs during I. multifiliis infection and recovery, the co-expressed pathways and DEGs were profiled. For the Ich and R-Ich groups, pathways with positive Z-scores in both the head kidney and liver were classified as PO-PO, while pathways with negative Z-scores in both organs were considered NE-NE. Likewise, pathways with a negative Z-score in the head kidney but a positive Z-score in the liver were classified as NE-PO, and pathways where Z-score > 0 and Z-score < 0 in the head kidney and liver, respectively, were called PO-NE ( Figure 4 and Tables S6, S7 ). During I. multifiliis infection, the metabolic (ko01100) and proteasome (ko03050) pathways were most highly represented in the PO-PO direction. On the other hand, the focal adhesion (ko04510) and circadian rhythm (ko04710) pathways had the lowest Z-scores in the head kidney and liver, and the featured pathways in the NE-PO and PO-NE directions were cell adhesion molecules (CAMs) (ko04514) and lysosome (ko04142), and cell cycle (ko04110) and MAPK signaling pathway (ko04013), respectively ( Tables S6, S7 ). Approximately 28, 25, 17, and 10% of pathways in the Con vs. Ich group were included in the PO-PO, NE-NE, NE-PO, and PO-NE directions, while 14, 7, 23, and 8% of pathways in the Con vs. R-Ich group belonged to the PO-PO, NE-NE, NE-PO, and PO-NE directions ( Table S8 ). The Z-score and the similarity between groups in the top 10% featured co-expressed pathways were visualized using an enrichment-network map ( Figure 4 ). In general, pathways related to cell replication, transcription, and cell cycles were highly expressed in the head-kidney, while protein folding and metabolism were mutually activated in both organs. But, the cellular signaling pathways such as focal adhesion, ECM-receptor interaction, and Rap1 signaling pathway were down-regulated in the head kidney. Noticeably, only the head kidney, not the liver, showed overall down-regulation of the immune system in the top 10% pathways ( Figure 4 ). In the figures, the dotted line connecting the DEGs and relevant pathways, along with the red and blue colors in nodes and edges, indicate up- and down-regulated DEGs in the head kidney and liver ( Figures S1–S4 ). Only a few DEGs were shared among organs, and most DEGs belonged only to either the head kidney or liver.

Enrichment network analysis based on top 10% co-expressed pathways in the head kidney and liver during I. multifiliis infection and recovery. Z-scores range from blue (lower Z-score) to red (higher Z-score). Edges and edge thickness, which represent biological similarity and the degree of overlap between two pathways, are shown in the light blue line. More detailed information is available in Tables S6, S7 .

 Featured Genes and Network Analysis Based on PLDA Analysis

One hundred featured genes among the total of 46,585 genes were selected after applying the Poisson linear discriminant (PLDA) algorithm, and all the featured genes were standardized by the gene expression level of the control group ( Table S9 ). The fold changes in the head kidney and liver were indicated in the node and edge parts of each featured gene, and network analysis was implemented with clustering by shared Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways ( Figure 5 ). The results showed that most featured genes were involved in cholesterol metabolism (Ko04979), ferroptosis (Ko04216), and the complement and coagulation cascades (Ko04610); also, the levels of several acute-phase proteins such as haptoglobin and serum amyloid A were increased. In addition, there was a reduction in the expression of genes associated with ribosomes and cellular nutrients ( Figure 5 ).

Classified genes from PLDA calculation. Nodes and edges in the circles color are the areas for the gene expression in the head-kidney and liver, respectively, and red and blue color show up- and down-regulated genes. The green circle represents GO pathway, and its size and letters indicate the number of genes belonging to each GO pathway. The gene(s) belonging to GO pathway(s) is marked by gray lines. Each gene expression is described in Table S9 .

 Featured Modules Based on WGCNA

The module in WGCNA is defined as highly interconnected gene clusters, and a total of 19 modules were identified in the module classification, and the red and gray modules out of 19 modules were selected as featured modules ( Figures 6 and S5 ). The red module is highly positively correlated with the degree of I. multifiliis infection (correlation coefficient = 0.46, P-value = 0.05), while the gray module is strongly negatively correlated (correlation coefficient = −0.68, P-value = 0.002). In the red module, the focal adhesion assembly, lipid phosphorylation, embryonic neurocranium morphogenesis, embryotic pattern specification, positive regulation of cell development, nucleobase-containing compound transport cofactor transport, cholesterol biosynthesis, cofactor transport, oxidoreductase activity, and heparin-binding pathways were the major functional groups that were positively correlated with the degree of I. multifiliis infection. On the other hand, epithelial tube morphogenesis, epithelial tube formation, nucleobase-containing compound transmembrane transport activity, regulation of post-transcriptional gene silencing, intraciliary transport particle B, oxidoreductase activity, response to fungus, carboxylic acid biosynthesis, endosomal membrane, endoplasmic reticulum parts, transition metal ion homeostasis, protein localization to the plasma membrane, regulation of BMP signaling pathway, and negative regulation of Wnt signaling pathway were featured functional groups in the gray module, which were down-regulated in accordance with the number of parasites in the gills ( Figure 6 ). Detailed information on each pathway, such as GO ID, term, ontology source, and P-value, is shown in Table S10 .

Correlations between each module’s eigen-gene (the first principal component among the genes in the module) and degree of parasite infection (A). Scatter plots for module membership and gene significance in the red and gray module (B, C). The size of circles and edge thickness represents the number of enriched genes in individual functional groups and the degree of overlap (biological similarity) between two pathways in the red and gray module, respectively (D, E).

 RNA-Seq Validation Using qPCR

In order to validate RNA-seq, eight DEGs were selected in consideration of the biomarkers known to be meaningful for the infectious disease and/or stress environment in former studies and the pattern of transcriptomic expression in this study (40, 42, 43). The expression levels of some DEGs were validated using qPCR. Relative fold changes compared to control were calculated with the ΔCt or ΔΔCt method using an internal control (Ef-1α). A total of eight genes (SAA; haptoglobin, Hap; heat shock protein, HSP70; hemopexin, Hmx; CD99 antigen-like protein, CD99; IgT; warm temperature acclimation protein, Wap65-1; collagen alpha-1 chain, Col) were selected to validate DEGs in the liver or head kidney for the Ich and R-Ich groups. Although the gene expression levels differed slightly between the qPCR and RNA-seq results, the direction and tendency were quite consistent. In the qPCR results, SAA and Col were the most up- and down-regulated gene; the same is true in the RNA-seq results ( Figure 7 ).

Comparison of gene expression between the RNA-seq and qPCR results in the liver (A) and head kidney (B) targeting SAA, Hap, HSP70, Hmx, CD99, IgT, Wap65-1, and Col. The qPCR results in the bar graph are shown as mean ± standard deviation.

 Discussion

Since fish gills and skin are the major portals of entry of I. multifiliis, RNA-seq analysis targeting the primary infection site was conducted accordantly in the former study (19). Syahputra et al. (4) investigated global transcriptomic expression in the gills of rainbow trout infected with I. multifiliis, and found the most significantly altered genes were related to immune pathways such as NOD-like receptor signaling, toll-like receptor signaling, platelet activation, chemokine signaling, leukocyte transendothelial migration, antigen processing and presentation, and T and B cell receptor signaling. However, not many studies have been done on transcriptomics analysis of the internal immunogenic organs despite the significant changes in genes related to the innate and adaptive immune system that occur following I. multifiliis infection (44). The gene expression of cytokines, chemokines, and complement factors in multiple organs (e.g., gill, spleen, liver, and head kidney) were compared using qPCR in previous studies (12, 19), and the authors showed notable different transcriptomic responses after an I. multifiliis infection. But, they profiled only a few selected genes out of more than tens of thousands, which are still insufficient to understand overall responses happening during the infection. This problem has been a major obstacle to understanding the systemic host responses against white spot diseases up until now. In this study, we successfully profiled major transcriptomics in the internal immunogenic organs, and the host biological interaction against I. multifiliis was investigated through the integration of the differential expression in major organs. Among immune-relevant organs, the fish head kidney and liver are well-known primary lymphoid tissues responsible for the immune system, central metabolism (carbohydrate, protein, and lipid metabolism) and detoxification (42, 45, 46). In particular, the liver is a central part of not only the primary immune system (responsible for phagocytosis and complement activation) but also the adaptive immune system (lymphocyte proliferation and antigen presentation) (47). Nevertheless, transcriptomic studies targeting these organs in the rainbow trout under I. multifiliis infection have not been well-described. Accordantly, we focused on profiling and interpreting dual-organ transcriptomic responses (head kidney and liver) to I. multifiliis infection through several types of multi-dimensional data analysis, such as DEGs, co-expression, and machine learning methods, and investigated the systemic interactions among the major organs while considering previous transcriptomic results in the gills (19). Although this study has limitations in that only three trout for each group were used and it was very difficult to obtain accurate disease history of the fish sampled from the fish farm, the results presented in this study greatly contribute to the systems biology of I. multifiliis infection in the fish.

During the course of I. multifiliis infection, some biological pathways were mutually expressed in the same direction in both the head kidney and liver, but other pathways were significantly altered in only one organ or had opposite expression patterns in the two organs. For example, pathways relevant to the endoplasmic reticulum (ER), mitochondria and proteasome, protein processing in the ER, oxidative phosphorylation, and proteasome were up-regulated in both organs. In general, the major function of the ER is to fold and process proteins (48). However, when many unfolded and misfolded proteins are present due to abnormal conditions and diseases, the ER increases its folding capacity to maintain homeostasis (49, 50). Increased expression of ER-related pathways may represent a host strategy to control the numerous inappropriate proteins caused by homeostasis imbalance and/or relevant indirect biological stress because of I. multifiliis infection as shown in some parasitic and bacterial infections (50, 51). Inácio et al. (52) demonstrated the protozoa infection (Plasmodium berghei) stimulated the ER-resident unfolded protein responses, which greatly contributed to parasite development by showing the decrease parasite load under knocking-down the ER transcription factor (cAMP responsive element-binding protein). This implied that the expression of ER-related pathways would be closely linked to the parasite infection in teleost as well. However, since I. multifiliis is an extracellular parasite unlike an intracellular parasite, Plasmodium berghei, further studies are needed to verify the direct relationship between ER and I. multifiliis infection.

In particular, when the system is overwhelmed by an excessive amount of inappropriate proteins, remaining abnormal proteins will be processed through the ER-associated protein degradation (ERAD) pathway and subsequently degraded by proteolytic activities in the proteasome (53). Hence, the high expression of the proteasomal genes seen in this study may be attributable to ERAD. This implies that rainbow trout under physiologically unstable conditions including I. multifiliis infection can increase expression of genes involved in the ER, mitochondria, proteasome, protein processing in the ER, and oxidative phosphorylation in different types of tissues. Also, genes linked to oxidative phosphorylation, the major metabolic pathway through which ATP is produced by cellular respiration, and the metabolism were strongly linked to significant responses to the infection in both organs (53, 54). Abdel-Hafez et al. (55) showed that liver glycogen can be depleted due to loss of osmoregulation during I. multifiliis infection, leading to a significant reduction in serum glucose levels, despite excessively high expression of metabolism-related and oxidative phosphorylation pathways. These transcriptomic responses with low glucose levels are the consequences of rapid energy depletion during I. multifiliis infection, and large expenditures of energy and metabolites could help to repair any damage and counteract the pathogen (56). In addition, most pathways related to cell replication, such as the cell cycle, DNA replication, meiosis, mRNA surveillance, and spliceosome pathways, were among the top 10% most up-regulated KEGG pathways in the head kidney, though they did not among the top 10% of pathways in the liver during I. multifiliis infection. Severe anemia, pale gills, and low numbers of erythrocytes are known as the major clinical signs of I. multifiliis infection (55, 57, 58). Likewise, we observed not only pale gills during infection but significantly lower Ht in the Ich and/or R-Ich groups, which implies that the host was in a state of considerable erythrocyte deficiency. Accordantly, infected trout should produce more erythrocytes than normal, and up-regulation of cell replication pathways provide evidence of activation of hematopoiesis in an attempt to counteract erythrocyte destruction.

Since the head kidney is the primary hematopoietic organ in the teleost, researchers generally think that high expression of cell replication- and hemoglobin-related genes occurs in the head kidney but not the liver (45, 46, 59). An increase in immunological pathways was a representative biological response that showed opposite directions of transcriptomic expression between the head kidney and liver. Overall, in infected trout, greater activation of immunological pathways in the liver was correlated with reduced activation of such pathways in the head kidney. Notably, major chemo-attractive and immune-relevant pathways were down-regulated in the head kidney, but slightly increased in the liver. Likewise, the expression of leukocyte cell-derived chemotaxis-2, one of the featured genes identified from PLDA analysis, was greatly increased in the liver but decreased in the head kidney, strongly indicating movement of leukocytes from the head kidney to the gills and other sites. Unlike other major bacterial and viral diseases that infect internal organs, the major I. multifiliis infection sites are the gills and skin, not the head kidney, and thus this difference is thought to cause leukocyte migration from hematopoietic tissues to the gills and skin (60). To allow for easy leukocyte migration to other tissues, the host needs to decrease the chemo-attractive and immune responses in the head kidney. Since platelet activation and T cell receptor signaling pathways are involved in the stimulation of cytokine expression and the mediation of leukocyte movement, down-regulation of these pathways in only the head kidney could help to maintain a low level of chemo-attractive responses, thereby releasing leukocytes from the head kidney (4, 61, 62). Syahputra et al. (4) also showed changes in the expression of pathways related to chemokine signaling and leukocyte activation, such as the Toll-like and NOD-like signaling pathways, in the gills, as in the liver in the current study. These patterns were also observed in the liver, which showed significant changes in the expression levels of genes related to major immunological responses such as antigen processing and presentation, NOD-like signaling, Fc gamma R-mediated phagocytosis, and natural killer cell-mediated cytotoxicity. Also, Jørgensen et al. (12) also observed overall up-regulation of cytokines (IL-1β, IL-6, INF-γ, and TNF-α) in the liver of rainbow trout at 1–4 week(s) post challenge with I. multifiliis, but no changes or down-regulation of INF- γ were found in the head-kidney of fish at the same sampling time point. Although the liver is not a primary infection site, Castro et al. (42) emphasized that it is an important immunocompetent organ in that numerous pathogens can be ingested orally and detected in the gastrointestinal tract via the portal vein. Since I. multifiliis have free living stages in the water, trout could accidentally uptake many live or dead I. multifiliis including their protein and DNA molecules. Also, the possibility that macrophages in the infection sites would move to the liver through the blood, and then present the I. multifiliis antigens to intraparenchymal macrophages should be considered (42). The activation of uptake of extracellular materials and molecules, including phagocytosis and endocytosis indicates that the liver has very important immunological roles against I. multifiliis infection. But, further studies are necessary to clarify the mode of action of activated phagocytosis and antigen processing and presentation in the liver under I. multifiliis infection. In the recovery stage, when compared to the Ich group, thousands of DEGs were found, but only several hundred DEGs (less than 10% of the number of DEGs between the Con and Ich groups) were found in the comparison between the Con and R-Ich groups. Also, most genes that were significantly altered during I. multifiliis infection returned to similar expression levels in the Con group. For serological results, however, the differences between the Con and Ich groups were much greater than between the Con and R-Ich groups, which is the opposite of the transcriptomic results. This would imply that although transcriptome responses were recovered in accordance with the recovery of I. multifiliis infection, the accumulated damages to the host in R-Ich group had not yet been recovered.

To invade the mucosal and epidermal layer, I. multifiliis can cause cellular destruction and histo-pathological changes such as necrosis, inflammation, and epidermal proliferation (7, 55, 63). Numerous proteases and lytic enzymes from parasites are well-known virulence factors that can destroy host tissues, such as the extracellular matrix and interstitial tissues, through catabolism of host membranes and cytoadherence (64). Several studies have reported on the existence of multiple types of cysteine proteases in I. multifiliis according to gelatin-precast zymography and genomic information from the macro-nucleus (7, 64–66). These virulence factors can directly affect erythrocytes, thus inducing hemolysis and anemia, as shown in this study (decrease in Ht from 39 to 21% after I. multifiliis infection) as well as other studies (55, 57, 58). When erythrocytes are ruptured, excessive Fe2+ and K+ flow into the blood. Although ferrous (Fe2+) is necessary for the function of NADPH oxidase, cytochrome P450, and the electron transport system, excessive Fe2+ catalyzes hydrogen peroxidase (H2O2) to hydroxyl radical (•OH) in “Fenton’s reaction”, the latter of which can induce greater oxidative stress than H2O2 (67, 68) ( Figure 8 ). Hence, massive hemolysis caused by I. multifiliis infection can trigger an influx of Fe2+ in the serum and cause systemic DNA damage, lipid peroxidation, oxidative stress, and ferroptosis. Ferroptosis is an iron-dependent cell death mechanism first defined by Dixon et al. (69). Interestingly, we observed that cell death pathways, such as the ferroptosis and NF-Kappa B signaling pathways, had high Z-scores in both organs. In the context of dysfunction and destruction of both organs, pathways relevant to cellular homeostasis, such as the functional cytoskeleton (ECM-receptor interaction, regulation of actin cytoskeleton), cell adhesion (focal adhesion, gap junction, cell adhesion molecules), and cellular signaling (Ras signaling, cGMP-PKG signaling, cAMP signaling, adrenergic signaling, calcium signaling, Rap1 signaling, and PI3K-Akt signaling), were simultaneously down-regulated in the Ich group, indicating that the function of the kidney was compromised. Additionally, impairment of the kidney and liver was verified by serological evidence (GOT, GPT, and BUN), electrolytes (K+), and osmotic imbalances in the blood. Potassium is a major intracellular electrolyte. Approximately, 98% of K+ is located in intracellular fluid (ICF). Destruction of cells (e.g., cell lysis, hemolysis) can lead to severe hyperkalemia and result in broad physiological disorders (70). In general, excessive K+ in extracellular fluid (ECF) is processed by the kidney. However, compromised renal function such as renal failure can lead to low K+ excretion through urine and failure to recovery K+ in the serum (71, 72). Likewise, the observed increase in BUN, which is readily eliminated by normally functioning kidneys and gills, is another indication of kidney and gill dysfunction in the teleost (73–75). Additionally, the levels of GOT and GPT, which are well-known biomarkers of hepatocyte damage, were approximately three to four times higher in the R-Ich group than the Con group, which provides clear evidence of liver dysfunction and damage. Even though I. multifiliis is normally infected in the skin and gills, we found that its virulence mechanisms and toxins systemically influence the host. The gills and skin are the same in that they are the site of primary infection, but their transcript responses under I. multifiliis infection may be different because of biological function, species, and etc. (19, 76). However, even though I. multifiliis was not directly infected with internal organs (e.g., live, head-kidney), it could greatly influence on systemic damages in rainbow trout as shown in Figure 8 .

Schematic diagram of systemic changes in the liver, head kidney, blood, and primary infection sites (e.g., gill and skin) during I. multifillis infection. The size of red and blue circles with letters represents the number of enriched genes in each KEGG pathway, and edge thickness between two pathways indicates the degree of overlap (biological similarity) during I. multifiliis infection.

One hundred featured genes were selected based on a machine learning method using the PLDA algorithm, which is a widely known nearest shrunken centroids algorithm for discrimination between different groups in micro-array and RNA-Seq data (77). This method is particularly useful for finding genes that are sensitive markers of the difference between control and treatment groups, and is specialized for the discovery of potential biomarkers or featured genes (23, 78). Also, we analyzed the transcriptomic results without considering the target organ to identify major systemic expression regardless of tissue type. Since this is a good approach to finding the key genes that are expressed, the pathways that are highly linked to key genes represent the major host responses to I. multifiliis infection in both organs. As a result, ferroptosis, complement and coagulation cascades, and cholesterol metabolism can be regarded as featured pathways that are co-expressed in multiple organs rather than confined to a single organ under I. multifiliis infection. Although some pathways (e.g., ferroptosis, phagosome) and genes were also investigated in the DEG analysis, others such as complement and coagulation cascades and hemoglobin relevant genes, not classified as DEGs, were more notably identified in PLDA. These results indicated that the application of different approaches for RNA-seq analysis, including machine learning, can broaden the perspectives to investigate biological meaningful pathways that could potentially be missed in a single method. In WGCNA, co-expression networks were generated using the expression tendency of all individual genes and the degree of infection, and red and gray modules were used to represent significant gene clusters that most positively and negatively correlated with parasite burden in this study. More severe I. multifiliis infection led to higher expression of the focal adhesion assembly, positive cell development regulation, and cholesterol biosynthesis pathways, which reflect the process of recovery. Also, heparin-binding growth factors including fibroblast growth factor and vascular endothelial growth factor in the red module are known to play roles in wound healing and vasculature development, which are important modes of action for recovery from an infection (79). In addition, the many pathways related to cell structures and signaling, such as the epithelial tube, endoplasmic reticulum, and plasma membrane, observed in the gray module can be interpreted as representing the cellular damage caused by I. multifiliis infection. This implies that the degree of I. multifiliis infection would mainly affect cellular signaling, development, and maintenance in both the head kidney and liver, indicating that the infection causes systemic responses, not only local responses in specific organ(s) as illustrated in Figure 8 .

In this study, we profiled dual-organ transcriptomics during I. multifiliis infection and recovery through several strategies for multi-dimensional data interpretation. Most highly expressed biological pathways include protein processing in the ER, oxidative stress, and proteasome in both organs. However, the majority of physiological responses showed marked differences between organs in response to I. multifiliis infection. The pathways involved in cell production and movement were most highly affected only in the head kidney, while endocytosis and some crucial cellular signaling responses such as NF-kB signaling were more notably expressed in the liver. Moreover, major immune pathways in the head kidney, especially leukocyte signaling and trans-endothelial migration, as well as chemo-attractive pathways were primarily down-regulated, but strong activation of antigen processing and presentation, leukocyte-mediated cytotoxicity, and phagocytosis was observed in the liver. Significant changes in the levels of serological factors (GOT, GPT, BUN, and K+) along with lower Ht reflect the destruction of tissues and hemolysis. Massive amounts of K+ and Fe2+ were released to the serum by hemolysis. Free Fe2+ stimulates toxic oxidative stress (•OH) through Fenton’s reaction, which induces systemic ferroptosis and leads to multiple organ dysfunction syndromes. Following I. multifiliis infection, control of the expression of various immune systems, restoration of damaged tissues and organs, side effects, and responses to infection simultaneously occur inside the body. The dual-organ transcriptomic strategy used in this study provides a robust tool for the simultaneous interpretation of a variety of biological processes.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

 Ethics Statement

The animal study was reviewed and approved by the Animal Research Ethics Committee at Pukyong National University.

 Author Contributions

HR: Writing—Original draft, Methodology, Conceptualization, Formal analysis, Visualization, Investigation, Methodology. NK, YL, JP, and BK: Methodology, Formal analysis. ML and C-IP: Project administration, Methodology. D-HK: Writing—Review and Editing, Project administration, Funding acquisition, Conceptualization. All authors contributed to the article and approved the submitted version.

 Funding

Ministry of Oceans and Fisheries, Korea.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Acknowledgments

This research was a part of the project titled ‘Development of rapid and sensitive diagnostic methods for the quarantine of aquatic animals and their products’, funded by the Ministry of Oceans and Fisheries, Korea.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.677730/full#supplementary-material

References Matthews R . Ichthyophthirius Multifiliis Fouquet and Ichthyophthiriosis in Freshwater Teleosts. Adv Parasitol (2005) 59:159–241. doi: 10.1016/S0065-308X(05)59003-1 Xu D Zhang Q Shoemaker CA Zhang D Moreira GS . Immune Responses of Channel Catfish Against Parasite Ichthyophthirius Multifiliis Following Theront Vaccination. The American Association of Immunologists, Inc. (2018). Christoffersen T Kania P von Gersdorff Jørgensen L Buchmann K . Zebrafish Danio Rerio as a Model to Study the Immune Response Against Infection With Ichthyophthirius Multifiliis. J Fish Dis (2017) 40:847–52. doi: 10.1111/jfd.12543 Syahputra K Kania PW Al-Jubury A Jafaar RM Dirks RP Buchmann K . Transcriptomic Analysis of Immunity in Rainbow Trout (Oncorhynchus Mykiss) Gills Infected by Ichthyophthirius Multifiliis. Fish Shellfish Immunol (2019) 86:486–96. doi: 10.1016/j.fsi.2018.11.075 Matthews R . Ichthyophthirius Multifiliis Fouquet, 1876: Infection and Protective Response Within the Fish Host. Parasitic Dis Fish (1994), 17–42. Jousson O Di Bello D Donadio E Felicioli A Pretti C . Differential Expression of Cysteine Proteases in Developmental Stages of the Parasitic Ciliate Ichthyophthirius Multifiliis. FEMS Microbiol Lett (2007) 269:77–84. doi: 10.1111/j.1574-6968.2006.00611.x von Gersdorff Jørgensen L . The Fish Parasite Ichthyophthirius Multifiliis–Host Immunology, Vaccines and Novel Treatments. Fish Shellfish Immunol (2017) 67:586–95. doi: 10.1016/j.fsi.2017.06.044 von Gersdorff Jørgensen L Heinecke RD Skjødt K Rasmussen K Buchmann K . Experimental Evidence for Direct In Situ Binding of IgM and IgT to Early Trophonts of Ichthyophthirius Multifiliis (Fouquet) in the Gills of Rainbow Trout, Oncorhynchus Mykiss (Walbaum). J Fish Dis (2011) 34:749–55. doi: 10.1111/j.1365-2761.2011.01291.x Sigh J Buchmann K . Comparison of Immobilization Assays and Enzyme-Linked Immunosorbent Assays for Detection of Rainbow Trout Antibody-Titres Against Ichthyophthirius Multifiliis Fouquet, 1876. J Fish Dis (2001) 24:49–51. doi: 10.1046/j.1365-2761.2001.00258.x Sigh J Lindenstrøm T Buchmann K . Expression of Pro-Inflammatory Cytokines in Rainbow Trout (Oncorhynchus Mykiss) During An Infection With Ichthyophthirius Multifiliis. Fish Shellfish Immunol (2004) 17:75–86. doi: 10.1016/j.fsi.2003.12.005 Sigh J Lindenstrøm T Buchmann K . The Parasitic Ciliate Ichthyophthirius Multifiliis Induces Expression of Immune Relevant Genes in Rainbow Trout, Oncorhynchus Mykiss (Walbaum). J Fish Dis (2004) 27:409–17. doi: 10.1111/j.1365-2761.2004.00558.x von Gersdorff Jørgensen L Nemli E Heinecke RD Raida MK Buchmann K . Immune-Relevant Genes Expressed in Rainbow Trout Following Immunisation With a Live Vaccine Against Ichthyophthirius Multifiliis. Dis Aquat Org (2008) 80:189–97. doi: 10.3354/dao01935 Buchmann K Lindenstrøm T . Sigh J. Partial Cross Protection Against Ichthyophthirius Multifiliis in Gyrodactylus Derjavini Immunized Rainbow Trout. J Helminthol (1999) 73:189–95. doi: 10.1017/S0022149X9900030X Heinecke RD Buchmann K . Inflammatory Response of Rainbow Trout Oncorhynchus Mykiss (Walbaum, 1792) Larvae Against Ichthyophthirius Multifiliis. Fish Shellfish Immunol (2013) 34:521–8. doi: 10.1016/j.fsi.2012.11.036 Olsen MM Kania PW Heinecke RD Skjoedt K Rasmussen KJ Buchmann K . Cellular and Humoral Factors Involved in the Response of Rainbow Trout Gills to Ichthyophthirius Multifiliis Infections: Molecular and Immunohistochemical Studies. Fish Shellfish Immunol (2011) 30:859–69. doi: 10.1016/j.fsi.2011.01.010 Long M Zhao J Li T Tafalla C Zhang Q Wang X . Transcriptomic and Proteomic Analyses of Splenic Immune Mechanisms of Rainbow Trout (Oncorhynchus Mykiss) Infected by Aeromonas Salmonicida Subsp. Salmonicida. J Proteomics (2015) 122:41–54. doi: 10.1016/j.jprot.2015.03.031 Sudhagar A Kumar G El-Matbouli M . Transcriptome Analysis Based on RNA-Seq in Understanding Pathogenic Mechanisms of Diseases and the Immune System of Fish: A Comprehensive Review. Int J Mol Sci (2018) 19:245. doi: 10.3390/ijms19010245 Gonzalez SF Buchmann K Nielsen ME . Complement Expression in Common Carp (Cyprinus Carpio L.) During Infection With Ichthyophthirius Multifiliis. Dev Comp Immunol (2007) 31:576–86. doi: 10.1016/j.dci.2006.08.010 Syahputra K Kania PW Al-Jubury A Marnis H Setyawan AC Buchmann K . Differential Immune Gene Response in Gills, Skin, and Spleen of Rainbow Trout Oncorhynchus Mykiss Infected by Ichthyophthirius Multifiliis. PloS One (2019) 14:e0218630. doi: 10.1371/journal.pone.0218630 Anders S Huber W . Differential Expression of RNA-Seq Data at the Gene Level–the DESeq Package Vol. 10. Heidelberg, Germany: European Molecular Biology Laboratory (EMBL) (2012). p. f1000research. Iancu OD Colville A Oberbeck D Darakjian P McWeeney SK Hitzemann R . Cosplicing Network Analysis of Mammalian Brain RNA-Seq Data Utilizing WGCNA and Mantel Correlations. Front Genet (2015) 6:174. doi: 10.3389/fgene.2015.00174 Merico D Isserlin R Stueker O Emili A Bader GD . Enrichment Map: A Network-Based Method for Gene-Set Enrichment Visualization and Interpretation. PloS One (2010) 5:e13984. doi: 10.1371/journal.pone.0013984 Goksuluk D Zararsiz G Korkmaz S Eldem V Zararsiz GE Ozcetin E . Mlseq: Machine Learning Interface for RNA-Sequencing Data. Comput Methods Programs BioMed (2019) 175:223–31. doi: 10.1016/j.cmpb.2019.04.007 Huang SS Makhlouf M AbouMoussa EH Ruiz Tejada Segura ML Mathew LS Wang K . Differential Regulation of the Immune System in a Brain-Liver-Fats Organ Network During Short-Term Fasting. Mol Metab (2020) 40:101038. doi: 10.1016/j.molmet.2020.101038 Libbrecht MW Noble WS . Machine Learning Applications in Genetics and Genomics. Nat Rev Genet (2015) 16:321–32. doi: 10.1038/nrg3920 Kim A Yoon D Lim Y Roh HJ Kim S Park C . Co-Expression Network Analysis of Spleen Transcriptome in Rock Bream (Oplegnathus Fasciatus) Naturally Infected With Rock Bream Iridovirus (Rbiv). Int J Mol Sci (2020) 21:1707. doi: 10.3390/ijms21051707 Bandoy DDR Weimer BC . Biological Machine Learning Combined With Bacterial Population Genomics Reveals Common and Rare Allelic Variants of Genes to Cause Disease. bioRxiv (2019) 739540. doi: 10.1101/739540 Deng J Kong W Mou X Wang S Zeng W . Identifying Novel Candidate Biomarkers of RCC Based on WGCNA Analysis. Personalized Med (2018) 15:381–94. doi: 10.2217/pme-2017-0091 Purssell RA Pudek M Brubacher J Abu-Laban RB . Derivation and Validation of a Formula to Calculate the Contribution of Ethanol to the Osmolal Gap. Ann Emerg Med (2001) 38:653–9. doi: 10.1067/mem.2001.119455 Berthelot C Brunet F Chalopin D Juanchich A Bernard M Noël B . The Rainbow Trout Genome Provides Novel Insights Into Evolution After Whole-Genome Duplication in Vertebrates. Nat Commun (2014) 5:1–10. doi: 10.1038/ncomms4657 Trapnell C Pachter L Salzberg SL . TopHat: Discovering Splice Junctions With RNA-Seq. Bioinformatics (2009) 25:1105–11. doi: 10.1093/bioinformatics/btp120 Trapnell C Williams BA Pertea G Mortazavi A Kwan G Van Baren MJ . Transcript Assembly and Quantification by RNA-Seq Reveals Unannotated Transcripts and Isoform Switching During Cell Differentiation. Nat Biotechnol (2010) 28:511–5. doi: 10.1038/nbt.1621 Trapnell C Roberts A Goff L Pertea G Kim D Kelley DR . Differential Gene and Transcript Expression Analysis of RNA-Seq Experiments With TopHat and Cufflinks. Nat Protoc (2012) 7:562. doi: 10.1038/nprot.2012.016 Kanehisa M Sato Y Morishima K . BlastKOALA and GhostKOALA: KEGG Tools for Functional Characterization of Genome and Metagenome Sequences. J Mol Biol (2016) 428:726–31. doi: 10.1016/j.jmb.2015.11.006 Cline MS Smoot M Cerami E Kuchinsky A Landys N Workman C . Integration of Biological Networks and Gene Expression Data Using Cytoscape. Nat Protoc (2007) 2:2366. doi: 10.1038/nprot.2007.324 Love MI Huber W Anders S . Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol (2014) 15:1–21. doi: 10.1186/s13059-014-0550-8 Conesa A Götz S García-Gómez JM Terol J Talón M Robles M . Blast2GO: A Universal Tool for Annotation, Visualization and Analysis in Functional Genomics Research. Bioinformatics (2005) 21:3674–6. doi: 10.1093/bioinformatics/bti610 Langfelder P Horvath S . WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:1–13. doi: 10.1186/1471-2105-9-559 Bindea G Mlecnik B Hackl H Charoentong P Tosolini M Kirilovsky A . ClueGO: A Cytoscape Plug-in to Decipher Functionally Grouped Gene Ontology and Pathway Annotation Networks. Bioinformatics (2009) 25:1091–3. doi: 10.1093/bioinformatics/btp101 Roh H Kim A Kim N Lee Y Kim D . Multi-Omics Analysis Provides Novel Insight Into Immuno-Physiological Pathways and Development of Thermal Resistance in Rainbow Trout Exposed to Acute Thermal Stress. Int J Mol Sci (2020) 21:9198. doi: 10.3390/ijms21239198 Benjamini Y Hochberg Y . Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat society: Ser B (Methodol) (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x Castro R Abós B Pignatelli J von Gersdorff Jørgensen L Granja AG Buchmann K . Early Immune Responses in Rainbow Trout Liver Upon Viral Hemorrhagic Septicemia Virus (VHSV) Infection. PloS One (2014) 9:e111084. doi: 10.1371/journal.pone.0111084 Kania PW Chettri JK Buchmann K . Characterization of Serum Amyloid A (SAA) in Rainbow Trout Using a New Monoclonal Antibody. Fish Shellfish Immunol (2014) 40:648–58. doi: 10.1016/j.fsi.2014.08.018 Wang Q Yu Y Zhang X Xu Z . Immune Responses of Fish to Ichthyophthirius Multifiliis (Ich): A Model for Understanding Immunity Against Protozoan Parasites. Dev Comp Immunol (2019) 93:93–102. doi: 10.1016/j.dci.2019.01.002 Whyte SK . The Innate Immune Response of Finfish–A Review of Current Knowledge. Fish Shellfish Immunol (2007) 23:1127–51. doi: 10.1016/j.fsi.2007.06.005 Bracamonte SE Johnston PR Knopf K Monaghan MT . Experimental Infection With Anguillicola Crassus Alters Immune Gene Expression in Both Spleen and Head Kidney of the European Eel (Anguilla Anguilla). Marine Genomics (2019) 45:28–37. doi: 10.1016/j.margen.2018.12.002 Parker GA Picut CA . Immune Functioning in Non Lymphoid Organs: The Liver. Toxicol Pathol (2012) 40:237–47. doi: 10.1177/0192623311428475 Oslowski CM Urano F . Measuring ER Stress and the Unfolded Protein Response Using Mammalian Tissue Culture System. Meth Enzymol (2011) 490:71–92. doi: 10.1016/B978-0-12-385114-7.00004-0 Ron D Walter P . Signal Integration in the Endoplasmic Reticulum Unfolded Protein Response. Nat Rev Mol Cell Biol (2007) 8:519–29. doi: 10.1038/nrm2199 Chaubey S Grover M Tatu U . Endoplasmic Reticulum Stress Triggers Gametocytogenesis in the Malaria Parasite. J Biol Chem (2014) 289:16662–74. doi: 10.1074/jbc.M114.551549 Ni M Lee AS . ER Chaperones in Mammalian Development and Human Diseases. FEBS Lett (2007) 581:3641–51. doi: 10.1016/j.febslet.2007.04.045 Inácio P Zuzarte-Luís V Ruivo MT Falkard B Nagaraj N Rooijers K . Parasite-Induced ER Stress Response in Hepatocytes Facilitates Plasmodium Liver Stage Infection. EMBO Rep (2015) 16:955–64. doi: 10.15252/embr.201439979 Martin S Blaney S Bowman A Houlihan D . Ubiquitin-Proteasome-Dependent Proteolysis in Rainbow Trout (Oncorhynchus Mykiss): Effect of Food Deprivation. Pflügers Archiv (2002) 445:257–66. doi: 10.1007/s00424-002-0916-8 Lee I Hüttemann M . Energy Crisis: The Role of Oxidative Phosphorylation in Acute Inflammation and Sepsis. Biochim Biophys Acta (BBA)-Molecular Basis Dis (2014) 1842:1579–86. doi: 10.1016/j.bbadis.2014.05.031 Abdel-Hafez G Lahnsteiner F Mansour N Licek E . Pathophysiology of Ichthyophthirius Multifiliis Infection in Rainbow Trout (Oncorhynchus Mykiss) and Chub (Leuciscus Cephalus). J Comp Pathol (2014) 151:394–9. doi: 10.1016/j.jcpa.2014.08.003 Antunes LC Arena ET Menendez A Han J Ferreira RB Buckner MM . Impact of Salmonella Infection on Host Hormone Metabolism Revealed by Metabolomics. Infect Immun (2011) 79:1759–69. doi: 10.1128/IAI.01373-10 Witeska M Kondera E Lugowska K . The Effects of Ichthyophthiriasis on Some Haematological Parameters in Common Carp. Turkish J Veterinary Anim Sci (2010) 34:267–71. doi: 10.3906/vet-0903-26 Tavares-Dias M Moraes FD Martins M Santana A . Haematological Changes in Oreochromis Niloticus (Osteichthyes: Cichlidae) With Gill Ichthyophthiriasis and Saprolegniosis. Bol Inst Pesca (2002) 28:1–9. Kondera E . Haematopoiesis in the Head Kidney of Common Carp (Cyprinus Carpio L.): A Morphological Study. Fish Physiol Biochem (2011) 37:355–62. doi: 10.1007/s10695-010-9432-5 Muller WA . Mechanisms of Leukocyte Transendothelial Migration. Annu Rev Pathol: Mech Dis (2011) 6:323–44. doi: 10.1146/annurev-pathol-011110-130224 Yun S Sim E Goh R Park J Han J . Platelet Activation: The Mechanisms and Potential Biomarkers. BioMed Res Int (2016) 2016:1–7. doi: 10.1155/2016/9060143 Lam WA Chaudhuri O Crow A Webster KD Kita A Huang J . Mechanics and Contraction Dynamics of Single Platelets and Implications for Clot Stiffening. Nat Mater (2011) 10:61–6. doi: 10.1038/nmat2903 Hines RS Spira DT . Ichthyophthiriasis in the Mirror Carp Cyprinus Carpio (L.) V. Acquired Immunity. J Fish Biol (1974) 6:373–8. doi: 10.1111/j.1095-8649.1974.tb04554.x Piña-Vázquez C Reyes-López M Ortíz-Estrada G de la Garza M Serrano-Luna J . Host-Parasite Interaction: Parasite-Derived and -Induced Proteases That Degrade Human Extracellular Matrix. J Parasitol Res (2012) 2012:1–24. doi: 10.1155/2012/748206 Klemba M Goldberg DE . Biological Roles of Proteases in Parasitic Protozoa. Annu Rev Biochem (2002) 71:275–305. doi: 10.1146/annurev.biochem.71.090501.145453 Coyne RS Hannick L Shanmugam D Hostetler JB Brami D Joardar VS . Comparative Genomics of the Pathogenic Ciliate Ichthyophthirius Multifiliis, Its Free-Living Relatives and a Host Species Provide Insights Into Adoption of a Parasitic Lifestyle and Prospects for Disease Control. Genome Biol (2011) 12:1–26. doi: 10.1186/gb-2011-12-10-r100 Crichton R . Iron Metabolism: From Molecular Mechanisms to Clinical Consequences. The Atrium, Southern Gate, Chichester, West Sussex, PO19 8SQ, United Kingdom: John Wiley & Sons (2016). Latunde-Dada GO . Ferroptosis: Role of Lipid Peroxidation, Iron and Ferritinophagy. Biochim Biophys Acta (BBA)-General Subj (2017) 1861:1893–900. doi: 10.1016/j.bbagen.2017.05.019 Dixon SJ Lemberg KM Lamprecht MR Skouta R Zaitsev EM Gleason CE . Ferroptosis: An Iron-Dependent Form of Nonapoptotic Cell Death. Cell (2012) 149:1060–72. doi: 10.1016/j.cell.2012.03.042 Hoskote SS Joshi SR Ghosh AK . Disorders of Potassium Homeostasis: Pathophysiology and Management. JAPI (2008) 56:685–725. Singh NN Das VK Srivastava AK . Insecticides and Ionic Regulation in Teleosts: A Review. Zoologica Poloniae (2002) 47:21–36. Lee JW . Fluid and Electrolyte Disturbances in Critically Ill Patients. Electrolyte Blood Press (2010) 8:72–81. doi: 10.5049/EBP.2010.8.2.72 Nelson K Jones J Jacobson S Reimschuessel R . Elevated Blood Urea Nitrogen (BUN) Levels in Goldfish as an Indicator of Gill Dysfunction. J Aquat Anim Health (1999) 11:52–60. doi: 10.1577/1548-8667(1999)011<0052:EBUNBL>2.0.CO;2 Bernet D Schmidt H Wahli T Burkhardt-Holm P . Effluent From a Sewage Treatment Works Causes Changes in Serum Chemistry of Brown Trout (Salmo Trutta L.). Ecotoxicol Environ Saf (2001) 48:140–7. doi: 10.1006/eesa.2000.2012 Yang J Chen H . Effects of Gallium on Common Carp (Cyprinus Carpio): Acute Test, Serum Biochemistry, and Erythrocyte Morphology. Chemosphere (2003) 53:877–82. doi: 10.1016/S0045-6535(03)00657-X Jørgensen L Korbut R Jeberg S Kania PW Buchmann K . Association Between Adaptive Immunity and Neutrophil Dynamics in Zebrafish (Danio Rerio) Infected by a Parasitic Ciliate. PloS One (2018) 13:e0203297. doi: 10.1371/journal.pone.0203297 Zararsız G Goksuluk D Korkmaz S Eldem V Zararsiz GE Duru IP . A Comprehensive Simulation Study on Classification of RNA-Seq Data. PloS One (2017) 12:e0182507. doi: 10.1371/journal.pone.0182507 Witten DM . Classification and Clustering of Sequencing Data Using a Poisson Model. Ann Appl Stat (2011) 5:2493–518. doi: 10.1214/11-AOAS493 Nanba D Higashiyama S . Dual Intracellular Signaling by Proteolytic Cleavage of Membrane-Anchored Heparin-Binding EGF-Like Growth Factor. Cytokine Growth Factor Rev (2004) 15:13–9. doi: 10.1016/j.cytogfr.2003.10.002