Original Research ARTICLE
Integration of Transcriptome, Gross Morphology and Histopathology in the Gill of Sea Farmed Atlantic Salmon (Salmo salar): Lessons From Multi-Site Sampling
- 1School of Biological Sciences, Institute of Biological and Environmental Sciences, University of Aberdeen, Aberdeen, United Kingdom
- 2Fish Health and Welfare, Marine Scotland Science, Aberdeen, United Kingdom
- 3Centre for Genome-Enabled Biology and Medicine, University of Aberdeen, Aberdeen, United Kingdom
- 4BioMar AS, Trondheim, Norway
- 5Scottish Sea Farms, Stirling, United Kingdom
The gill of teleost fish is a multifunctional organ involved in many physiological processes such as gas exchange, osmotic and ionic regulation, acid-base balance and excretion of nitrogenous waste. Due to its extensive interface with the environment, the gill plays a key role as a primary mucosal defense tissue against pathogens, as manifested by the presence of the gill-associated lymphoid tissue (GIALT). In recent years, the prevalence of multifactorial gill pathologies has increased significantly, causing substantial losses in Atlantic salmon aquaculture. The transition from healthy to unhealthy gill phenotypes and the progression of multifactorial gill pathologies, such as proliferative gill disease (PGD), proliferative gill inflammation (PGI) and complex gill disorder (CGD), are commonly characterized by epithelial hyperplasia, lamellar fusion and inflammation. Routine monitoring for PGD relies on visual inspection and non-invasive scoring of the gill tissue (gross morphology), coupled with histopathological examination of gill sections. To explore the underlying molecular events that are associated with the progression of PGD, we sampled Atlantic salmon from three different marine production sites in Scotland and examined the gill tissue at three different levels of organization: gross morphology with the use of PGD scores (macroscopic examination), whole transcriptome (gene expression by RNA-seq) and histopathology (microscopic examination). Our results strongly suggested that the changes in PGD scores of the gill tissue were not associated with the changes in gene expression or histopathology. In contrast, integration of the gill RNA-seq data with the gill histopathology enabled us to identify common gene expression patterns associated with multifactorial gill disease, independently from the origin of samples. We demonstrated that the gene expression patterns associated with multifactorial gill disease were dominated by two processes: a range of immune responses driven by pro-inflammatory cytokines and the events associated with tissue damage and repair, driven by caspases and angiogenin.
The gill of teleost fish is a multifunctional organ involved in many physiological processes such as gas exchange, osmotic and ionic regulation, acid-base balance and excretion of nitrogenous waste (Evans et al., 2005). To facilitate these functions, the gill tissue has evolved into a highly complex system of branching vascular structures (lamellae), separated from the external milieu only by a thin layer of gill epithelium and mucosa (Koppang et al., 2015; Salinas, 2015). The densely packed lamellar structure of the gill is highly advantageous because it provides a large surface area for oxygen transfer, amounting to approximately 0.1–0.4 m2 of lamellar surface per kg of body mass (Maina, 2011; Park et al., 2014). However, having such an extensive interface with the environment comes at a high price (reviewed in Nilsson, 2007). Firstly, the large respiratory surface area may contribute to the increased water and ion fluxes that need to be counteracted by energetically expensive ion pumping. That so-called osmorespiratory compromise, which is a trade-off between high gill permeability (to promote respiratory gas exchange) and low gill permeability (to limit unfavorable water and ion fluxes), has been demonstrated in freshwater fish during exercise and in hypoxia-intolerant species (e.g., salmonids) during the periods of exposure to low-oxygen water (Onukwufor and Wood, 2018; Morgenroth et al., 2019). Secondly, the highly complex vasculature of the gill is prone to mechanical injuries, increasing the risk of hemorrhage. Because the gill is the only organ in fish that receives the entire cardiac output, large or frequent gill hemorrhages can be fatal. Thirdly, the large lamellar surface area may facilitate the uptake of toxic substances, both those that occur naturally (such as ammonia, algal toxins and metal ions) as well as man-made pollutants (e.g., pesticides, herbicides, fertilizers, detergents, industrial chemicals and microplastics) (McIntyre et al., 2018; Kintner and Brierley, 2019). Finally, the close proximity of blood and environmental water across the large respiratory surface area provides (1) major ports of entry for pathogens via transepithelial transport and (2) attachment sites for water-born parasites. To combat pathogens and parasites, the gill is equipped with the gill-associated lymphoid tissue (GIALT), adding another dimension to the importance of this organ to fish functioning and survival (Koppang et al., 2015; Salinas, 2015; Xu et al., 2020).
To avoid ‘the large interface’ penalties, fish constantly adjust the functional size of the gill tissue to match their current oxygen demands (Nilsson, 2007). When faced with changes in their need for oxygen uptake (due to changes in metabolic rate or water oxygen content), fish are known to (1) regulate the water flow over the gill by adjusting the volume and frequency of buccal pumping and/or (2) regulate the blood flow inside the gill to alter the perfusion levels of lamellae. New evidence suggests that many teleost species may also undergo extensive gill tissue remodeling, during which lamellae are either completely embedded in a cell mass (to minimize the respiratory surface area) or fully protruded (to maximize the functional size of the tissue), with many potential intermediate stages between these extremes (Sollid et al., 2003; Anttila et al., 2015). Whether these regulatory and compensatory mechanisms aiming to optimize gill performance are sufficient for fish to cope with climate change, especially in the areas of high human impacts, remains unknown.
Fish are predominantly ectothermic and any increases in water temperature resulting from climate change will increase their metabolic rate and thus demands for oxygen. The Q10 for metabolism in fish (the fold increase in metabolic rate with a 10% increase in water temperature) varies from 1.5 to 2 (White et al., 2006). But the elevated oxygen demand is not the only challenge that the gill tissue is facing. Warmer water holds less oxygen and rising global temperatures alter oceanic circulation, which altogether may decrease the oxygen supply available for fish to support their increased metabolism (McBryan et al., 2013). It has been estimated that the open ocean lost ∼2% (77 billion metric tons) of its oxygen over the past 50 years (reviewed in Breitburg et al., 2018). The extent of deoxygenation of coastal waters may even be greater, due to eutrophication from agriculture and sewage runoff (Diaz and Rosenberg, 2008). Furthermore, environmental changes, including temperature increases, have been linked to enhanced expression of marine infectious diseases (Burge et al., 2014). Many marine organisms, including marine pathogens, are shifting their distributions poleward or to deeper waters as ocean temperatures warm (Nye et al., 2009; Pinsky et al., 2013). As a result, many species of fish are at risk of being exposed to novel pathogens and parasites that were not part of their evolutionary history (LeBlanc et al., 2019). Climate change may also affect the fish-pathogen interactions that are already well established, for example by altering pathogen virulence, overall pathogenicity and/or host immune response (Benedicenti et al., 2019; Laurin et al., 2019). Given the fact that in fish most environmental stressors converge at the level of the gill, performance of this tissue and its capacity to accommodate environmental change are central to understand short- and long-term impacts of warming and hypoxia on fish species (Akbarzadeh et al., 2018; McCormick and Regish, 2018; Houde et al., 2019). Monitoring gill health is especially important in farmed marine fish, which are restricted to coastal waters and have limited microhabitat choice.
In recent years, the prevalence of gill damage and disease in farmed marine fish has significantly increased, leading to substantial losses in the Scottish and global Atlantic salmon (Salmo salar) aquaculture industry (Young et al., 2007; Shinn et al., 2015; Nowak and Archibald, 2018). Although gill-related mortality varies from site to site and between years, July to December in Northern Hemisphere is considered the higher risk period because the sea during these months is warmer. With the projected increase in sea water temperature associated with climate change, gill pathologies are now regarded as the greatest challenge to the farmed salmon sector due to direct and indirect impacts on productivity. Losses may occur through direct mortalities, poor growth rates linked to decreased feed efficiency and because of the increased risk of co-infection (Rodger, 2007; Downes et al., 2018). The transition from healthy to unhealthy gill tissue happens through the progression of various gill pathologies such as amoebic gill disease (AGD), proliferative gill disease (PGD), proliferative gill inflammation (PGI) and complex gill disorder (CGD), which are commonly characterized by epithelial hyperplasia, lamellar fusion and inflammation (Steinum et al., 2010; Matthews et al., 2013; Mitchell et al., 2013; Herrero et al., 2018; Gjessing et al., 2019). The diagnostic of gill pathologies is often difficult because single-cause diseases (such as AGD) are compounded by multi-cause (multifactorial) disorders such as PGD, PGI or CGD, with complex etiology and potentially synergistic effects of co-infections (Gjessing et al., 2017; Gunnarsson et al., 2017; English et al., 2019). Due to overlapping symptoms of multifactorial gill inflammation, the terms PGD, PGI and CGD are often used interchangeably, with the first two historically referring to salmon farmed in Scotland (PGD) and Norway (PGI), while the last one highlighting the growing complexity of gill diseases and encompassing the syndromes referred to as PGD and PGI (Rozas-Serri, 2019). The multifactorial gill diseases are associated with a number of causative agents, including phytoplankton, zooplankton, viruses, fungi, bacteria and larger parasites (Gjessing et al., 2017; Herrero et al., 2018), with climate change contributing to emerging new aquaculture pathogens (Bayliss et al., 2017). Furthermore, gill diseases not only directly affect fish health and performance, but also creates a challenge for farmed salmon producers as they are required to treat compromised fish, increasing the risk of mechanical injuries to the gill vasculature (Powell et al., 2005). To meet the challenge of salmon gill health, the sector needs non-invasive monitoring tools for (1) early detection of gill pathologies, (2) diagnostics that is sufficient to establish appropriate treatment and (3) evaluation of effectiveness of treatment and recovery of gill function. It is important that such initiatives target multifactorial diseases, as they are becoming progressively more prevalent (Rahel and Olden, 2008).
The current practice is to regularly check sea farmed salmon for any signs of gill inflammation by in situ gross morphology examination of gill arches, using the PGD scoring system from 0 (no macroscopic pathology) to 5 (severe macroscopic pathology) (Gill Health Initiative, 2015; Bloecher et al., 2018). This is a non-invasive procedure that requires only light anesthesia and can be performed on a weekly basis in fish from across a facility. However, the applicability of the PGD scoring system to grade inflammation has never been tested in a systematic way. Because PGD is a multifactorial disease, progression from low to high PGD scores is likely to have a large component of case specificity, depending on the combination of pathogens involved, the overall health status of the fish and the environmental conditions. Despite complex etiology, the phenotypes generated by PGD are macroscopically similar (inflammation and hyperplasia of respiratory epithelium), which suggests that at least some of the underlying mechanisms associated with this pathology may be common (Mitchell and Rodger, 2011; Bruno et al., 2013). Unlocking these shared mechanisms is essential for early detection of gill disease, developing strategies to improve gill health and for fine-tuning salmon husbandry practices to the challenging conditions of a rapidly changing marine environment.
To identify the common gene expression patterns associated with gill inflammation, we examined Atlantic salmon from three marine production sites in Scotland. The sites and time of sampling were chosen to ensure (1) differences in fish cohort (hatchery origin, on-site treatments and overall health status), (2) diversity of pathogen and (3) differences in local environment. For each site, we selected fish with low and high PGD scores, analyzed their gill transcriptome (RNA-seq) and gill histopathology (microscopic examination), and then integrated these data to explore gene expression patterns associated with multifactorial gill disease.
Materials and Methods
Sampling and Gross Morphology
Fish (Atlantic salmon, Salmo salar) were sampled at three marine production sites belonging to Scottish Sea Farms (A on Isle of Mull and B and C in Shetland) between October 2017 and March 2018 (for details see Table 1 and Supplementary Figure 1). All fish were of strain Fanad and originated from the same egg fertilization batch. They were reared in different hatcheries (Couldoran, Pettigo-Damph and Knock-Frisa for sites A, B and C, respectively) for one year and entered the sea in spring 2017.
During sampling, randomly selected fish (76 in total) were placed in an anesthetic bath (∼20 g of MS-222/150 L) for 5–10 min, followed by macroscopic examination (gross morphology) of gill tissue as per routine procedure to monitor fish welfare indicators. This procedure is performed by fish health professionals, using semi-quantitative 6-grade scoring systems for PGD, from 0 (no visual pathology) to 5 (severe visual pathology) (Gill Health Initiative, 2015; Bloecher et al., 2018). Each of the 8 gill arches (4 on left and 4 on right side) was scored separately, generating in total 8 PGD scores per fish. Although scoring for PGD is typically non-destructive and requires only light anesthesia (e.g., 10 g of MS-222/150 L), fish in our study were subjected to terminal anesthesia for subsequent harvesting of tissue samples. Immediately after PGD scoring, fish were bled and the gill arch with the highest PGD score was excised for both transcriptome profiling and histopathological examination. For transcriptomic profiling, 3 transverse sections (gill arch ∼3 mm long with 3–4 pairs of filaments) from dorsal, medial and ventral regions of the gill (Figure 1) were transferred to RNAlater® (Sigma-Aldrich, St. Louis, MO, United States), kept at 4°C overnight for equilibration and then stored at −80°C prior to RNA extraction. All remaining tissue from the same gill arch was transferred to freshly prepared seawater Davidson’s fixative (Cadoret et al., 2013) for 24 h, followed by a short-term exposure to 10% neutral buffered formalin before tissue processing for histopathological examination.
Figure 1. Approach to representatively divide the gill arch with the highest PGD score between transcriptome profiling and histopathological examination. For transcriptomic profiling, 3 transverse sections from dorsal, medial and ventral regions of the gill were subjected to total RNA extraction to generate individual RNA samples, which were subsequently pooled. All remaining tissue was used for histopathology.
Selection of Fish With Low and High PGD Scores
Grouping all 76 sampled fish by their highest PGD score and site demonstrated that each site had fish with PGD scores 1, 2 or 3, but not 0 or 5, and only rarely 4 (Table 2). Thus, the comparison between fish with low and high PGD scores was restricted to PGD scores 1 and 3, respectively. In total, 45 gill arches (27 with PGD score 1 and 18 with PGD score 3) were subjected to total RNA extraction, but the final RNA-seq sample size dropped to 43 gill arches (from 43 fish), due to losses at the both pre- and post-sequencing stage (Table 2).
Table 2. Results of gross gill scoring for proliferative gill disease (PGD), performed on all 8 gill arches of individual fish to identify the gill arch with the highest PGD score for each fish.
Total RNA extraction was performed on individual gill transverse sections (45 gills × 3 sections) after removing the arch tissue and leaving only full-length filaments for further processing (Figure 1). Briefly, 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, Germany). Following isolation, the RNA was quantified by spectrophotometry (NanoDrop Technologies, Wilmington, DE, United States) and its integrity was confirmed by electrophoresis (Agilent Technologies, Santa Clara, CA, United States). The three individual RNA samples that originated from the same gill were then pooled to generate a single RNA sample per gill tissue (n = 45 RNA pools in total), with an equimolar contribution of RNA from dorsal, medial and ventral regions of the gill to each pool. All but one pooled gill RNA samples had a 260/280 ratio > 1.8 and RIN number > 9.3, thus meeting the criteria for RNA-sequencing. The sample with degraded RNA was eliminated from further processing.
RNA-seq Library Preparation and Sequencing
RNA-seq library preparation and sequencing were carried out by Edinburgh Genomics at the University of Edinburgh (United Kingdom). The libraries for each of the 44 samples were constructed using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina, San Diego, CA, United States), according to the manufacturer’s instructions. The paired-end sequencing (50 bp from each end) was performed on the NovaSeq 6000 system with S2 flow cell (Illumina, San Diego, CA, United States) at a sequencing depth of ∼50 million read pairs per library. The raw reads in BCL format were converted to FastQ format with bcl2fastq2 Conversion Software v2.19.1 (Illumina, San Diego, CA, United States). All raw sequences have been deposited in the ArrayExpress repository1 under accession number E-MTAB-8855.
To assess the quality of the sequencing data, reads were analyzed with FastQC v0.11.8 (Andrews, 2010). Sequencing adaptors and sequences shorter than 20 bp were removed using Flexbar v3.4.0 (Dodt et al., 2012). Filtered reads were then mapped to the Atlantic salmon reference genome ICSASG_v2 (GenBank: GCF_000233375.1, Lien et al., 2016) using HISAT2 v2.1.0. (Kim et al., 2015) with the stranded library preparation parameter. Overall, alignment rates ranged from 93.2 to 95.7%. Aligned reads were counted at gene locations using featureCounts v1.6.4 (Liao et al., 2014). For multi-mapping reads, a fractional count (1/n) was generated for each reported alignment of the multi-mapping read, with n reflecting the total number of alignments reported for that read.
Differential Gene Expression
Differential expression analysis was performed using the Bioconductor package edgeR v3.22.5 (Robinson et al., 2010) in R v3.5.1 (R Core Team, 2018). Genes with a CPM (count per million) < 1 in three or more samples were removed, resulting in 35996 genes for analysis. Filtered counts were subsequently normalized using a trimmed mean of M-values (TMM) between each pair of samples. Based on exploratory data analysis, one library was identified as an outlier and removed from subsequent analysis (for details on sample size see Table 2). Data for the remaining 43 fish were modeled using a negative binomial generalized log-linear mixed model that included both group (PGD scores 1 and 3) and site (A, B and C) as fixed effects. In total, four contrasts were generated from the same model: PGD 3 vs PGD 1 (with the site effect blocked) and then sites A vs C, B vs C and A vs B (with the group effect blocked). Differentially expressed genes (DEGs) were identified at false discovery rate (FDR) < 0.01 and absolute Log2 fold change (FC) > 1.
Functional Analysis of Gene Expression
Salmon DEGs were mapped to human orthologs to generate HGNC (HUGO Gene Nomenclature Committee) gene identifiers for functional analysis of the RNA-seq results. This approach has been demonstrated to improve biological interpretation of the salmon gene expression profiles by providing access to well-annotated databases and tools for mammalian model organisms, despite limitations of the mapping due to the extra genome duplication events in teleost fish and species-specific differences in gene function and molecular pathways (Song et al., 2014; Król et al., 2016). Mapping was done by aligning salmon transcript sequences to the protein sequences from the human genome (release 88, downloaded from Ensembl at2) using BLASTX (version 2.2.31) (Camacho et al., 2009) with an E-value cut off of 0.00001 and a maximum of one target sequence for every transcript. As one transcript can have multiple hits against one target sequence, custom Python scripts were used to filter the blast results to contain only the hit with the highest identity for each transcript. Although the majority of the salmon genes mapped to a unique human ortholog, some salmon genes mapped to the same human ortholog (for details see Results). To obtain a single expression value (mean Log2 FC) per HGNC gene identifier, the expression of the salmon genes mapped to the same human ortholog was either averaged (if contributing salmon genes had similar expression profiles) or based on the expression of the salmon gene that was more abundant (if contributing salmon genes had contrasting expression profiles). For biological interpretation of the RNA-seq results, human orthologs of the salmon DEGs along with their Log2 FC values were analyzed using Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity) to explore (1) enrichment of canonical pathways, (2) upstream regulators and (3) downstream effects associated with these genes. The same set of genes was also submitted to PANTHER Classification System (Mi et al., 2013) to perform Gene Ontology (GO) enrichment analysis.
Gill tissue was routinely dehydrated in ethanol, equilibrated in xylene and embedded in paraffin wax according to standard histological techniques (Bancroft and Gamble, 2007). Sagittal sections (3 μm) of the gill arch were cut with a microtome and mounted onto microscope slides. These sections were then subjected to haematoxylin and eosin (H&E) staining. All sections were digitized at 40 × magnification, using the Olympus dotSlide 2.1 Virtual Slide System (Olympus Corporation, Tokyo, Japan). The resultant images were randomized to ensure blinded examination and then scored using a system developed to assess gill histopathology in sea farmed Atlantic salmon (Mitchell et al., 2012), with minor amendments. The details of the scoring system are presented in Table 3.
Table 3. Semi-quantitative scoring system used for histopathological examination of gill tissue in sea farmed Atlantic salmon (adapted and modified from Mitchell et al., 2012).
Prior to the analysis, three variables were excluded (bE, bT and ob) as scores were invariant across fish (see Results). To compare the gill histopathology scores between PGD groups (1 and 3) and also between fish from different sites (A, B and C), non-metric multidimensional scaling (NMDS) was performed using the metaMDS function from the vegan package (Oksanen et al., 2019) in R. For the NMDS analysis, a Jaccards dissimilarity index was used to calculate the dissimilarity matrix, two dimensions were specified and 100 random starts used. Similarities between fish histopathology scores were visualized using a biplot with 95% confidence ellipses around the group centroids and variable vectors included.
Phenotypic Characteristics of Fish
Body mass of 43 fish subjected to the RNA-seq experiment varied from 1.5 to 4.0 kg, with their body length ranging from 47 to 62 cm (Figures 2A,B). The effects of group (PGD 1 and PGD 3) and site (A, B and C) on body mass and length were not significant (P > 0.05, 2-way ANOVA). As a result, the fish did not differ in their Fulton’s condition factor K (P > 0.05, 2-way ANOVA), which varied from 1.1 to 1.7 (Figure 2C). Given the small sample size (Table 2), these results should be treated with caution.
Figure 2. Individual body mass (A), body length (B) and Fulton’s condition factor K (C) of 43 fish subjected to RNA-seq experiment, plotted by group (PGD 1 and PGD 3) and site (A, B and C). For formula to calculate K, see Table 1.
The NMDS analysis of 43 gill transcriptomes revealed no differences between 26 fish with PGD score 1 and 17 fish with PGD score 3, as indicated by their overlapping 95% confidence intervals (Figure 3A). This finding was reinforced by differential gene expression analysis, which showed 0 DEGs for comparison PGD 3 vs PGD 1 (Table 4 and Supplementary Table 1). In contrast, gill transcriptomes were clearly separated by site, with no overlap in 95% confidence intervals for 20 fish from site A, 10 fish from site B and 13 fish from site C (Figure 3B). The site-specific differences in the gill transcriptomes are reflected in the number of DEGs identified between the sites: 1360 for A vs C, 708 for B vs C and 240 for A vs B (Table 4, for the lists of DEGs see Supplementary Tables 2–4).
Table 4. Results of differential gene expression analysis performed on gill transcriptome of 43 fish to elucidate the differences between groups (PGD 3 vs PGD 1) and sites (A vs C, B vs C and A vs B).
Figure 3. Gill transcriptome of 43 fish grouped by PGD score (A) and site (B). Each panel shows a NMDS plot of gene expression profiles between different fish (numbers are fish IDs). The distances on the plots correspond to the leading fold change (FC), which is the average (root-mean-square) Log2 FC for the 500 genes most divergent between each pair of fish. Ellipses indicate 95% confidence intervals, overlapping for fish grouped by PGD score (A) but not for fish grouped by site (B). The stress value of the NMDS ordination is 0.162.
Histopathological examination of the gill tissue generated 15 scores (ranging from 0 to 3) per fish, focusing on 4 index and 11 ancillary criteria (Supplementary Table 5). Distribution of the histopathological scores by PGD score (1 and 3) and site (A, B and C) for each of 15 criteria is presented in Supplementary Figure 2. The NMDS analysis of the gill histopathology results for 43 fish showed no differences between 26 fish with PGD score 1 and 17 fish with PGD score 3, as indicated by their overlapping 95% confidence intervals (Figure 4A). In contrast, gill histopathology differed between sites, with 1) fish from site C being clearly separated from other fish and 2) fish from site A overlapping with fish from site B (Figure 4B). The position of the site-associated clusters along the NMDS1 axis indicates that fish from site C (with lower scores) show relatively low level of gill histopathology, while fish from sites A and B (with higher scores) have gill tissue with relatively moderate changes in histopathology. The main drivers of these differences are lamellar hyperplasia (LH), lamellar fusion (LF), cellular anomalies (CA) and Neoparamoeba-like protist parasites (pp), as evidenced by the vectors paralleled to the NMDS1 axis in Figure 4B. Overall, fish from site A and independently from site B developed significantly higher grade inflammation than fish from site C, thus facilitating the search for the molecular bases of multifactorial gill disease.
Figure 4. Gill histopathology of 43 fish grouped by PGD score (A) and site (B). Each panel shows a NMDS plot of histopathological profiles between different fish (numbers are fish IDs). The distances on the plots are calculated from the scores of 12 histopathological parameters (LH, LF, LO, CA, in, eg, cc, cd, ib, ch, pp and op), with 3 parameters (bE, bT and ob) removed from the analysis because all scores were 0 (for parameter abbreviations and details see Table 3 and Supplementary Figure 2). Ellipses indicate 95% confidence intervals, overlapping for fish grouped by PGD score (A) and for fish from sites A and B but not from C (B). Vectors are fitted to visualize the contribution of histopathological parameters to the ordination. The stress value of the NMDS ordination 0.190.
Identification of Gene Expression Patterns Associated With Multifactorial Gill Disease
Our approach to identify the common gene expression patterns of non-specific gill inflammation is explained in Figure 5. Specifically, all gill transcriptomes from sites A and B (with moderate gill histopathology) were independently compared to all gill transcriptomes from site C (with low gill histopathology), which resulted in 1360 and 708 DEGs (at FDR < 0.01 and absolute Log2 FC > 1) for A vs C and B vs C contrasts, respectively (Table 4, for the lists of DEGs see Supplementary Tables 2, 3). A relatively large number of genes were common between these two comparisons (462 DEGs in total), including 354 protein-coding genes, 35 immunoglobulin gene segments, 15 pseudogenes and 58 non-coding RNAs (Table 5, for the list of common DEGs see Supplementary Table 6). The Log2 FC values for these genes were averaged between A vs C and B vs C contrasts to provide one expression value per common gene. Importantly, all common DEGs show the same direction of change in both the moderate-to-low histopathology contrasts, i.e., they are either upregulated or downregulated in relation to site C, as indicated in Figure 5. That commonality of the response suggests that the identified DEGs are part of the gene expression profile indicative of progression from lower to higher grade inflammation of gill tissue. To understand the underlying mechanisms, we focused on the 354 protein-coding DEGs and predicted their functionality through mapping to the human orthologs.
Figure 5. Approach to identify gene expression patterns associated with multifactorial gill disease. Based on the histopathological examination, A and B were classified as two independent sites with relatively moderate gill histopathology, while C was classified as a site with relatively low gill histopathology. Firstly, the lists of differentially expressed genes (DEGs) for A vs C and B vs C contrasts were generated (FDR < 0.01 and absolute Log2 FC > 1), yielding 1360 and 708 transcripts, respectively. Secondly, both lists of DEGs were checked for common genes (462 in total). Finally, the list of common genes was screened for a commonality of response to ensure that only genes either upregulated or downregulated in both sets of DEGs were considered to constitute a common gene expression profile of non-specific gill inflammation.
Mapping Salmon Genes to Human Orthologs
The results of blasting all 462 salmon DEGs associated with multifactorial gill disease against the protein sequences from the human genome (BLASTX, E-value < 0.00001, top hit) are presented in Table 5 and Supplementary Table 6. Among them, 311 of 354 protein-coding salmon genes were matched to HGNC gene identifiers, but not all of them were unique. Specifically, 191 protein-coding salmon transcripts mapped uniquely to one human ortholog, while the remaining 120 protein-coding salmon transcripts mapped to 44 human orthologs, with > 1 salmon gene mapping to the same HGNC gene identifier (Supplementary Table 7). The expression of salmon genes mapped to the same human ortholog was averaged to provide one expression value (mean Log2 FC) per HGNC gene identifier, if the response of the salmon genes was consistent (i.e., all contributing salmon genes were either upregulated or downregulated, but not both). In two cases with the contrasting expression of salmon genes (PRF1 and CTSV, for details see Supplementary Table 7), the mean Log2 FC was based on the expression of the salmon gene that was more abundant (higher CPM). As a result, the expression patterns of 311 protein-coding salmon DEGs were represented by 235 human orthologs and their corresponding 235 expression values (Log2 FC for 191 unique mapping and mean Log2 FC for 44 multiple mapping), which were then used to predict functionality.
Functional Analysis of Gene Expression Patterns Associated With Multifactorial Gill Disease
Biological interpretation of the transcriptomic changes in gills with moderate histopathology (sites A and B) vs low histopathology (site C) was performed on the expression profiles of 311 protein-coding salmon DEGs converted into 235 human orthologs (Supplementary Table 7), using (1) IPA for enrichment of canonical pathways, upstream regulators and downstream effects and (2) PANTHER Classification System for enrichment of GO terms (Biological Process). The predictions made by IPA also require the expression values, while GO enrichment analysis is based on the list of DEGs.
IPA identified 13 canonical pathways that were altered in gills with moderate histopathology (sites A and B) vs low histopathology (site C) at P-value < 0.01, with contribution of 35 DEGs in total (Figure 6 and Supplementary Table 8). Most of these pathways are associated with 1) cellular immune response (IL-17 Signaling, IL-6 Signaling, Granzyme A Signaling, Crosstalk between Dendritic Cells and Natural Killer Cells, Agranulocyte Adhesion and Diapedesis and HMGB1 Signaling), 2) cytokine signaling (IL-17 Signaling, IL-6 Signaling, Acute Phase Response Signaling, Role of JAK family kinases in IL-6-type Cytokine Signaling, TNFR2 Signaling and HMGB1 Signaling) and (3) tissue damage and repair. The tissue damage and repair were evidenced by alterations of pathways related to cellular stress and injury (Autophagy and HMGB1 Signaling), apoptosis (TNFR2 Signaling and JAK/Stat Signaling), cellular growth (STAT3 Pathway and JAK/Stat Signaling) and proliferation and development (STAT3 Pathway and JAK/Stat Signaling).
Figure 6. Top canonical pathways altered in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA) and P-value < 0.01. For details and list of corresponding genes see Supplementary Table 8.
IPA analysis of upstream regulators is based on prior knowledge of predictable effects between transcriptional regulators (e.g., transcription factors, cytokines, microRNAs, receptors, kinases, chemicals and drugs) and their target genes, stored in the Ingenuity Knowledge Base. When such analysis was performed on the gene expression profiles indicative of multifactorial gill disease, IPA identified 15 top upstream regulators associated with 101 of the 235 submitted DEGs, all highly significant (overlap P-value < 0.0001) and activated (z-score > 2) (Figure 7 and Supplementary Table 9). Among them were one endotoxin (LPS), one pattern recognition receptor (NOD2), 10 cytokines (IL1B, IFNG, TNF, IL2, LIF, IL12B, IL12A, IL1, IFNA2 and IL6), two growth factors (PDGF and AGT) and one transcription factor (NFKB).
Figure 7. Top upstream regulators predicted from gene expression in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA). For details and list of corresponding genes see Supplementary Table 9.
IPA analysis of downstream effects predicts potential outcomes from gene expression data, using the Ingenuity Knowledge Base of differential gene expression in varying disease and functional states. The predictions made for gills with moderate histopathology (sites A and B) vs low histopathology (site C) included 13 downstream disease/functions at P-value < 9.4E-08, with nearly all DEGs (230 of 235) contributing to the effects (Figure 8 and Supplementary Table 10). The dominant features of these predictions are (1) immune and inflammatory diseases (Inflammatory Response, Infectious Diseases, Inflammatory Disease, Immune Cell Trafficking, Immunological Disease and Humoral Immune Response), (2) tissue damage and repair (Organismal Injury and Abnormalities, Tissue Morphology, Connective Tissue Disorders, Cellular Movement, Cell Death and Survival, Cellular Function and Maintenance) and (3) intra-tissue communication (Cell-To-Cell Signaling and Interaction).
Figure 8. Top downstream effects predicted from gene expression in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA). For details and list of corresponding genes see Supplementary Table 10.
GO enrichment analysis identified 15 GO terms (Biological Process) that were overrepresented with the majority of DEGs (209 of 235) associated with multifactorial gill inflammation, generating fold enrichments from 1.2 to 2.6 at Bonferroni-corrected P-value < 0.05 (Figure 9 and Supplementary Table 11). The enriched GO terms pointed towards increased (1) intra-tissue communication (Cellular response to cytokine stimulus and Cell communication), (2) presence of external stimuli (Response to external stimulus and Response to oxygen-containing compound), (3) activated immune response (Immune system process) and 4) ongoing tissue remodeling (Animal organ development and Regulation of biological quality).
Figure 9. Top GO terms (Biological Process) associated with differentially expressed genes (DEGs) in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using PANTHER overrepresentation test and Bonferroni-corrected P-value < 0.05. The GO terms are sorted by the fold enrichment of the most specific categories, with their parent terms indented directly below. For details and list of corresponding genes see Supplementary Table 11.
Top Genes Associated With Multifactorial Gill Disease
Top genes were defined here as the DEGs with the largest magnitude of change in expression (absolute Log2 FC > 2), thus including the protein-coding genes that were at least 4 times higher or 4 times lower expressed in gills with moderate histopathology (sites A and B) vs low histopathology (site C). The list of these top genes is presented in Table 6, with 43 upregulated salmon genes mapped to 25 unique human orthologs and 14 downregulated salmon genes mapped to 12 unique human orthologs.
Our study is the first to focus on the gill transcriptome of Atlantic salmon in a highly variable and complex environment of the marine production sites. We specifically chose to work on fish farmed in the open sea-based system (such as net pens), because (1) the use of net pens is currently the most common form of salmon production (Philis et al., 2019) and (2) in the open pens, only nets separate the fish from the environment, making them prone to developing gill pathologies but also assisting with the tissue recovery through high water exchange and oxygenation (Shinn et al., 2015; Nowak and Archibald, 2018). Farming in net pens is typically associated with the multifactorial gill diseases (PGD, PGI or CGD) rather than the single-cause gill pathologies, such as AGD (Herrero et al., 2018; Gjessing et al., 2019; Laurin et al., 2019). So far, the extensive transcriptomic profiling of the gill tissue has been performed only in salmon exposed to the controlled environment of the closed-system tanks, following a single-pathogen challenge to induce either AGD (Morrison et al., 2006; Wynne et al., 2008a; Bloecher et al., 2018; Boison et al., 2019; Robledo et al., 2019) or infectious salmon anemia (Valenzuela-Miranda et al., 2015). In most cases, fish were subjected to a single-dose challenge of the infectious organisms, with only few studies evaluating the effects of re-infection (Bloecher et al., 2018; Boison et al., 2019) or co-exposure to other infectious agents, such as hydroids (Bloecher et al., 2018) and Yersinia ruckeri (Valdenegro-Vega et al., 2015). To closely resemble the dynamics of multifactorial gill diseases observed in the production systems, fish need to be continuously exposed to a variety of infectious agents and environmental stressors that impact their health.
To evaluate the robustness and reliability of the PGD scores (gross morphology) in conveying information about gill health, we sampled salmon at three geographically distant locations in autumn and spring, without prior knowledge of the biotic and abiotic conditions of these sites. By doing this, we included in our data set the gill samples that were diverse in terms of origin, time of sampling, fish background and the overall environmental milieu (Table 1). The multi-site sampling approach to identify the common aspects of the multifactorial gill inflammation has been used before (Kvellestad et al., 2005; Steinum et al., 2010; Gjessing et al., 2019), but not in the context of the transcriptome analysis or underlying molecular mechanisms. Our study is the first to investigate the gill tissue from salmon at three different levels of organization: gross morphology (macroscopic examination), histopathology (microscopic examination) and whole transcriptome (gene expression) (Figure 1). We have focused on the gill arch with the highest PGD score per fish, as this arch may provide more accurate information about the general state of the gill health than the arch sampled randomly or the arch sampled because of its position, as practiced elsewhere (Wise et al., 2008; Steinum et al., 2010; Benedicenti et al., 2019; Gjessing et al., 2019). For better integration of the transcriptomic data with the gill gross morphology and histopathology, the tissue fragments used for the total RNA extraction aimed to represent the whole surface area of the gill rather than some specific areas of the tissue (Figure 1). This is in contrast to many gill gene expression (qPCR) studies, which tend to focus on either interbranchial lymphoid tissue or visible pathologies of the gill tissue, such as mucoid patches, hyperplastic lesions and lamellar fusions (Pennacchi et al., 2016; Marcos-López et al., 2018; Benedicenti et al., 2019).
PGD Scores as Measures of Gill Health
Effective treatment of salmon gill diseases requires development of diagnostic and prognostic tools that are non-invasive and suitable for a frequent use in the rapidly changing marine environment (Herrero et al., 2018; Gjessing et al., 2019; Rozas-Serri, 2019). We have specifically focused on the applicability of the PGD scores (gross morphology) to reflect gill health and underlying pathology, as this scoring system is already in use as part of routine fish welfare indicator assessment of salmon health in seawater farms worldwide (Gill Health Initiative, 2015; Bloecher et al., 2018). By comparing 43 gill samples with low (1) and high (3) PGD scores across three marine production sites in Scotland, we found that these two groups of gill tissue classified as different at the macroscopic level (PGD 1 and PGD 3) were in fact indistinguishable at the level of whole-transcriptome gene expression (Figure 3A) and also indistinguishable at the level of microscopic histopathology (Figure 4A). Furthermore, we could not identify any single gene that was expressed differently between the two groups (Table 4). Our results strongly suggest that the changes in gross morphology were not supported by the changes in gene expression or histopathology. The lack of detectable transcriptomic and/or histopathological changes associated with the progression of the PGD scores questions the suitability of the gross scoring system as diagnostic and prognostic tools to monitor and control both existing and emerging gill diseases. Instead, the PGD scores are good proxies for monitoring changes in macroscopic gill surface area available for respiration, the knowledge of which may be essential for fish health management and husbandry practices in aquaculture settings.
Table 6. Top genes altered in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C), based on FDR < 0.01 and absolute Log2 FC > 2.
The PGD scoring system of salmon gill is based on the 6 macroscopic grades, from 0 (no visual pathology) to 5 (severe visual pathology) (Gill Health Initiative, 2015; Bloecher et al., 2018). It is important to realize that our comparison was done on the groups of tissue that differed only by two grades (PGD 1 and PGD 3) and were in the middle range of the spectrum, with a bias towards less damaged gills. Although our intention was to compare the transcriptome and histopathology of gill tissues with the full spectrum of PGD scores (e.g., PGD 0 and PGD 5) from the same farm location, fish with such a large range in gill gross morphology were not found during the sampling events, not only within the three production sites described in the current study (Table 2), but also at the four other locations (Scottish Sea Farms) visited by us in years 2017–2019 (R. Bickerdike and S. A. M. Martin, unpublished data). It is unusual for the salmon farmed in the sea to have all 8 gill arches scored as PGD 0 (apart from the first few weeks following the transfer from the freshwater facility), while the macroscopic gill damage classified as PGD 5 is very rare and would not be expected to be found until later in the production cycle, and as a result of a specific acute event or from cumulative significant environmental insults over a period of time (R. Bickerdike, personal communication).
Discrepancies between macroscopic and microscopic examination of the gill tissue have been reported by previous studies, mainly in the context of AGD (Pennacchi et al., 2016). For example, it has been shown that the AGD scoring system (gross morphology) may not be a reliable method of confirming the disease in cases of light severity of AGD (Clark and Nowak, 1999; Zilberg et al., 2001). This is because small AGD-associated lesions, which affect < 10 lamellae and are easy to detect under the microscope, are typically overlooked during the visual inspection of the gill (Adams et al., 2004). As a result, fish exposed to Neoparamoeba perurans and then classified as clear of the AGD symptoms during the macroscopic gill scoring for AGD may in fact need to be re-classified at the level of histopathological examination (Wynne et al., 2008b). Further complexity to the diagnostic problems is added when the gill disease is multifactorial (Wise et al., 2008; Gjessing et al., 2019; Noguera et al., 2019). More broadly, the poor diagnostic and prognostic value of the PGD scores demonstrated in our study is consistent with the limited applicability of the gross morphology to diagnose complex diseases in livestock and humans, which typically need extensive histopathology and molecular profiling to confirm and prognosticate (Hoffmann et al., 2009; Ahmed and Abedalthagafi, 2016; Mobadersany et al., 2018).
Histopathology-Directed Analysis of Gill Transcriptome Between Sites
By sampling sea farmed Atlantic salmon at three production sites (A, B and C) in Scotland, we demonstrated that the gill samples from different sites had different histopathology (Figure 4B) and different transcriptomic profiles (Figure 3B), with 240 to 1360 genes identified as differentially expressed between the sites (Table 4). The drivers of this pronounced site-to-site variability in gill histopathology and transcriptome are unknown and may reflect the differences in (1) fish cohort, including their hatchery origin, past and present husbandry practices and overall health status, (2) local biotic conditions (e.g., diversity of pathogen), (3) local abiotic conditions (e.g., sea water temperature) and (4) interplay and interactions between all the above factors. Identification of these drivers requires more rigorous spatial and temporal sampling regime (Jokinen et al., 2012; Maestrini and Basso, 2018) and is beyond the scope of the current study.
The NMDS plots revealed that the gill samples from site C were different (i.e., clearly separated from sites A and B) not only at the level of their transcriptome (Figure 3B), but also in terms of their histopathology (Figure 4B). Because the NMDS ordination of the histopathological data is based on the same set of parameters (i.e., LH, LF, LO, CA, in, eg, cc, cd, ib, ch, pp and op) across all fish, lower NMDS1 values of gill samples from site C reflect lower scores of the contributing parameters (lower grade inflammation) than higher NMDS1 values of gill samples from sites A and B (Figure 4B, Supplementary Table 5, and Supplementary Figure 2). In contrast, the NMDS ordination of the transcriptomic data (Figure 3B) is based on the expression of 500 genes that are most divergent between each pair of fish (amounting to 903 sets of genes for all combinations of 43 fish in total) rather than being performed on the same set of genes across all fish (Robinson et al., 2010). Thus, the higher leading Log2 FC dim 1 values of gill samples from site C refer to their generally higher magnitude of transcriptomic changes rather than to the specific set of genes that separates them from sites A and B, with latter requiring more detailed analysis of the existing data to perform.
Because gill samples were more different between sites than within sites (Figures 3B, 4B), we performed the histopathology-directed analysis of gill transcriptome on all fish belonging to each site rather than on the individual fish within or between sites (Figure 5). Specifically, we classified all fish from site C as having gills with low histopathology, while all fish from site A and independently all fish from site B as having gills with moderate histopathology. Comparing the two independent sets of the moderately inflamed gill transcriptomes (A and B) with the low inflammation gill transcriptome (C), in the A vs C and B vs C contrasts, led to the identification of the common DEGs that are likely to constitute the gene expression patterns of multifactorial gill inflammation. Similar approaches utilizing independent sampling and site-to-site comparisons to establish common transcriptomic responses to a range of environmental factors have been successfully applied in studies ranging from aquatic toxicology (Oleksiak, 2008; Defo et al., 2018) to human medicine (Mueller et al., 2009; Rohart et al., 2017).
Gene Expression Patterns of Multifactorial Gill Disease
Transcriptome profiling by RNA-seq has become a powerful tool for identification of genes and molecular pathways involved in the progression from health to disease in farmed fish (Sun et al., 2012; Martin and Król, 2017; Houston and Macqueen, 2019; Ronza et al., 2019). Combining RNA-seq assays with an independent sampling regime is especially important for understanding multifactorial diseases, because it allows for the potential identification of the common gene expression patterns of the disease that is essential for developing treatment and prevention strategies (Cookson et al., 2009; Bhuju et al., 2012). This approach has for example been used in the context of non-specific inflammatory diseases in fish intestine (reviewed in Martin et al., 2016). The present study is the first to focus on the transcriptomic patterns of multifactorial gill disease, a condition increasingly prevalent in salmon sea farms.
Transcriptomic comparison of the two independent sets of moderate histopathology gill samples (sites A and B) with the low histopathology gill samples (site C) led to the identification of 462 common DEGs in total, 311 of which were protein-coding transcripts mapped to 235 human orthologs (Table 5). The identification of these genes as part of the gene expression patterns of multifactorial gill disease was further supported by their uniform direction of changes in the expression levels (Figure 5). Subsequent functional analysis of these genes by IPA and GO pointed towards a range of immune responses as the most dominant feature of the moderate histopathology gills.
At the gene level, the immune responses were driven by pro-inflammatory cytokines (IL17F, CXCL9, CXCL10, CCL4L1 and TNF superfamily members: LTA and TNFSF14), cytokine receptors (TNFRSF1B, TNFRSF6B and CXCR1) and regulators of cytokine expression and signaling (C1QTNF3, SOCS1 and SOCS3), all of which were upregulated apart from TNFSF14 (Supplementary Table 7). The proteins encoded by IL17 family genes have been shown to stimulate the production of several other cytokines, including IL6 and IL8 in splenocytes of rainbow trout using recombinant IL17A/F2a (Monte et al., 2013) as well as IL1β, IL6, IL8 and TNF-α in head kidney leukocytes of grass carp using recombinant IL17A/F1 (Du et al., 2015). Furthermore, the emerging role of IL17 cytokines in gill mucosal immunity is supported by reports showing differential expression of these genes following the exposure to Aeromonas hydrophila in common carp (Dong et al., 2019) and after challenge with Ichthyophthirius multifiliis in rainbow trout (Syahputra et al., 2019). According to our BLASTX results, CXCL9 represents salmon IL8, whose high expression in our data set (Log2 FC = 3.4, Supplementary Table 7) is consistent with upregulation of IL17F (Zou and Secombes, 2016). The implication of IL8 in gill immune responses has been previously demonstrated in Atlantic cod (Caipang et al., 2010) and rainbow trout (Olsen et al., 2011; Santana et al., 2016). The protein encoded by LTA (salmon TNF-α2) mediates a large variety of inflammatory, immunostimulatory and antiviral responses in mammals (Upadhyay and Fu, 2013), showing induction at the level of gene expression in gills of Atlantic bluefin tuna during natural infection with Digenea (Pleić et al., 2015), but not in Atlantic salmon following the infection with Neoparamoeba spp. to induce AGD (Morrison et al., 2007). Cytokines were also overwhelmingly present in our functional analysis of gene expression patterns in the moderately inflamed gills. Specifically, 6 of the top 13 IPA canonical pathways were related to cytokine signaling (Figure 6 and Supplementary Table 8) and 10 of the top 15 IPA upstream regulators were predicted to be cytokines (Figure 7 and Supplementary Table 9). The presence of cytokines among DEGs and also as upstream regulators is consistent with a cytokine signaling cascade, during which one cytokine stimulates its target cells to produce another cytokines (Tisoncik et al., 2012). As a result of cytokine presence among DEGs, 6 of the top 13 IPA downstream effects were predicted to be immunological and inflammatory diseases (Figure 8 and Supplementary Table 10). Altogether, our results identified cytokine-driven immune response as a hallmark of multifactorial gill disease.
Besides the immune response, the multifactorial gill disease from sites A and B had a common transcriptomic profile indicative of tissue damage and repair. Part of the tissue damage could be inflicted by high levels of nitric oxide (NO), suggested by upregulation of NOS2 (Log2 FC = 2.4, Supplementary Table 7). This gene encodes a NO synthase that is inducible by a combination of lipopolysaccharide (LPS) and certain cytokines (Okamoto et al., 2004). Although NO is produced as a first-line defense against invading pathogens, its strong cytotoxic effects may also damage the tissue of the host (Abramson et al., 2001). Furthermore, high levels of NO may contribute to oxidative stress through production of reactive oxygen species (ROS) (Girouard et al., 2009), which could potentially explain the upregulation of GPX2 in our data set (Log2 FC = 3.9, Supplementary Table 7). The protein encoded by this gene belongs to the glutathione peroxidase family, members of which catalyze the reduction of organic hydroperoxides and hydrogen peroxide (H2O2) by glutathione, and thereby protect cells against oxidative damage (Brigelius-Flohé and Maiorino, 2013). The potential role of GPX2 in the antioxidant defense against increased ROS at gill mucosal surfaces has been recently discussed in the context of exposing Atlantic salmon to peracetic acid (Soleng et al., 2019). It is well established that the repair of inflamed tissue requires efficient removal of damaged cells through controlled cell death (apoptosis) and concurrent cell proliferation to regenerate damaged structures and build up lost tissue, with both processes closely linked to the activity of cysteine-dependent aspartate-directed proteases (caspases) (Fogarty and Bergmann, 2017). In our study, both CASP14 (salmon caspase-14-like) and other genes associated with apoptosis and autophagy (RNF152, CTSL, CTSG, CTSV, RAB32 and TGM2) were upregulated, pointing towards ongoing tissue repair in the moderately inflamed gills (Supplementary Table 7). The increased activity of caspases during gill inflammation has been previously documented in rainbow trout challenged with Aeromonas salmonicida (Rojas et al., 2015) and grass carp after infection with Flavobacterium columnare (Chen et al., 2019). Because of the complexity and size of gill vasculature, the repair of the gill tissue requires extensive vasculogenesis (formation of new blood vessels from vascular precursor cells), angiogenesis (process of outgrowing vessels from the existing vasculature) and arteriogenesis (remodeling of arteries where collateral arterial anastomoses undergo abluminal expansion) (Shi et al., 2017). One of the most potent mediators of new blood vessel formation is angiogenin (Hoang and Raines, 2017), represented in our data set by ANG (salmon angiogenin-1 and ribonuclease-like 3), which was the gene with the highest level of upregulation (Log2 FC = 4.6, fold change ∼24) in the moderate histopathology vs low histopathology gills (Table 6 and Supplementary Table 7). This result is consistent with the earlier study conducted by Valdenegro-Vega et al. (2014), showing increased protein levels of angiogenin (fold change ∼12) in the gills of Atlantic salmon following four successive infections with Neoparamoeba perurans. Among other genes involved in angiogenesis were ADM, BMP10 and VCAN, all upregulated in the moderately inflamed gills (Supplementary Table 7). Overall, 5 of top 13 IPA canonical pathways (Figure 6 and Supplementary Table 8) and 5 of top 13 IPA downstream effects (Figure 8 and Supplementary Table 10) were associated with cell death and proliferation related to tissue damage and repair, highlighting their importance in multifactorial gill disease.
Conclusion and Future Directions
The gill is central for understanding the impacts of climate change on fish health. Recent increases in gill pathologies in sea farmed Atlantic salmon highlight the need to establish the molecular basis of multifactorial gill disease (frequently referred to as PGD, PGI or CGD) to improve diagnosis and preventive measures of this condition. To ensure the multifactorial etiology of gill disease, we sampled Atlantic salmon from three different production sites in Scotland and then examined the gill tissue at three different levels of organization: gross morphology with the use of PGD scores (macroscopic examination), whole transcriptome (gene expression by RNA-seq) and histopathology (microscopic examination). By exploring the association between gill transcriptome and gill gross morphology (PGD scores), we found that the PGD scores were less effective in providing a graded assessment of gill health status than expected, and they did not convey any information about the underlying pathology and/or tissue deterioration. In contrast, integration of the gill RNA-seq data with the gill histopathology enabled us to identify common gene expression patterns associated with multifactorial gill disease. We demonstrated that the gene expression patterns associated with multifactorial gill disease were dominated by two processes: a range of immune responses driven by pro-inflammatory cytokines and the events associated with tissue damage and repair, driven by caspases and angiogenin.
Previous studies of the gill inflammation have typically focused on single-cause pathologies (such as AGD), in the experiments performed on fish exposed to the controlled environment of the closed-system tanks. Although these studies are very important from a mechanistic point of view of specific pathologies, they may have limited relevance for the fish facing a changing environment at farm locations, e.g., rising sea water temperature, deoxygenation and surges of existing and emerging pathogens with the complexity of co-infections, coupled with husbandry practices. Performing more transcriptomic studies in the field rather than in the lab would benefit both academia and industry. The multi-site approach is also important because in this study, we demonstrated large and significant effects of sampling site (sensu lato) on gill transcriptome and histopathology. The drivers of this site-to-site variability are presently unknown and require more specific sampling regime.
Our histopathology-directed analysis of gill transcriptomes identified in total 462 genes that we claim to constitute the gene expression profile of multifactorial gill inflammation, including 354 protein-coding genes, 35 immunoglobulin gene segments, 15 pseudogenes and 58 non-coding RNAs. It is, however, important to remember that our analysis was based on the gill samples from three production sites in one specific part of the world (coastal waters of Scotland), with sampling events covering October, November and March. Substantial amount of work is therefore needed to test the association of these genes with multifactorial gill diseases at different times of year in Scotland and worldwide. The diagnostic and therapeutic value of these transcripts is currently unknown and require further studies.
Data Availability Statement
The RNA-seq data have been deposited in the ArrayExpress repository (http://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-8855. The R and Python scripts are available on request to AD and SS, respectively.
Ethical review and approval was not required for the animal study because the sampling of fish was carried out under established protocols for routine health assessments in accordance with RSPCA Assured Welfare Standards for farmed Atlantic salmon and under supervision of the company veterinarian.
SM, AD, EK, VV, KG, and RB conceived and designed the study, interpreted results, and gave final approval of the manuscript. All authors (apart from SS and AD) were involved in sampling. EK extracted RNA, performed IPA analysis, and drafted the manuscript with AD, who oversaw bioinformatics and statistical analysis. SS mapped salmon genes to human orthologs to enable functional analysis of gene expression. PN performed histopathological examination of samples, generated scores, and interpreted results. EC provided support in data exploration. All authors discussed and commented on the manuscript.
The study was supported by the Scottish Aquaculture Innovation Centre (SAIC grant SL 2017 08, ‘Nutritional Aspects of Gill Disease in Atlantic Salmon’). The authors declare that this study received funding from BioMar AS and Scottish Sea Farms (SSF). The funders had the following involvement with the study: contribution to the experimental design, selection of sampling sites, participation in sampling, contribution to the interpretation of the results, and approval of the manuscript.
Conflict of Interest
KG and VV were employed by the company BioMar AS. RB was employed by the company Scottish Sea Farms (SSF).
The remaining 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.
We thank SSF farm personnel for accommodating our research, performing gross morphology scoring and helping with sampling. Edinburgh Genomics is partly supported through core grants from NERC (R8/H10/56), MRC (MR/K001744/1), and BBSRC (BB/J004243/1). We also thank two reviewers for their comments on the earlier draft of the article. EK dedicates this paper to her late mother, Irena Król.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00610/full#supplementary-material
Adams, M. B., Ellard, K., and Nowak, B. F. (2004). Gross pathology and its relationship with histopathology of amoebic gill disease (AGD) in farmed Atlantic salmon, Salmo salar L. J. Fish Dis. 27, 151–161. doi: 10.1111/j.1365-2761.2004.00526.x
Akbarzadeh, A., Günther, O. P., Houde, A. L., Li, S., Ming, T. J., Jeffries, K. M., et al. (2018). Developing specific molecular biomarkers for thermal stress in salmonids. BMC Genomics 19:749. doi: 10.1186/s12864-018-5108-9
Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed February 11, 2020).
Anttila, K., Lewis, M., Prokkola, J. M., Kanerva, M., Seppänen, E., Kolari, I., et al. (2015). Warm acclimation and oxygen depletion induce species-specific responses in salmonids. J. Exp. Biol. 218, 1471–1477. doi: 10.1242/jeb.119115
Bayliss, S. C., Verner-Jeffreys, D. W., Bartie, K. L., Aanensen, D. M., Sheppard, S. K., Adams, A., et al. (2017). The promise of whole genome pathogen sequencing for the molecular epidemiology of emerging aquaculture pathogens. Front. Microbiol. 8:121. doi: 10.3389/fmicb.2017.00121
Benedicenti, O., Pottinger, T. G., Collins, C., and Secombes, C. J. (2019). Effects of temperature on amoebic gill disease development: does it play a role? J. Fish Dis. 42, 1241–1258. doi: 10.1111/jfd.13047
Bhuju, S., Aranday-Cortes, E., Villarreal-Ramos, B., Xing, Z., Singh, M., and Vordermeier, H. M. (2012). Global gene transcriptome analysis in vaccinated cattle revealed a dominant role of IL-22 for protection against bovine tuberculosis. PLoS Pathog. 8:e1003077. doi: 10.1371/journal.ppat.1003077
Bloecher, N., Powell, M., Hytterød, S., Gjessing, M., Wiik-Nielsen, J., Mohammad, S. N., et al. (2018). Effects of cnidarian biofouling on salmon gill health and development of amoebic gill disease. PLoS One 13:e0199842. doi: 10.1371/journal.pone.0199842
Boison, S. A., Gjerde, B., Hillestad, B., Makvandi-Nejad, S., and Moghadam, H. K. (2019). Genomic and transcriptomic analysis of amoebic gill disease resistance in Atlantic salmon (Salmo salar L.). Front. Genet. 10:68. doi: 10.3389/fgene.2019.00068
Breitburg, D., Levin, L. A., Oschlies, A., Grégoire, M., Chavez, F. P., Conley, D. J., et al. (2018). Declining oxygen in the global ocean and coastal waters. Science 359:eaam7240. doi: 10.1126/science.aam7240
Burge, C. A., Eakin, C. M., Friedman, C. S., Froelich, B., Hershberger, P. K., Hofmann, E. E., et al. (2014). Climate change influences on marine infectious diseases: implications for management and society. Annu. Rev. Mar. Sci. 6, 249–277. doi: 10.1146/annurev-marine-010213-135029
Cadoret, K., Bridle, A. R., Leef, M. J., and Nowak, B. F. (2013). Evaluation of fixation methods for demonstration of Neoparamoeba perurans infection in Atlantic salmon, Salmo salar L., gills. J. Fish Dis. 36, 831–839.
Caipang, C. M., Lazado, C. C., Brinchmann, M. F., and Kiron, V. (2010). Infection-induced changes in expression of antibacterial and cytokine genes in the gill epithelial cells of Atlantic cod, Gadus morhua during incubation with bacterial pathogens. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 156, 319–325. doi: 10.1016/j.cbpb.2010.04.009
Chen, K., Zhou, X. Q., Jiang, W. D., Wu, P., Liu, Y., Jiang, J., et al. (2019). Dietary phosphorus deficiency caused alteration of gill immune and physical barrier function in the grass carp (Ctenopharyngodon idella) after infection with Flavobacterium columnare. Aquaculture 506, 1–13. doi: 10.1016/j.aquaculture.2019.03.018
Defo, M. A., Douville, M., Giraudo, M., Brodeur, P., Boily, M., and Houde, M. (2018). RNA-sequencing to assess the health of wild yellow perch (Perca flavescens) populations from the St. Lawrence River, Canada. Environ. Pollut. 243, 1657–1668. doi: 10.1016/j.envpol.2018.09.133
Dong, C., Kong, S., Zheng, X., Zhang, J., Nie, G., Li, X., et al. (2019). Genome-wide identification of interleukin-17 (IL17) in common carp (Cyprinus carpio) and its expression following Aeromonas hydrophila infection. Gene 686, 68–75. doi: 10.1016/j.gene.2018.10.038
Downes, J. K., Yatabe, T., Marcos-Lopez, M., Rodger, H. D., MacCarthy, E., O’Connor, I., et al. (2018). Investigation of co-infections with pathogens associated with gill disease in Atlantic salmon during an amoebic gill disease outbreak. J. Fish Dis. 41, 1217–1227. doi: 10.1111/jfd.12814
Du, L., Feng, S., Yin, L., Wang, X., Zhang, A., Yang, K., et al. (2015). Identification and functional characterization of grass carp IL-17A/F1: an evaluation of the immunoregulatory role of teleost IL-17A/F1. Dev. Comp. Immunol. 51, 202–211. doi: 10.1016/j.dci.2015.03.014
English, C. J., Tyml, T., Botwright, N. A., Barnes, A. C., Wynne, J. W., Lima, P. C., et al. (2019). A diversity of amoebae colonise the gills of farmed Atlantic salmon (Salmo salar) with amoebic gill disease (AGD). Eur. J. Protistol. 67, 27–45. doi: 10.1016/j.ejop.2018.10.003
Evans, D. H., Piermarini, P. M., and Choe, K. P. (2005). The multifunctional fish gill: dominant site of gas exchange, osmoregulation, acid-base regulation, and excretion of nitrogenous waste. Physiol Rev. 85, 97–177. doi: 10.1152/physrev.00050.2003
Fogarty, E., and Bergmann, A. (2017). Killers creating new life: caspases drive apoptosis-induced proliferation in tissue repair and disease. Cell Death Differ. 24, 1390–1400. doi: 10.1038/cdd.2017.47
Gill Health Initiative (2015). Gill Diseases in Scotland by Iain Berrill (SSPO). Available online at: https://www.fishhealth.ie/fhu/health-surveillance/aquaplan-fish-health-management-ireland/gill-health-initiative (accessed February 11, 2020).
Girouard, H., Wang, G., Gallo, E. F., Anrather, J., Zhou, P., Pickel, V. M., et al. (2009). NMDA receptor activation increases free radical production through nitric oxide and NOX2. J. Neurosci. 29, 2545–2552. doi: 10.1523/jneurosci.0133-09.2009
Gjessing, M. C., Steinum, T., Olsen, A. B., Lie, K. I., Tavornpanich, S., Colquhoun, D. J., et al. (2019). Histopathological investigation of complex gill disease in sea farmed Atlantic salmon. PLoS One 14:e0222926. doi: 10.1371/journal.pone.0222926
Gjessing, M. C., Thoen, E., Tengs, T., Skotheim, S. A., and Dale, O. B. (2017). Salmon gill poxvirus, a recently characterized infectious agent of multifactorial gill disease in freshwater- and seawater-reared Atlantic salmon. J. Fish Dis. 40, 1253–1265. doi: 10.1111/jfd.12608
Gunnarsson, G. S., Karlsbakk, E., Blindheim, S., Plarre, H., Imsland, A. K., Handeland, S., et al. (2017). Temporal changes in infections with some pathogens associated with gill disease in farmed Atlantic salmon (Salmo salar L). Aquaculture 468, 126–134. doi: 10.1016/j.aquaculture.2016.10.011
Herrero, A., Thompson, K. D., Ashby, A., Rodger, H. D., and Dagleish, M. P. (2018). Complex gill disease: an emerging syndrome in farmed Atlantic salmon (Salmo salar L.). J. Comp. Pathol. 163, 23–28. doi: 10.1016/j.jcpa.2018.07.004
Hoffmann, B., Beer, M., Reid, S. M., Mertens, P., Oura, C. A. L., van Rijn, P. A., et al. (2009). A review of RT-PCR technologies used in veterinary virology and disease control: sensitive and specific diagnosis of five livestock diseases notifiable to the World Organisation for Animal Health. Vet. Microbiol. 139, 1–23. doi: 10.1016/j.vetmic.2009.04.034
Houde, A. L. S., Akbarzadeh, A., Günther, O. P., Li, S., Patterson, D. A., Farrell, A. P., et al. (2019). Salmonid gene expression biomarkers indicative of physiological responses to changes in salinity and temperature, but not dissolved oxygen. J. Exp. Biol. 222:jeb198036. doi: 10.1242/jeb.198036
Houston, R. D., and Macqueen, D. J. (2019). Atlantic salmon (Salmo salar L.) genetics in the 21st century: taking leaps forward in aquaculture and biological understanding. Anim. Genet. 50, 3–14. doi: 10.1111/age.12748
Jokinen, C. C., Edge, T. A., Koning, W., Laing, C. R., Lapen, D. R., Miller, J., et al. (2012). Spatial and temporal drivers of zoonotic pathogen contamination of an agricultural watershed. J. Environ. Qual. 41, 242–252. doi: 10.2134/jeq2011.0203
Kintner, A., and Brierley, A. S. (2019). Cryptic hydrozoan blooms pose risks to gill health in farmed North Atlantic salmon (Salmo salar). J. Mar. Biol. Assoc. U.K. 99, 539–550. doi: 10.1017/s002531541800022x
Koppang, E. O., Kvellestad, A., and Fischer, U. (2015). “Fish mucosal immunity: gill,” in Mucosal Health in Aquaculture, eds B. H. Beck and E. Peatman (Amsterdam: Academic Press), 93–133. doi: 10.1016/b978-0-12-417186-2.00005-4
Król, E., Douglas, A., Tocher, D. R., Crampton, V. O., Speakman, J. R., Secombes, C. J., et al. (2016). Differential responses of the gut transcriptome to plant protein diets in farmed Atlantic salmon. BMC Genomics 17:156. doi: 10.1186/s12864-016-2473-0
Kvellestad, A., Falk, K., Nygaard, S. M., Flesja, K., and Holm, J. A. (2005). Atlantic salmon paramyxovirus (ASPV) infection contributes to proliferative gill inflammation (PGI) in seawater-reared Salmo salar. Dis. Aquat. Organ. 67, 47–54. doi: 10.3354/dao067047
Laurin, E., Jaramillo, D., Vanderstichel, R., Ferguson, H., Kaukinen, K. H., Schulze, A. D., et al. (2019). Histopathological and novel high-throughput molecular monitoring data from farmed salmon (Salmo salar and Oncorhynchus spp.) in British Columbia, Canada, from 2011–2013. Aquaculture 499, 220–234. doi: 10.1016/j.aquaculture.2018.08.072
LeBlanc, F., Ditlecadet, D., Arseneau, J. R., Steeves, R., Boston, L., Boudreau, P., et al. (2019). Isolation and identification of a novel salmon gill poxvirus variant from Atlantic salmon in Eastern Canada. J. Fish Dis. 42, 315–318. doi: 10.1111/jfd.12922
Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656
Lien, S., Koop, B. F., Sandve, S. R., Miller, J. R., Kent, M. P., Nome, T., et al. (2016). The Atlantic salmon genome provides insights into rediploidization. Nature 533, 200–205. doi: 10.1038/nature17164
Marcos-López, M., Calduch-Giner, J. A., Mirimin, L., MacCarthy, E., Rodger, H. D., O’Connor, I., et al. (2018). Gene expression analysis of Atlantic salmon gills reveals mucin 5 and interleukin 4/13 as key molecules during amoebic gill disease. Sci. Rep. 8:13689. doi: 10.1038/s41598-018-32019-8
Matthews, C. G. G., Richards, R. H., Shinn, A. P., and Cox, D. I. (2013). Gill pathology in Scottish farmed Atlantic salmon, Salmo salar L., associated with the microsporidian Desmozoon lepeophtherii Freeman et Sommerville, 2009. J. Fish Dis. 36, 861–869.
McBryan, T. L., Anttila, K., Healy, T. M., and Schulte, P. M. (2013). Responses to temperature and hypoxia as interacting stressors in fish: implications for adaptation to environmental change. Integr. Comp. Biol. 53, 648–659. doi: 10.1093/icb/ict066
McCormick, S. D., and Regish, A. M. (2018). Effects of ocean acidification on salinity tolerance and seawater growth of Atlantic salmon Salmo salar smolts. J. Fish Biol. 93, 560–566. doi: 10.1111/jfb.13656
McIntyre, J. K., Lundin, J. I., Cameron, J. R., Chow, M. I., Davis, J. W., Incardona, J. P., et al. (2018). Interspecies variation in the susceptibility of adult Pacific salmon to toxic urban stormwater runoff. Environ. Pollut. 238, 196–203. doi: 10.1016/j.envpol.2018.03.012
Mitchell, S. O., Baxter, E. J., Holland, C., and Rodger, H. D. (2012). Development of a novel histopathological gill scoring protocol for assessment of gill health during a longitudinal study in marine-farmed Atlantic salmon (Salmo salar). Aquacult. Int. 20, 813–825. doi: 10.1007/s10499-012-9504-x
Mitchell, S. O., Steinum, T. M., Toenshoff, E. R., Kvellestad, A., Falk, K., Horn, M., et al. (2013). ‘Candidatus Branchiomonas cysticola’ is a common agent of epitheliocysts in seawater-farmed Atlantic salmon Salmo salar in Norway and Ireland. Dis. Aquat. Organ. 103, 35–43. doi: 10.3354/dao02563
Mobadersany, P., Yousefi, S., Amgad, M., Gutman, D. A., Barnholtz-Sloan, J. S., Velázquez Vega, J. E., et al. (2018). Predicting cancer outcomes from histology and genomics using convolutional networks. Proc. Natl. Acad. Sci. U.S.A. 15, E2970–E2979. doi: 10.1073/pnas.1717139115
Monte, M. M., Wang, T., Holland, J. W., Zou, J., and Secombes, C. J. (2013). Cloning and characterisation of rainbow trout interleukin-17A/F2 (IL-17A/F2) and IL-17 receptor A: expression during infection and bioactivity of recombinant IL-17A/F2. Infect. Immun. 81, 340–353. doi: 10.1128/iai.00599-12
Morgenroth, D., Ekström, A., Hjelmstedt, P., Gräns, A., Axelsson, M., and Sandblom, E. (2019). Hemodynamic responses to warming in euryhaline rainbow trout: implications of the osmo-respiratory compromise. J. Exp. Biol. 222:jeb207522. doi: 10.1242/jeb.207522
Morrison, R. N., Cooper, G. A., Koop, B. F., Rise, M. L., Bridle, A. R., Adams, M. B., et al. (2006). Transcriptome profiling the gills of amoebic gill disease (AGD)-affected Atlantic salmon (Salmo salar L.): a role for tumor suppressor p53 in AGD pathogenesis? Physiol. Genomics 26, 15–34. doi: 10.1152/physiolgenomics.00320.2005
Morrison, R. N., Zou, J., Secombes, C. J., Scapigliati, G., Adams, M. B., and Nowak, B. F. (2007). Molecular cloning and expression analysis of tumour necrosis factor-α in amoebic gill disease (AGD)-affected Atlantic salmon (Salmo salar L.). Fish Shellfish Immunol. 23, 1015–1031. doi: 10.1016/j.fsi.2007.04.003
Mueller, P. P., Drynda, A., Goltz, D., Hoehn, R., Hauser, H., and Peuster, M. (2009). Common signatures for gene expression in postnatal patients with patent arterial ducts and stented arteries. Cardiol. Young 19, 352–359. doi: 10.1017/S1047951109004260
Noguera, P., Olsen, A. B., Hoare, J., Lie, K. I., Marcos-López, M., Poppe, T. T., et al. (2019). Complex gill disorder (CGD): a histopathology workshop report. Bull. Eur. Assoc. Fish Pathol. 39, 172–176.
Nye, J. A., Link, J. S., Hare, J. A., and Overholtz, W. J. (2009). Changing spatial distribution of fish stocks in relation to climate and population size on the Northeast United States continental shelf. Mar. Ecol. Prog. Ser. 393, 111–129. doi: 10.3354/meps08220
Okamoto, T., Gohil, K., Finkelstein, E. I., Bove, P., Akaike, T., and van der Vliet, A. (2004). Multiple contributing roles for NOS2 in LPS-induced acute airway inflammation in mice. Am. J. Physiol. Lung Cell. Mol. Physiol. 286, L198–L209.
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). Vegan: Community Ecology Package. R package version 2.5-6. Available online at: https://CRAN.R-project.org/package=vegan (accessed February 11, 2020).
Olsen, M. M., Kania, P. W., Heinecke, R. D., Skjoedt, K., Rasmussen, K. J., and Buchmann, K. (2011). Cellular and humoral factors involved in the response of rainbow trout gills to Ichthyophthirius multifiliis infections: molecular and immunohistochemical studies. Fish Shellfish Immunol. 30, 859–869. doi: 10.1016/j.fsi.2011.01.010
Onukwufor, J. O., and Wood, C. M. (2018). The osmorespiratory compromise in rainbow trout (Oncorhynchus mykiss): the effects of fish size, hypoxia, temperature and strenuous exercise on gill diffusive water ?uxes and sodium net loss rates. Comp. Biochem. Physiol. A 219–202, 10–18. doi: 10.1016/j.cbpa.2018.02.002
Pennacchi, Y., Adams, M. B., Nowak, B. F., and Bridle, A. R. (2016). Immune gene expression in the gills of Atlantic salmon (Salmo salar L.) following experimental reinfection with Neoparamoeba perurans. Aquaculture 464, 410–419. doi: 10.1016/j.aquaculture.2016.07.025
Philis, G., Ziegler, F., Gansel, L. C., Jansen, M. D., Gracey, E. O., and Stene, A. (2019). Comparing life cycle assessment (LCA) of salmonid aquaculture production systems: status and perspectives. Sustainability 11:2517. doi: 10.3390/su11092517
Pleić, I. L., Bušelić, I., Trumbić, Ž., Bočina, I., Šprung, M., and Mladineo, I. (2015). Expression analysis of the Atlantic bluefin tuna (Thunnus thynnus) pro-inflammatory cytokines, IL-1β, TNFα1 and TNFα2 in response to parasites Pseudocycnus appendiculatus (Copepoda) and Didymosulcus katsuwonicola (Digenea). Fish Shellfish Immunol. 45, 946–954. doi: 10.1016/j.fsi.2015.06.008
Powell, M. D., Harris, J. O., Carson, J., and Hill, J. V. (2005). Effects of gill abrasion and experimental infection with Tenacibaculum maritimum on the respiratory physiology of Atlantic salmon Salmo salar affected by amoebic gill disease. Dis. Aquat. Org. 63, 169–174. doi: 10.3354/dao063169
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Robledo, D., Hamilton, A., Gutiérrez, A. P., Bron, J. E., and Houston, R. D. (2019). Characterising the mechanisms underlying genetic resistance to amoebic gill disease in Atlantic salmon using RNA sequencing. bioRxiv [Preprint]. doi: 10.1101/699561
Rohart, F., Eslami, A., Matigian, N., Bougeard, S., and Lê Cao, K. A. (2017). MINT: a multivariate integrative method to identify reproducible molecular signatures across independent experiments and platforms. BMC Bioinformatics 18:128. doi: 10.1186/s12859-017-1553-8
Rojas, V., Camus-Guerra, H., Guzman, F., and Mercado, L. (2015). Pro-inflammatory caspase-1 activation during the immune response in cells from rainbow trout Oncorhynchus mykiss (Walbaum 1792) challenged with pathogen-associated molecular patterns. J. Fish Dis. 38, 993–1003. doi: 10.1111/jfd.12315
Ronza, P., Robledo, D., Bermúdez, R., Losada, A. P., Pardo, B. G., Martínez, P., et al. (2019). Integrating genomic and morphological approaches in fish pathology research: the case of turbot (Scophthalmus maximus) enteromyxosis. Front. Genet. 10:26. doi: 10.3389/fgene.2019.00026
Santana, P. A., Guzmán, F., Forero, J. C., Luna, O. F., and Mercado, L. (2016). Hepcidin, cathelicidin-1 and IL-8 as immunological markers of responsiveness in early developmental stages of rainbow trout. Dev. Comp. Immunol. 62, 48–57. doi: 10.1016/j.dci.2016.04.014
Shinn, A. P., Pratoomyot, J., Bron, J. E., Paladini, G., Brooker, E. E., and Brooker, A. J. (2015). Economic costs of protistan and metazoan parasites to global mariculture. Parasitology 142, 196–270. doi: 10.1017/s0031182014001437
Soleng, M., Johansen, L. H., Johnsen, H., Johansson, G. S., Breiland, M. W., Rormark, L., et al. (2019). Atlantic salmon (Salmo salar) mounts systemic and mucosal stress responses to peracetic acid. Fish Shellfish Immunol. 93, 895–903. doi: 10.1016/j.fsi.2019.08.048
Sollid, J., De Angelis, P., Gundersen, K., and Nilsson, G. (2003). Hypoxia induced adaptive and reversible gross morphological changes in crucian carp gills. J. Exp. Biol. 206, 3667–3673. doi: 10.1242/jeb.00594
Song, Y., Salbu, B., Teien, H. C., Sørlie Heier, L., Rosseland, B. O., Høgåsen, T., et al. (2014). Hepatic transcriptomic profiling reveals early toxicological mechanisms of uranium in Atlantic salmon (Salmo salar). BMC Genomics 15:694. doi: 10.1186/1471-2164-15-694
Steinum, T., Kvellestad, A., Colguhoun, D. J., Heum, M., Mohammad, S., Grontvedt, R. N., et al. (2010). Microbial and pathological findings in farmed Atlantic salmon Salmo salar with proliferative gill inflammation. Dis. Aquat. Org. 91, 201–211. doi: 10.3354/dao02266
Sun, F. Y., Peatman, E., Li, C., Liu, S. K., Jiang, Y. L., Zhou, Z. C., et al. (2012). Transcriptomic signatures of attachment, NF-kappa B suppression and IFN stimulation in the catfish gill following columnaris bacterial infection. Dev. Comp. Immunol. 38, 169–180. doi: 10.1016/j.dci.2012.05.006
Syahputra, K., Kania, P. W., Al-Jubury, A., Marnis, H., Setyawan, A. C., and Buchmann, K. (2019). Differential immune gene response in gills, skin, and spleen of rainbow trout Oncorhynchus mykiss infected by Ichthyophthirius multifiliis. PLoS One 14:e0218630. doi: 10.1371/journal.pone.0218630
Tisoncik, J. R., Korth, M. J., Simmons, C. P., Farrar, J., Martin, T. R., and Katze, M. G. (2012). Molecular innate immunity in teleost fish: review and future perspectives. Into the eye of the cytokine storm. Microbiol. Mol. Biol. Rev. 76, 16–32. doi: 10.1128/MMBR.05015-11
Valdenegro-Vega, V. A., Cook, M., Crosbie, P., Bridle, A. R., and Nowak, B. F. (2015). Vaccination with recombinant protein (r 22C03), a putative attachment factor of Neoparamoeba perurans, against AGD in Atlantic salmon (Salmo salar) and implications of a co-infection with Yersinia ruckeri. Fish Shellfish Immunol. 44, 592–602. doi: 10.1016/j.fsi.2015.03.016
Valdenegro-Vega, V. A., Crosbie, P., Bridle, A., Leef, M., Wilson, R., and Nowak, B. F. (2014). Differentially expressed proteins in gill and skin mucus of Atlantic salmon (Salmo salar) affected by amoebic gill disease. Fish Shellfish Immunol. 40, 69–77. doi: 10.1016/j.fsi.2014.06.025
Valenzuela-Miranda, D., Boltaña, S., Cabrejos, M. E., Yáñez, J. M., and Gallardo-Escárate, C. (2015). High-throughput transcriptome analysis of ISAV-infected Atlantic salmon Salmo salar unravels divergent immune responses associated to head-kidney, liver and gills tissues. Fish Shellfish Immunol. 45, 367–377. doi: 10.1016/j.fsi.2015.04.003
Wise, D. J., Griffin, M. J., Terhune, J. S., Pote, L. M., and Khoo, L. H. (2008). Induction and evaluation of proliferative gill disease in channel catfish fingerlings. J. Aquat. Anim. Health 20, 236–244. doi: 10.1577/H08-023.1
Wynne, J. W., O’Sullivan, M. G., Cook, M. T., Stone, G., Nowak, B. F., Lovell, D. R., et al. (2008a). Transcriptome analyses of amoebic gill disease-affected Atlantic salmon (Salmo salar) tissues reveal localized host gene suppression. Mar. Biotechnol. 10, 388–403. doi: 10.1007/s10126-007-9075-4
Wynne, J. W., O’Sullivan, M. G., Stone, G., Cook, M. T., Nowak, B. F., Lovell, D. R., et al. (2008b). Resistance to amoebic gill disease (AGD) is characterised by the transcriptional dysregulation of immune and cell cycle pathways. Dev. Comp. Immunol. 32, 1539–1560. doi: 10.1016/j.dci.2008.05.013
Xu, Z., Takizawa, F., Casadei, E., Shibasaki, Y., Ding, Y., Sauters, T. J. C., et al. (2020). Specialization of mucosal immunoglobulins in pathogen control and microbiota homeostasis occurred early in vertebrate evolution. Sci. Immunol. 5:eaay3254. doi: 10.1126/sciimmunol.aay3254
Young, N. D., Crosbie, P. B., Adams, M. B., Nowak, B. F., and Morrison, R. N. (2007). Neoparamoeba perurans n. sp., an agent of amoebic gill disease of Atlantic salmon (Salmo salar). Int. J. Parasitol. 37, 1469–1481. doi: 10.1016/j.ijpara.2007.04.018
Zilberg, D., Gross, A., and Munday, B. L. (2001). Production of salmonid amoebic gill disease by exposure to Paramoeba sp. harvested from the gills of infected fish. J. Fish Dis. 24, 79–82. doi: 10.1046/j.1365-2761.2001.00271.x
Keywords: proliferative gill disease, gene expression, RNA-seq, immune response, gill inflammation, aquaculture, climate change
Citation: Król E, Noguera P, Shaw S, Costelloe E, Gajardo K, Valdenegro V, Bickerdike R, Douglas A and Martin SAM (2020) Integration of Transcriptome, Gross Morphology and Histopathology in the Gill of Sea Farmed Atlantic Salmon (Salmo salar): Lessons From Multi-Site Sampling. Front. Genet. 11:610. doi: 10.3389/fgene.2020.00610
Received: 29 February 2020; Accepted: 19 May 2020;
Published: 19 June 2020.
Edited by:Dong Xia, Royal Veterinary College (RVC), United Kingdom
Reviewed by:Aleksei Krasnov, Nofima, Norway
Mark D. Fast, University of Prince Edward Island, Canada
Copyright © 2020 Król, Noguera, Shaw, Costelloe, Gajardo, Valdenegro, Bickerdike, Douglas and Martin. 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.