Novel Genomic Regions Associated with Intramuscular Fatty Acid Composition in Rabbits

Simple Summary A divergent selection experiment on intramuscular fat (IMF) content was carried out during nine generations in rabbits. The IMF content was successfully improved through generations. Besides, selection for IMF content generated a correlated response on its composition. Association analyses were performed to understand the genetic background of IMF composition using two rabbit lines divergently selected for IMF content. Several genomic regions and genes were identified, revealing the polygenic nature of the intramuscular fatty acid composition in rabbits. Abstract Intramuscular fat (IMF) content and its composition affect the quality of meat. Selection for IMF generated a correlated response on its fatty acid composition. The increase of IMF content is associated with an increase of its saturated (SFA) and monounsaturated (MUFA) fatty acids, and consequently a decrease of polyunsaturated fatty acids (PUFA). We carried out a genome wide association study (GWAS) for IMF composition on two rabbit lines divergently selected for IMF content, using a Bayes B procedure. Association analyses were performed using 475 individuals and 90,235 Single Nucleotide Polymorphisms (SNPs). The main objectives were to identify genomic regions associated with the IMF composition and to generate a list of candidate genes. Genomic regions associated with the intramuscular fatty acid composition were spread across different rabbit chromosomes (OCU). An important region at 34.0–37.9 Mb on OCU1 was associated with C14:0, C16:0, SFA, and C18:2n6, explaining 3.5%, 11.2%, 11.3%, and 3.2% of the genomic variance, respectively. Another relevant genomic region was found to be associated at 46.0–48.9 Mb on OCU18, explaining up to 8% of the genomic variance of MUFA/SFA. The associated regions harbor several genes related to lipid metabolism, such as SCD, PLIN2, and ERLIN1. The main genomic regions associated with the fatty acids were not previously associated with IMF content in rabbits. Nonetheless, MTMR2 is the only gene that was associated with both the IMF content and composition in rabbits. Our study highlighted the polygenic nature of the fatty acids in rabbits and elucidated its genetic background.


Introduction
Intramuscular fat (IMF) content and its fatty acid composition are key traits influencing meat quality [1]. The increase of IMF improves the juiciness, tenderness, and flavor of the meat [2]. IMF content and its fatty acid composition determine also the nutritional value of the meat. Individual saturated (SFA) and monounsaturated (MUFA) fatty acids could be biosynthesized. However, linoleic (precursor of n-6 fatty acids) and alpha linoleic (precursor of n-3 fatty acids) are both polyunsaturated fatty acids (PUFA) that cannot be biosynthesized, thus they must be obtained from the diet [3]. High intake of SFA increases low-density lipoprotein (LDL) cholesterol which is associated with cardiovascular diseases. In contrast, PUFA decrease LDL cholesterol [4]. Meat consumers are interested in the fatty acid composition of the meat, demanding healthy products rich in PUFA respect to SFA.
The high heritability of IMF and its substantial variability allowed its improvement through selection in cattle [5], chickens [6,7], and pigs [8]. In rabbits, a divergent selection experiment on IMF was carried out at the Universitat Politècnica de València during nine generations revealing a correlated response on its fatty acid composition [9,10]. The high-IMF line had greater percentages of MUFA than the low-IMF line. In contrast, the low-IMF line had greater percentages of PUFA [9].
Several association analyses were performed in order to unravel the genetic background of the IMF fatty acid composition. Previous Genome-Wide Association Studies (GWAS) reported genetic markers and genes associated with fatty acid composition in pigs [11][12][13] and cattle [14][15][16]. The main candidate genes were Stearoyl-CoA desaturase-1 (SCD), Leptin receptor (LEPR), Fatty Acid Synthase (FASN), and Fatty Acid Elongase (ELOVL). These studies underlined the complexity and the polygenic nature of the fatty acids in pigs [11][12][13] and cattle [14][15][16]. To our knowledge, there are no previous association analyses for the intramuscular fatty acid composition in rabbits. A previous study identified four genomic regions related to IMF content in rabbits, using two rabbit lines divergently selected for IMF [10]. IMF content and its fatty acid composition can share genomic regions with pleiotropic effects, since they have moderate to high genetic correlations [9]. In this study, we carried out a GWAS for intramuscular fatty acid composition, using the same previous rabbit lines divergently selected for IMF, in order to identify genomic regions associated with this trait and to detect putative candidate genes.

Animals
The present study was carried out on 478 individuals from the 9th generation of two lines divergently selected for IMF (239 rabbits each line). The base population had 13 sires and 83 does [17]. High (high-IMF) and low (low-IMF) lines had approximately 8 sires and 40 does per generation within each line, from the 1st to the 7th generation [9]. The latest generations (8th and 9th) comprised approximately 10 sires and 60 does per generation for each line [10]. Rabbits were weaned at 4 weeks of age and housed collectively until slaughter. Both lines were contemporarily reared in the same building having the same management and fed ad libitum with a standard commercial diet, containing 15.1% of crude protein, 14.5% of crude fiber and 2.8% of fat. The housing ventilation was controlled, and the photoperiod was constant (16:8 h). Two full sibs of the first parity of each doe were slaughtered at 9 weeks of age by exsanguination prior electrical stunning. Carcasses were chilled for 24 h at 4 • C. From each animal, Longissimus thoracis et lumborum muscle (LM) was excised and freeze-dried. IMF content was measured in LM by near-infrared reflectance spectroscopy (NIRS) and expressed in g/100 g of LM [18]. Twenty percent of the samples were chemically analyzed to test NIRS results. More details of the experiment can be found in previous works [9,17,18].

Phenotypes
Fatty acids were quantified by gas chromatography. First, we prepared fatty acid methyl esters (FAME) [19]. Further, FAME were analyzed by gas chromatography (FOCUS, Thermo, Milan, Italy), using a split/splitless injection, flame ionization detection, and a fused silica capillary column SPTM 2560 (Supelco, Bellefonte, PA, USA) (100 m × 0.25 mm × 0.2 µm film thickness). We injected FAME using a split ratio of 1/100 and Helium as a carrier gas (20 cm/s). The oven temperature was set at 140 • C during 5 min, then increased (4 • C/min) to 240 • C and maintained at 240 • C during 30 min. Detector and injector were at a temperature of 260 • C. We determined the fatty acids by comparing their retention times with the standards supplied by Supelco (Bellefonte, PA, USA), after that individual Animals 2020, 10, 2090 3 of 17 fatty acids were quantified using C13:0 as internal standard. The fatty acids contents were expressed as percentage of total fatty acids.

Genotypes and Quality Control
Rabbits were genotyped using the Affymetrix Axiom OrcunSNP Array (Affymetrix Inc., Santa Clara, CA, USA). The Single Nucleotide Polymorphisms (SNP) array contained 199,692 molecular markers. Quality control was performed using Axiom Analysis Suite v.4.0.3 of Thermo Fisher Scientific (Santa Clara, CA, USA). SNPs with a call rate greater than 0.95, a minor allele frequency greater than 0.05, and a known chromosome position were retained for the association analyses. In addition, individuals with a missing genotype frequency greater than 0.05 were excluded from the dataset. The missing genotypes were imputed by Beagle software v.4.1 [20] using 96 ancestors. A score of R 2 > 0.75 was used as quality parameter for imputed genotypes. After quality control, 475 animals (low-IMF: 236; high-IMF: 239) and 90,235 SNPs remained in the dataset.

Genome-Wide Association Study (GWAS)
The association analyses were performed using Bayesian procedures [21]. A Multiple-Marker Regression was performed under the following Bayes B model as in [21,22]: where y is the phenotypes vector. X is the incidence matrix for fixed effects and b is the vector of the fixed effects of month (five levels), sex (two levels), and parity order (three levels), with bounded flat prior density. z j is the genotypes vector for a SNP at locus j (j = 1, . . . , k, where k is the number of SNPs after quality control). α j is the random substitution effect for SNP j, distributed as t ν 0, σ 2 α , where σ 2 α was derived as in [21,22] from the genetic variance of each trait estimated using data of the fatty acids from all generations of the divergent selection experiment for IMF [9]. δ j is a random 0/1 variable (δ j = 1 represents the presence of SNP j with probability 1−π and δ j = 0 represents the absence of SNP j with probability π). e is the residual term, distributed as N 0, Iσ 2 e , where σ 2 e has a bounded flat prior density. It was assumed π = 0.9988, according to a previous GWAS for IMF in rabbits and to exploratory analyses for the studied traits [10]. Marginal posterior distributions of the model unknowns were estimated using Monte Carlo Markov chains. A total of 500,000 iterations were performed with a burn-in of 100,000 and a lag of 40. In this study, 1982 non-overlapping genomic windows of 1 Mb were a priori allocated to the 21 autosomes using the SNPs that remained after quality control. The windows have an average of 45 SNPs per 1 Mb window. The contributions of the windows to the genomic variance were computed as the posterior distribution of the genomic variance explained by the SNPs within each 1 Mb genomic window. All GWAS were performed using the GenSel software [22].
Bayes factors (BF) were estimated to evaluate the statistical relevance of the association between SNPs and traits. Bayes factors were calculated as follows [23]: wherep j is the posterior probability of SNP j.

Genomic Regions Associated with the Fatty Acid Composition
The genomic windows exceeding 1.0% of the genomic variance of the trait were considered to be associated with the trait. In addition, the genomic windows exceeding 0.5% of the genomic variance of the trait and having SNPs with a BF greater than 10 were also considered as associated with the trait. These thresholds of 1% and 0.5% represent 20 and 10 times the expected percentage of the genomic variance explained by each genomic window, respectively. SNPs with a BF greater than 10 were considered as associated with the trait [24]. The associated genomic windows were extended to ± 500 Kb from the first and last associated SNP in order to consider nearby associated SNPs, taking into account the linkage disequilibrium (LD). These new genomic windows are termed as "extended regions".

Identification of Candidate Genes
Genes were retrieved from the extended regions. Candidate genes were obtained from the Ensembl Genes 98 database using the Oryctolagus cuniculus as the reference genome [25]. The biological functions were retrieved from the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v.6.8 [26]. The genes related to lipid metabolism have been also investigated using the gene ontology (GO) [27].

Linkage Disequilibrium Analysis
In addition, we examined the LD in the most relevant extended regions to visualize the relationships between the SNPs. LD analyses were performed using plink software [28] and LDheatmap function from R [29].

Permutation Test
Results should be validated in order to overcome the problem of multiple testing. The phenotypes were resampled 100,000 times to remove any possible relationship between the phenotypes and genotypes, using the following covariates: sex, month, parity order, and the three first components of the principal component analysis of the genotypes (explaining 36% of the variance) in order to avoid population structures. Permutation testing was performed on the new data corresponding to the null hypothesis. Empirical p-values (EMP1) were provided for each SNP. EMP1 was estimated as K/N, where N is the total number of permutations and K is the number of times the SNP comes out associated in the simulated data. SNPs with EMP1 values close to 10 −5 were deemed true positives [30]. Permutation testing was performed using Plink software [28].

Genomic Data
The association analyses were performed using the 475 individuals and 90,235 SNPs remaining in the dataset after quality control. Allele frequencies differ among lines due to the random genetic drift and selection. As the animals came from a divergent selection experiment for IMF, genes with large causal variants will show extreme frequencies for alternate alleles between lines due to selection. We expect that high-IMF line will show genes with high frequencies given a reference allele, while the same genes will have lower frequencies in low-IMF line. Therefore, these genes will show intermediate frequencies when both lines are used. We did not consider the line as a fixed effect to take advantage from these extreme differences caused by selection and detect the genes with large causal variants. We also performed GWAS with the fixed effect of line in the model. The same associated genomic regions were obtained with lower percentages of the genomic variances of the traits. Divergent selection experiments can improve the detection power, according to simulation studies, especially for large causal variants and their nearby SNPs [31,32]. This statement is also certain for correlated traits. Thus, GWAS would detect in our experiment important genomic regions associated with the correlated traits, despite the small population size.
Many other genes were identified in the associated genomic regions; however, their functional annotations did not show a direct relationship with lipids metabolism. The current study identified several genomic regions associated with the fatty acid composition of IMF in rabbits, showing its polygenic nature, since the associated genomic regions did not explain large percentages of the genomic variances. We did not detect genomic regions with mayor effects like those found in pigs [11].   Many other genes were identified in the associated genomic regions; however, their functional annotations did not show a direct relationship with lipids metabolism. The current study identified several genomic regions associated with the fatty acid composition of IMF in rabbits, showing its polygenic nature, since the associated genomic regions did not explain large percentages of the genomic variances. We did not detect genomic regions with mayor effects like those found in pigs [11].

Linkage Disequilibrium Analysis
Linkage disequilibrium (LD) was assessed throughout the most relevant extended regions (at 34-37.9 Mb on OCU1, at 148.5-150.2 Mb on OCU3, and at 46.0-49.2 Mb on OCU18). Different LD patterns were detected (Figures 4-6). The extended region at OCU1 displayed several short blocks for both high-IMF and low-IMF lines, indicating the presence of several causal variants (QTL). Correlations between SNPs were greater in high-IMF line than in low-IMF line ( Figure 4). Besides, the genomic region at 148.5-150.2 Mb on OCU3 displayed different LD patterns between lines. Correlations between SNPs were greater in low-IMF line than in high-IMF line ( Figure 5). The QTL at 46.0-49.2 Mb on OCU18 formed a strong block (r 2 = 1) for both lines, while correlations between SNPs were greater for high-IMF line ( Figure 6). Thus, the putative QTLs were exposed to different selection pressures for IMF content.  On OCU1, a large genomic region at 34.0-37.9 Mb was associated with C14:0, C16:0, C18:0, and SFA, explaining 3.54%, 11.2%, 0.87%, and 11.3% of the genomic variances, respectively ( Table 2). The same region was associated with C18:2n6, C20:4n6, PUFA, PUFA/SFA, and n-6/n-3, explaining 3.18%, 1.57%, 2.41%, 2.98%, and 1.45% of the genomic variances, respectively. This genomic region harbors two genes related to lipids metabolism, alkaline ceramidase 2 (ACER2) and perilipin 2 (PLIN2). ACER2 promotes the hydrolysis of ceramides, the generation of sphingosine, and consequently the biosynthesis of sphingolipids [33,34]. Besides, the protein encoded by PLIN2 covers the intracellular lipid droplets and plays a major role in their stabilization [35]. In pigs, PLIN2 was considered as a useful marker for lean growth and its expression was positively associated with IMF content [36,37]. In rabbits, a selection signature for IMF was located in this genomic region, and presented high values for cross population-composite likelihood ratio (XP-CLR) and cross population-extended haplotype homozygosity (XP-EHH) [38]. Hence, this region harbors genes of pleiotropic effect given its association with the IMF content and its composition. Moreover, this region disclosed two validated SNPs with high BF. The first SNP at position 34.70 Mb was highly associated with C16:0 (BF = 177) and SFA (BF = 115), while the second SNP at 35.18 Mb was associated with SFA (BF = 59) (Table S1, Figure 1).
On the same chromosome (OCU1), the genomic window at 121.0-121.9 Mb was associated with C14:0, C18:0, C18:1n9, and MUFA. This region harbors the MTMR2 gene involved in lipid metabolic processes according to DAVID database. MTMR2 was also associated with the IMF content in rabbits, the genomic region at 120.6-121.9 Mb on OCU1 accounted for 2.03% of the genomic variance of IMF content [10]. Besides, MTMR2 gene was located at a signature of selection for IMF in pigs [39].
The genomic region at 25.0-26.9 Mb on OCU3 accounted for 6.05% of the genomic variance of linoleic acid (C18:2n6). Two candidate genes were retrieved, FGF1 and NR3C1. FGF1 was also considered as a candidate gene for IMF in chickens [40], while an expression Genome-Wide Association Study (eGWAS) reported the implication of NR3C1 gene in lipid metabolism in pigs [41]. On the same chromosome, the genomic region at 149-149.9 Mb was associated with almost all the studied traits except C16:1n7, C18:3n3, MUFA/SFA, and n-6/n-3. This region includes ENSOCUG00000000157 novel gene related to lipid biosynthesis. The novel gene is also known as ST3GAL1 in humans, pigs, and mice. The validated SNPs at this region showed high values of the BF (up to 141 for the ratio PUFA/SFA, 114 for C18:0 and 101 for C14:0) (Table S1). Moreover, these associations were corroborated by high values of the posterior probability of association (PPA). The PPA of C14:0, C16:0, and C18:0 were 0.96, 0.98, and 0.99, respectively ( Table 2).
Three different regions on OCU6 were associated with the fatty acid composition. The first region at the start of the chromosome was associated with C16:0, PUFA, and PUFA/SFA. However, it did not map to putative candidate genes related to lipid metabolism. The second genomic window at 5-5.9 Mb was associated with C18:2n6, and included two genes related to lipid metabolism (LITAF and SNX29). LITAF regulates the response to lipopolysaccharide according to gene oncology (GO: 0032496), whereas SNX29 affects phosphatidylinositol (family of lipids) binding (GO: 0035091). SNX29 was also a candidate gene for IMF in pigs [42]. The third region at 8-9.9 Mb was associated with C16:1n7, C18:1n9, and MUFA. This region harbors two genes related to lipid metabolism (SMG1 and GDE1) according to DAVID database.
The genomic region at 12.0-12.5 Mb on OCU9 accounted for 1.05% of the genomic variance of C16:1n7. This region harbors PPARG gene involved in adipocyte differentiation [44]. PPARG was also associated with fat deposition in sheep [45]. In pigs, two polymorphisms in the transcriptional regulatory region of this gene were related to IMF content [46]. Moreover, a correlation (0.64) was detected between mRNA PPARG and IMF content in Korean cattle [47]. On the same chromosome, the region at 64-65.9 Mb was found to be associated with PUFA and PUFA/SFA. This region harbors NPC1 and OSBPL1 genes. NPC1 is associated with lipid transport [48], while OSBPL1 is involved in both lipid transport (GO: 0006869) and binding (GO: 0008289).
On OCU10, the genomic region at 5-5.6 Mb was associated with the ratio MUFA/SFA. This region harbors the ITGB8 gene related to lipid metabolic process. This gene was also located at a signature of selection for IMF content in pigs [39].
Two genomic regions on OCU13 were found to be associated with the fatty acids. The first region at 33-34.9 Mb accounted for 3.42% of the genomic variance of linoleic acid (C18:2n6). This region includes five genes related to lipid metabolism according to DAVID database (CRP, DCAF8, FCER1A, PIGM, and ATP1A2). The second region at 126.0-126.9 Mb was associated with C18:0, C18:1n9, and MUFA, harboring HEYL, MFSD2A, and PPT1 genes. MFSD2A regulates the docosahexaenoic acid (C22:6n3) transport [49]. Validated associated SNPs at this region showed high values of the BF (96 for C18:0, 84 for C18:1n9, and 72 for MUFA; Table S1; Figure 3). Stearic (C18:0) and oleic (C18:1n9) acids were associated with the genomic window at 82-82.9 Mb on OCU14. This region harbors the ADIPOQ gene that encodes for the adiponectin involved in fatty acids metabolism [50]. ADIPOQ was also associated with IMF in cattle [51]. Moreover, a previous study reported a high correlation between its expression and IMF in cattle [52]. On the same chromosome, the genomic region at 78.2-78.9 Mb was associated with the ratio MUFA/SFA, presenting ATP11B gene involved in phospholipid transport (GO: 0015914).
The genomic region at 10-10.9 Mb on OCU15 was associated with C16:0, SFA, C18:2n6, PUFA, and the ratio PUFA/SFA. Two genes related to lipid metabolism mapped to this region (LRAT and TLR2). On the same chromosome, the genomic region at 70-71.9 Mb accounted for 3.59% of the genomic variance of linoleic acid (C18:2n6), harboring ANXA3 and GK2 genes related to lipid metabolism.
Two regions on OCU16 were associated with the fatty acids. The first region at 62.0-62.9 Mb was associated with C18:2n6 and MUFA/SFA, harboring HSD11B1 gene related to lipid metabolic process. The second region at 71-72.9 Mb was associated with C18:3n3 and n-6/n-3, explaining 1.67% and 1.04% of the genomic variances, respectively. This region harbors ENSOCUG00000006240 novel gene known as NR5A2 in humans. NR5A2 regulates cholesterol homeostasis (GO: 0042632).
The genomic region at 11-11.9 Mb on OCU18 was associated with the ratio MUFA/SFA, harboring SAMD8 gene related to ceramides biosynthetic process (GO: 0046513). Besides, a large genomic region at 46-48.9 Mb on OCU18 was associated with the majority of the traits, explaining 7.9% of the genomic variance of the ratio MUFA/SFA. This region harbors the following genes: GOT1, ERLIN1, ENSOCUG00000001375, and ENSOCUG00000014801. ERLIN1 promotes lipid binding; it encodes for proteins that bind to cholesterol and interacts with the sterol regulatory element binding protein (SREBP). This latter stimulates the cholesterol and fatty acid synthesis in liver [53,54]. Moreover, a previous GWAS reported that ERLIN1 was associated with lipid metabolism in cattle [55]. The novel genes ENSOCUG00000001375 and ENSOCUG00000014801 are known in pigs and humans as Stearoyl-CoA Desaturase (SCD). SCD gene is responsible for the biosynthesis of MUFA from SFA. A strong association between SCD and the fatty acid composition was reported in pigs [11,56], cattle [16], and goats [57]. In addition, an eGWAS study showed a high correlation (0.78) between SCD and PPARG expressions [58]. The latter gene was located on the associated genomic region at OCU9.
The genomic region at 21.0-21.9 Mb on OCU19 was associated with linoleic acid (C18:2n6) and disclosed the ADAP2 gene related to phosphatidylinositol bisphosphate binding (GO: 1902936).
Many other genes were identified in the associated genomic regions; however, their functional annotations did not show a direct relationship with lipids metabolism. The current study identified several genomic regions associated with the fatty acid composition of IMF in rabbits, showing its polygenic nature, since the associated genomic regions did not explain large percentages of the genomic variances. We did not detect genomic regions with mayor effects like those found in pigs [11].

Linkage Disequilibrium Analysis
Linkage disequilibrium (LD) was assessed throughout the most relevant extended regions (at 34-37.9 Mb on OCU1, at 148.5-150.2 Mb on OCU3, and at 46.0-49.2 Mb on OCU18). Different LD patterns were detected (Figures 4-6). The extended region at OCU1 displayed several short blocks for both high-IMF and low-IMF lines, indicating the presence of several causal variants (QTL). Correlations between SNPs were greater in high-IMF line than in low-IMF line ( Figure 4). Besides, the genomic region at 148.5-150.2 Mb on OCU3 displayed different LD patterns between lines. Correlations between SNPs were greater in low-IMF line than in high-IMF line ( Figure 5). The QTL at 46.0-49.2 Mb on OCU18 formed a strong block (r 2 = 1) for both lines, while correlations between SNPs were greater for high-IMF line ( Figure 6). Thus, the putative QTLs were exposed to different selection pressures for IMF content.
for both high-IMF and low-IMF lines, indicating the presence of several causal variants (QTL). Correlations between SNPs were greater in high-IMF line than in low-IMF line ( Figure 4). Besides, the genomic region at 148.5-150.2 Mb on OCU3 displayed different LD patterns between lines. Correlations between SNPs were greater in low-IMF line than in high-IMF line ( Figure 5). The QTL at 46.0-49.2 Mb on OCU18 formed a strong block (r 2 = 1) for both lines, while correlations between SNPs were greater for high-IMF line ( Figure 6). Thus, the putative QTLs were exposed to different selection pressures for IMF content.

Polygenic Nature of the Intramuscular Fatty Acid Composition
Taken together, the present study underlined different genomic regions associated with the fatty acid composition of IMF in rabbits. The results were corroborated having high values of Bayes factors (Figures 1-3). The present study highlighted also the polygenic nature of the fatty acid composition and its complexity in rabbits. Similar findings were reported by previous studies for the fatty acid composition in pigs [11][12][13] and cattle [14,16]. To our knowledge, this is the first GWAS for fatty acid composition in rabbits. A pervious study was performed on the IMF content [10]. Comparing the results, MTMR2 is the only gene that was associated with both IMF content and its composition. In addition, ACER2 and PLIN2 genes were associated with the fatty acid composition and were also detected at a selection signature for IMF content [38].

Conclusions
Association analyses for the fatty acid composition were carried out on two rabbit lines divergently selected for IMF content. Several associated regions were detected. The most relevant regions were found at 34-37.9 Mb on OCU1, at 149.0-149.9 Mb on OCU3, and at 46.0-48.9 Mb on OCU18. The associated regions disclose several genes related to lipid metabolism such as SCD,

Polygenic Nature of the Intramuscular Fatty Acid Composition
Taken together, the present study underlined different genomic regions associated with the fatty acid composition of IMF in rabbits. The results were corroborated having high values of Bayes factors (Figures 1-3). The present study highlighted also the polygenic nature of the fatty acid composition and its complexity in rabbits. Similar findings were reported by previous studies for the fatty acid composition in pigs [11][12][13] and cattle [14,16]. To our knowledge, this is the first GWAS for fatty acid composition in rabbits. A pervious study was performed on the IMF content [10]. Comparing the results, MTMR2 is the only gene that was associated with both IMF content and its composition. In addition, ACER2 and PLIN2 genes were associated with the fatty acid composition and were also detected at a selection signature for IMF content [38].

Conclusions
Association analyses for the fatty acid composition were carried out on two rabbit lines divergently selected for IMF content. Several associated regions were detected. The most relevant regions were found at 34-37.9 Mb on OCU1, at 149.0-149.9 Mb on OCU3, and at 46.0-48.9 Mb on OCU18. The associated regions disclose several genes related to lipid metabolism such as SCD, PLIN2, ERLIN1, and PPARG. The main genomic regions in which we found genes related to lipid Taken together, the present study underlined different genomic regions associated with the fatty acid composition of IMF in rabbits. The results were corroborated having high values of Bayes factors (Figures 1-3). The present study highlighted also the polygenic nature of the fatty acid composition and its complexity in rabbits. Similar findings were reported by previous studies for the fatty acid composition in pigs [11][12][13] and cattle [14,16]. To our knowledge, this is the first GWAS for fatty acid composition in rabbits. A pervious study was performed on the IMF content [10]. Comparing the results, MTMR2 is the only gene that was associated with both IMF content and its composition. In addition, ACER2 and PLIN2 genes were associated with the fatty acid composition and were also detected at a selection signature for IMF content [38].

Conclusions
Association analyses for the fatty acid composition were carried out on two rabbit lines divergently selected for IMF content. Several associated regions were detected. The most relevant regions were found at 34-37.9 Mb on OCU1, at 149.0-149.9 Mb on OCU3, and at 46.0-48.9 Mb on OCU18. The associated regions disclose several genes related to lipid metabolism such as SCD, PLIN2, ERLIN1, and PPARG. The main genomic regions in which we found genes related to lipid metabolism were not detected in our previous experiment for IMF. MTMR2 is the only gene that was associated with both the IMF content and its composition. To our knowledge, this is the first GWAS for fatty acid composition in rabbits. Our study highlighted the polygenic nature of the fatty acid composition as no large percentages of the genomic variances were explained by the associated genomic regions. Further analyses throughout the associated regions would be needed in order to validate their importance.

Conflicts of Interest:
The authors declare no conflict of interest.