Gill Transcriptomic Responses to Toxin-producing Alga Prymnesium parvum in Rainbow Trout

The gill of teleost fish is a multifunctional organ involved in many physiological processes, including protection of the mucosal gill surface against pathogens and other environmental antigens by the gill-associated lymphoid tissue (GIALT). Climate change associated phenomena, such as increasing frequency and magnitude of harmful algal blooms (HABs) put extra strain on gill function, contributing to enhanced fish mortality and fish kills. However, the molecular basis of the HAB-induced gill injury remains largely unknown due to the lack of high-throughput transcriptomic studies performed on teleost fish in laboratory conditions. We used juvenile rainbow trout (Oncorhynchus mykiss) to investigate the transcriptomic responses of the gill tissue to two (high and low) sublethal densities of the toxin-producing alga Prymnesium parvum, in relation to non-exposed control fish. The exposure time to P. parvum (4–5 h) was sufficient to identify three different phenotypic responses among the exposed fish, enabling us to focus on the common gill transcriptomic responses to P. parvum that were independent of dose and phenotype. The inspection of common differentially expressed genes (DEGs), canonical pathways, upstream regulators and downstream effects pointed towards P. parvum-induced inflammatory response and gill inflammation driven by alterations of Acute Phase Response Signalling, IL-6 Signalling, IL-10 Signalling, Role of PKR in Interferon Induction and Antiviral Response, IL-8 Signalling and IL-17 Signalling pathways. While we could not determine if the inferred gill inflammation was progressing or resolving, our study clearly suggests that P. parvum blooms may contribute to the serious gill disorders in fish. By providing insights into the gill transcriptomic responses to toxin-producing P. parvum in teleost fish, our research opens new avenues for investigating how to monitor and mitigate toxicity of HABs before they become lethal.


INTRODUCTION
The gill of teleost fish is a complex organ whose function goes beyond extracting oxygen from water and excreting carbon dioxide. Apart from the respiratory gas exchange, the gill plays a key role in osmotic and ionic regulation, acid-base balance, and excretion of nitrogenous waste (1). Because these processes are predominantly surface-dependent, the gill tissue consists of a highly complex system of branching vascular structures (primary and secondary lamellae) that are separated from the environment only by a thin layer of gill epithelium and mucosa (2,3). As a result, the epithelial surface area of the gill is typically larger than the total surface of skin (4). A substantial complexity of the teleost gill has also been demonstrated at the cellular resolution, with recent identification of 20 distinct cell clusters in the gill of Atlantic salmon (Salmo salar), using a single-nuclei RNA-seq approach (5).
Having the large epithelial surface of the gill system open to the external milieu comes with some disadvantages. Among them is the risk of mechanical injuries and gill abrasion, which may cause haemorrhage and contribute to gill inflammation (6). Furthermore, the large surface area of the gill may facilitate the uptake of toxic substances, including those occurring naturally (e.g., algal toxins, metal ions and ammonia) and the whole range of man-made pollutants (e.g., industrial chemicals, pesticides and microplastics) (7)(8)(9). Last, but not least, the gill surface provides major ports of entry for pathogens (via transepithelial transport) or site of interaction with other harmful invertebrates, such as water-born parasites and cnidarian jellyfish (10). To meet these challenges, the gill is strategically equipped with its own immune system, called the gill-associated lymphoid tissue (GIALT), thereby substantially contributing to overall fish health and survival (2,3). The GIALT consists of a panel of resident immune cells, including B and T cells, monocytes, macrophages, neutrophils, thrombocytes, dendritic-like cells, natural killer-like cells, eosinophilic granule cells, rodlet cells, and melanin-containing cells among the others (11,12).
Evidence is growing that climate change is putting extra strain on gill health and performance in many environments, including aquaculture (13). Rising water temperatures elevate metabolism of fish and their demands for oxygen (14), but at the same time decrease the oxygen content of water (15), intensifying the osmorespiratory conflict between the functionally large gill surface area (to promote respiratory gas exchange) and the need for the reduced gill surface area (to limit water and ion fluxes) (16,17). To avoid unfavourable changes in the environment, many marine organisms (including pathogens and parasites) are shifting their distributions as ocean temperatures warm (18). Among the hallmarks of the rapidly changing aquatic ecosystems are the occurrences of harmful algal blooms (HABs), whose increasing frequency, magnitude, and duration as well as migration poleward have been linked to ocean warming, marine heat waves, oxygen loss, eutrophication, and pollution (19,20). The climate-related range shifts and expansion in marine and freshwater environments put many species of fish at risk of being exposed to a variety of novel algal toxins, pathogens, and parasites, which were never before part of their own evolutionary history (21). While wild fish can relocate, farmed fish are heavily restricted in their movements, which makes their gills particularly vulnerable to the environmental insults. Thus, investigating the link between the climate change phenomena and gill health has wide implications not only for biodiversity and sustainability of aquatic ecosystems but also for food security at the global scale (13).
One of the HAB species that is particularly toxic to fish and other gill-breathing animals is the haptophyte alga Prymnesium parvum, undergoing a rapid range expansion in coastal and inland waters worldwide. Commonly referred to as the golden alga, P. parvum has been documented to kill~135 metric tons of farmed Atlantic salmon (Salmo salar) in Norway during a bloom in 2007 (22). In inland water bodies of Texas,~34 million fish valued at 13 million US dollars were lost due to P. parvum-related fish kills between 1981 and 2008 (23). As proposed decades ago (24), there is consensus that P. parvum acts on the gill tissue, but the exact mechanisms of its action are under investigation. Some research suggests that P. parvum cells need to be in direct physical contact with the gill surface to release harmful toxins, as a part of the toxin-assisted micropredation (25,26). Others propose that the presence of P. parvum toxins in the water is sufficient to kill fish (27,28). Although the complete suite of potential toxic compounds produced by P. parvum may not have been fully characterized, increasing evidence points towards prymnesins (a class of ladder-frame polyether phycotoxins) being responsible for enhanced fish mortality and massive fish kills (29)(30)(31). Yet, the experimental manipulations of fish with the use of toxin-producing P. parvum are relatively scarce.
Some studies use fish only for testing the acute toxicity of P. parvum cultures to establish their LC 50 values (mortality assay), with little focus on the fish themselves (28,(32)(33)(34)(35)(36)(37). Likewise, red cells extracted from fish blood have been used ex vivo to evaluate the haemolytic activity (a proxy for toxicity) of P. parvum isolates (38,39). In contrast, the sublethal doses of P. parvum were employed to investigate the effects of HABs on the whole-animal respiratory physiology (e.g., oxygen consumption, ventilation volume and frequency) in rainbow trout (Oncorhynchus mykiss) and European plaice (Pleuronectes platessa), with both studies pointing towards substantially reduced capacity of gills to extract oxygen from the environment (8,40). Both fish larvae (fathead minnow Pimephales promelas and zebrafish Danio rerio) as well as fish liver (Hepa-E1 and PLHC-1) and gill (G1B and RTgill-W1) cell lines were exposed to the sublethal doses of P. parvum to explore various aspects of the toxin-induced oxidative stress and antioxidant defence, including lipid peroxidation, oxidative DNA damage as well as gene expression and activity of antioxidant enzymes (41)(42)(43). It has also been demonstrated that the sublethal exposure of rainbow trout to P. parvum modulates susceptibility of fish to infectious agents such as viral haemorrhagic septicaemia virus (VHSV) (44). Surprisingly, little research has been done on the gill transcriptome, despite the gill being considered a key target tissue for P. parvum action.
In the current study, we used juvenile rainbow trout to investigate the transcriptomic responses of the gill tissue to two (high and low) sublethal densities of toxin-producing P. parvum, in relation to non-exposed control fish. The exposure time to P. parvum (4-5 h) was sufficient to identify three different phenotypic responses among the exposed fish, enabling us to focus on the common gill transcriptomic responses to P. parvum that were independent of dose and phenotype. Instead of profiling individual genes, we performed a microarray experiment to evaluate the gene expression changes at the level of whole tissue transcriptome (45). The functional analysis of the differentially expressed genes (DEGs) induced by P. parvum exposure included identification of associated canonical pathways, upstream regulators, and downstream effects.

Fish and Housing
All animal work was performed at the University of Aarhus (Denmark) in summer 2016. We used juvenile females of outbred rainbow trout (Oncorhynchus mykiss, body mass 10-15 g) from an all-female stock, hatched and reared in fresh water (temperature 10-11°C) under pathogen-free laboratory conditions. Fish were allocated to 8 experimental tanks (16 fish per 10-L tank, water temperature 15-17°C, water salinity 1.1%, water oxygen saturation > 90%) and then allowed to acclimate to the experimental conditions for 5 days. During the acclimation and experiment, fish were fed a standard commercial diet.

Exposure to Prymnesium parvum
The haptophyte P. parvum (Kalmar University Culture Collection, strain KAC 39) were cultured in F/2 medium (temperature 15°C, salinity 0.9%, photoperiod 14 h:10 h light: dark) as described previously (44). The exposure of fish to P. parvum was performed by adding exponentially growing cultures of P. parvum to the water to create environments with high (~4 x 10 4 cells per mL of water) and low (~1.5 x 10 4 cells per mL of water) densities of algae. During the transfer, the viability of the algal cells was confirmed by microscopy, after which the cells were pipetted into the water column and gently mixed to ensure their homogenous distribution within the tank. The P. parvum densities were chosen to mimic natural blooms (46), with both doses expected to have sublethal effects on fish (44). The effects of high and low doses of algae were evaluated using 3 replicate tanks. The remaining 2 tanks had no P. parvum added and served as negative controls.

Fish Sampling
Immediately after the exposure to P. parvum, fish were closely and continuously monitored for any behavioral and physiological abnormalities until the end of experiment (5 h post-exposure). This led to the identification of three different phenotypic responses that gradually developed among the exposed fish, referred to as high response, moderate response and low response (for details see Table 1). The high and moderate phenotypic responses were expressed by fish exposed to the high dose of P. parvum, while low phenotypic responses were observed in fish exposed to the low dose of P. parvum. Based on the dose and phenotype, fish were classified into 4 groups with high exposure/high response (HH), high exposure/ moderate response (HM), low exposure/low response (LL) and control group (C) with no exposure/no response ( Table 1).
Control fish (n = 16) were sampled 3-4 h after starting the experiment, followed by the experimental fish (n = 48), which were sampled 4-5 h after the exposure to P. parvum. Sampling was alternated between the HH, HM and LL groups for experimental balance. The experiment was terminated at 5 h post-exposure to reduce the risk of fish suffering from the algal toxicity (44).
Before sampling, fish were first euthanized by immersion in 0.01% benzocaine, then bled by removing the caudal fin and finally subjected to excision of gill arches. The collected gill tissue (2 separate gill arches per fish) was quickly blotted with absorbing paper to remove residual blood, followed by immediate transfer to RNAlater ® (Sigma-Aldrich, St. Louis, MO, United States). The gill samples (8 gill arches from 4 fish from the same group per tube) were kept at 4°C overnight for equilibration and then stored at −20°C prior to RNA extraction.

RNA Extraction and Sample Pooling
Total RNA extraction was performed on individual gill samples (n = 64), including the gill arch and full-length filaments. The RNA was isolated by homogenization of~100 mg of gill tissue in TRIzol ® Reagent (Ambion by Life Technologies, Carlsbad, CA, United States), using 3 mm tungsten carbide beads and a TissueLyser II Disruption System (Qiagen GmbH, Hilden,

Microarray Experiment
We used a custom designed Agilent oligonucleotide microarray platform Trout_imm_v1 (Agilent design ID: 028918) with 4 x 44 K probes per slide, developed for rainbow trout and validated using RTqPCR (47). The experiment consisted of 16 hybridisations (4 groups × 4 RNA pools), which were performed on 4 microarray slides in a semi-randomised order (each slide with 4 different groups). All RNA pools (n = 16) were subsampled to generate a common control, with contribution of the gill total RNA from 32-64 fish (4 groups × 4 RNA pools x 2-4 fish per pool). RNA amplification and labelling, followed by microarray hybridisation, scanning and feature extraction were performed as described previously (47)(48)(49). Briefly, antisense amplified RNA (aRNA) was generated from~2 mg total RNA per experimental or common control pool, using Amino Allyl MessageAmp ™ II aRNA Amplification Kit (Ambion by Life Technologies, Carlsbad, CA, USA). Pools of aRNA were then coupled with amine reactive Cy fluorescent dyes (Amersham ™ Cy ™ 3 and Cy ™ 5 Mono-reactive Dye Packs; GE Healthcare UK Limited, Little Chalfont, UK). All experimental samples were labelled with Cy3, while Cy5 was used to label the common control. The reaction products were purified using a DyeEx ® 2.0 Spin Kit (Qiagen GmbH, Hilden, Germany). Dye incorporation and post-labelling aRNA yield were quantified by spectrophotometry (NanoDrop Technologies, Wilmington, DE, USA). Each hybridisation was performed using 825 ng of Cy3-labelled experimental sample and 825 ng of Cy5-labelled common control. The aRNA was first fragmented and then hybridised at 65°C for 17 h in an Agilent hybridisation oven. Following hybridisation, slides were subjected to washing steps, after which they were air dried in the dark and scanned within 2 h. Scanning was carried out at 5 mm resolution on a Gene-Pix Personal 4100A scanner (Axon Instruments, Molecular Devices Corp., Sunnyvale, CA, USA), with the PMT values adjusted manually to ensure the mean intensity ratio of Cy3:Cy5 signal was close to 1. Agilent Feature Extraction Software (version 9.5.3) was used to identify features and to extract the raw intensity values for subsequent statistical analysis.

Microarray Data Analysis
Feature intensities were pre-processed and analysed for differential gene expression using the Bioconductor package limma (50) in R v3.5.0 (51). Briefly, loess normalisation was performed within each array to account for intensity-dependent variation in dye bias along with quantile normalisation used to stabilize experimental variances across arrays. Normalised data were subsequently filtered to remove control features and nonresponsive RNA targets with equal expression across 4 groups. Differential expression of RNA targets between groups was assessed using linear modelling, with the contrasts set up to compare each group of the fish exposed to P. parvum with the non-exposed control group (i.e., HH vs C, HM vs C and LL vs C). Multiple testing was accounted for by controlling the false discovery rate at 5% using the Benjamini-Hochberg procedure. The RNA targets with adjusted p-value < 0.05 and absolute Log 2 FC > 1 were considered as differentially expressed. All RNA targets meeting these criteria were 1) checked for the number of unique fish genes (with approximately half of the probes being redundant) and 2) mapped to human orthologs to enhance the functional analysis of gene expression ( Table 2, Supplementary  Tables 1-3).

Fish-to-Human Orthologs
Differentially expressed RNA targets were mapped to human orthologs to generate HGNC (HUGO Gene Nomenclature Committee) gene identifiers needed for functional analysis. This approach was demonstrated to enhance functional analysis of fish genes by providing access to well-annotated databases and tools for mammalian model organisms (mice, rats, and humans), despite differences between fish and mammals in gene function and molecular pathways (49,52,53). Mapping was done by aligning the microarray probe TABLE 2 | Results of differential gene expression analysis performed on the gill transcriptome of rainbow trout exposed to the toxin-producing alga Prymnesium parvum.

Group comparison
Number of differentially expressed genes 1 Supplementary RNA targets (representing fish mRNAs bound to the microarray probes) were mapped to human orthologs (BLASTX, E-value < 0.00001, top hit) to generate HGNC gene identifiers (HGNC, HUGO Gene Nomenclature Committee). Numbers refer to the total number of RNA targets (including replicate and redundant probes), and to the unique numbers of fish genes and human orthologs mapped to these targets. 2 for volcano plots of differential expression of RNA targets see Supplementary Figure 3.
The contrasts were set up to compare fish from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups vs control group (no exposure/no response). Genes were considered differentially expressed at adjusted p-value < 0.05 and absolute Log 2 FC > 1.

Functional Analysis of Gene Expression
Human orthologs of the differentially expressed fish genes were analysed using Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity). In total, 3 sets of differentially expressed genes (DEGs) were submitted to IPA (along with their Log 2 FC values), representing 3 comparisons between the fish groups: HH vs control (1436 DEGs), HM vs control (1069 DEGs) and LL vs control (567 DEGs) ( Table 2,  Supplementary Tables 1-3). These gene sets were analysed using the Ingenuity Knowledge Base (genes only) as a reference set, including in the analysis all species (mice, rats, and humans) as well as all tissues and cell lines (default settings).
The main focus of the functional analysis were canonical pathways, upstream regulators and downstream effects, the significance of which was based on the Benjamini-Hochberg (B-H) multiple testing correction p-value, with the overall activation/inhibition states predicted by the IPA z-score algorithm (z-score ≥ 2 predicts an increase in activity while zscore ≤ −2 predicts a decrease in activity). Gene ratios for each canonical pathway were calculated as the number of DEGs contributing to the pathway divided by the total number of genes that constitute the pathway and are present in the reference set.

Differentially Expressed Genes (DEGs)
The principal component analysis (PCA) of the gill transcriptome profiles in rainbow trout from HH (high exposure/high response), HM (high exposure/moderate response), LL (low exposure/low response) and C (control, no exposure/no response) groups demonstrated a partial overlap between groups of fish exposed to P. parvum (HH, HM and LL) and their clear separation from the non-exposed control group (C) ( Figure 1). This finding was reinforced by the differential gene expression analysis, which identified 1436, 1069 and 567 DEGs for comparisons of fish from HH vs control groups, HM vs control groups and LL vs control groups, respectively, with all DEGs referring to the human orthologs of fish genes ( Table 2,  Supplementary Tables 1-3, Supplementary Figure 3). The HH fish exposed to the high dose of P. parvum and displaying the strongest clinical response had the highest number of transcripts altered (1436 DEGs) followed by the similarly exposed but less affected HM fish (1069 DEGs), while the LL fish exposed to the low dose of P. parvum had approximately half of the number of transcripts altered (567 DEGs) relative to HH and HM groups. Inspection of DEGs from the three groups of fish exposed to P. parvum (HH, HM and LL) revealed a relatively large number of common DEGs (382), despite two different doses of P. parvum used in the experiment (high and low) and three different phenotypic responses of fish to the toxic algae (high, moderate and low) (Figure 2A, Supplementary Table 4). The majority of the common genes (375 of 382) showed the same direction of change (288 genes upregulated and 87 genes downregulated) across the three groups of fish, with the magnitude of change (Log 2 FC values) being highly correlated (correlation between HH and HM, r = 0.96, p < 0.001, Figure 2B; correlation between HH and LL, r = 0.91, p < 0.001, data not shown; correlation between HM and LL, r = 0.95, p < 0.001, data not shown). In addition, the mean expression of the 382 common genes was not significantly different between the groups, averaging 1.0 ± 1.7, 0.9 ± 1.6 and 0.8 ± 1.4 Log 2 FC (mean ± standard deviation) for HH, HM and LL groups, respectively (one-way ANOVA, p = 0.090). The common genes with the highest levels of upregulation and downregulation are presented in Table 3. . Ellipses indicate 95% confidence intervals. The percentage of total variance explained by PC1 and PC2 is given in parentheses.

Canonical Pathways
Based on the three sets of DEGs with 1436, 1069 and 567 transcripts, IPA identified 34, 95 and 33 canonical pathways that were significantly altered in the gill transcriptome of fish from HH, HM and LL groups, respectively, at B-H p-value < 0.001 (Supplementary Tables 5-7). Inspection of these pathways revealed a relatively large number of pathways (18) that were common between the three groups of fish ( Figure 3A), consistent with the large number of common DEGs (Figure 2A).  Figure 4) were both significantly activated, based on their IPA z-score > 2 for fish from HH, HM and LL groups. According to the predictions made by IPA, binding of cytokines to their specific receptors induced significant changes in intracellular and second messenger signalling pathways (represented by Glucocorticoid Receptor Signalling) and nuclear receptor signalling pathways (represented by PPAR Signalling and LXR/RXR Activation). The gills of fish exposed to P. parvum were likely affected by hypoxia, as indicated by nearly fully activated (z-scores from 0.9 to 1.9) HIF1a Signalling pathway, which is consistent with the increased production of gill mucus observed in all exposed fish ( Table 1 Table 4). The P. parvum-induced cellular stress and injury were also supported by alterations of IGF-1 Signalling and Prolactin Signalling pathways, both promoting cellular proliferation and differentiation of a variety of cell types in the attempts to repair the damaged tissue. Finally, 3 of the 18 common pathways were related to human inflammatory diseases (Role of Macrophages, Fibroblasts and Endothelial Cells in Rheumatoid Arthritis, Systemic Lupus Erythematosus in B Cell Signalling Pathway and Hepatic Fibrosis Signalling Pathway) and 3 others were associated with cancer (Tumour Microenvironment Pathway, Role of Tissue Factor in Cancer and PI3K/AKT Signalling).
Comparison of the gene ratios (the number of DEGs contributing to the pathway divided by the total number of genes that constitute the pathway) for the 18 common pathways across the three groups of fish demonstrated that the fish exposed to the high dose of P. parvum (HH and HM) had consistently similar ratios, and thus a similar number of DEGs contributing to a given pathway, while the gene ratios observed in the fish exposed to the low dose of the toxic algae were substantially lower ( Figure 4, Supplementary Table 8).

Upstream Regulators
IPA upstream regulator analysis identified 26, 28 and 14 upstream regulators (excluding chemicals and drugs) that could explain changes in the gene expression patterns observed in the gill tissue of HH, HM and LL fish, respectively, following the exposure to P. parvum (Supplementary Table 9). These upstream regulators were considered significant at Benjamini-Hochberg multiple testing correction p-value < 1E−15 and absolute z-score ≥ 2. Among the identified upstream regulators, 10 of them were shared between the three groups of fish ( Figure 3B).
The 10 common upstream regulators (z-scores from 2.7 to 5.5) included 4 cytokines (TNF, IL-1B, IL-6 and IFNG), with TNF predicted to be responsible for expression of 284, 223 and 143 target genes (DEGs) in the gill transcriptome of HH, HM and LL fish, respectively ( Figure 6, Supplementary Table 10). Among the common upstream regulators were also 3 growth factors (TGFB1, HGF and EGF), 2 complexes (PDGF BB and NFkB) and a peptidase (F2). One common upstream regulator (IL-6) was consistently upregulated across the three groups of fish (Log 2 FC values from 2.1 to 2.7), while others were not significantly altered in our experimental dataset (apart from F2 transcript downregulated in HM fish).

Downstream Effects
IPA downstream effect analysis predicted 28, 31 and 17 downstream effects based on the gene expression changes observed in the gill tissue of HH, HM and LL fish, respectively, following the exposure to P. parvum (Supplementary Table 11). The downstream effects were considered significant at Benjamini-Hochberg multiple testing correction p-value < 1E −15. A relatively large number of downstream effects (16) were common between the three groups of fish ( Figure 3C).
The predicted 16 common downstream effects covered cellular (7 effects), tissue and organ (5 effects) and organismal

DISCUSSION
Climate change is altering aquatic ecosystems worldwide, including patterns, distribution, and intensity of HABs in marine, brackish, and freshwater environments (19,20). These effects have been linked to 1) HAB species becoming more competitive relative to non-HAB species within plankton communities, 2) increased toxin production by toxic HAB species, and 3) HAB species reaching higher biomass due to changes in hydrology (55). Blooms present a challenge to fish health particularly in aquaculture, where containment of fish prevents avoidance behaviours. In anticipation of growing HAB problems, intensified research is needed to uncover molecular mechanisms by which toxin-producing algae affect wild and farmed fish, with the use of controlled laboratory settings, cultured HAB species and high-throughput technologies. Furthermore, there is also a need to perform some experiments at the level of whole animals. Although valuable information can be gained from experiments on fish cell lines, they lack physiological milieu and behavioural responses that are important in toxicological studies (56). We are the first to report on the gill transcriptomic responses to high and low sublethal doses of the toxin-producing P. parvum in commercially important rainbow trout, taking into account the fish phenotype.
Visual inspection of the fish following exposure to P. parvum allowed for identification of two distinct phenotypes among the fish originating from the same tanks and exposed to the high dose of the toxin-producing alga ( Table 1). Fish with a more advanced phenotype (HH) showed advanced lethargy, loss of balance and dark skin colour, while fish with a less advanced phenotype (HM) showed only mild lethargy, with both phenotypes having increased respiratory effort and increased production of gill mucus. Our interpretation of these differences is that fish with the more advanced phenotype were likely metabolically more active than fish with the less advanced phenotype. Fish with higher metabolic rate and thus higher oxygen uptake have been demonstrated to 1) increase water flow over the gill by adjusting the volume and frequency of buccal pumping, 2) increase blood flow inside the gill to alter the perfusion levels of lamellae, and 3) initiating remodelling of gill tissue towards more protruded lamellae (16,57). All these adjustments make the gill tissue potentially more available for toxin uptake (58). Salmonids are particularly known for their substantial intraspecific variation in metabolic rate, which has been linked to the differences in individual's behaviour (e.g., dominance, foraging or stress avoidance) and performance (e.g., growth) (59)(60)(61). With our single sampling approach at 4-5 h TABLE 3 | Top common genes altered in the gill transcriptome of rainbow trout from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups in relation to C (control, no exposure/no response) group.

Gene (HGNC symbol and name)
Gene post-exposure, we were unable to establish whether the phenotypic differences persisted or converged with time. In contrast to the high dose exposure, no phenotypic differences were observed among the fish exposed to the low dose of P. parvum (Table 1). Similarly, we cannot exclude the possibility that different metabolic phenotypes might manifest themselves at the lower dose over a more prolonged period of time. Despite the differences in the phenotype, fish exposed to the high dose of P. parvum (HH and HM) had a relatively similar number of genes altered (1436 and 1069) in the gill tissue, while the fish exposed to approximately half of the dose of P. parvum (LL) had approximately half of the genes altered (567) relative to HH and HM groups (Figure 2A, Table 2). The positive relationship between the dose of P. parvum and the number of DEGs points towards an overall sensitivity of the gill transcriptome to toxic algae dose. Dose-dependent transcriptomic effects have also been demonstrated in the gut of Atlantic salmon exposed to plant proteins (49) and recently in human nasal airway epithelium cultures exposed to air pollution particulate matter < 2.5 mm (PM 2.5 ) (62). Furthermore, fish from HH and HM groups had not only twice more DEGs identified in their gills relative to LL fish, but also twice as many upstream regulators (26, 28 and 14, respectively) ( Figure 3B) and twice as many downstream effects (28, 31 and 17, respectively) ( Figure 3C). In contrast, fish exposed to the high dose of P. parvum differed in the number of canonical pathways, with only 34 pathways identified in the gill of fish showing the more advanced phenotype (HH) and 95 pathways identified in the gill of fish with the less advanced phenotype (HM) ( Figure 3A). Fewer DEGs contributing to more pathways (1069 and 95, respectively) may reflect more coordinated patterns of gene expression in HM fish, while more DEGs contributing to fewer pathways (1436 and 34, respectively) may suggest the onset of dysregulated gene expression patterns in HH fish. The transition from coordinated to dysregulated gene expression patterns has been linked to tissue and organ malfunction, leading to disease and ultimately death (63,64). Thus, the high dose of P. parvum used in our experiment may be considered as the borderline of what the fish could handle, with the HH phenotype representing individuals that were pushed beyond their physiological control and the HM phenotype reflecting fish that were still able to cope with the algal exposure by mobilizing a broader range of defence mechanisms initiated at the level of gill transcriptome. The fact that the LL phenotype shared many of the altered canonical pathways as well as upstream regulators and downstream effects with the HH and HM fish may reflect a more general and balanced response of the LL fish to the P. parvum exposure. The mucus production and release from gills of all exposed fish likely represent a physical defence mechanism for clearing the gills from algae/algal toxin. Accordingly, few algae were found in direct contact with the gills while high numbers of the motile algal cells were found trapped in the secreted mucus (Supplementary Figure 5).
The gill transcriptomic responses to two different doses of P. parvum in fish expressing three different phenotypes were to a relatively large extent overlapping. The total number of unique DEGs for HH, HM and LL fish was 1854, with 382 (20.6%) of them being common ( Figure 2A) and having highly correlated gene expression patterns ( Figure 2B). Likewise, the three groups of fish shared 18 (17.6%) of 102 canonical pathways ( Figure 3A), 10 (30.3%) of 33 upstream regulators ( Figure 3B) and 16 (50.0%) of 32 downstream effects ( Figure 3C), according to the functional analysis of gene expression performed by IPA. These overlaps are important to investigate because they likely represent a core transcriptomic response characterizing the gill's attempt to maintain homeostasis when responding to P. parvum, regardless of the algal dose and fish phenotype. A similar approach has been recently used to gain insights into the mechanisms of multifactorial gill disease in Atlantic salmon (53), adaptation of the foodborne pathogen Listeria monocytogenes to desiccation (65) and resilience to heat stress in laying hens (66).
DEGs identified as common for HH, HM and LL fish (382 transcripts in total) represent a diverse group of genes, encoding 23 cell surface receptors (including 7 G-protein coupled receptors), 30 transporters, 2 ion channels, 11 cytokines (including CSF3 and CXCL2 with the highest Log 2 FC values of all DEGs), 128 enzymes (including 17 peptidases, 13 kinases and 10 phosphatases), 3 growth factors, 54 transcription factors (including 2 ligand-dependent nuclear receptors) and 4 translation regulators among others ( Table 3, Supplementary  Table 4). Such diversity is consistent with the cellular and molecular mechanisms of action by which most toxins disrupt eukaryotic cells, first interacting with the host cell surface and then compromising the host intracellular processes associated FIGURE 4 | Common canonical pathways altered by P. parvum exposure in the gill transcriptome of rainbow trout from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups in relation to C (control, no exposure/no response) group. Pathways were identified as significantly altered by Ingenuity Pathway Analysis (IPA) at Benjamini-Hochberg multiple testing correction p-value < 0.001. Gene ratios refer to the number of DEGs contributing to each pathway (our experimental dataset) divided by the total number of genes that constitute the pathway (Ingenuity Knowledge Base). For details and list of corresponding genes see Supplementary Tables 5-8. with energy metabolism, cytoskeleton stability, gene expression, posttranslational modifications, motility, secretion, cell division or other more specific functions (67,68). Some toxins are known to specifically disrupt the ionic equilibrium maintained by the cell membrane barrier, which can be done either directly (by forming pores or causing an enzymatic degradation of the lipid bilayer) or indirectlyby acting on ion pumps or ion-gated channels that are responsible for maintaining the ion gradients (69). A classic example of the pore-forming toxins is pardaxin secreted by the finless sole (Pardachirus marmoratus), which is used as a defensive mechanism against predators including sharks via targeting their gills (70). Other toxins interfere with the signalling cascades of the eukaryotic cells by acting on a variety of cell membrane targets, including ligand-gated ionotropic receptors, G-protein coupled receptors, tyrosine kinase receptors, integrin receptors and certain lipid species present in the bilayer plasma membrane of the cell. Examples here are domoic acid produced by the red alga Chondria armata that acts on ligand-gated ionotropic glutamate receptors (71) and okadaic acid produced by the marine algae Halichondria okadai and Halichondria melandocia that affects protein kinases and phosphatases (69). In our study, some of the DEGs identified as common for HH, HM and LL fish (e.g., ion channel and transporter transcripts) clearly point towards P. parvum disrupting the ionic regulation of the gill tissue, while other common DEGs (e.g., kinase and phosphatase transcripts) implicate the dysregulation of the signalling cascades in the gill cells. The wide range of gill transcriptomic responses to P. parvum likely reflect the heterogeneity of effects of the algal cells themselves and their released toxic compounds along with counteractive reactions in the gills such as mucus production and secretion, tissue repair and innate immune reactions. Signalling pathway in the gill transcriptome of rainbow trout exposed to the high dose of P. parvum with high phenotypic response (HH group) in relation to the non-exposed control fish (C group). The pathway was identified as significant (Benjamini-Hochberg multiple testing correction p-value < 0.001) and activated (z-score ≥ 2) by Ingenuity Pathway Analysis (IPA). Among 126 genes that constitute the pathway, 26 were significantly altered by P. parvum exposure (yielding the gene ratio of 0.206), including 23 genes upregulated (in red) and 3 genes downregulated (in green). The upregulated genes were CEBPB, CCAAT enhancer binding protein beta; FOS, Fos proto-oncogene; AP-1 transcription factor subunit; HRAS, HRas proto-oncogene; GTPase; IL1B, interleukin 1 beta; IL1R2, interleukin 1 receptor type 2; IL1RAPL1, interleukin 1 receptor accessory protein like 1; IL1RAPL2, interleukin 1 receptor accessory protein like 2; IL6R, interleukin 6 receptor; IL8/CXCL8, C-X-C motif chemokine ligand 8; JUN, Jun proto-oncogene; AP-1 transcription factor subunit; LBP, lipopolysaccharide binding protein; MAP4K4, mitogen-activated protein kinase kinase kinase kinase 4; MCL1, MCL1 apoptosis regulator; BCL2 family member; NFKBIA, NFKB inhibitor alpha; PIK3R5, phosphoinositide-3-kinase regulatory subunit 5; RAP2B, RAP2B; member of RAS oncogene family; RASD1, ras related dexamethasone induced 1; SOCS1, suppressor of cytokine signalling 1; SOCS3, suppressor of cytokine signalling 3; SOS1, SOS RasRac guanine nucleotide exchange factor 1; SRF, serum response factor; TRAF6, TNF receptor associated factor 6 and VEGFA, vascular endothelial growth factor A. The downregulated genes were IKBKAP/ELP1, elongator acetyltransferase complex subunit 1; MAPK13, mitogen-activated protein kinase 13 and TNFAIP6, TNF alpha induced protein 6. For details see Supplementary Tables 1, 5 and 8.
Inspection of 18 canonical pathways identified as commonly altered in HH, HM and LL fish suggests that exposure to P. parvum promotes a strong inflammatory response within the gill tissue and thus induces gill inflammation, as evidenced by Acute Phase Response Signalling, IL-6 Signalling, IL-8 Signalling and IL-17 Signalling pathways (Figure 4, Supplementary Table 8). Whether these changes were associated with the resident immune cells (GIALT) or the recruitment of new immune cells to the gill tissue remains unknown. The acute phase response is a rapid inflammatory reaction that provides protection against microorganisms by non-specific defence mechanisms (72). It is typically manifested by an increase in inflammatory factors (such as pro-inflammatory cytokines) and a change in concentration of several plasma proteins (the acute phase proteins) that can become measurable as early as 4-5 h after a single inflammatory stimulus. One of the key regulators of the acute phase response is IL-6 (a pleiotropic cytokine with roles in inflammation, immune response, haematopoiesis, and the endocrine and nervous systems), which signals through JAK-STAT and RAS-MAPK pathways to regulate transcription of target genes (73,74). In our study, 26 (HH), 23 (HM) and 17 (LL) of the 126 genes that constitute the IL-6 Signalling pathway were significantly altered, with the overall direction of change pointing towards a significant activation of this pathway in the three groups of fish ( Figure 5, Supplementary Tables 5-8). Among the target genes of IL-6 is IL-8, a member of the C-X-C family of chemokines that plays a central role in inflammation, angiogenesis, and tumour growth, with the IL-8 receptors expressed on several cell types like neutrophils, monocytes, endothelial cells, and tumour cells (75). Signalling by IL-8 or other ligands of the IL-8 receptors can trigger inflammation in cells like neutrophils leading to chemotaxis, the respiratory burst, granule release, and increased cell adhesion (76). All fish exposed to P. parvum showed a robust upregulation of IL-8 transcription (with Log 2 FC values from 1.3 to 2.3), along with highly significant alteration of IL-8 Signalling consistent with activation (z-scores from 1.8 to 3.4) (Supplementary Tables 4 and 8). Another highly versatile cytokine with a strong pro-FIGURE 6 | Common upstream regulators of gene expression changes in the gills of rainbow trout exposed to P. parvum from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups in relation to C (control, no exposure/no response) group, predicted by Ingenuity Pathway Analysis (IPA). Upstream regulators were considered significant at Benjamini-Hochberg multiple testing correction p-value < 1E−15 and absolute z-score ≥ 2, excluding chemicals and drugs (for details see Supplementary Tables 9 and 10).
inflammatory profile is IL-17, which promotes expansion and recruitment of innate immune cells such as neutrophils and stimulates production of beta-defensins and other anti-microbial peptides (77). Because IL-17 is primarily secreted by T cells, it plays a central role in integrating adaptive and innate immune responses (78). The IL-17 Signalling pathway was fully activated in HH, HM and LL fish (z-scores > 3.0), with 28 (HH), 28 (HM) and 24 (LL) of the 187 genes that make up the pathway being significantly altered, including CSF3 with the highest Log 2 FC values of all 382 common DEGs ( Table 3, Supplementary  Figure 4). Although the interpretation of our results by IPA in the context of human inflammatory diseases may be irrelevant, the presence of highly affected pathways related to autoimmune diseases (i.e., Role of Macrophages, Fibroblasts and Endothelial Cells in Rheumatoid Arthritis, Systemic Lupus Erythematosus in B Cell Signalling Pathway and Hepatic Fibrosis Signalling Pathway) highlight the severity and scope of the gill inflammatory response towards P. parvum (Figure 4). A similar argument can be made for the three cancer-related pathways (i.e., Tumour Microenvironment Pathway, Role of Tissue Factor in Cancer and PI3K/AKT Signalling). The exposure to P. parvum does not induce cancer, but some of the key features associated with cancer (such as local inflammation, cell proliferation, angiogenesis, and hypoxia) are also the functions of the DEGs identified in the current study (79).
If not resolved, the acute inflammatory response may contribute to tissue injury and chronic inflammation, an underlying cause of human chronic inflammatory diseases such as arthritis, diabetes, metabolic syndrome, neurodegenerative diseases, asthma, allergy, tissue fibrosis, certain types of cancer, cardiovascular and periodontal diseases (80). Resolution of inflammation involves highly coordinated actions of various immune and non-immune cells and pathways to first clear damaged cells and proinflammatory cytokines (by apoptosis and efferocytosis) and then FIGURE 7 | Common downstream effects of gene expression changes in the gills of rainbow trout exposed to P. parvum from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups in relation to C (control, no exposure/no response) group, predicted by Ingenuity Pathway Analysis (IPA). Downstream effects were considered significant at Benjamini-Hochberg multiple testing correction p-value < 1E−15 (for details see Supplementary Tables 11 and 12).
at least partially reconstitute tissue integrity and function (81). In our study, the attempts to resolve inflammation in the gill tissue following exposure to P. parvum are evidenced by Glucocorticoid Receptor Signalling and IL-10 Signalling pathways, commonly altered in HH, HM and LL fish ( Figure 5, Supplementary  Tables 5-8). Glucocorticoids are a major subclass of steroid hormones that produce their effects on responsive cells by activating the glucocorticoid receptor to modulate the transcription of target genes, thereby regulating a large number of immune, metabolic, cardiovascular and behavioural functions (82). Despite the ability of glucocorticoids to induce transcription of antiinflammatory genes, the main anti-inflammatory effects of glucocorticoids are through repression of pro-inflammatory cytokine genes (83). These anti-inflammatory actions are also complemented by the ability of glucocorticoids to induce apoptosis of immune cells, including monocytes and T cells. Among the 462 genes that constitute the Glucocorticoid Receptor Signalling pathway, 65, 60 and 37 genes were altered in HH, HM and LL fish, respectively, including IL-10 transcription with Log 2 FC values from 2.3 to 2.5 (Supplementary Table 4). IL-10 is a cytokine with diverse effects on hematopoietic cells, which regulates the growth and differentiation of B, T, natural killer and dendritic cells (84). One of its best-known functions is to limit and terminate the inflammatory response. The mechanism behind this action is signalling through the IL-10 receptor that causes inhibition of JAK/STAT-dependent signalling, leading to the transcription inhibition of pro-inflammatory genes like IL-1 and TNF-a (85). In our study, both pro-and anti-inflammatory cytokine transcripts (Supplementary Table 4) and signalling pathways (Supplementary Table 8) were altered, but because we used a single sampling approach (with no temporal resolution of gene expression data), it is impossible to conclude whether the inflammation of gill tissue following the 4-5 h exposure to P. parvum was progressing or resolving.
The presence of gill inflammation inferred from the gene expression changes in the fish exposed to P. parvum is also supported by the results of the IPA upstream regulator analysis ( Figure 6, Supplementary Table 10). Specifically, 6 of the 10 common upstream regulators (all identified as activated) are known for their pro-inflammatory action, including TNF, TGFB1, IL-1B, IL-6, IFNG and NFkB complex, with TGFB1 and NFkB complex also contributing to mounting the anti-inflammatory response (73,(86)(87)(88)(89). The remaining 4 common upstream regulators have been linked to inhibition of inflammatory responses (PDGF BB), inflammation-induced coagulation (F2) and the growth, proliferation and differentiation of numerous cell types that are necessary for tissue repair and regeneration (HGF and EGF) (90)(91)(92). It is important to clarify here that the IPA upstream analysis does not take into account the gene expression observed for the predicted upstream regulator itself, because the gene expression for the upstream regulator may not differ between the experimental and control groups. In our study, only one upstream regulator transcript (IL-6) was commonly upregulated in the fish exposed to P. parvum, which suggests that the remaining 9 common upstream regulators were likely activated by other means rather than by increased gene expression.
Finally, the results of the IPA downstream effect analysis are consistent with gill inflammation inferred from the transcriptomic data, as evidenced by the common whole-animal downstream effects such as Inflammatory Response, Inflammatory Disease, and Infectious Diseases (Figure 7, Supplementary Table 12). These predictions are in good agreement with the secretion of gill mucus observed in all fish exposed to P. parvum ( Table 1, Supplementary Figure 5). Furthermore, gill inflammation in salmonids is characterised by excessive cellular death by apoptosis, followed by proliferation and migration of epithelial cells to replace the damaged mucosal surface and to provide a defensive barrier that isolates the gill tissue from the environmental insults (93,94). These processes were in our study represented by a number of predicted cellular downstream effects that were common for all fish exposed to P. parvum, including Cell Death and Survival, Cellular Movement, Cellular Development, Cellular Growth and Proliferation, Cellular Function and Maintenance, Cell-To-Cell Signalling and Interaction, and Immune Cell Trafficking. The cellular effects pointing towards the inflammatory processes were further supported by the tissue and organ downstream effects such as Tissue Morphology, Haematological System Development and Function, Lymphoid Tissue Structure and Development, Connective Tissue Disorders, and Tissue Development, providing a strong support for harmful algal gill pathology in rainbow trout exposed to P. parvum (95,96).
Our results have important implications for understanding the molecular basis of HAB-induced gill injury, which has remained largely unknown due to the lack of high-throughput transcriptomic studies performed on teleost fish under controlled laboratory conditions. Firstly, we demonstrated that the gill transcriptomic responses are sensitive and proportional to the sublethal concentrations of P. parvum. However, the whole-animal phenotypic and transcriptomic responses to near-lethal algal levels may vary, depending on whether the individual fish are able to maintain physiological homeostasis as observed for HH vs HM phenotypes. Thus, the phenotype variabilities for a given algal exposure should be taken into account when evaluating the effects of sublethal doses of HABs on fish. Secondly, we identified more coordinated patterns of gill gene expression in fish with relatively healthy phenotypes (LL and HM fish), and less coordinated (more dysregulated) patterns of gill gene expression in HH fish with advanced clinical presentation. The distinction between the coordinated and dysregulated gene expression patterns was based on the number of DEGs in relation to the number of canonical pathways identified as significantly altered, highlighting the importance of functional analysis of gene expression in aquatic toxicology studies. Thirdly, we demonstrated that the transcriptomic changes in the gill tissue of fish exposed to different doses of P. parvum (high and low) and expressing three different phenotypes (high, moderate and low responses) showed a high degree of overlap, pointing towards common molecular defence mechanisms against P. parvum. Finally, the inspection of the common transcriptomic features (i.e., DEGs, canonical pathways, upstream regulators, and downstream effects) in the gills of HH, HM and LL fish led us to the conclusion that P. parvum caused a strong inflammatory response, associated with gill inflammation and the attempts to calm the inflammation down, as well as potential downstream tissue damage and defects in gill function. While we could not determine if the inferred gill inflammation was progressing or resolving, our study clearly suggests that even HABs with the sublethal levels of toxicity may contribute to the serious gill disorders in rainbow trout. By providing insights into the gill transcriptomic responses to toxinproducing P. parvum in teleost fish, our research opens new avenues for investigating how to monitor and mitigate toxicity of HABs before they become lethal.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available in the ArrayExpress repository (http://www.ebi.ac.uk/arrayexpress/) under accession numbers A-MEXP-2315 (Trout_imm_v1 array design) and E-MTAB-10541 (microarray raw data).

AUTHOR CONTRIBUTIONS
NL, PH, NA, and SM conceived and designed the study, with the input from AB and DF. DS and NL performed the animal work. SM oversaw molecular work, done by MC and EK. MC carried out differential expression analysis in R. EK performed functional analysis of gene expression and deposited raw microarray data at ArrayExpress. EK, MC, and SM drafted and prepared the manuscript for publication. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank two reviewers for their comments on the earlier draft of the article.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021. 794593/full#supplementary-material Supplementary Figure 1 | 3D PCA (principal component analysis) plot of the gill transcriptome profiles in rainbow trout from HH (high exposure/high response), HM (high exposure/moderate response), LL (low exposure/low response) and C (control, no exposure/no response) groups. Each circle refers to all genes from one microarray hybridisation performed on 2-4 fish, with 4 hybridisations (biological replicates) per group (16 hybridisations in total). The percentage of total variance explained by PC1, PC2 and PC3 is given in parentheses.
Supplementary Figure 2 | Principal component analysis (PCA) of the gill transcriptome profiles in rainbow trout exposed to P. parvum (HH, HM and LL groups) and non-exposed control fish. Each circle refers to all genes from one microarray hybridisation performed on 2-4 fish, with 12 hybridisations representing the exposed fish and 4 hybridisations representing non-exposed fish. Ellipses indicate 95% confidence intervals. The percentage of total variance explained by PC1 and PC2 is given in parentheses.
Supplementary Figure 3 | Volcano plots of differential expression of RNA targets in the gill transcriptome of rainbow trout from HH (high exposure/high response), HM (high exposure/moderate response) and LL (low exposure/low response) groups in relation to C (control, no exposure/no response) group. Genes were considered differentially expressed at adjusted p-value < 0.05 and absolute Log 2 FC > 1. Upregulated genes are in green, while downregulated genes are in red.
Supplementary Figure 5 | Mucus secretion in the gills of rainbow trout exposed to the toxin-producing alga Prymnesium parvum. When exposed to UV-light, the two chloroplasts of P. parvum become fluorescent (A, B). Only few algal cells tended to reach the gill lamellae (C, D), while a relatively large number of algal cells was trapped in the secreted mucus (E, F). Photos were taken with a Leica DM4008B-M microscope equipped with epifluorescence, using either transmitted white light (A, C, E) or UV-light of the same field (B, D, F). Bars correspond to 10 µM (A, B) or 25 µM (C-F).