GWAS identifies genetic variants associated with omega-3 fatty acid composition of Atlantic salmon fillets

The objective of this study was to identify genetic variants associated with long-chain omega-3 fatty acid content in Atlantic salmon muscle, in order to identify genes underlying the genetic variation in these traits, by performing a genome-wide association study. Fatty acid composition, including eicosapentaenoic acid (EPA), docosahexaenoic acid (DHA) content and the ratios DHA/DPA (docosapentaenoic acid) and DHA/ ALA (α-linolenic acid) of skeletal muscle was determined in 642 Atlantic salmon individuals from the Salmobreed broodstock program. Further, a 57 K singlenucleotide polymorphism (SNP) array was used to genotype the 642 individuals to search for SNPs associated with skeletal muscle content of individual omega-3 fatty acids as well as ratios between fatty acids, using a mixed model approach. We identified markers showing a significant association with the ratio of DHA to DPA located on chromosome 19, close to the candidate gene elovl2, which is directly involved in the conversion of DPA to DHA. The results further suggested that genetic variation affecting fillet EPA and DHA content is present on Atlantic salmon chromosome 21, as the GWAS analysis of EPA, DHA and DHA/ALA ratio all pointed to this chromosome. No known genes of direct relevance to the omega-3 bioconversion pathway were found here. We discuss the relevance of other interesting genes related to lipid metabolism and health located in this region.


Background
Atlantic salmon (Salmo salar L.) is an important farmed fish species, known for its high content of the health-promoting omega-3 long-chain polyunsaturated fatty acids (n-3 LC PUFA) eicosapentaenoic (EPA; 20:5n-3) and docosahexaenoic acid (DHA; 22:6n-3) in muscle. Levels of n-3 LC PUFA have decreased significantly the last few decades because of the substitution of a large portion of marine ingredients with more sustainable plant ingredients in the diet of farmed Atlantic salmon (FAO, 2016;Sprague et al., 2016;Ytrestoyl et al., 2015). Quantitative genetic analyses have estimated genetic parameters of Atlantic salmon omega-3 fatty acids (FA) and showed the potential of selective breeding to increase n-3 LC PUFA levels in salmon fillets, as they are heritable traits (Horn et al., 2018;Leaver et al., 2011).
In addition to the omega-3 FA bioconversion pathway, other metabolic processes including selective uptake, β-oxidation and deposition of FAs are known to affect muscle FA composition in Atlantic salmon. Vegusdal et al. (2004) found for instance that Atlantic salmon muscle cells take up EPA more rapidly than 18:1n-9. Several studies have observed a higher selective deposition of DHA in muscle when compared to other dietary FAs Torstensen et al., 2004), which is probably due to DHA being a relatively poor substrate for β-oxidation (Sargent et al., 2002). EPA, on the other hand, seems to be metabolized to a larger extent than DHA in muscle (Torstensen et al., 2004). A recent study of Atlantic salmon by our group found that muscle EPA content was negatively associated with expression of genes involved in lipid catabolism and hypothesized that individual differences in EPA content is caused by selective oxidation of FAs (Horn et al., 2019). There is, however, limited knowledge on which genes and metabolic pathways that are most important for determining EPA and DHA content in Atlantic salmon muscle.
Genome-wide association studies (GWAS) are a powerful tool that can identify Single Nucleotide Polymorphisms (SNPs) linked to genes associated with phenotypic traits. GWAS of muscle FA composition have been performed in cattle (Chen et al., 2015;Saatchi et al., 2013;Zhu et al., 2017) and pigs (Ramayo-Caldas et al., 2012), where several significantly associated regions and loci were identified. In cattle, FA composition is influenced by a few host genes with major effects and many genes of smaller effects (Chen et al., 2015). In fish, one study in Asian seabass (Lates calcarifer) identified multiple quantitative trait loci (QTL) for omega-3 FA in muscle and different QTL were associated with different omega-3 FA (Xia et al., 2014), showing the complexity of the genetic architecture of FAs in fish. To increase the understanding of the genetics of muscle n-3 LC PUFA content in Atlantic salmon, GWAS is a relevant approach that has not yet been explored in Atlantic salmon.
The aim of the current study was to identify genetic variants associated with long-chain omega-3 fatty acid content in Atlantic salmon muscle, by performing a genome-wide association study.

Fish population and recordings
A total of 642 fish from 92 sires and 194 dams were included in the analysis. All sires had four or more offspring from more than one dam. The fish originated from one year-class of the Atlantic salmon breeding population of SalmoBreed AS. The fish were transferred to sea at a mean weight of 0.1 kg, and slaughtered approximately 12 months later, at a mean weight of 3.6 kg (ranging from 1.2 to 6.4 kg). The fish were fed a commercial broodstock feed from Skretting (https://www. skretting.com/en/products/atlantic-salmon/?lifephase=474980) containing 70% fish oil and 3.1% each of EPA and DHA. All the fish were reared under the same conditions, and were fasted 13-14 days prior to slaughter.
At slaughter, sex was determined visually by inspection of the gonads, body weight (g) was recorded, and skeletal muscle samples for lipid and FA analysis were taken from Norwegian Quality Cut were collected, frozen and stored at −20°C.

Muscle fat and fatty acid analysis
Muscle fat content was measured by extracting total lipids from homogenized skeletal muscle samples of individual fish, according to the Folch method (Folch et al., 1957). Using one milliliter from the chloroform-methanol phase, FA composition of total lipids was analyzed following the method described by Mason and Waller (1964).The extract was dried briefly under nitrogen gas and residual lipid extract was trans-methylated overnight with 2′,2′-dimethoxypropane, methanolic-HCl, and benzene at room temperature. The methyl esters formed were separated in a gas chromatograph (Hewlett Packard 6890; HP, Wilmington, DE, USA) with a split injector, using an SGE BPX70 capillary column (length 60 m, internal diameter 0.25 mm, and film thickness 0.25 μm; SGE Analytical Science, Milton Keynes, UK) and a flame ionization detector. The results were analyzed using HP Chem Station software. The carrier gas was helium, and the injector and detector temperatures were both 270°C. The oven temperature was raised from 50 to 170°C at the rate of 4°C / min, and then raised to 200°C at a rate of 0.5°C / min and finally to 240°C at 10°C / min. Individual FA methyl esters were identified by reference to well-characterized standards. The content of each FA was expressed as a percentage of the total amount of FAs in the analyzed sample. FA composition for the most abundant FA in the muscle (16:0, 18:1n-9, 18:2n-6) and FA that are part of the omega-3 bioconversion pathway are given in Supplementary file 1.

Genotyping
The fish were genotyped using a customized 57 k axiom Affymetrix SNP Genotyping Array (NOFSAL02).
From the initial 57 k SNPs on the NOFSAL02 SNP chip, we retained those with call rate > 0.8, minor allele frequencies > 0.05, and Hardy-Weinberg equilibrium correlation p-value > .001. A total of 49,726 SNPs passed quality control filtering and were used in the GWAS.

Statistical analyses
The GWAS was carried out with the SNP and Variation Suite v8.x (SVS; Golden Helix, Inc., Bozeman, MT, www.goldenhelix.com).
A single-locus mixed linear model (EMMAX; Kang et al., 2010) was used with body weight and sex as covariates (Horn et al., 2018), and a genomic relationship matrix constructed in SVS: where y is a vector of FA trait records, μ is the overall mean, b 1 is a fixed regression coefficient of body weight as a covariable, b 2 is a (2 × 1) vector of fixed effects of sex (Horn et al., 2018), X 1 and X 2 are incidence matrices for the effects contained in b 1 and b 2 , a i is the regression coefficient of the genotype of SNP i; S i is a vector of genotypes of SNP i (coded as 0, 1, or 2 for the homozygote, heterozygote, and alternative homozygote genotype, respectively); u is a vector of additive genetic effects distributed as u~N(0,G σ 2 u ), where σ 2 u is the additive genetic variance, G is the genomic relationship matrix; Z is the corresponding incidence matrix to additive genetic effects, and e is the vector of random residual effects with e~N(0,I σ 2 e ). Genomic heritability (narrow sense) was estimated as the ratio of additive genetic variance to total phenotypic variance. Variance components were estimated by ASReml 4.0 (Gilmour et al., 2015) by applying the following linear mixed model: where y is a vector of the FA trait, μ is the overall mean, b 1 is a vector of fixed effect of body weight as a covariable, b 2 is a vector of fixed effect of sex (Horn et al., 2018), X 1 and X 2 are incidence matrices for the effects contained in b 1 and b 2 , u 1 is a vector of additive genetic effects distributed as u 1~N (0, Gσ u 2 ), where σ u 2 is the additive genetic variance, G is the genomic relationship matrix; Z is the corresponding incidence matrix to additive genetic effects, and e is the vector of random residual effects with e~N(0, Iσ e 2 ). The p-values were Bonferroni corrected to account for multiple testing, i.e. genome wide significance was reported when a nominal pvalue of α/TotSNPs = 0.05 / 49,726 = 1 × 10 −6 was obtained, where α is the desired genome wide significance level. The Bonferroni criterion is considered very strict (Johnson et al., 2010). Therefore, a threshold for suggestive association was also applied, which is less stringent and imply that across the genome one false positive result is expected (Lander and Kruglyak, 1995). When testing for suggestive association nominal p-values of 1 / 49,726 = 2 × 10 −5 were used.
Manhattan plots were created to visualize the association analyses. Quantile-quantile (Q-Q) plots were used to analyze the extent to which the observed distribution of the test statistic followed the expected (null) distribution in order to assess potential systematic bias due to population structure or to the analytical approach.
QTL regions were defined to include SNPs positioned within one million bp distance from a significant SNP, and whose likelihood was not more than a factor 10 lower than that of the significant SNP (1-LOD-dropp-off intervals) (Lander and Botstein, 1989). If regions for the same trait overlapped, the regions were combined.
The GWAS resulted in 10 and 11 SNPs surpassing the suggestive association threshold for the ratio traits DHA/ALA and DHA/DPA, respectively, but only one and two SNPs were significant for DHA and EPA, respectively (Table 2). Quantile-quantile plots presenting the distribution of observed vs expected p-values are presented in Fig. 1. These showed no indication of generally inflated p-values of the GWAS results.
The GWAS results for the ratio of DHA to DPA showed a clear signal with a genome-wide significant peak on chromosome 19 (Fig. 2). DHA/ DPA ratio can be an indicator of conversion of DPA to DHA. There are several processes involved in this conversion. Firstly, 22:5n-3 is elongated to 24:5n-3 by elolv2 or elovl4 (Morais et al., 2009), then Δ6fad forms 24:6n-3, and lastly peroxisomal β-oxidation forms DHA (Sprecher et al., 1995). The strongest association for DHA/DPA ratio mapped to AX-87328720 on chromosome 19 (p-value: 4.28 × e −8 ) ( Table 2). Five of the top 10 SNPs were located in the region~70-74,000 K ( Table 2). One of the genes in this interval is fatty acid elongase 2 (elovl2) (located at~73,300 K), which is directly involved in the conversion of DPA to DHA by elongating DPA to 24:5n-3, and is therefore a strong candidate to explain this association.
The third most significant SNP for DHA/DPA ratio, AX-87250961, was located on chromosome 19 at~57,323 K, and had a p-value of 2.06e-6 (Table 2). Located < 600 K base pairs (bp) from this SNP is long-chain acyl-CoA synthetase family member 5 (acsl5). The different long chain acyl-CoA synthases are important in channelling FAs to either β-oxidation or storage (Digel et al., 2009), and may in this way affect the DHA/DPA ratio. A recent study with primary human muscle cells concluded that acsl5 in human skeletal muscle functions to increase mitochondrial β-oxidation (Kwak et al., 2019).
For DHA, the GWAS resulted in a peak on chromosome 21 with one SNP surpassing the suggestive significance threshold (Fig. 3, Table 2).

Table 2
Markers surpassing the suggestive association threshold (P-value = 2 × 10 −5 ) for the traits DHA, EPA, ratio DHA/ALA and ratio DHA/DPA and the size of the QTL region for top SNPs indicated by 1-LOD-drop-off intervals. Genome-wide significant markers are indicated in bold. S.S. Horn, et al. Aquaculture 514 (2020) 734494 For EPA, two SNPs surpassed the suggestive significance threshold and they were located on the same chromosome as the DHA peak, although in a different region (Table 2). EPA had a very low genomic heritability (0.05; Table 1), and the pattern looked more noisy compared to DHA (Fig. 3).
The GWAS results for DHA to ALA ratio further strengthened the significance of chromosome 21 in relation to n-3 LC PUFA content in Atlantic salmon muscle. The analysis of this trait resulted in a clear signal on chromosome 21, with 10 SNPs surpassing the suggestive significance threshold, and four SNPs surpassing the Bonferroni genome-wide threshold (Fig. 3, Table 2). Three QTL regions were identified for this trait. The top SNP for DHA was among the significant SNPs also for DHA/ALA.
To identify candidate genes, we searched annotated genes with functional relevance to FA or lipid metabolism. No known genes of direct relevance to omega-3 FA metabolism were found in the region of interest mentioned above. However, a few potential candidate genes linked to EPA and DHA in a less direct way were annotated to this region. Three of these genes are related to T-cell signaling; T-cell surface glycoprotein CD3 zeta chain-like (LOC106581970), Toll-like receptor 7 (LOC106581875) and Toll-like receptor 8 (LOC106581917). Toll like receptors (TLR) are a family of transmembrane proteins that plays an essential role in the innate immune system through controlling inflammatory and immunological responses (Moresco et al., 2011;Rogero and Calder, 2018). TLRs controls T-cell responses both directly and indirectly (Rahman et al., 2009). N-3 LC PUFAs, such as EPA and DHA  are shown to exert anti-inflammatory actions through the TLR signaling pathways, and modulation of T-cell activation (Denys et al., 2005;Lee et al., 2003;Rogero and Calder, 2018). In Atlantic salmon, it has been shown that EPA and DHA affect TLR activation (Arnemo et al., 2017). Although the effect of EPA and DHA on expression of these genes have been documented, it is unknown whether the genes can influence EPA and DHA content in muscle.
A few other genes that could be related to EPA and DHA were also found in the region of the most significant SNPs on chromosome 21; such as renin receptor (renr), which was located close to the top significant SNP for ratio DHA/ALA (~150 K bp apart). Renin receptor has in mammalian studies been shown to be involved in the renin-angiotensin system, which regulates blood pressure (Fyhrquist and Saijonmaa, 2008). It is documented that EPA and DHA reduce blood pressure in mammals (Guo et al., 2019;Miller et al., 2014), and it has been suggested that this effect is via the renin-angiotensin system (Francois and Coffman, 2004;Jayasooriya et al., 2008). However, it is not known if EPA and DHA have a cardiovascular effect through renin-angiotensin system in Atlantic salmon. Another role of renin receptor is that it may reduce excessive autophagy in skeletal muscle during starvation (Mizuguchi et al., 2018). Because the fish in this study were starved for 14 days prior to slaughter, it can be hypothesized that this role of renin receptor may influence muscle FA composition in the fish of this trial.
Based on the results of this study, we cannot pinpoint which genes are determining EPA and DHA content in skeletal muscle of Atlantic salmon. The positional candidate genes discussed above are only indirectly linked to EPA and DHA, and have no documented direct effect on FA metabolism or content. Therefore, they are not considered strong candidate genes. There may be several reasons why we do not find any strong candidate genes for the n-3 LC PUFA traits on chromosome 21: • The SNP-chip does not contain markers in several of the genes known to influence lipid metabolism in mammals and/or fish, such as fatty acid synthase, peroxisome proliferator-activated receptors pparb2b and pparb2a, carnitine palmitoyl transferases cpt1b and cpt2, hormone-sensitive lipase, and acetyl-CoA carboxylase. Including SNPs in these genes could potentially identify QTLs not in linkage disequilibrium with the SNPs in this study.
• The region of the most significant SNPs on chromosome 21 contained several uncharacterized genes. Some of these might prove to be involved in omega-3 FA metabolism in the future. Additional studies including fine mapping and functional validation are needed to identify the causal variants for further elucidation of the variants underlying the FA composition traits • It is also possible that the GWAS could not directly identify the causal mutations in our study since all causal mutations are unlikely to be included in the 57 K SNP-chip, which is relatively small compared to the number of genes.
• The signals detected may be due to regulatory regions rather than genes. Genetic variation does not lead directly to phenotype, but instead alters intermediate molecular pathways that in turn affects the traits. Transcriptional regulatory sequences are as important for gene function as the coding sequences that determine the linear array of amino acids in a protein (Wray et al. 2003).
• The power of GWAS to detect genetic variation affecting Atlantic salmon FA traits is influenced by the number of samples, sampling methods and the procedure chosen for phenotypic recording. In this study, phenotypic recordings for FA traits were available for only 642 Atlantic salmon. However, they were all reared and slaughtered under the same conditions, fed the same feed and slaughtered at the same age, reducing the influences of environmental factors to a minimum. FA composition was determined by gas chromatography, which is highly accurate. The latter may explain why we found some genome-wide significant results, despite the relatively small number of records.
The current study found a candidate gene for DHA/DPA ratio directly involved in omega-3 bioconversion activity (elovl2). However, we did not find a signal in or near genes known to be involved in omega-3 bioconversion for the traits EPA, DHA or DHA/ALA ratio. The genes of this pathway should supposedly have been identified in these analyses if this pathway was important for determining EPA and DHA content in muscle. This study therefore suggests that the genes of the omega-3 bioconversion pathway are not clearly associated with EPA and DHA content in muscle of Atlantic salmon. The lack of association with omega-3 bioconversion pathway is in line with the results from a gene expression study of the same fish material that rather suggested selective beta-oxidation as an explanatory factor for variation in EPA and DHA content in muscle (Horn et al., 2019). A potential major factor behind the lack of association with omega-3 bioconversion pathway genes is that the fish was fed a high fish oil diet. High dietary levels of EPA and DHA suppress the activity of the omega-3 bioconversion pathway Ruyter and Thomassen, 1999;Tocher et al., 2003;Zheng et al., 2004). The starvation period may also have affected the metabolic activity and fatty acid composition. The differences observed in omega-3 content are therefore more likely caused by differences in other metabolic processes rather than omega-3 bioconversion. Another possible explanation is that although omega-3 bioconversion is important for FA composition in liver , it might not be as important for muscle FA composition. The role of this pathway in determining muscle FA composition has not been shown, and longchain omega-3 fatty acid content in liver and muscle appear to be two different traits, as they display a very low correlation (Horn et al., 2019).

Conclusions
Our results suggest that skeletal muscle n-3 LC PUFA content traits have a polygenic architecture with potentially a few QTL explaining moderate levels of the genetic variation in Atlantic salmon.
Markers showing a significant association with the DHA/DPA ratio were located on chromosome 19, close to the candidate gene elovl2, which is directly involved in the conversion of DPA to DHA.
The results suggest that genetic variation affecting fillet EPA and DHA content is present on Atlantic salmon chromosome 21.
A larger dataset is needed to validate these QTLs. Including SNPs in genes known to be involved in fatty acid metabolism could potentially identify additional QTLs not in linkage disequilibrium with the SNPs in this study.

Declaration of Competing Interest
The authors declare no conflict of interest.