Characterising the mechanisms underlying genetic resistance to amoebic gill disease in Atlantic salmon using RNA sequencing

Abstract

Background

Gill health is one of the main concerns for Atlantic salmon aquaculture, and Amoebic Gill Disease (AGD), attributable to infection by the amoeba Neoparamoeba perurans, is a frequent cause of morbidity. In the absence of preventive measures, increasing genetic resistance of salmon to AGD via selective breeding can reduce the incidence of the disease and mitigate gill damage. Understanding the mechanisms leading to AGD resistance and the underlying causative genomic features can aid in this effort, while also providing critical information for the development of other control strategies. AGD resistance is considered to be moderately heritable, and several putative QTL have been identified. The aim of the current study was to improve understanding of the mechanisms underlying AGD resistance, and to identify putative causative genomic factors underlying the QTL. To achieve this, RNA was extracted from the gill and head kidney of AGD resistant and susceptible animals following a challenge with N. perurans, and sequenced.

Results

Comparison between resistant and susceptible animals primarily highlighted differences mainly in the local immune response in the gill, involving red blood cell genes and genes related to immune function and cell adhesion. Differentially expressed immune genes pointed to a contrast in Th2 and Th17 responses, which is consistent with the increased heritability observed after successive challenges with the amoeba. Five QTL-region candidate genes showed differential expression, including a gene connected to interferon responses (GVINP1), a gene involved in systemic inflammation (MAP4K4), and a positive regulator of apoptosis (TRIM39). Analyses of allele-specific expression highlighted a gene in the QTL region on chromosome 17, cellular repressor of E1A-stimulated genes 1 (CREG1), showing allelic differential expression suggestive of a cis-acting regulatory variant.

Conclusions

In summary, this study provides new insights into the mechanisms of resistance to AGD in Atlantic salmon, and highlights candidate genes for further functional studies that can further elucidate the genomic mechanisms leading to resistance and contribute to enhancing salmon health via improved genomic selection.

Background

Gill health is currently one of the major concerns for Atlantic salmon farming worldwide. Fish gills are multifunctional organs fundamental for gas exchange, ionoregulation, osmoregulation, acid-base balance and ammonia excretion, but also play an important role in hormone production and immune defence [1]. Gills are constantly exposed to the marine environment, and are often the first line of defence against pathogens. Gill damage is often observed in Atlantic salmon under farming conditions, and can pose a significant welfare, management and economic burden. While the aetiology of gill disorders is complex, Amoebic Gill Disease (AGD) is currently regarded as a key threat to gill health [2, 3]. This disease adversely affects the gill, and can result in respiratory distress, and ultimately mortality if left untreated. Initially limited to Tasmania, AGD is currently causing major economic and fish welfare burden to Norwegian, Scottish and Australian salmon aquaculture [4]. The causative agent of this disease is the amoeba Neoparamoeba perurans, an opportunistic pathogen that typically requires with expensive and laborious fresh water or hydrogen peroxide treatments [5], and there are currently very limited opportunities for prevention.

A promising avenue to decrease the incidence of AGD in farmed Atlantic salmon is to increase genetic resistance of aquaculture stocks to N. perurans. There is significant genetic variation in resistance to AGD in commercial Atlantic salmon populations [6,7,8,9], therefore selective breeding has potential to improve gill health via a reduction in amoebic load and associated gill damage. The use of genetic markers through genomic selection can expedite genetic gain in aquaculture breeding programmes (e.g. [8, 10,11,12]), however, the need to genotype a large number of animals and to perform disease challenges in every generation involves a relatively high cost. The discovery of the mechanisms leading to resistance and the underlying causative genetic variants has the potential to reduce this cost via incorporation of functional SNPs into the genomic prediction models.

Discovering the genes and pathways that lead to successful immune responses to pathogens is a major goal in genetics and immunology research. Understanding disease resistance can aid selective breeding via incorporation of putative causative variants with greater weighting in genomic prediction models, which can improve selection accuracy and reduce the need for routine trait recording [13, 14]. Such information can also inform the development of improved disease challenge models, and more successful prevention or treatment strategies through an increased knowledge of host-pathogen interactions. Finally, with the potential role for targeted genome editing (e.g. using CRISPR/Cas9) in future aquaculture breeding programmes, understanding the functional mechanisms underlying disease resistance traits is key to identifying target genes and variants for editing [15]. Previous studies into AGD-infected Atlantic salmon have suggested that the amoebae might elicit an immunosuppressive effect on the innate response of the host, with a concurrent up-regulation of the adaptive Th2-mediated response [16,17,18]. Th2 cytokines were also found consistently up-regulated when comparing AGD infected and non-infected samples, and lesion and non-lesion areas [18]. The heritability of resistance to AGD has been shown to increase after successive cycles of disease challenge / treatment [7], which could suggest that the ability to elicit a successful adaptive immune response is partly under genetic control. Finally, a higher expression of genes related to adaptive immunity has been previously reported in more AGD-resistant salmon compared to their more susceptible counterparts using a microarray approach to measure gene expression [19].

In a previous study by our team, several QTL regions with a significant contribution to genetic AGD resistance were identified in Atlantic salmon derived from a commercial breeding programme [8]. In the current study, the gill and head kidney transcriptomes of AGD resistant and susceptible Atlantic salmon from the same population were sequenced and compared. The main goals of the study were a) to assess the differences in local and systemic immune responses between AGD resistant and susceptible Atlantic salmon, and b) to use gene expression data to identify positional and functional candidate genes underlying the previously detected resistance QTL.

Results

Sampling and sequencing

Fish were classified into resistant or susceptible based on their mean gill score and their gill amoebic load. A previous study by our group has shown a high positive genetic correlation between these two traits (higher gill score associated with higher amoebic load), and both are considered indicator traits for resistance to AGD. RNA sequencing (RNA-Seq) was performed on the gill and head kidney of 12 resistant and 12 susceptible fish. Resistant animals had a mean gill score of 2.92 ± 0.13, mean amoebic load (qPCR ct value, high ct value corresponds to low amoebic loads) of 37.12 ± 3.63 and mean weight of 543 ± 116 g at the point of sampling; susceptible animals had a mean gill score of 4.12 ± 0.20, mean amoebic load of 25.99 ± 1.80 and mean weight of 409 ± 96 g. Sequencing of one of the gill samples rendered an extremely low number of reads and therefore was discarded. The remaining samples had an average of 24 M filtered paired-end reads, which were pseudoaligned to transcripts (determine, for each read, which transcript it is compatible with) in the Atlantic salmon genome (ICSASG_v2; Genbank accession GCF_000233375.1 [20]). Exploratory analyses based on distance measures revealed two head kidney samples as outliers and they were removed (Additional file 1). Therefore, the final dataset comprised of 23 gill and 22 head kidney samples from 24 individuals. The two organs showed clearly distinct patterns of gene expression, as would be expected. The difference in global gene expression pattern between resistant and susceptible samples in both tissues was much less pronounced, but still evident in the gill in particular (Fig. 1). Similar results were described in a Norwegian commercial population [9].

Fig. 1

Principal component analysis. RNA-Seq samples clustered according to their gene expression. The larger symbols represent group means, and ellipses represent 95% confidence intervals for the groups

Differential expression

A total of 115 and 42 differentially expressed transcripts (following multiple-testing correction, false discovery rate (FDR) p-value < 0.05) were detected between resistant and susceptible samples in gill and head kidney tissues respectively (Fig. 2, Additional file 2). The clearest evidence for differential immune responses was found in gill, where several differentially expressed immune-related transcripts were detected. Most differentially expressed transcripts in head kidney were not obviously related to AGD or disease resistance. To gain an overall view of the results, a Gene Ontology (GO) enrichment test was performed in both gill and head kidney for sets of differentially expressed transcripts according to three different significance criteria (p-value < 0.01, 0.05 and 0.1) (Fig. 3). In the gill, various relevant GO terms were observed, such as “Response to stress”, “Cytoskeleton” and “Circulatory system process”. A larger number of enriched GO terms were observed in head kidney. While most of them cannot be directly connected to AGD related responses (i.e. “cell proliferation” or “kinase activity”), terms such as “response to stress” or “protein modification process” were observed. For instance, of 22 genes showing p-values < 0.01, 15 of them were assigned to “Response to stress”. Similar analyses for KEGG pathways did not reveal any significant enrichment.

Fig. 2

Heatmap of differentially expressed genes between resistant and susceptible samples. Heatmaps of all differentially expressed genes in gill (a) and head kidney (b). Samples and genes were clustered according to gene expression (mean centered and scaled normalized counts)

Fig. 3

Gene Ontology enrichment for differentially expressed genes. GO enrichment is shown for all differentially expressed genes in gill and head kidney according to three different significant criteria (FDR p-value < 0.1, 0.05 and 0.01). The height of the bars represents fold enrichment (percentage of genes assigned to the GO term in the set of differentially expressed genes compared to the percentage assigned to that GO term in the transcriptome of that tissue)

Detailed inspection of the differentially expressed transcripts in the gill revealed that they can be grouped into three broad categories concordant with GO enrichment results: 1) immune response (“Response to stress”), 2) red blood cells and coagulation (“Circulatory system process”) and 3) cell adhesion or cell shape (“Cytoskeleton”).

Amongst the immune-related transcripts showing differential expression in the gill was interleukin-17 receptor E (IL17RE), which was highly expressed in resistant animals (Log₂ fold change value - logFC = 1.1). In mice IL17RE is the receptor for IL-17C, which has an essential role in host mucosal defense against infection and is critical for a successful immune response against bacterial infection [21]. The IL-17C – IL17RE pair also stimulates T-helper cell 17 responses, which has a proinflammatory effect [22]. IL-17C expression was also shown to have a negative correlation with amoebic load in a previous study of Atlantic salmon, and the Th17 pathway in general was found to be significantly down-regulated in response to AGD [16]. This could be a mechanism of immune evasion elicited by the parasite, which might be more effective in susceptible fish than resistant. Another highly expressed transcript is involved in T-cell function, T-cell specific surface glycoprotein CD28 (CD28; logFC = 1.60). CD28 promotes T-cell survival and proliferation, and enhances the production of multiple cytokines including IL4 [23] IL4 has been found to be up-regulated in response to AGD [15], and this gene induces differentiation of naïve helper T cells to Th2 cells. The Th2 pathway was found to be up-regulated in late stages of AGD [16]. This pathway is linked to humoral immune responses against extracellular parasites and to tissue repair [24], and therefore is an expected response to AGD. A higher prevalence of this type of response in resistant animals would also be consistent with the observed increase of the heritability of resistance after successive cycles of disease challenge / treatment [7], reflecting genetic variability in the effectiveness of the adaptive response, and / or variation in immune memory.

Several genes connected to red blood cells were found to be differentially expressed, including five different haemoglobin subunit transcripts, which were highly expressed and clearly up-regulated in resistant samples in the gill (logFC ~ 2). Haemoglobin α and β subunits have been previously found down-regulated in AGD lesions at the transcript [25] and protein level [26], and reduced hematocrit has been described in AGD infected Atlantic salmon, linked mainly to gill damage [27]. However, it has also been suggested that this haemoglobin dysregulation could be related to antimicrobial peptides derived from haemoglobin β [26], which have been described to have parasiticidal properties in channel catfish [28, 29]. The plasma protease C1 inhibitor gene (SERPING1) was also up-regulated (logFC = 1.2). This gene inhibits the complement system and also has anti-inflammatory functions [30]. Complement proteins have been found in gill mucus of AGD infected Atlantic salmon [31]. The lower expression of SERPING1 in susceptible samples might simply be a reflection of the higher extent of gill damage in these animals, requiring activation of the complement system and increase of local inflammatory responses.

There are also a few differentially expressed transcripts connected to cell adhesion and cell shape, including a cadherin gene (cadherin-related family member 5; logFC = 4.5) and an actin related gene (actin filament associated protein 1-like 1; logFC = 1.3). Another cadherin gene (Cadherin 1) was previously found dysregulated in response to AGD, along with two additional cell adhesion related genes [25]. The Cdc42 effector protein 2 (CDC4EP2; logFC = 0.6) was also up-regulated in resistant fish, and has been associated with roles in actin filament assembly and control of cell shape [32]. A previous study identified an enrichment of cell-adhesion genes in severely affected animals compared to others with healthier gills infected by AGD [9]. These changes are consistent with the epithelial hyperplasia and other structural changes caused by the parasite in the gill of infected animals [26].

In head kidney, most of the differentially expressed (DE) transcripts are seemingly unconnected to biological processes that have previously been related to AGD. Tumor necrosis factor alpha-induced protein 8-like protein 1 (TNFAIP8L1; logFC = − 0.9) was more highly expressed in susceptible samples. This gene inhibits apoptosis by suppressing the activity of caspase-8 [33]. The down-regulation of pro-apoptotic genes has been connected to AGD severity [18], however previous studies have not found up-regulation of tumour necrosis factor-alpha (TNA-α), which could potentially suggest some immunomodulation mechanism from the parasite [34]. Nonetheless, the lack of a clear picture in head kidney might reflect the relative importance of the local and systemic immune responses in response to AGD. Previous studies have found that the transcriptomic differences between affected and unaffected gills of AGD infected salmon are similar to those between affected gills of infected salmon and the gills of healthy salmon, suggesting indeed a localized response to AGD [25].

The regulation of transcripts upon infection is a strong indication of the involvement of the gene product in the immune and physiological response of the host to the pathogen, but a comparison between resistant and susceptible animals can offer insight into the mechanisms determining the success of the immune response against the pathogen. The main caveat of this approach is that it is difficult to distinguish cause and consequence, i.e. is the gene differentially expressed because it confers resistance or due to differential disease progression? Additional evidence, such as the co-localization of differentially expressed genes with QTL or the identification of cis regulatory variants in the QTL regions can further contribute to understanding of the mechanisms of disease resistance, and discover underlying candidate genes.

Integration with previous QTL

The overlap between previously identified putative QTL regions in this population [8] and the differentially expressed genes was explored (Fig. 4). A differentially expressed gene, interferon-induced very large GTPase 1 (GVINP1), was found in one of the QTL regions of chromosome 18, which explained ~ 20% of the genetic variance in resistance to AGD (second largest QTL). Very little is known about the function of this gene, but it has been shown to respond to both type I and type II interferon response in mammals [35]. The genes showing FDR corrected p-values < 0.1 (a total of 268 genes) were also investigated, and four additional genes were found in these QTL regions. MAP4K4, located in a putative QTL region of chromosome 17, surpassed this threshold, and is involved in systemic inflammation in mammals [36], and TRIM39 in the second QTL region in chromosome 18, a positive regulator of apoptosis [37].

Fig. 4

Differentially expressed genes located in resistance QTL. The location of the QTL regions in the chromosomes are shown in grey. Genes with significance values < 0.05 are in red, those with significance values < 0.1 are in orange. Positive fold changes correspond to higher expression in resistant samples

Allele specific expression

To explore potential cis-acting variation underlying the resistance QTL, an allele specific expression (ASE) test was performed for the SNPs in transcripts within the QTL regions (Fig. 5), finding a significant ASE event in a gene in chromosome 17; cellular repressor of E1A-stimulated genes 1 (CREG1). In humans this protein is connected to the regulation of cellular proliferation and differentiation [38], and antagonizes the proliferative effects of adenovirus E1A protein [39]. This gene showed a log fold change of 0.75 between resistant and susceptible samples in gill (FDR p-value = 0.25). A second significant ASE event was found in an uncharacterised gene in chromosome 18 (LOC106576659; Fig. 5), however this gene showed no differences in fold change between resistant and susceptible samples.

Fig. 5

Allele specific expression CREG1. Barplot showing the read counts for each allele for those SNPs in the QTL regions showing allele specific expression. The two SNPs are located in CREG1 (Chromosome 17–24,545,527 bp) and the uncharacterized gene LOC106576659 (Chromosome 18–57,163,493 bp)

The polygenic nature of resistance to AGD means that different resistance mechanisms might be operating in each different family. The connection between genotypes and expression, through expression QTL (eQTL) or ASE, can provide strong evidence for functional candidate genes underlying QTL. While eQTL studies require a relatively large number of animals, the advantage of ASE is that the statistical test in performed separately in each heterozygous individual. It is well known that most causative variants are part of regulatory elements and affect gene expression [40, 41], therefore the detection of ASE in a QTL can provide strong evidence linking the function of a gene to the QTL and the phenotype of interest.

Discussion

The potential benefits of the identification of causative variation impacting on complex traits are substantial, ranging from fundamental knowledge of the biology underlying the traits of interest to their application for enhancing these traits in farmed populations. However, even with the addition of various layers of information such as RNA sequencing, determining the causative gene underlying a QTL is not straightforward, especially because the QTL regions tend to be large and contain a large number of genes, as previously described for sea lice resistance QTL [42].

Eventually, functional assays are necessary to provide actual evidence of its causality. The advent of CRISPR-CAS9 has made this much more feasible in non-model species. Likewise, this technology now provides the opportunity of using this information to introduce or fix favorable alleles in farmed populations [15]. The genetic architecture of quantitative traits usually varies across populations, and indeed AGD resistance QTL seem to vary across different Atlantic salmon commercial populations [6, 7]. While the use of genome editing in farmed animals requires societal and regulatory changes, the transference of causative variants across populations can lead to a rapid increase in disease resistance [15], with long-lasting effects on animal welfare and food security. Nevertheless, the discovery of causative variants and genes can be used to increase the weight of causative variants in genomic selection, increasing its accuracy and therefore speeding up genetic gain in each generation [13]. More widely, basic knowledge about the pathways leading to resistance to disease can inform drug development or preventive measures such as functional feeds. To summarise, finding the underlying cause(s) of resistance to disease can provide large benefits for aquaculture and society in the form of healthier animals, increased food security and sustainable economic gain, directly through their implementation in breeding schemes in the present and through genome editing in the future. Overlaying genome-wide association studies with gene expression differences between genetically distinct individuals, as performed for resistance to AGD in the current study, is an important step towards identifying these causative mechanisms.

Conclusions

The transcriptomic differences between AGD resistant and AGD susceptible Atlantic salmon are limited, which might not be surprising considering the polygenic nature of the trait. These differences were more evident in the gill than in head kidney, potentially highlighting the importance of the local immune response. Genes involved in immune response (Th2 and Th17 pathways), red blood cells and cell adhesion could be part of the mechanisms leading to AGD resistance, albeit it is difficult to discriminate cause and consequence. The integration of previously discovered QTL and expression data pointed to potential candidate genes of interest, such as GVINP1, MAP4K4 or TRIM39. An additional candidate gene, CREG1, showed allele specific expression in one of the QTL regions. Follow-up studies to investigate the functional role of these genes in the response to AGD could help improve understanding of the molecular mechanism of resistance to this parasite, and contribute to improving gill health in farmed populations through incorporation of functional data to improve genomic prediction, or potentially via genome editing in the longer term.

Methods

Experimental design

The AGD challenge experiment was performed using 797 Atlantic post-smolt salmon from 132 nuclear families (~ 18 months, mean weight after challenge ~ 464 g) originating from a commercial breeding programme (Landcatch, UK). The challenge experiment was performed as described in [8]. In brief, seeder fish with a similar level of AGD infection (similar gill damage) were produced specifically for this study by cohabitation with fish infected from an ongoing in vivo culture. The experimental challenge consisted of three separate cycles of infection, established by cohabitation at a ratio of 15% seeder to naïve fish, with a recovery period after the first two infections [7], using a 4 m³ seawater tank in the experimental facilities of University of Stirling’s Marine Environmental Research Laboratory, Machrihanish (Scotland, UK). The fish were kept under a 16-h light and 8-h dark photoperiod, starting at 05:00; the fish were fed Biomar organic salmon feed, automatic every 20 min to approximately 1% biomass; water supply was ambient flow-through filtered to approximately 90 μm, for the duration of the trial, water temperature was between 13 & 14 °C and salinity was 33–35 ppt. Fish were treated with fresh water 21 days after the start of the two first challenges, and fresh water treatment was followed by a week of recovery and the addition of infected seeder fish after that week. The fish were checked visually four times daily during this period. In the third cycle of infection, the disease was allowed to progress until the sampling point, when fish were terminated by an overdose of anaesthetic (Phenoxyethanol, 0.5 mg/L) followed by destruction of the brain according to UK Home Office regulations. All fish were sampled and phenotyped during three consecutive days. Gill damage was recorded for both gills, and scored from 0 to 5 according to the severity of the lesions [43]. A single operator recorded all the gill lesion scores, and the classification was guided by pictures. Additional operators scored some fish, and the scores never differed by > 0.5. Further, one of the gills was stored in ethanol for amoebic load estimation using qPCR with N. perurans specific primers. Amoebic load has been used as a suitable indicator of resistance to AGD in Atlantic salmon [7]. All fish in this study had phenotypes for mean gill score (mean of the left gill and right gill scores) and amoebic load (qPCR values using N. perurans specific primers, amplified from one of the gills). Twenty-four fish, each from a different full-sib family, were selected for RNA sequencing (Additional file 3) based on high or low levels of resistance according to the measured traits (mean gill score and amoebic load as measured by qPCR). The number of samples was decided based on the availability of animals with extreme phenotypes from different families. The estimated power of the RNA-Seq experiment for this sample size is > 0.8 for genes showing log fold change > |1| and a dispersion of 0.5 [44]. Gill and head kidney samples were obtained from each animal and stored in RNAlater at 4 °C for 24 h, and then at − 20 °C until RNA extraction.

RNA extraction and sequencing

RNA from all the 48 samples was obtained using a standard TRI Reagent RNA extraction protocol. Briefly, approximately 50 mg of tissue was homogenized in 1 ml of TRI Reagent (Sigma, St. Louis, MO) by shaking using 1.4 mm silica beads, followed by the addition of 100 μl of 1-bromo-3-chloropropane (BCP) for phase separation. Afterwards, 500 μl of isopropanol were added for precipitation, which was followed by subsequent washes with 65–75% ethanol. RNA was resuspended in RNAse-free water and treated with Turbo DNAse (Ambion). Qiagen RNeasy Mini kit columns were used to clean up the samples, and RNA integrity was checked on Agilent 2200 Bioanalyzer (Agilent Technologies, USA). RNA-Seq libraries were prepared using the Illumina Truseq mRNA stranded RNA-Seq Library Prep Kit following standard protocols. Library quantity and quality were quantified using the Bioanalyzer 2100 (Agilent), and then sequenced on three lanes of an Illumina Hiseq 4000 as 75 base paired-end reads at the facilities of Edinburgh Genomics (UK). Raw reads were deposited in NCBI’s Sequence Read Archive (SRA) under BioProject accession number PRJNA552604.

Read mapping

Bioinformatic analyses were performed as previously described [45]. Briefly, the quality of the sequencing output was assessed using FastQC v.0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and residual adaptors and low quality (< 20) bases were trimmed using Trimmomatic v.0.38 [46]. Only reads where both pairs were longer than 36 bp post-filtering were retained. STAR v.2.6.1a [47] was used to map the filtered reads to the most recent Atlantic salmon genome assembly (ICSASG_v2; Genbank accession GCF_000233375.1 [20]). Kallisto v0.44.0 [48] and the latest Atlantic salmon genome annotation (NCBI Salmo salar Annotation Release 100) were used to quantify transcript expression.

Differential expression

Differential expression analyses were performed using R v.3.5.2 (https://www.r-project.org/). Kallisto HDF5 binary files were used obtain differential gene expression estimates with the Bioconductor package DESeq2 v.3.4 [49]. Briefly, the ‘median of ratios’ method was used to calculate size factors for each sample, and normalization of count data was carried out to account for differences in library depth. After dispersion estimates were fitted to the mean intensity reduced towards the expected dispersion values, the expression of gene was fitted to a negative binomial model and significance assessed using the Wald test. Genes showing Benjamini-Hochberg false discovery rate (FDR) p-values < 0.05 and log₂ fold change values (logFC) > 0.5 were considered differentially expressed. Prior to differential expression the whole dataset was evaluated using hierarchical clustering and principal component analyses, and outlier samples were removed for downstream analyses. PCA plots were created using the R package factoextra (http://www.sthda.com/english/rpkgs/factoextra/).

Gene Ontology (GO) enrichment analyses were performed in R v.3.5.2 (https://www.r-project.org/) using Bioconductor packages GOstats v.2.48.0 [50] and GSEABase v.1.44.0 [51]. GO term annotation for the Atlantic salmon transcriptome was obtained using the R package Ssa.RefSeq.db v1.3 (https://gitlab.com/cigene/R/Ssa.RefSeq.db). The over-representation of GO terms in differentially expressed gene lists compared to the corresponding transcriptomes (gill or head kidney) was assessed with a hypergeometric test. A GO term was considered enriched if it showed ≥5 DE genes assigned and a p-value < 0.05. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out using KOBAS v3.0.3 [52]. Briefly, salmon genes were annotated against KEGG protein database [53] to determine KEGG Orthology (KO). KEGG enrichment for differentially expressed gene lists was tested by comparison to the whole set of expressed genes in the corresponding tissue using Fisher’s Exact Test. KEGG pathways with ≥5 DE genes assigned and showing a Benjamini-Hochberg FDR corrected p-value < 0.05 were considered enriched for differential expression. The reference tissue transcriptome for both GO and KEGG enrichment comprised only those genes with mean normalized counts value > 5.

Allele specific expression

Gene expression estimates and genotypes obtained from the RNA sequencing were used to investigate allele specific expression. The samtools v1.6 software [54] was used to identify SNPs, and call genotypes for those SNPs in individual samples. PCR duplicates, reads with mapping quality < 20 and bases with phred quality scores < 20 were excluded. SNPs within 5 bp of an indel, with quality < 20, MAF < 0.05 or less than 4 reads supporting the alternative allele were discarded. The putative effect of the SNPs was assessed using the official salmon genome annotation (NCBI Salmo salar Annotation Release 100) and the SnpEff v.4.2 software [55]. Allelic specific expression was assessed using the R package AllelicImbalance [56]. For every SNP in the regions of interest, read counts were obtained for each allele in heterozygous animals, those with less than 10 reads were filtered, and a binomial test was performed to assess the significance of the allelic differences. Only those genomic positions called as heterozygotes in a minimum of 4 and a maximum of 36 (75%) samples (75%) were considered. An allele specific expression event was considered significant if the mean p-value of all heterozygotes was < 0.05. All significant events were manually inspected.