Experiments in this study were conducted in accordance with Guidelines of Animal Research and Ethics Committees of Ocean University of China (Permit Number: 20141201), U.K. Animals (Scientific Procedures) Act, 1986 and associated guidelines, EU Directive 2010/63/EU for animal experiments, and use of laboratory animals (NIH Publications No. 8023, revised 1978) National Institutes of Health Guide for the Care and Use of Laboratory Animals (NIH publication no. 8023, revised 1978). No endangered or protected animal species were used. The effects of sex were not considered because trout juveniles are sexually immature.

Rainbow trout juveniles were obtained from Linqu Salmon and Trout Aquatic Breeding LLC (Weifang, Shandong, China). These juveniles were from the same full-sibling family batch and spawned on the same day with synchronized development. Trout were acclimatized for 14 days in indoor cuboidal tanks equipped with a water pump, chiller system, sand filter, and biofilter at the Experimental Fish Facility in Key Laboratory of Mariculture, Ocean University of China. According to the Standards of Linxia Salmon and National Trout Elite Breeding and Protection Farm (Linxia, Gansu, China, Approved by Department of Agriculture, China, 2009), trout were cultured at ~16°C and ~7 mg/L of dissolved oxygen. Trout were fed a commercial pellet twice a day at 7% of total body weight.

The V. anguillarum strain was obtained from the Laboratory of Pathology and Immunology of Aquatic Animals, Ocean University of China (41, 42). The bacteria were grown overnight at 28°C in 2216E medium. The bacterial suspension was then centrifuged and resuspended with 0.01 M phosphate-buffered saline (PBS, pH = 7.2). V. anguillarum suspension density was adjusted to serial dilutions for preliminary testing: 10⁹, 10⁸, or 10⁷ colony forming units (CFU)/ml (33).

This manuscript used the same RNA-Seq samples previously described in two papers evaluating the growth hormone and insulin-like growth factor axes, as well as the caspase gene family in rainbow trout (33, 34). Previous studies showed 10⁷ to 10⁹ CFU/ml of V. anguillarum could cause vibriosis in rainbow trout and other teleosts (5, 41–43). Our published paper further showed that V. anguillarum of 10⁷ CFU/ml at 20°C exhibited mild to moderate symptoms of vibriosis disease with a relatively lower mortality (33). Therefore, trout were challenged by 10⁷ CFU/ml of V. anguillarum at 20°C. In the control group, 90 trout were randomly distributed into three tanks, with 30 trout in each tank. The control trout (CT) were intraperitoneally injected with 200 μl physiological saline (saline-challenged, 0.9% NaCl). In the challenged group, 90 trout were equally and randomly distributed into three tanks. Trout of the challenged group were challenged by intraperitoneal injection of 200 μl V. anguillarum (10⁷ CFU/ml). In challenged groups, the first three erratically swimming moribund trout showing severe symptoms, such as hemorrhage in fins, in tank #1 were pooled as sample #1 of the symptomatic trout (ST). After 120 h post-challenge, the three surviving trout with minor or no symptoms were pooled as sample #1 of the asymptomatic trout (AT) ( Figure 1 ). Likewise, sample #2 of ST and AT, as well as sample #3 of ST and AT, were collected from tank #2 and tank #3, respectively ( Figure 1 ). Trout were anesthetized by MS-222 (35–45 mg/L [ppm]) before sampling. Biological samples of organs and tissues (brain, spleen, kidney, liver, and gill) were isolated and stored at −80°C for further analysis.

A total of 27 libraries [3 tissues (brain, kidney, spleen) × 3 replicated samples (each sample contained three pooled individuals × 3 treatment groups] was constructed via the TruSeq™ RNA Sample Prep Kit (Illumina, CA, USA). This study used the same RNA-Seq data with our previously published paper (33, 34), but we focused on different functional genes and used various analyses. The sequence reads are available from the NCBI sequence read archive (SRA) with the accession number of PRJNA667799.

The amino acid sequences of trout novel C3-1 proteins, and zebrafish (Danio rerio), southern catfish (Silurus meridionalis), rat (Rattus norvegicus), and human (Homo sapiens) C3 proteins were used for the phylogenetic analysis and sequence alignment. Phylogenetic analyses were plotted using the Neighbor-joining (N-J) method via MEGA 7, with 1000 bootstrap replications for phylogeny. The SWISS-MODEL between trout and mammalian C3 proteins was generated using the SWISS-MODEL (https://swissmodel.expasy.oAT/) (44, 45). The mammalian C3 with an intact thioester at 3Å resolution [PDB ID: 2B39 (46)] was used as the template. Comparison of the domains between trout and mammalian C3 and the cartoon, stick, and sphere structures of the proteins were generated with the PyMOL software package (47, 48).

Based on published papers on biomedical and fishery studies (49–51), the peak intensity tables of selected genes were uploaded to the websites of MetaboAnalyst and NetworkAnalyst (ATtps://www.xialab.ca/tools.xATml) for data processing and analyses (52). The uploaded data (count normalized by DESeq2 package in the R software (53)) were performed by sum normalization, thus obtaining the belt data (Poisson) distribution for further statistical analysis ( Figure S1 ). In the multivariate analysis module of MetaboAnalyst, the normalized data were then subject to principal component analysis (PCA) and partial least squares discriminant analysis (PLSDA) for pattern discovery ( Figure S1 ). Genes of each pairwise comparison (ST/CT, AT/CT, or AT/ST) were selected to create a heatmap (Based on log₁₀(normalized count+1)) and correlation analysis (with Pearson’s correlation) (51).

The heatmap displayed the transcriptional profile of genes associated with the stress response, cytokines and cellular functions, and the complement system between ST and CT ( Figures 2A–C ). The overall transcriptional profiles of target genes in ST and CT in response to V. anguillarum infection were summarized by PCA ( Figure 2D ). Red dots show the vector containing overall gene expression in CT, and green dots showed the vector containing overall gene expression in ST. Separated PCA vectors were present, indicating that the V. anguillarum infection resulted in different profiles of genes associated with the stress response, cytokines and cellular functions, and the complement system between ST and CT ( Figure 2D ). The loading plot of PCA shows the genes exerting stronger influences on PCA analysis ( Figure 2E , points far away from the zero point, Table S1 ).

The volcano plots showed that, compared to CT, the ST showed significantly downregulated kidney mrα, mrβ, c7-2, and cd93, and spleen grα, grβ, c7-2, and c1qa, and brain c7-2 and c3-4 ( Figures 2F, G ). Compared to CT, the kidney il11, mbl-h2, and c3-1b1, and spleen il1β1, il1β2, il8, tnfα2, c3-1a, c3-1b1, and c3-3 were significantly upregulated in ST ( Figures 2F, G ). The genes showed in volcano plots were labeled in the loading plot ( Figure 2E ).

The correlation analysis of all target genes is depicted using a heatmap ( Figure 2H and Figure S2 ). The Pearson correlation coefficients showed that the kidney mrα or mrβ exhibited strong negative relationships with the cytokines of il1β3, il4, il8, and tnfα3 ( Figures 2H–J ). The spleen grβ showed negative relationships with il1β1, il1β2, il8, and tnfα2 ( Figures 2H, K ).

The transcriptional profiles of genes involved in the stress response, cytokines, cellular functions, and the complement system between AT and CT were shown by heatmap ( Figures 3A–C ). Separated PCA plots indicate that genes related to cytokines, the stress response, cellular functions, and the complement system were differently expressed in AT and CT ( Figure 3D ). The loading plot showed the genes significantly involved in the separated PCA plots ( Figure 3E , points far away from the zero point, Table S2 ).

The volcano plots showed that, compared to CT, the AT showed downregulated kidney c7-2, cd93, and spleen mrα, mrβ, grβ, hsd11β2, c7-2, c6, and c8g, and brain c7-2 ( Figures 3F, G ). The AT exhibited significantly upregulated kidney il1β2, il11, c4, and mbl-h2, and spleen il1β1, il1β2, il1β3, il6, tnfα3, cfb, cfp1, c3-1b2, c3-3, and bcf2-b ( Figures 3F, G ). These genes were highlighted in the loading plot ( Figure 3E ).

Heatmap showing the Pearson correlation coefficients of genes ( Figure 3H and Figure S3 ). Pearson correlation coefficients showed that the spleen mrα or mrβ exhibited strong negative relationships with the cytokines of il1β1, il1β2, il1β3 and tnfα3 ( Figures 3H–J ). The spleen grβ exhibited negative Pearson correlation coefficients with il1β3, and tnfα3 ( Figures 3H, L ), while the kidney grα showed negative relationships with il1β1, il1β3, tnfα1 and tnfα3 ( Figures 3H, K ).

We compared the overall gene expression between ST and AT by heatmap ( Figures 4A–C ). In PCA plots, vectors showing gene expression in ST were separated from those showing gene expression in AT, demonstrating differential gene expression between ST and AT in response to V. anguillarum infection ( Figure 4D ). The loading plot showed key genes resulting in discrimination and stronger influences on PCA vectors ( Figure 4E , points far away from the zero point, Table S3 ).

Volcano plots showed the expression of genes (kidney hsd11β2, sod3, mrα, and mrβ, and spleen c1qa, and brain pomcβ, cat, c3-4, and c7-2) in ST were significantly lower than those of AT ( Figures 4F, G ). Compared to AT, ST showed upregulated gene expression of spleen c7-1 and c3-1b1 and brain c8g ( Figures 4F, G ). The Pearson correlation coefficients of target genes are indicated by a heatmap ( Figure 4H and Figure S4 ). The kidney grα and grβ showed strong negative relationships with the cytokines of il10b and il4, respectively ( Figures 4H, I, K ).

We identified three novel c3 gene subtypes in RNA-Seq data. Based on the alignment of the amino acid sequences, these three C3 proteins showed the conserved functional domains, including the ANATO domain, thioester domain, and C3-convertase cleavage site ( Figure 5A and the whole sequences alignment are shown in Figure S7 ). Based on mammalian C3 (PDB ID: 2B39), the SWISS-MODEL illustrated conserved motifs between trout and mammalian C3 with blue cartoons ( Figure 5B and Figures S7, S8, S9A, S9B, S9C ). Red and green boxes mark the ANATO and thioester domains, respectively ( Figures 5B, C ). The comparison of the thioester (green) and ANATO (red) domains between trout and mammalian C3 are shown in cartoons ( Figure 5C and Figures S9B, S9C ). The conserved amino acid sequences of GCGEQ in thioester domain were labeled ( Figure 5C , top figure). The locations were adjacent, and the identities were identical for both GCGEQ sequences of mammalian and trout C3 ( Figure 5C , top figure; Figure S9C ). Likewise, the ANATO domains of both mammalian and trout C3 are similarly organized, and their amino acid sequences were highly identical ( Figure 5C , bottom figure; Figure S9B ). The gene expression levels of three c3 were shown ( Figures S9D–F ). Compared to CT, the ST showed significantly upregulated c3-1a and c3-1b1 expression in the kidney and spleen ( Figures S9D, E ).

Compared to CT, the AT showed upregulated Ko04610 (complement and coagulation cascades) in the kidney and spleen ( Figure 5D ). In contrast, the ST showed a downregulated Ko04610 pathway in the brain, kidney, and spleen ( Figure 5D ). No significant changes in the Ko04610 pathway were observed in the kidney and spleen between AT and ST. The overlapping genes in the Ko04610 pathway are shown in Venn diagrams ( Figures S9G, S9H ), and their expression levels among CT, AT, or ST were shown by heatmap ( Figures S9I–S9L ).

We showed grα and grβ were shared in the list of DEGs between groups of CT and ST or CT and AT. In the gene ontologies (GO) terms involved in grα and grβ, 8 GO terms were shared between the comparisons of CT and ST or CT and AT ( Figure 6A , details in Table 2 ). Three GO terms were specifically enriched in the comparison of CT and ST ( Figure 6B , details in Table 2 ), and five GO terms were specifically enriched in the comparison of CT and AT ( Figure 6B , details in Table 2 ). Likewise, tnfα subtypes were shared in the DEGs list between CT and ST or CT and AT. Five GO terms associated with tnfα were both identified between the comparisons of CT and ST or CT and AT ( Figure 6C , details in Table 2 ). Three GO terms were specifically enriched in the comparison of CT and ST ( Figure 6D , Details in Table 2 ), and 14 GO terms were specifically enriched in the comparison of CT and AT ( Figure 6D , details in Table 2 ).

Genes of mrα and mrβ were identified in the list of DEGs between the comparison of ST and AT. Based on the KEGG database, four pathways that are associated with functions of steroid hormones were enriched ( Figure 7A ), including ko04960 (aldosterone-regulated sodium reabsorption, Figure 7B ), ko04978 (mineral absorption), ko00140 (steroid hormone biosynthesis, Figure 7C ) and ko04913 (ovarian steroidogenesis, Figure 7D ).

The C3 serves as a major acute-phase protein (55). The expression of c3 gene subtypes is significantly upregulated in response to bacterial or LPS stimulation in multiple teleosts, including the dojo loach (Misturnus anguillicaudatus), rainbow trout, southern catfish (Silurus meridionalis), and grass carp (Ctenopharyngodon idella) (56–59). Consistently, our study found that the trout c3 gene subtypes showed upregulation in responses to V. anguillarum infection. Salmonidae species, such as trout and salmon, experienced four rounds of genome duplication. Consequently, the genetic expansions are characterized by duplicated functional gene copies (paralogs) in Salmonidae fishes (60, 61). Previous studies identified multiple trout c3 subtypes (c3-1, c3-3, and c3-4) with functional diversity (62, 63). Our study identified three novel subtypes within c3-1 (c3-1a, c3-1b1, and c3-1b2) ( Figure 5 and Figures S7–S9 ). These genes exhibited conserved sequence identity but specific expression patterns in responses to V. anguillarum infection ( Figure 5 and Figures S7–S9 ), indicating that these genes can encode bioactive proteins with diversity in functions.

The complement system served as a major governor of inflammatory responses (64). The homeostasis of inflammatory reactions plays a vital role in modulating health balance. Either inefficient or overactive activation of the complement system could disturb the homeostasis, which is detrimental for health balance (64, 65). Compared to CT, the kidney and spleen of ST exhibited downregulated complement cascades (Ko04610). Previous studies in mice indicated that the inefficient activation of complement cascades might be associated with increased susceptibility to infectious diseases (64, 66). Therefore, the ST showed severe symptoms in response to V. anguillarum infection. The complement and coagulation cascades belong to a complex inflammation regulatory network (67, 68). In most of the pathophysiological processes, both the complement and coagulation cascades are activated simultaneously (69). Consistent with the downregulated complement cascades, key genes in coagulation cascades were downregulated in the kidney and spleen of ST, including vwf subtypes (von Willebrand factor), α2m (α-2-macroglobulin), and f13a (Coagulation factor XIII A chain). The ST also showed downregulated platelet activation (ko04611, Figure S10 ). The downregulated coagulation cascades and platelet activation probably caused severe hemorrhages in the fins, kidneys, and other visceral masses in ST, all of which were lethal to the trout. Studies in biomedical and fishery sciences showed healthy individuals could efficiently regulate the complement system, thus not only preventing the complement(s) exhaustion but also enabling the complement(s) to restore (21, 70). However, the moribund trout might fail to efficiently regulate the complement system. The complement exhaustion further reduced the defense to pathogen infection and eventually caused the worse outcomes (death) (71, 72).

The complement system can activate the innate immune system and thus play an essential role in linking the innate and adaptive systems in mammals and teleosts (18, 35, 55, 73). The AT showed upregulated complement and coagulation cascades, enabling the AT to fight the inflammatory pathogenesis and prevents life-threatening bleeding (69). Consistently, AT had higher fga and fgb expression (fga and fgb: fibrinogen α/β chain, which has a significant function in hemostasis, Figure S10 ). Based on these pieces of evidence, we propose that the different responses of complement and coagulation cascades might link to varying phenotypes of trout in response to V. anguillarum infection. A recent study showed that complement cascades serve as a bridge between immunomodulation in trout in response to bacterial infection (74), consistent with what we found in our study.

The cytokine networks govern the normal development and physiology in animals, and dysregulations of cytokine networks are involved in pathophysiological alternations (75). In humans, the IL1 serves as the most potent endogenous pyrogens in organisms affected by infectious diseases (76, 77). Likewise, IL1 plays an apical role in initiating inflammatory responses in teleosts (78), and V. anguillarum infection results in significantly upregulated il1β in teleosts, including Atlantic cod (Gadus morhua), sea bream, and European sea bass (79–82). Our study showed that ST and AT exhibited significantly upregulated il1β subtypes ( Figures 2 and 3 , Table 1 ). ST and AT also showed upregulated tnfα subtypes ( Figures 2 and 3 , Table 1 ), which was consistent with a previous study showing that the functions of IL1 and TNF largely overlap in teleosts (83). Indeed, the IL1 and TNF work synergistically, and the TNF usually acts as the first cytokine to follow an IL1 surge in an inflammatory response (83). Like IL1β, IL11 could regulate a series of important immunomodulatory effects by affecting proliferation and differentiation of hematopoietic progenitors, thus serving as a multifunctional modulator (84, 85). Studies showed kidney il11 was significantly upregulated in response to bacterial pathogens in golden pompano (Trachinotus ovatus) (86), which is in line with our results ( Figures 2 and 3 , Table 1 ).

In addition to upregulated cytokine genes (il1β subtypes, tnfα subtypes, and il11), ST and AT showed specifically upregulated il8 and il6, respectively ( Figures 2 and 3 ). IL6 and IL8 are two important proinflammatory cytokines that play an important role in regulating local or systemic inflammation (87). Studies showed both IL1α and IL1β subtypes could initiate the signal transduction and trigger the expression of IL6 and IL8 (12, 88). Consistently, this study revealed strong positive relationships between the expression of il6/il8 and il1β subtypes ( Figures 2 and 3 , Figures S2, S3, S6 ). For example, the il1β3 and il6 were both upregulated in AT rather than ST ( Figures 3 ). During evolution, the IL1α is evolving faster than IL1β, thus resulting in decreased sequence and functional homology between trout and mammalian IL1α orthologs (89, 90). Our further studies will investigate whether the evolutionally conserved IL1β exhibits subtype-specific IL6/IL8 expression regulation.

Compared to trout in ST, trout in AT exhibited more upregulated GO terms associated with immune defenses and the resulting intracellular signaling ( Figure 6 and Table 2 ), including GO:0051607, defense response to virus; GO:0035631, CD40 receptor complex; GO:0002768, immune response-regulating cell surface receptor signaling pathway; GO:2000353, positive regulation of endothelial cell apoptotic process; GO:0043123, positive regulation of I-κB kinase/NF-κB signaling; GO:0051092, positive regulation of NF-κB transcription factor activity; and GO:0042531, positive regulation of tyrosine phosphorylation of STAT protein. Despite limited studies on TNF-regulated intercellular and intracellular signaling transduction in teleosts, the in vivo studies on humans and rodents provide a potential model that could describe the immune mechanisms specifically activated in AT. Relevant to the GO terms of GO:2000353, GO:0043123, GO:0051092, and GO:0042531, previous biomedical studies showed TNFα activates the intracellular NF-κB signaling, while the cytoplasmic STAT serves as a negative regulator of TNFα-triggered NF-κB activation (91). The activation of NF-κB signaling and NF-κB transcriptional factors maintains an evolutionarily conserved and important role in initiating and coordinating the innate and adaptive immune responses (92). The phosphorylated STAT dimer will translate and localize to the nucleus, where it cannot interact with the TNFα-receptor complex. STAT localization to the nucleus allows a more robust TNFα-triggered NF-κB activation (91), enabling the trout to activate the immune defenses in response to V. anguillarum infection ( Figure 6E ).

In addition to the GR, the teleost MR also serves as a receptor for stress perception. Our results showed the asymptomatic trout showed upregulated kidney mrα and mrβ expression. Consistently, previous studies showed MR and/or GR are expressed in immune tissues and regulate the immunomodulation (93–95). Moreover, increased stress hormone levels are observed in trout and zebrafish treated with the V. anguillarum vaccine (1, 96). Indeed, bidirectional communication exists between stress and immune responses, and low levels of stress (eustress) may result in enhanced immune competence (97). The slightly upregulated mrα and mrβ could act as an alarm and stimulate the immune system to fight against V. anguillarum infection, consistent with previous studies (10, 11).

Studies on humans, rodents, and other mammals showed that cytokines could affect the genes associated with the stress response through cytokine-specific mechanisms. For example, IL1 and IL6 exhibit positive effects, while the TNFα exhibits the opponent manners (16, 98, 99). These cytokines have also been reported to dysregulate and/or block the functions of GR subtypes (17). In teleosts, the immune responses regulated by the interactions between the genes in the stress response and cytokine networks are not homogeneous. Previous studies showed it is greatly affected by specific characteristics of challenges (environmental stressors or disease pathogens), target tissues (in vitro or in vivo; peripheral tissues or mucosal surfaces), and the adaptive life story of each species (bream, bass, or trout) (1, 5, 13, 100). In this study, the results showed that downregulated mr and gr subtypes exhibited strong negative relationships with cytokine genes of il1β and tnfα subtypes in AT and ST ( Figures 2 and 3 ), which were partially consistent with the results in mammalian studies. Previous studies in sea bream showed that the stress response can suppress the gene expression of cytokines (13). These results indicated that the stress response and cytokine networks are intimately and bidirectionally linked, enabling teleosts to cope with challengers from environmental stimuli and/or pathogen invasion (8, 10, 83, 97).

After V. anguillarum infection, ST and AT both exhibited significantly downregulated GO terms associated with functions of grα and grβ (GO:0004883, glucocorticoid receptor activity; GO:0031963, cortisol receptor activity; GO:0005496, steroid binding), and significantly upregulated GO terms that are involved in tnfα-regulated immune responses (GO:0042832, defense response to protozoan; GO:0030890, positive regulation of B cell proliferation) ( Figure 6 and Table 2 ). The upregulated il1β, il6, il8, and tnfα genes are markers of M1 macrophage polarization, which activates the proinflammatory cytokine cascade against the pathogen invasion (101). The M1 macrophage-triggered proinflammatory cytokine cascade is suppressed by glucocorticoids and GR in basal conditions, but is activated by downregulated glucocorticoids and GR in an active infection (102). In this study, both ST and AT exhibited downregulated GO terms associated with cortisol and cortisol receptor functions and upregulated M1 macrophage polarization markers. These results suggest that activation of the proinflammatory cytokine cascade by M1 macrophage polarization is a general response for trout to fight pathogen invasion.

Compared to trout of AT, four KEGG pathways involved in steroid hormone biosynthesis and functions were downregulated in ST ( Figure 7 ). Steroid hormones, which include corticosteroids and sex steroids, play an important role in regulating homeostasis via modulating metabolism, immunomodulation, salt balance, water balance, and reproduction. The KEGG analysis revealed that the genes associated with the biosynthesis of corticosteroids and sex steroids (ko04913 and ko00140) were significantly downregulated, suggesting that V. anguillarum infection severely dysregulated the homeostasis of the steroid hormone network in trout of ST ( Figure 7E ). The dyshomeostasis of steroid hormone might disturb the bidirectional link between stress and immune responses. Thus, steroid hormone receptors (such as kidney mr subtypes) might fail to transduce the “alarm” of pathogen infection to immune systems in symptomatic trout. Based on previous studies, the well-orchestrated stress response can be divided into three phases: alarm, resistance, and exhaustion (103–105). Downregulated steroid hormone biosynthesis might indicate that the ST was in an exhaustion phase, which is consistent with the human study showing death may be associated with an exhausted adrenal cortex (106–108).

Previously published chapters in the book of Fish Physiology (Biology of Stress in Fish, Volume 35) showed that, with the perception of health challenges, the induction of neuroendocrine cascades serves as the primary responses. The secondary response to stressors includes the physiological adjustments of hydromineral balance and immune function (3, 7, 109). Hydromineral dysfunction is a typical stress response because the altered adrenaline, which is induced by stressors, can change the gill blood flow and gill permeability and dysregulate the cardiovascular and respiratory functions (7, 109, 110). Consistently, our studies showed significantly downregulated KEGG pathways associated with aldosterone-regulated salt and water balance (ko04960 and ko04978) ( Figure 7E ), indicating that ST trout show hydromineral dyshomeostasis. Previous studies in Chinook salmon (Oncorhynchus tshawytscha) showed the hydromineral balance is changed during euthanasia (111), which is consistent with our KEGG results. Based on this data, we propose positive feedback between severe infection and imminent death: (1) infection and its resulting stress response disturb the hydromineral homeostasis, thus resulting in a moribund condition. (2) The moribund condition further exacerbated the dyshomeostasis of hydromineral functions, leading to death.