Influence of Varying Dietary ω6 to ω3 Fatty Acid Ratios on the Hepatic Transcriptome, and Association with Phenotypic Traits (Growth, Somatic Indices, and Tissue Lipid Composition), in Atlantic Salmon (Salmo salar)

Simple Summary Plant oils are routinely used in fish feeds as a fish oil replacement. However, these terrestrial alternatives typically contain high levels of ω6 fatty acids (FA) and, thus, high ω6 to ω3 (ω6:ω3) FA ratios, which influence farmed fish and their consumers. The ω6:ω3 ratio is known to affect many biological processes (e.g., inflammation, FA metabolism) and human diseases; however, its impacts on fish physiology and the underlying molecular mechanisms are less well understood. In this study, we used 44 K microarrays to examine which genes and molecular pathways are altered by variation in dietary ω6:ω3 in Atlantic salmon. Our microarray study showed that several genes related to immune response, lipid metabolism, cell proliferation, and translation were differentially expressed between the two extreme ω6:ω3 dietary treatments. We also revealed that the PPARα activation-related transcript helz2 is a potential novel molecular biomarker of tissue variation in ω6:ω3. Further, correlation analyses illustrated the relationships between liver transcript expression and tissue (liver, muscle) lipid composition, and other phenotypic traits in salmon fed low levels of fish oil. This nutrigenomic study enhanced the current understanding of Atlantic salmon gene expression response to varying dietary ω6:ω3. Abstract The importance of dietary omega-6 to omega-3 (ω6:ω3) fatty acid (FA) ratios for human health has been extensively examined. However, its impact on fish physiology, and the underlying molecular mechanisms, are less well understood. This study investigated the influence of plant-based diets (12-week exposure) with varying ω6:ω3 (0.4–2.7) on the hepatic transcriptome of Atlantic salmon. Using 44 K microarray analysis, genes involved in immune and inflammatory response (lect2a, itgb5, helz2a, p43), lipid metabolism (helz2a), cell proliferation (htra1b), control of muscle and neuronal development (mef2d) and translation (eif2a, eif4b1, p43) were identified; these were differentially expressed between the two extreme ω6:ω3 dietary treatments (high ω6 vs. high ω3) at week 12. Eight out of 10 microarray-identified transcripts showed an agreement in the direction of expression fold-change between the microarray and qPCR studies. The PPARα activation-related transcript helz2a was confirmed by qPCR to be down-regulated by high ω6 diet compared with high ω3 diet. The transcript expression of two helz2 paralogues was positively correlated with ω3, and negatively with ω6 FA in both liver and muscle, thus indicating their potential as biomarkers of tissue ω6:ω3 variation. Mef2d expression in liver was suppressed in the high ω6 compared to the balanced diet (ω6:ω3 of 2.7 and 0.9, respectively) fed fish, and showed negative correlations with ω6:ω3 in both tissues. The hepatic expression of two lect2 paralogues was negatively correlated with viscerosomatic index, while htra1b correlated negatively with salmon weight gain and condition factor. Finally, p43 and eif2a were positively correlated with liver Σω3, while these transcripts and eif4b2 showed negative correlations with 18:2ω6 in the liver. This suggested that some aspects of protein synthesis were influenced by dietary ω6:ω3. In summary, this nutrigenomic study identified hepatic transcripts responsive to dietary variation in ω6:ω3, and relationships of transcript expression with tissue (liver, muscle) lipid composition and other phenotypic traits.


Introduction
Plant-based oils are commonly used in aquafeeds to replace fish oil (FO), due to decreasing global availability, rising market price, and concerns regarding the ecological sustainability of the finite fishery resources upon which FO production depends [1,2]. Indeed, plant oils (PO) were shown to be more economical and environmentally sustainable [3], and their inclusion as alternatives to FO in many experimental diets did not affect the growth and survival of farmed Atlantic salmon (Salmo salar) [4][5][6]. However, terrestrial oils are devoid of long-chain polyunsaturated fatty acids (LC-PUFA), such as eicosapentaenoic acid (EPA, 20:5ω3), docosahexaenoic acid (DHA, 22:6ω3), and arachidonic acid (ARA, 20:4ω6), which are abundant in FO. These LC-PUFA have important functions in vertebrate health, reproduction, neural development, and growth, among other biological processes [5,7]. This has resulted in decreased fillet EPA and DHA levels in farmed fish that were fed with PO as a partial or full replacement for FO, compromising their nutritional quality for human consumers [8][9][10]. Further, previous studies reported impacts on fish health and physiology with the dietary replacement of FO by PO (e.g., liver steatosis, altered complement pathway and phagocytic activity, and modulated expression of genes involved in immune response) [11][12][13][14][15]. Another concern is that most terrestrial oils used in aquafeeds, and the farmed seafood consuming them, may not provide adequate ratios of ω6 to ω3 (ω6:ω3) fatty acids (FA) due to high ω6 FA content [10,[16][17][18]. Previous human nutrition studies reported that high dietary ω6:ω3 promotes the pathogenesis of many diseases, including cardiovascular, inflammatory, autoimmune, and cognitive diseases, as well as obesity and cancer [19][20][21][22]. An optimal ratio of ω6:ω3 is important for maintaining the homeostasis of many biological processes such as cell apoptosis, inflammation, fatty acid and cholesterol metabolism, and others [23][24][25]. However, the underlying molecular mechanisms are still poorly understood in fish, and it is not known which genes are involved in variation in dietary and tissue ω6:ω3 in salmon fed high levels of terrestrial-based oils.
A feeding trial was performed to examine the impact of five plant-based diets with varying ω6:ω3 on salmon growth, tissue (i.e., muscle, liver) lipid composition, liver LC-PUFA synthesis, and transcript expression (targeted qPCR) of lipid metabolism and eicosanoid synthesis-related genes [26]. The objective of our current study was to utilize a 44 K salmonid oligonucleotide microarray [27][28][29] for the examination of the impact of the two extreme ω6:ω3 diets (i.e., high ω6 and high ω3) on the hepatic transcriptome at week 12. We hypothesized that salmon fed the two diets with the most extreme lipid compositions (i.e., High ω3 and High ω6) would show the most extensive transcriptomic differences. The current study used the same fish as in Katan et al. [26]. The aim was to identify novel biomarker genes and molecular pathways that are altered by variation in ω6:ω3. To aid in the elucidation of the relationships between liver transcripts and phenotypic traits (i.e., growth parameters, somatic indices, tissue FA and lipid class composition), correlation analyses were also performed.

Fish and Experimental Diets
Five experimental diets with varying ratios of ω6:ω3 were formulated and manufactured by Cargill Innovation Center (Dirdal, Norway). The diets had ω6:ω3 of 1:3 (high ω3), 1:2 (medium ω3), 1:1 (balanced), 2:1 (medium ω6) and 3:1 (high ω6). Dietary formulations and their lipid composition were published previously [26]. However, as they pertain to the current study, the formulation and lipid profiles of the relevant diets (i.e., high ω3, balanced, and high ω6) are also included as supplementary material herein (Tables S1 and  S2). All diets contained the same sources and equal levels of marine and plant proteins, but had different mixes of plant-based oils (i.e., linseed (flax), soy, and palm). All diets were formulated to be isonitrogenous and isoenergetic (Table S1), and to meet the nutritional requirements of salmonids [30].
Atlantic salmon pre-smolts were transported from Northern Harvest Sea Farms (Stephenville, NL, Canada) in October 2015, and held in the Dr. Joe Brown Aquatic Research Building (Ocean Sciences Centre, Memorial University of Newfoundland, St. John's, NL, Canada) in 3800-L tanks. After their arrival, fish were graded in order to select the most uniform population, and this was followed by PIT (Passive Integrated Transponder; Easy AV, Avid Identification Systems, Norco, CA, USA)-tagging for individual identification. Then, post-smolts (203 ± 24 g mean initial weight ± SE) were randomly distributed into twenty 620-L tanks (40 fish tank −1 ), and subjected to a 2.5-week acclimation period. After the completion of the acclimation period, fish were switched from the commercial diet (Nutra Transfer NP, 3 mm, Skretting Canada, St. Andrews, NB, Canada), and fed with the experimental diets (4 tanks diet −1 ) for 12 weeks. The photoperiod was maintained at 24 h light. Fish were fed overnight using automatic feeders, and apparent feed intake was recorded throughout the trial. Mortalities were also recorded during the trial. For additional details regarding the rearing conditions and recordings, refer to Katan et al. [26].

Sample Collection
Growth performance parameters (e.g., fork-length, weight, organ indices) were measured at the beginning and the end of the 12-week feeding trial [26]. At the end of the trial, salmon were starved for 24 h, and then 5 fish per tank were euthanized with an overdose of MS-222 (400 mg L −1 ; Syndel Laboratories, Vancouver, BC, Canada) and dissected for tissue collection. For gene expression analyses, liver samples (50-100 mg) were collected in 1.5 mL nuclease-free tubes, flash-frozen in liquid nitrogen, and stored at −80 • C until RNA extractions were performed. Liver and muscle samples, for lipid analyses, were collected, processed, and stored as described in Katan et al. [26]. Only liver samples from fish that showed weight gains within one standard deviation below and above the mean value of each tank were utilized for this study, in order to reduce biological variability in the gene expression data among fish. Tank means rather than dietary treatment means, were chosen for sample selection, so that variability between tanks could be included in the statistical analysis.

RNA Extraction, DNase Treatment, Column Purification and cDNA Synthesis
The TissueLyser II system (at 25 Hz for 2.5 min) with 5 mm stainless steel beads, QIA-GEN, Mississauga, ON, Canada) was used to homogenize liver samples in TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Samples were subjected to RNA extraction according to manufacturer's instructions. Due to low 260/230 ratios (i.e., 1.0-1.6) following TRIzol extraction, all RNA samples were then re-extracted (phenol-chloroform) and precipitated following standard methods [31]. This was followed by DNaseI treatment and column purification using RNase-free DNase Set and RNeasy Mini Kit (QIAGEN). All procedures were conducted according to manufacturer instructions, and as described in Xue et al. [29]. RNA integrity was verified by 1% agarose gel electrophoresis, and RNA purity and quantity were assessed by NanoDrop UV spectrophotometry (NanoDrop, Thermo Scientific, Mississauga, ON, Canada). DNased and column-purified RNA samples had A260/280 and A260/230 ratios of 1.8-2.2. All cDNAs were synthesized by reverse transcription of 1 µg of DNaseI-treated, column-purified total RNA from each sample, with 1 µL of random primers (250 ng; Invitrogen), 1 µL of dNTPs (0.5 mM final concentration; Invitrogen), 4 µL of 5× first-strand buffer (1× final concentration; Invitrogen), 2 µL of DTT (10 mM final concentration; Invitrogen) and 1 µL of Moloney murine leukemia virus (M-MLV) reverse transcriptase (RT) (200 U; Invitrogen) at 37 • C for 50 min, following the manufacturer's instructions, and as described in Xue et al. [29]. The total reaction volume was 20 µL. Finally, all cDNAs were diluted 40× with nuclease-free water (Invitrogen) prior to the qPCR.

Microarray Hybridization and Data Acquisition
Eight fish (2 from each of the 4 dietary tanks) from each of the 2 extreme ω6:ω3 treatments (high ω6 or high ω3) were used in the microarray analysis (i.e., 16 fish total), using a common reference design. Four array slides were used in the current study, and each array contained 2 fish per treatment, which were randomly selected. The common reference was made by an equal quantity of each DNase I-treated, column-purified total RNA liver sample. The microarray experiment was performed as described in Xue et al. [29]. Briefly, anti-sense amplified RNA (aRNA) was in vitro transcribed from 1 µg of each columnpurified RNA or reference pooled RNA using Ambion's Amino Allyl MessageAmp II aRNA Amplification kit (Life Technologies, Burlington, ON, Canada), following the manufacturer's instructions. The quantity and quality of aRNA were assessed using NanoDrop spectrophotometry and 1% agarose gel electrophoresis, respectively. Then, 20 µg of each aRNA were precipitated overnight, following standard molecular biology procedures, and re-suspended in coupling buffer. Each individual aRNA sample was labeled with Cy5 (i.e., experimental samples), whereas the reference pool was labeled with Cy3 (i.e., common reference) fluor (GE HealthCare, Mississauga, ON, Canada), following the manufacturer's instructions. The "microarray" function of the NanoDrop spectrophotometer was used to measure the labeling efficiency of the aRNA. The labeled aRNA (825 ng) from each experimental sample (i.e., Cy5) was mixed with an equal quantity of labeled aRNA from the common reference (i.e., Cy3), for each array, and the resulting pool was fragmented, following the manufacturer's instructions (Agilent, Mississauga, ON, Canada). Each pool was co-hybridized to a consortium for Genomic Research on All Salmonids Project (cGRASP)designed 4 × 44 K salmonid oligonucleotide microarray (GEO accession # GPL11299) [27] (Agilent). Finally, the arrays were hybridized at 65 • C for 17 h with rotation (10 rpm), using an Agilent hybridization oven. The microarray slides were washed immediately after hybridization as per the manufacturer's instructions.
Each microarray slide was scanned at 5 µm resolution with 90% laser power using a ScanArray Gx Plus scanner and ScanExpress v4.0 software (Perkin Elmer, Waltham, MA, USA), and the Cy3 and Cy5 channel photomultiplier tube (PMT) settings were adjusted to balance the fluorescence signal between channels. The resulting raw data were saved as TIFF images, and the signal intensity data were extracted using Imagene 9.0 (BioDiscovery, El Segundo, CA, USA). Removal of low-quality or flagged spots on the microarray, as well as log 2 -transformation and Loess-normalization of the data, were performed using R and the Bioconductor package mArray [32]. Features absent in more than 25% (i.e., 4 out of 16 arrays) of the arrays were omitted, and the missing values were imputed using the EM_array method and the LSimpute package [33,34]. The final dataset used for statistical analyses consisted of 10,264 probes for all arrays (GEO accession number: GSE139418; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE139418.

Microarray Data Analysis
The Significance Analysis of Microarrays (SAM) algorithm [35] was performed to identify genes that were significantly differentially expressed between the two extreme ω6:ω3 treatments. A false discovery rate (FDR) threshold of 10% was used with the Bioconductor package siggenes [36] in R. For the identification of additional transcripts that were differentially expressed between the two dietary treatments, the Rank Products (RP) method was also used, as this method is less sensitive to high biological variability [37,38]. The latter analysis was performed at a percentage of false-positives (PFP) threshold of 10%, using the Bioconductor package RankProd [39]. Due to high background signal in the first slide (i.e., slide # 11,502), no genes were initially identified as significantly differentially expressed; therefore, this slide was removed from the analyses. In order to maximize our capacity to identify differentially expressed genes, gene lists were obtained with 2 and 3 of the remaining slides (consisting of 4 and 6 fish per treatment, respectively). Each slide is composed of 4 arrays, i.e., 4 biological replicates analyzed per slide.
The resulting gene lists were annotated using the contiguous sequences (contigs) that were used for the design of the 60 mer oligonucleotide probes of the array [27]. Annotations were performed manually with a BLASTx alignment against the NCBI non-redundant (nr) amino acid database using an E-value threshold of 10 −5 . The best BLASTx hits corresponding to putative Homo sapiens orthologues were used to obtain gene ontology (GO) terms manually from the UniProt Knowledgebase (http://www.uniprot.org/, accessed on 4 November 2020).

qPCR Study and Data Analysis
Transcript expression levels of 10 genes of interest (GOI) ( Table 1), identified as differentially expressed in the microarray analyses, were assayed by qPCR. In addition to the high ω6 and high ω3 treatments, the qPCR analysis also included liver samples from fish fed the balanced diet. In addition to the microarray-identified GOI, BLASTn searches using publicly available Atlantic salmon cDNA sequences (i.e., in NCBI's nonredundant nucleotide (nt) and expressed sequence tags (EST) databases) were used to identify paralogues for each GOI, as described in Caballero-Solares et al. [40].  [40]. d The Salmo salar sequence of mtco2 used in the qPCR assay showed 87% identity with the 60 mer microarray probe (C060R108) affiliated with rainbow trout (Oncorhynchus mykiss). e Primers that were previously published in Caballero-Solares et al. [15] (annotated as VHSV-induced protein in that study). Alias pric285 stands for peroxisomal proliferator-activated receptor A interacting complex 285. f The Atlantic salmon sequences of itgb5 used in the qPCR assay showed 86% identity with the 60mer microarray probe (C002R106) affiliated with rainbow trout. Primers were designed based on common regions between DW540995 and XM014165323. g Primers that were previously published in Katan et al. [26].
Paralogue-specific primers were used for eif4b, htra1, lect2 and helz2 (Table S3 and Figures S1-S4). The sequences of the primer pairs used in qPCR, GenBank accession number of sequences used for primer design, and other details are presented in Table 1. Notably, primers for the transcript lhpl4 (GenBank accession number NM_001146670) failed quality testing, and thus, this transcript was not included in the qPCR study. In addition, the 60 mer microarray probe for mtco1 (C188R069) is affiliated with a rainbow trout sequence, and had relatively low identity (i.e., <85%) with available Salmo salar sequences (using NCBI's EST and nt databases) and, therefore, was excluded from the qPCR study. The program Primer 3 (http://frodo.wi.mit.edu, accessed on 19 October 2019) was used for primer design. Each primer pair was quality-tested, including standard curve and dissociation curve to ensure that a single product was amplified with no primer dimers [32,41]. Primer pairs were quality-tested using the 7500 Fast Real Time PCR system (Applied Biosystems/Life Technologies, Foster City, CA, USA). The amplification efficiency [42] of each primer pair was determined using a 5-point 1:3 dilution series starting with cDNA representing 10 ng of input total RNA. Two pools were generated (i.e., high ω3 pool and high ω6 pool), with each pool consisting of 8 fish (and each fish contributing an equal quantity to the pool). The reported primer pair amplification efficiencies are an average of the two pools, except if one pool showed poor efficiency or spacing (i.e., p43, eif2a, htra1a, helz2a and helz2b, where one pool was used due to low expression levels). A 5-point 1:2 dilution series was used for the primers mtco2 and helz2b as these transcripts had lower expression levels (fluorescence threshold cycle (C T ) values of~30 and 31, respectively). Furthermore, amplicons were checked by 1.5% agarose gel electrophoresis and compared with the 1 kb plus DNA Ladder (Invitrogen) to ensure that the correct size fragment was amplified.
All PCR amplifications were performed in a total reaction volume of 13 µL and consisted of 4 µL of cDNA (5 ng input total RNA), 50 nM each of forward and reverse primer and 1× Power SYBR Green PCR Master Mix (Applied Biosystems), and nuclease-free water (Invitrogen). The qPCR reactions, including no-template controls, were performed in technical triplicates using the ViiA 7 Real-Time PCR System (384-well format) (Applied Biosystems, Foster City, CA, USA) and the Power SYBR Green I dye chemistry. The Real-Time analysis program consisted of 1 cycle of 50 • C for 2 min, 1 cycle of 95 • C for 10 min, followed by 40 cycles (of 95 • C for 15 s and 60 • C for 1 min), with the fluorescence signal data collection after each 60 • C step. When a C T value within a triplicate was greater than 0.5 cycles from the other two values, it was considered to be an outlier, discarded and the average C T of the remaining two values was calculated. The relative quantity (RQ) of each transcript was calculated using ViiA 7 Software v1.2 (Applied Biosystems) for Comparative CT (∆∆C T ) analysis [45], with primer amplification efficiencies incorporated ( Table 1). The expression levels of each GOI were normalized to both normalizer genes, and the sample with the lowest normalized expression was used as the calibrator sample (i.e., RQ = 1.0) for each GOI, as in [46]. Transcript expression data are presented as RQ values (mean ± SE) relative to the calibrator.

qPCR Data
A general linear model with tank nested in diet, followed by a Tukey pairwise comparison (p < 0.05), was used to identify significant differences among dietary treatments at week 12. In cases where significant tank effect was observed (p < 0.05), a one-way ANOVA followed by a Tukey pairwise comparison post-hoc test was performed (Minitab 17 Statistical Software, State College, PA, USA). The RQ data are presented as mean ± SE. Each dietary treatment group was tested for outliers using Grubb's test (p < 0.05). In total, 9 RQ values were identified as statistical outliers in the entire dataset (i.e., out of 322 values), and excluded from the study. Each GOI had a minimum of 6 samples per dietary treatment, while most GOI had a sample size of 7-8 per dietary treatment. The qPCR fold-changes were calculated by dividing the mean RQ value of the high ω6 fish by that of the high ω3 fish. Finally, residuals were tested to verify normality, independence, and homogeneity of variance. Normality was examined using the Anderson-Darling test. If the test failed (p < 0.05), a one-way ANOVA on ranks was performed, which was followed by the Kruskal-Wallis test (SigmaPlot, Systat Software, Inc., Version 13, San Jose, CA, USA). In all cases, differences were considered statistically significant when p < 0.05.

Correlation Analyses of qPCR and Lipid Composition Data
Tissue lipid composition (muscle and liver) and growth performance of salmon fed varying ω6:ω3 diets were previously published by Katan et al. [26]. Pearson correlation analyses were performed in the current study to identify the relationships between hepatic transcript expression (i.e., qPCR data), tissue composition (i.e., % FA and lipid classes), and growth parameters (i.e., weight gain (WG), condition factor (CF)), using individual fish. All GOI in the qPCR study were used in the correlation analysis in order to identify differences between the liver and muscle tissue. Only ω3 and ω6 FA that accounted for >0.5% of the total FA in the tissue (average of each treatment) were included in the analyses. Furthermore, hierarchical clustering was used to group transcripts and lipid composition (using group average in PRIMER (Version 6.1.15, Ivybridge, UK)). IBM SPSS Statistics was used for the correlation analyses.
SAM analysis identified p43, eif2a, eif4b1, and itgb5 as differentially expressed genes (FDR < 10%) between the high ω6 and high ω3 fed salmon, using three slides (slides # 11,503-11,505) ( Table 2). These genes were down-regulated in the high ω6 compared with the high ω3 fed fish (fold-change values ranged from −2.79 to −5.12). One feature (itgb5) was represented in both SAM and RP analysis, and was down-regulated in the high ω6 compared to the high ω3 fed fish, in both analyses (−5.12 and −5.25foldchange, respectively).  a Refers to the identity of the probe on the 44 K array. Probes that are shown in bold font are features that were identified by SAM (FDR < 10%), and the remaining features were identified by RP analysis (PFP < 10%). The probe with an asterisk represents a feature that was identified in both SAM and RP analysis. Two 4 × 44 K array slides (slides # 11,504-11,505; representing 4 fish per treatment) were used in the RP analysis. However, the RP-identified mef2d was obtained using 3 slides (slides # 11,503-11,505; representing 6 fish per treatment). SAM-identified features were obtained using three slides (slides # 11,503-11,505). b Genes were identified by BLASTx, using the contig from which the microarray probe was designed against the NCBI non-redundant database. The best BLASTx hit with E-value < 10 −5 and an informative protein name was used, and presented with species name, GenBank accession number, E-value and % amino acid (AA) identity. c All microarray-identified genes, with the exception of lhpl4 and mtco1, were quantified by qPCR (see Materials and Methods). d Gene Ontology (GO) terms were selected from putative Homo sapiens orthologues (i.e., best BLASTx hit). Representative GO terms were identified (i.e., redundancies were not included), and divided into the categories: biological process (BP), molecular function ( Putative identities were determined for the 12 microarray-identified features, and functional annotations (i.e., GO terms) were collected for them ( Table 2). The microarrayidentified gene lhpl4 (4.78-fold up-regulated) is involved in the nervous system, with GO annotations "regulation of inhibitory synapse assembly" and "GABA receptor binding". The feature htra1b (3.57-fold up-regulated) was classified as a gene involved in cell proliferation and showed the functional annotations "positive regulation of epithelial cell proliferation", "proteolysis" and "extracellular space" ( Table 2). Several informative microarray features represented genes involved in translation, such as p43, eif2a, eif4b1, and rpl18 (−2.79 to −4.37-fold down-regulated), with the associated GO terms "tRNA binding", "aminoacyl-tRNA synthetase multienzyme complex", "translational initiation", "ribosome assembly", and "structural constituent of ribosome". Furthermore, the GO terms "defense response to virus", "inflammatory response" and "response to wounding" were also identified with p43. The features mtco2 and mtco1 (−3.27-and −7.11-fold down-regulated, respectively) were classified as mitochondrion respiratory chain components, and showed the GO terms "electron transport chain", "oxidation-reduction process" and "cytochrome-c oxidase activity". Other microarray-identified features corresponded to immune-and inflammation-related genes such as lect2a (with the GO terms "chemotaxis" and "metal ion binding") and itgb5 ("antigen processing and presentation" and "phagocytic vesicle"), and showed down-regulation in the high ω6 fed fish (−3.48 to −5.25-fold-change, respectively). The gene mef2d (−4.54-fold down-regulated) is involved in muscle cell proliferation, and in neuronal cell differentiation and survival, with the associated GO terms "muscle organ development", "skeletal muscle cell differentiation" "nervous system development", "apoptotic process" and "DNA-binding transcription factor activity" (Table 2). Further, the microarray-identified feature itgb5 was associated with the GO term "muscle contraction". Finally, the gene helz2a (−4.71-fold down-regulated) was classified as a gene involved in lipid metabolism regulation by peroxisome proliferator-activated receptor alpha (PPARα), and showed the functional annotations "regulation of lipid metabolic process", "nuclear receptor transcription activity", "ATP binding", "metal ion binding", "hydrolase activity" and "ribonuclease activity".

qPCR Study
Ten microarray-identified genes were used in the qPCR study. All genes, with the exception of mtco2 and rpl18, showed an agreement in the direction of expression foldchange (i.e., up-or down-regulation) between the microarray and qPCR studies ( Table 3). The microarray-identified helz2a showed significantly lower transcript expression in the high ω6 compared to the high ω3 fed fish (−1.49-fold; p = 0.04). The paralogue helz2b showed significantly lower expression in both the high ω6 and balanced groups compared to the high ω3 fed fish (−1.61-fold; p = 0.002). The transcript mef2d showed significantly lower expression in the high ω6 compared to the balanced fed fish (−1.27-fold; p = 0.03), and a lower expression trend in the high ω6 compared to the high ω3 fish (−1.22-fold; p = 0.06). Both paralogues of htra1 were numerically higher (although not statistically significant) in the high ω6 compared to the balanced and high ω3 fish (1.34-1.57-fold and 3.75-2.09-fold; p = 0.25 and 0.07, respectively) ( Table 3).

Overlapping Lipid-Gene Correlations between the Liver and Muscle Analyses
Some significant correlations showed an overlap between the liver and muscle analyses ( Figure 3). The hepatic transcript expression of htra1a was correlated positively with 20:2ω6 in both tissues. The expression of eif2a was correlated positively with Σω3, and negatively with 18:2ω6, Σω6 and ω6:ω3, while that of mef2d showed negative correlations with ω6:ω3 in both liver and muscle. The transcript expression of helz2a was correlated positively with 18:3ω3, and negatively with 18:2ω6, 20:3ω6 and Σω6, while that of helz2b correlated positively with 18:3ω3, 20:3ω3 and 20:4ω3, and negatively with 18:2ω6, 20:2ω6, 20:3ω6, Σω6 and ω6:ω3 (Figure 3). The hepatic transcript expression of most genes (all except helz2b with 20:4ω3) showed stronger positive correlations with liver compared to muscle FA. However, negative correlations were mostly (all except eif2a with 18:2ω6, and both helz2 paralogues with 20:3ω6) more significant with muscle compared with the liver FA ( Figure 3).

Discussion
The microarray study indicated that dietary variation in ω6:ω3 resulted only in small changes in the liver transcriptome of salmon fed plant-based diets. This can partly be explained by the fact that growth performance and somatic indices were not significantly affected by diet [26]. It was previously shown that different replacements of FO with camelina oil had no impact on Atlantic cod (Gadus morhua) growth, and resulted in only one microarray-identified gene that showed a significant difference in spleen basal expression between treatments [47]. Furthermore, our results are in line with previous microarray studies, which demonstrated that dietary replacement of fish meal (FM) and FO with terrestrial ingredients resulted in subtle gene expression changes in Atlantic salmon distal intestine [48], head kidney [49], and liver [50]. However, Atlantic salmon fed soy and linseed oils showed large alterations in hepatic gene expression compared to those fed FO [51]. Differences in the numbers of responsive transcripts between Leaver et al. [51] and the current study could be related to dietary lipid sources and studied time points.
Several transcripts that play important roles in immune and inflammatory response (lect2a, itgb5, helz2a, p43), lipid metabolism (helz2a), cell proliferation (htra1b), muscle and neuronal cell development (mef2d), and translation (eif2a, eif4b1, p43) were identified by our microarray study as diet-responsive. All transcripts, with the exception of mtco2 and rpl18, showed an agreement in the direction of expression fold-change between the microarray and the qPCR analyses ( Table 3). The 60mer microarray probe representing the transcript mtco2, which was designed using a rainbow trout cDNA sequence, showed only 87% similarity (see Materials and Methods) with the Salmo salar cDNA sequence that was used in the qPCR study (and other S. salar sequences in NCBI databases). This fact, as well as other limitations (e.g., mRNA regions targeted by the qPCR primers and microarray probe may not be the same; possibility of contig misassembly) could have contributed to the disagreement between microarray and qPCR results [32].
Hepatic helz2a showed a significant differential expression between the high ω6 and high ω3 fed fish in the microarray experiment, and both paralogues of this transcript (i.e., helz2a and helz2b) were significantly down-regulated in the high ω6 compared to the high ω3 fed fish in the qPCR analysis. Interestingly, the transcript expression of helz2b was positively correlated with ω3 (i.e., 18:3ω3, 20:3ω3, 20:4ω3), and negatively with ω6 PUFA (i.e., 18:2ω6, 20:2ω6, 20:3ω6, 20:4ω6), Σω6 and ω6:ω3, in the liver tissue ( Figure 1). In the muscle tissue, these PUFA (with the exception of 20:4ω6) were also correlated with hepatic helz2b expression ( Figure 2). These data suggest that helz2 is a potential novel molecular biomarker of tissue variation in ω6:ω3. The protein encoded by this gene is a nuclear transcriptional co-activator for PPARα [52][53][54], which is a master regulator of numerous genes involved in lipid metabolism processes (e.g., FA oxidation, and metabolism of bile acids, triacylglycerols, and retinoids) [55]. Additionally, HELZ2 was shown to have an antiviral function in mammals [56,57], and its gene was identified as an ancestral (between mammals and fish) interferon stimulated gene (ISG) with conserved components of antiviral immunity [58]. Helz2 transcripts (referred to as VHSV-induced protein (vig1)) showed up-regulation in the head kidney of Atlantic salmon exposed to the viral mimic polyriboinosinic polyribocytidylic acid (pIC) [15]. The negative correlation between liver helz2b and HSI is not surprising, given the involvement of PPARα in hepatic FA β-oxidation, and in liver steatosis [59,60]. In addition, the observed positive correlations with ω3 PUFA are in line with the anti-inflammatory properties of PPARα [60]. In our previous study [26], it was observed that the fatty acid binding protein-encoding transcript fabp10 showed an upregulation trend (p = 0.06) in the high ω3 compared to the balanced and high ω6 fed fish. Thus, this suggests that the high ω3 diet may have influenced the transport of ω3 FA in liver cells, and played a role in the activation of PPARα. The interaction between liver fatty acid transport and PPARα activation has been shown in previous mammalian studies [61,62]. Another potential mechanism that could explain the positive correlation between helz2b expression with ω3 PUFA, is that ω3 PUFA bind to PPARα with higher affinity than ω6 PUFA [63]. However, there is still a lack of knowledge about the interaction between dietary ω3 and ω6 PUFA, and the mechanisms by which they regulate PPARα activators [64].
Mef2d was identified in the microarray as down-regulated by the high ω6 diet, with qPCR showing significantly lower expression in the high ω6 compared to the balanced diet fed fish. Furthermore, hepatic mef2d expression was correlated positively with liver 20:5ω3 and Σω3, and negatively with the ratio of ω6:ω3 in both liver and muscle tissues. The gene mef2, characterized in zebrafish (Danio rerio) [65] and common carp (Cyprinus carpio) [66], is involved in skeletal and cardiac muscle development and differentiation, as well as in neuronal cell development [67][68][69][70]. Wei et al. [71] reported a significant increase in skeletal muscle mef2c transcript expression in pigs fed with linseed-enriched (10%) as compared with a control diet (0%). Additionally, a study with Atlantic salmon revealed that feeding with a synthetic FA (i.e., tetradecylthioacetic acid (0.25%) compared to a control diet (0%) increased the cardiosomatic index, and the cardiac expression of mef2c [72]. In relation to the liver, previous studies showed that members of the MEF2 family regulate the activation of hepatic stellate cells -a type of cell involved in liver fibrosis-in mice [73] and rats [74], and ω3 PUFA inhibited the proliferation and activation of these cells in mouse liver [75]. Further, Wang et al. [74] reported that mef2d was induced during hepatic stellate cells activation.These data may support the idea that the transcript mef2 is responding to dietary FA (particularly ω3 PUFA) in vertebrates. However, as most studies examined the role of mef2 expression in skeletal [71,76] and cardiac [72,77] muscle development, the interactions between mef2d and liver physiology are less understood in fish. Future studies should investigate the influence of dietary ω6:ω3 on liver mef2d expression, and their interaction with hepatic stellate cells in fish.
Serine protease HTRA1-encoding transcript (htra1b) was up-regulated in the high ω6 compared to the high ω3 fed fish, in the microarray study. Htra1b showed a similar trend of higher expression in the high ω6 fed fish (p = 0.07) in the qPCR analysis, and a similar fold-change (i.e., high ω6/high ω3) in the microarray and qPCR studies (Table 3). Further, hepatic htra1b was positively correlated with ΣSFA, 22:6ω3 and ST, and negatively with ΣMUFA, 20:3ω3 and TAG in the liver, while, in the muscle, it showed positive and negative correlations with ω6:ω3 and 18:4ω3, respectively. The transcript htra1a was positively correlated with 20:2ω6 and ω6:ω3, and negatively with ω3 PUFA (i.e., 20:5ω3 and Σω3) in the liver, and showed positive correlation with 20:2ω6 in the muscle. Serine protease HTRA1 function was linked to cell growth and apoptosis, as well as immune and inflammatory responses (by inhibiting the TGF-beta pathway) in mammalian tissues (e.g., eye, bone and liver) [78][79][80][81]. It was previously shown that dietary FM replacement with terrestrial plant alternatives induced higher hepatic htra1 transcript expression in Atlantic salmon [82]. Conversely, replacing both dietary FM and FO by terrestrial plant alternatives down-regulated the transcription of htra1a and htra1b in Atlantic salmon liver [40]. Discrepancies between Caballero-Solares et al. [40] and the present study extend to the FA-transcript correlation analyses; while, in the present study, the transcript expression of htra1 paralogues correlated positively and negatively with tissue ω6 and ω3 FA levels, respectively, the opposite tendency was observed in Caballero-Solares et al. [40]. However, unlike the previous studies [40,82], our study tested different mixes of vegetable oils while keeping FM/FO inclusion levels equal across diets. Therefore, although the studies cannot be directly compared, such discrepancies suggest that the regulation of HTRA1-mediated processes in the liver of Atlantic salmon depends on the combination of protein and lipid sources included in the diet. Finally the negative correlation observed between htra1b and growth parameters (i.e., WG and CF) in the current study ( Figure 2) is interesting, as mammalian HTRA1 was negatively linked to skeletal muscle development and bone formation [80,83,84].
The immune related microarray features lect2a and itgb5 were down-regulated in the high ω6 compared to the high ω3 fed fish, and an agreement was observed in the direction of expression fold-change between the microarray and qPCR studies. LECT2 is a multifunctional protein that plays a role in cell growth, neutrophil chemotactic activity, and innate immune response against pathogens in fish [85][86][87][88]. LECT2 also functions as a hepatokine that modulates the inflammatory response in mammals [89,90]. Earlier microarray studies reported down-regulation of hepatic lect2 in Atlantic salmon fed terrestrial as compared with marine diets [29,40]. This may indicate that pro-inflammatory plant-based diets suppress the constitutive transcript expression of hepatic lect2. The observed negative correlation between both paralogues of lect2 and VSI suggests that lect2 suppression is related to higher lipid deposition. However, the interaction between fat deposition and immune response is very complex, and requires further investigation in fish [90,91]. Additionally, we showed a positive correlation between the transcript expression of both paralogues of lect2 and liver sterol content. Interestingly, a previous study reported down-regulation of hepatic lect2 in Atlantic salmon fed a cholesterol-supplemented diet as compared with a non-supplemented plant-based diet, and this coincided with reduced plasma phytosterols (i.e., sitosterol and campesterol) [92]. Although dietary sterol levels were not significantly different in our study, liver sterol concentration did vary among treatments [26]. Clearly, more studies are required in order to elucidate the impact of dietary and tissue cholesterol and phytosterols on the constitutive transcript expression of lect2 in fish. Finally, the correlations observed in our study between itgb5 and liver FA are in line with the notion that FA can regulate the mRNA and protein levels of integrins and other adhesion proteins in leukocytes and endothelial cells [93,94]. Interestingly, European seabass (Dicentrarchus labrax) fed a plant-based diet showed a down-regulation in hepatic Integrin beta-2 compared to those fed a marine diet [95], while a reduction in the ω6:ω3 ratio of human lung cancer cells resulted in a delayed adhesion, and down-regulation of integrin-α2 [96]. Taken together, these data suggest that the transcript expression of integrins may be impacted by dietary or tissue ω6:ω3.
Similar to lect2a and itgb5, the transcripts p43, eif2a, and eif4b1 were down-regulated in the high ω6 fed fish in the microarray experiment, and they showed an agreement in the direction of dietary modulation with the qPCR study. The transcript expression of p43, eif2a, and eif4b2 was negatively correlated with 18:2ω6, and both p43 and eif2a were positively correlated with liver Σω3. Further, eif2a expression was positively correlated with 18:3ω3 and Σω3, and negatively correlated with 18:2ω6, 20:2ω6, Σω6 and ω6:ω3 in the muscle. The protein p43 is associated with a multi-tRNA synthetase complex, and regulates tRNA channeling in mammals [97]. In addition, p43 also encodes an apoptosis-induced cytokine, which regulates inflammation, wound healing, and angiogenesis [98][99][100]. Phosphorylation of the protein eIF4B was shown to stimulate translation in zebrafish [101,102] and yeast [103]. However, phosphorylation of eIF2A repressed translation in response to accumulation of misfolded proteins in the ER of several fish species [104][105][106]. Thus, changes in the expression patterns of p43, eif2a and eif4b, and the correlations observed with tissue lipid composition, suggest that some aspects of protein synthesis were influenced by dietary and tissue ω6:ω3. Previous mammalian studies demonstrated that translation is inhibited in apoptotic cells, and this was correlated with enhanced cleavage of the eukaryotic translation initiation factors eIF4B, eIF2, and others [107,108]. Thus, the fact that the high ω6 fed fish showed up-regulation of htra1b ( Table 2) may suggest that apoptosis was associated with the observed modulation of translation-related transcripts (Table 2), and their correlations with tissue FA (Figures 1 and 2). However, as our microarray study did not identify other well-known apoptosis biomarkers (e.g., genes encoding caspases and Bcl-2 family members), this can only be postulated. Further, the stimulatory effects of ω3 PUFA on protein synthesis [109] could be another potential mechanism explaining the positive correlations observed between p43, eif2a, and liver ω3 FA. Research examining the impact of replacing FO/FM with plant-based diets on protein synthesis in salmonids has been contradictory. Some authors showed an induction [82], while others showed a suppression [110] of these and other translation-related transcripts. Indeed, protein synthesis regulation in fish is a dynamic process, and is influenced by dietary formulations, genetic [50] and abiotic factors, protein requirement, growth, and the tissues examined [111,112].

Conclusions
Our 44 K microarray study demonstrated that high ω6 and high ω3 plant-based diets with varying ratios of ω6:ω3 (i.e., 2.7 and 0.4, respectively) resulted in a relatively low number of differentially expressed transcripts in salmon liver. However, the identified transcripts and/or their functional annotations suggested important roles in lipid metabolism (helz2a), cell proliferation (htra1b), immune and inflammatory response (lect2a, itgb5, helz2a, p43), control of muscle and neuronal cell development (mef2d), and translation (eif2a, eif4b1, p43). Two paralogues of helz2 were down-regulated in the high ω6 compared to the high ω3 fed fish in the qPCR study. Significant positive correlations were observed between the hepatic transcript expression of helz2b and ω3 PUFA, while negative correlations were identified with ω6 PUFA and ω6:ω3, in both the liver and muscle tissues. This indicated that the PPARα activation-related transcript helz2 is a potential novel molecular biomarker of tissue variation in ω6:ω3. Given these data and the importance of helz2 as an ancestral vertebrate interferon stimulated gene, future studies should investigate the dietary ω6:ω3 impact on Atlantic salmon anti-viral response. The transcript mef2d was suppressed in the high ω6 compared to the balanced fed fish, and was negatively correlated with ω6:ω3 in both tissues. Our microarray study revealed that the upregulation of hepatic htra1b concurred with the suppression of immune-and inflammatory-related transcripts (i.e., lect2a, p43, helz2a, helz2b, and itgb5). This supported the idea proposed by other researchers [40,82] of a link between the dietary modulation of htra1 and that of immune-related transcripts. Finally, the transcripts p43, eif2a, and eif4b1 were significantly down-regulated in the high ω6 compared to the high ω3 fed fish in the microarray, and showed an agreement in the direction of expression fold-change between the microarray and qPCR studies. These data, along with the significant correlations observed between p43, eif2a and eif4b2 expression and tissue PUFA, suggested that the molecular regulation of protein synthesis in the liver may have been impacted by dietary ω6:ω3.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biology10070578/s1, Table S1: Formulation and nutrient composition (%) of experimental diets fed to Atlantic salmon. This diet information was published in Katan et al. [26]. However, it is included here as this information is pertinent to the current study as well, Table S2: Lipid and FA composition (%) of experimental diets fed to Atlantic salmon. This diet information was published in Katan et al. [26]. However, it is included here as this information is pertinent to the current study as well, Table S3: Paralogues of genes involved in the qPCR, and their identity (%) over the aligned nucleotide regions, Figure S1: Alignment of the nucleotide sequences of htra1a (GenBank accession number NM_001141717) and htra1b (GenBank accession number EG831192). Conserved regions are highlighted in yellow. Alignments were performed using AlignX (Vector NTI Advance 11). Forward primers are bolded and underlined, whereas reverse primers are bolded without an underline, Figure S2: Alignment of the nucleotide sequences of eif4b1 (GenBank accession number BT072661) and eif4b2 (accession number DY739566). Conserved regions are highlighted in yellow. Alignments were performed using AlignX (Vector NTI Advance 11). Forward primers are bolded and underlined, whereas reverse primers are bolded without an underline, Figure S3: Alignment of the nucleotide sequences of lect2a (GenBank accession number BT059281) and lect2b (GenBank accession number DV106130). Conserved regions are highlighted in yellow. Alignments were performed using AlignX (Vector NTI Advance 11). Forward primers are bolded and underlined, whereas reverse primers are bolded without an underline, Figure S4: Alignment of the nucleotide sequences of helz2a (GenBank accession number BT072427) and helz2b (GenBank accession number EG928625). Conserved regions are highlighted in yellow. Alignments were performed using AlignX (Vector NTI Advance 11). Forward primers are bolded and underlined, whereas reverse primers are bolded without an underline.