Muscle fat and moisture contents are interrelated attributes that affect the organoleptic quality and nutritional value of muscle foods [19]. In fish, high-fat content may influence fillet processing and reduce the firmness leading to fillet downgrading [1]; moreover, it significantly impacts texture, juiciness, and flavor. In mammals, increased marbling scores are positively related to beef tenderness, accounting for ~ 9% of the shear force variation [34]. The inability to retain moisture during postmortem storage, in both fish and mammals, is associated with a high drip loss and, in turn, reduces the industry profitability by influencing processing yield and palatability [35, 36]. In the pork industry, drip loss results in up to 10% product losses affecting profitability at wholesale and retail levels [35]. Increased knowledge of the genetic basis of muscle quality traits will facilitate to further the commercial breeding in salmonids. GWA studies are powerful tools to identify genetic variants associated with complex traits [19]. However, no GWA studies were previously conducted to dissect the genetic architecture of fillet fat and moisture contents in rainbow trout. The estimated heritability for fat and moisture content was 0.31 ± 0.07 and 0.42 ± 0.07, respectively, suggesting existence of adequate genetic variability in the NCCCWA fish population to allow genetic improvement through selective breeding. A higher rate of genetic gain is obtained when genomic information is used [37].
In this study, we used genomic windows of 50 SNPs of a 50K SNP chip, to perform GWA analyses, in addition to the single-marker analysis approach, to identify genomic regions associated with the traits. Given that the 50 K, SNP chip contains SNPs of potential association with intramuscular fat content [28; 49], all fish used to build the SNP-chip were excluded from the GWA analysis in the present study.
The fish population used for the current GWA analyses had an average muscle fat content of 9.2 ± 1.91 (%) and moisture content of 69.93 ± 1.75 (%). Variations in fat and moisture content are shown in Figure (1). Previous studies reported a significant correlation between changes in fat and moisture content in fish [10, 15]. Consistently, our data showed a significant negative correlation between fat and moisture content (R = -0.88; p-value = 6.3E-262).
QTL affecting muscle fat and moisture content using WssGBLUP
All 35,322 SNPs (70.6%) that passed QC were used in the WssGBLUP analysis. A complete list of proportions of additive genetic variance for fat content explained by all genomic windows is provided in Table S1. Of them, a total of 137 genomic sliding windows explaining at least 2% (arbitrary value) of the additive genetic variance for fat content are listed in Table S2. Most of the SNP sliding windows (n = 124; ~91%) were located within 62 protein-coding genes. Genomic loci affecting the additive variance for fat content were clustered in 5 chromosomes (1, 4, 5, 19, and 29) (Fig. 2).
Chromosome 19 harbored the highest number (n = 50) and the most significant peaks affecting fat content (up to 5.51%) (Table S2, Fig. 2). Many of the SNPs were located within the CDS of the SNP-harboring genes (n = 58) as well as their 3’UTR (n = 55). In order to understand the biological significance of the QTL associated with fat content, we annotated the SNP-harboring genes and searched their functions in the literature (described below).
Similarly, a complete list of the proportions of additive genetic variance for moisture content explained by all windows identified in this study is provided in Table S3. A total of 178 genomic sliding windows revealing at least 2% of the additive genetic variance for moisture content are listed in Table S4. Most of the SNP sliding windows (n = 165; ~93%) were located within 86 genes coding for proteins. Genomic loci affecting the additive variance for moisture content were clustered on 5 chromosomes (5, 14, 19, 25, and 29) (Fig. 3). Chromosome 29 harbored the highest number (n = 48), whereas the most significant peaks affecting moisture content (up to 4.46%) were identified on chromosome 19 (Table S4, Fig. 3). Many of the SNPs were located within CDS of the SNP-harboring genes (n = 68) as well as their 3’UTR (n = 72).
Common genes affecting muscle fat and moisture content
As shown above, a negative linear relationship has been established between fat and moisture content in this selectively bred rainbow trout population (YC 2010 and YC 2012), suggesting a common mechanism underlying the genetic variation in the two traits. This negative correlation was consistent with other studies in fish and mammals [6, 7, 10, 15]. The current WssGBLUP identified common SNPs affecting the additive genetic variance for fat and moisture content on chromosomes 19 and 29 (Tables S2 & S4). The majority of the common SNPs (n = 47) were located on chromosome 19. Thirty-two SNPs, out of 47, involved in lipid metabolism were identified in 16 protein-coding genes on chromosome 19 (Table 1). Briefly, cathepsin B had a single 3’UTR SNP. Cathepsin B regulates very-low-density lipoprotein (VLDL) secretion and free fatty acid uptake in response to oleic acid exposure in mice [38]. Thioredoxin-related transmembrane protein 1-like (TMX1) had three SNPs. Loss of TMX increases lipid peroxidation in TMX−/− mice, which, in turn, enhances oxidative stress [39]. Guanine nucleotide-binding protein GI/GS/GO gamma-2 subunit (GNG2) had a single 3’UTR SNP. GNG2 expression is positively correlated with adipocyte size [40]. SNPs in genes encoding beta-taxilin and Alpha-L-fucosidase 2 (FUCA2) were covering windows explaining the highest proportion of the additive genetic variation for fat and moisture content. Adipose tissue of obesity susceptible and resistant rats differentially expressed beta-taxilin under a high-fat diet [41]. FUCA2 is a glycolipid processing enzyme [42]. Two SNPs in F-box only protein 30 (FBXO30) and the microtubule-binding protein ensconsin were ranked next to beta-taxilin and FUCA2. An SNP in FBXO30 was located in a genomic region explaining 4.95% of the additive genetic variance for polyunsaturated fatty acids in cattle [43]. Knock-down of microtubule-binding or -associated proteins led to changes in fat accumulation during adipogenesis [44]. Dihydropyrimidinase-related protein 5-like (CRMP5) had a single synonymous SNP. CRMP5 has GO terms belong to lipid metabolic processes [45]. Five SNPs were identified in a gene encoding trifunctional enzyme subunit alpha, mitochondrial (HADHA). This gene is involved in fatty acid beta-oxidation [45]. mRNA decay activator protein ZFP36L1 had three SNPs in the 3’UTR. Knockdown of mammalian ZFP36L1 led to downregulation of ERK activation and inhibition of adipogenesis [46]. A single 3’UTR SNP was identified in ELM2 and SANT domain-containing protein 1 (ELMSAN1). Epigenome-wide association analysis showed DNA methylation changes in ELMSAN1 were associated with body mass index (a key measure of adiposity) [47]. Prostaglandin reductase 2 (PTGR2) had two nonsynonymous SNPs. This enzyme catalyzes reduction of the conjugated α,β-unsaturated double bond of 15-keto-PGE2 in an NADPH-dependent manner, which is a critical step in inhibition of PPARγ-mediated adipocyte differentiation [48]. Spectrin beta chain, erythrocytic (SPTB) gene had two SNPs. The SPTB interacts with phospholipids in natural membranes and has a role in controlling the fluidity of the inner lipid leaflet of the cell membrane [49].
Chromosome 29 had 14 SNPs in genomic windows explaining at least 2% of the additive genetic variance for intramuscular fat and moisture content (Tables S2 & S4). Of them, seven SNPs were involved in lipid metabolism (Table 1). A single SNP was identified in a gene encoding short-chain specific acyl-CoA dehydrogenase, mitochondrial (ACADS). This enzyme has a role in fatty acid beta-oxidation [50]. An intronic SNP was identified in a gene coding for arrestin domain-containing protein. The latter has GO terms belonging to fat pad and skin development and regulates the body mass [51]. Myocyte enhancer factor 2c (MEF2C) had the highest number of SNPs (n = 5) on chromosome 29. MEF2C is a transcription factor involved in skeletal muscle differentiation; however, it has been reported as a constituent of a mechanism that programs gene expression involved in development of brown adipocytes [52]. MEF2A and MEF2D isoforms exhibited in vivo differential expression in mammalian striated muscle and white adipose tissue of insulin-deficient diabetic mice [53]. To our knowledge, the role of MEF2C in white adipose tissue remains uncertain.
In addition, twelve SNPs in genes involved in transmembrane transport and cytoskeleton remodeling were identified in common QTL affecting additive variance for fat and moisture content (Table 2). The majority of these SNPs were identified on chromosome 19 (n = 11). Three synonymous SNPs were identified in a gene encoding intersectin-2 (ITSN2). This protein is necessary for the clathrin-mediated endocytosis and actin cytoskeleton remodeling [54]. Six SNPs were identified in 3 genes involved in vesicle-mediated transport (i.e., exocytosis); dnaJ homolog subfamily C member 5B, visinin-like protein 1, and syntaxin-binding protein 5. The actin cytoskeleton remodeling controls each step of exocytosis [55]. Three SNPs were identified in microtubule-associated protein RP/EB family member 3 (MAPRE3) and centrin-3. MAPRE3 and centrin-3 control the dynamics of the microtubule cytoskeleton [45, 56].
Unique genes affecting the additive genetic variance for moisture
The actin cytoskeleton interacts with the cell membrane to control water transport [57]. In the current study, thirty-six variants in genes (n = 14) involved in cytoskeleton remodeling were identified, affecting the additive variance for moisture content in rainbow trout (Table 3). Briefly, bone morphogenetic protein receptor type-2 (BMPR2) had a single synonymous SNP. BMPR2 is known to interact with the cytoskeleton, and BMPR2 mutant mice exhibited cytoskeletal defects [58]. A single SNP was identified in a gene encoding muscle associated receptor tyrosine kinase (MUSK). Activation of MUSK in myotubes regulates the reorganization of the actin cytoskeleton [59]. Two SNPs were identified in THAP domain containing 1 (THAP1), which has a role in regulation of the mitotic cell cycle [60]. The gene encoding asparaginyl-tRNA synthetase (NARS) had 3 SNPs in windows explaining the highest additive variance (up to 3.46%; Table 3). Mutations in NARS leads to cell cycle arrest in the S phase [61]. The actin cytoskeleton undergoes dramatic changes during the cell cycle [62]. Ten SNPs were identified in three genes coding for cyclin-I (CCNI), cyclin-G1 (CCNG1), and cyclin-G2 (CCNG2). Cyclins function as regulators of the cell cycle and actin cytoskeleton dynamics [63]. The serine/threonine-protein, phosphatase 2A (PP2A), had a 3’UTR SNP. This phosphatase is associated with microtubule stabilization, where it binds and dephosphorylates the microtubule-associated proteins [64]. Annexin A6 (ANXA6) had two synonymous SNPs. ANXA6 contributes to membrane and cytoskeleton organization in a Ca2+-dependent manner [65]. Tubulin beta-4B chain (TUBB4B) had four synonymous SNPs within 1Kb of chromosome 25. TUBB4B is a critical component of microtubules [45]. Five SNPs, clustered in ~ 2Kb, were identified in a gene coding for mid1-interacting protein 1 (MID1IP1) [66]. This protein enhances fatty acid biosynthesis [66] and stabilizes microtubule organization [45]. Two SNPs were identified in a gene encoding tubulin-specific chaperone A (TBCE). TBCE is a tubulin-folding protein required for proper microtubule cytoskeleton organization [67]. Additionally, mutations in TBCE drive muscle atrophy [67]. Proteinase-activated receptor 1 (PAR1) and PAR2 had four SNPs. PAR-mediated, RhoA activation is vital for cytoskeletal reorganization [68].
A strong positive correlation between moisture and protein content has been established in different species. A simultaneous decline in protein and moisture content was previously reported in mammals [7]. Moisture content in rainbow trout exhibited a bidirectional relationship with protein content depending on physiological/metabolic status. For example, a negative correlation between moisture and protein content were previously reported under muscle catabolic conditions associated with full sexual maturation (R2 = 0.994, p < 0.01) [8]; whereas, a positive correlation was reported in female trout, on a high plane of nutrition, that were approaching spawning [10]. This was explained by selective mobilization of either protein during spawning or fat before spawning; in either case, the depleted macromolecule was replaced by water. It is noteworthy that protein content variation of the current study was not statistically significant between the 4 high-ranked families versus 4 low-ranked families (data not shown). The current WssGBLUP analysis indicated that thirteen SNPs in genes involved in protein degradation were involved in the additive genetic variance of moisture content (Table 4). Briefly, E3 ubiquitin-protein ligase RNF170 is an E3 ubiquitin-protein ligase that plays an essential role in the ubiquitination and degradation of inositol 1,4,5-trisphosphate receptor type 1 (ITPR1). The latter controls the calcium release from the endoplasmic reticulum [69], which affects the muscle protein content in rainbow trout [20] and has a profound effect on the regulation of cytoskeleton [70]. Cystatin-1, which possesses a peptidase inhibitor activity, had a single 5’UTR SNP. Thioredoxin-like 1 (TXNL1) had two synonymous SNPs. The knockdown of TXNL1 moderately stabilizes the ubiquitin-protein conjugates suggesting a connection between protein reduction and proteolysis [71]. Pre-mRNA-processing factor 19 (PRPF19) and ubiquitin-conjugating enzyme E2 D2 (UBE2D2) had four SNPs. These ligases catalyze polyubiquitin chain assembly and play a role in proteasomal protein degradation [72, 73]. Nuclear factor NF-kappa-B p105 subunit (NFKB1) had two 3’UTR SNPs. NFKB1 is involved in the negative regulation of cellular protein metabolic process [74] and apoptotic process [75].
In addition to the ubiquitin-protein ligases, SNPs in three genes involved in lysosomal/phagosomal pathways were identified. Ras-related protein rab7 (RAB7A) harbored a 3’UTR SNP. RAB7A is a major regulator of endo-lysosomal maturation/ trafficking and protein targeting to lysosome inducing autophagosome formation [76]. Thus, RAB7A positively regulates the protein catabolic process [77]. V-type proton ATPase subunit B (ATP6V1B2) had a 3’UTR SNP. V-ATPase is responsible for acidifying the intracellular compartments, including lysosomes [78]. The gene encoding the β chain of the adaptor protein-3 (AP-3) complex had a single synonymous SNP. Deletion in AP3B1 perturbs assembly of AP-3 complex and, in turn, trafficking of transmembrane lysosomal proteins [79].
Taken together, most of the common genomic loci affecting the highest proportion of the additive variance were involved in lipid metabolism, suggesting a common mechanism underlying intramuscular fat and moisture content and, partially, explaining the strong negative correlation between the fat and moisture content in this selectively bred rainbow trout population. Unique loci affecting moisture content were primarily involved in cytoskeleton regulations and protein turnover. Inhibition protease activity, such as calpains, reduced degradation of proteins responsible for cell membrane-cytoskeleton attachments and postmortem drip channel formation in muscle. The presence of calcium enhances proteolysis by µ-Calpain, of myofibrillar and other cytoskeletal proteins during postmortem storage [36].
Single marker GWA analyses
To identify single SNP markers associated with variation in fat and moisture content, we analyzed SNPs that passed QC filtration (n = 29,451) using a generalized score test; this test incorporates multiple covariates in the analysis and accounts for family structure using a kinship matrix [32]. In this study, 8 and 24 significant SNPs, surpassing the genome-wide significance level, had a potential impact on the fat and moisture content (Bonferroni-corrected p < 1.69E-06; Figs. 4 & 5 and Tables S5 & S6), respectively.
SNPs associated with the fat content were mainly located on chromosome 5 (n = 7) and SNP-harboring genes have roles in lipid metabolism (Table 5). The list includes 78 kDa glucose-regulated protein (GRP78), spindle and kinetochore associated complex subunit 1 (SKA1), apelin receptor B (APLNR-B), desmoplakin, podocan, and calcium-binding mitochondrial carrier protein SCaMC-1 (SLC25A24). Briefly, two missense mutations were identified in genes coding for GRP78 and SKA1. GRP78 is essential for adipocyte differentiation and a balanced secretion of adipokines. Deletion of GRP78 causes lipoatrophy in mice observed as a dramatic reduction in gonadal and subcutaneous adipose tissue [80]. SKA1 was downregulated in adipose tissues between samples from obese and healthy control children and has been suggested as a candidate biomarker for childhood obesity [81]. Three synonymous mutations were identified in genes encoding APLNR, desmoplakin, and podocan. APLNR knockout mice demonstrated excess fatty acid accumulation in skeletal muscle [82]. Abnormalities in desmoplakin have been associated with changes in lipid metabolism [83]. Podocan belongs to the small leucine-rich proteoglycans (SLRPs) that bind to low-density lipoprotein receptor-related protein (LRP-1) [84]. A 3’UTR SNP was identified in a gene coding for SLC25A24. Mice fed a high-fat diet exhibited increased expression level of SLC25A24; whereas, adipocyte differentiation was suppressed in Slc25a24- knockout [85].
Puromycin-sensitive aminopeptidase (NPEPPS), on chromosome 17, had a nonsynonymous SNP explaining the highest variability in fat content (R2 = 3.2%) (Table 5). NPEPPS impacts different physiological processes, including protein turnover and cell cycle regulation. NPEPPS was upregulated in mitten crabs fed with a linseed oil rich in linoleic acid [86]. However, the effect NPEPPS on lipid metabolism in fish needs further investigations.
SNPs associated with moisture content (n = 24) were associated with protein turnover, calcium metabolism, and cytoskeleton integrity (Table S5). Most of these SNPs (n = 11; ~46%) were located on chromosome 17. An SNP in a gene coding for acylphosphatase 2 ranked at the top of the list (R2 = 7.4%) (Table 5); however, its physiological role is not clear. Eight SNPs associated with moisture content were identified in five genes engaged in protein metabolism (Table 5). These genes are NPEPPS, eukaryotic initiation factor 4A-III (EIF4A3), eukaryotic translation initiation factor 4B (EIF4B), ribosome binding protein 1 (RRPB1), and F-box only protein 46 (FBXO46). Briefly, the aminopeptidase, NPEPPS, was associated with variation in fat and moisture content, suggesting a correlation between moisture and fat content. Five SNPs were identified in two genes encoding EIF4A3 and EIF4B, suggesting a role for the protein translation machinery in determining variation in moisture content. RRPB1 is an ER integral membrane protein implicated in polysome assembly and, therefore, protein synthesis [87]. RRPB1 has been suggested as essential in regulation of UPR signaling molecules and autophagy [88]. Finally, the F-box family SCF-E3 ubiquitin ligase, FBXO46, had a single 3’UTR SNP.
Expression of genes involved in the actin cytoskeleton and cytoskeletal organization was positively correlated with drip loss [89]. In this study, a total of 10 SNPs associated with moisture content were identified in seven genes engaged in cytoskeleton regulation (Table 5). These genes are encoding serum response factor (SRF), kinesin-1 heavy chain (KIF5B), inositol-trisphosphate 3-kinase C (IP3KC), supervillain (SVIL), calcineurin subunit B type 1 (PPP3R1), eIF4A and EIF4B. Briefly, SRF is a master regulator of the actin cytoskeleton [90]. Mutations in KIF5A caused cytoskeletal defects in humans [91]. Two, 3’UTR SNPs were identified in IP3KC and PPP3R1. IP3K and calcium/calcineurin signaling play critical roles in maintaining Ca2+ homeostasis that has a profound effect on the cytoskeleton [70, 92, 93]. Supervillin (SVIL) is one of the first components of the costameric membrane skeleton to assemble during muscle formation. It establishes a high-affinity connection between the membranes and actin cytoskeleton [94]. Translation initiation factors, including eIF4A and EIF4B, associate with the actin cytoskeleton, which affects protein synthesis [95].
In our previous work, we profiled transcriptome expression of fish families (YC 2010) showing contrasting phenotypes in fat content, which revealed only 17 differentially expressed transcripts associated with fat content [14]. About 90% of the genetic variation among individuals comes from SNPs [96], and therefore, identifying SNP markers associated with complex traits is most suitable for genetic evaluation in selection programs. Two previous GWA studies identified a few SNPs responsible for the additive variance for fat content in Atlantic salmon and Common Carp [1, 19]. The current GWA analysis identified a total of 137 SNPs in windows explaining at least 2% of the additive genetic variance for fat content, suggesting a better characterization of the genetic basis underlying variation in fat content. The discrepancies among the different GWA studies might be due to; 1) usage of different algorithms in the GWA studies, 2) variation in population size, 3) substantial difference in the capacity of the SNP arrays, 4) polygenic nature of intramuscular fat content, 5) different thresholds in each study including sliding window size [13].
Compared to our WssGBLUP analysis, the single marker GWA analysis revealed a smaller number of SNP markers associated with variation in intramuscular fat and moisture content. In addition, these two GWA approaches revealed different significant peaks associated with traits of interest. Aguilar et al. (2019) showed that the highest peak based on the p-value was not the same based on the proportion of variance explained, and this is because the latter depends on allele frequency, i.e., high effect but low frequency decreases the variance explained[97]. This result is consistent with other studied traits, such as fillet firmness, protein content [20], and bodyweight gain [98] in rainbow trout. The potential factors associated with observed heterogeneity between the two approaches are different algorithms, thresholds, and windows size used in each approach. The WssGBLUB was more effective than the single marker GWA in examining the genetic architecture of studied traits and identifying common QTL between traits. This method has proven to be optimal for breeding populations given the data structure: phenotyped individuals may not have genotypes, and there is a long history of pedigree recording [97]. Common QTL identified in this study may explain the high negative correlation between fat and moisture content. The recombinational progression of QTL and nearby markers determines the information content of haplotypes [99]. However, SNP-harboring genes identified by the two approaches had similar biological functions and were involved in lipid metabolism, protein turnover, and cytoskeletal remodeling. Routine use of single-SNP and multi-makers for GWA analysis was previously recommended to take advantage of the complete information content of the genotypes [99].