# Bivariate GWAS performed on rabbits divergently selected for intramuscular fat content reveals pleiotropic genomic regions and genes related to meat and carcass quality traits

Bolívar Samuel Sosa-Madrid, Agostina Zubiri-Gaitán, Noelia Ibañez-Escriche, Agustín Blasco, Pilar Hernández

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bolívar Samuel Sosa-Madrid, Agostina Zubiri-Gaitán, Noelia Ibañez-Escriche, Agustín Blasco, Pilar Hernández. Bivariate GWAS performed on rabbits divergently selected for intramuscular fat content reveals pleiotropic genomic regions and genes related to meat and carcass quality traits. Genetics Selection Evolution, 2025, 57 (1), pp.36. 10.1186/s12711-025-00971-5. hal-05160423

HAL Id: hal-05160423

https://hal.science/hal-05160423v1

Submitted on 14 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Bivariate GWAS performed on rabbits divergently selected for intramuscular fat content reveals pleiotropic genomic regions and genes related to meat and carcass quality traits

Bolívar Samuel Sosa-Madrid¹*†, Agostina Zubiri-Gaitán¹†@, Noelia Ibañez-Escriche¹, Agustín Blasco¹ and Pilar Hernández¹*

## Abstract

**Background** Meat quality plays an important economic role in the meat industry and livestock breeding programmes. Intramuscular fat content (IMF) is one of the main meat quality parameters and its genetic improvement has led breeders to investigate its genomic architecture and correlation with other relevant traits. Genetic markers associated with causal variants for these traits can be identified by bivariate analyses. In this study, we used two rabbit lines divergently selected for IMF to perform bivariate GWAS with the aim of detecting pleiotropic genomic regions between IMF and several weight, fat, and meat quality traits. Additionally, whole-genome sequencing data from these lines were used to identify potential causal variants associated with the genetic markers.

**Results** The main pleiotropic region was found on Oryctolagus cuniculus chromosome (OCC) 1 between 35.4 Mb and 38.2 Mb, explaining up to 2.66% of the IMF genetic variance and being associated with all traits analysed, except muscle lightness. In this region, the potentially causal variants found pointed to *PLIN2*, *SH3GL2*, *CNTLN*, and *BNC2* as the main candidate genes affecting the different weight, fat depots and meat quality traits. Other relevant pleiotropic regions found were those on OCC3 (148.94–150.89 Mb) and on OCC7 (27.07–28.44 Mb). The first was associated with all fat depot traits and explained the highest percentage of genetic variance, up to 10.90% for scapular fat. Several allelic variants were found in this region, all located in the novel gene *ENSOCUG00000000157* (orthologous to *ST3GAL1* in other species), involved in lipid metabolism, suggesting it as the main candidate affecting fat deposition. The region on OCC7 was associated with most meat quality traits and explained 8.48% of the genetic variance for pH. No allele variants were found to segregate differently between the lines in this region; however, it remains a promising region for future functional studies.

<footnote>†Bolívar Samuel Sosa-Madrid and Agostina Zubiri-Gaitán share the first position of the current article.</footnote>

<footnote>*Correspondence:
Bolívar Samuel Sosa-Madrid
bosomad@upvnet.upv.es
Pilar Hernández
phernan@dca.upv.es</footnote>
**Conclusions** Our results showed that bivariate models assuming pleiotropic effects are valuable tools to identify genomic regions simultaneously associated with IMF and several weight, fat and meat quality traits. Overall, our results provided relevant insights into the correlations and relationships between traits at the genomic level, together with potential functional mutations, which would be relevant for exploration in rabbit and other livestock breeding programmes.

## Background

Meat quality has an increasingly important economic role in the meat industry and livestock breeding programmes [1], in which the ultimate goal is to meet consumers' demands [2, 3]. It is a complex trait determined by multiple characteristics, including nutritional properties, organoleptic characteristics, and technological factors. The intramuscular fat content (IMF) is considered as one of the main meat quality parameters, since it affects several of the aforementioned characteristics, mainly the sensory properties such as flavour, juiciness and tenderness [3, 4]. Even though the relationship between IMF and these sensory traits varies between studies, higher IMF content is generally associated with greater tenderness and juiciness, and better flavour [5, 6]. However, increasing the IMF content can also affect the nutritional quality of meat, which is relevant for consumers' acceptance, and it is also related to carcass quality traits, especially because of its positive correlation with carcass fat content [1, 2] which could lead to a deterioration of carcass quality.

The importance of IMF as a promising selection criterion has led to the development of several genetic experiments for example in beef cattle [7, 8], pigs [9, 10], chicken [11, 12], sheep [13], and rabbits [14, 15]. In fact, a successful divergent selection experiment for IMF was carried out in rabbits [14, 15], and genomic analyses have identified genomic regions and candidate genes associated with IMF [16, 17] and fatty acid composition [16]. Other traits were also measured in this selection experiment, such as carcass weight, fat depots, pH and colour parameter traits; and given their correlated responses and their moderate-to-high genetic correlations with IMF [14, 15], the selection for IMF could also affect causal variants and genes related to those traits. When a causal variant affects two or more traits in the same individual or species, this variant has a pleiotropic effect [18]. In breeding programmes, pleiotropic effects contribute to genetic covariances, together with linkage disequilibrium. Therefore, the discovery of genomic regions associated with favourable and unfavourable effects is important for understanding the genomic architecture of economically important traits, such as meat and carcass quality traits. Moreover, genetic markers linked to

causal variants associated with IMF and other quality traits can be identified by bivariate analyses, which are assumed to have a greater power of detection [19, 20].

The aim of this study was to perform bivariate GWAS analyses to identify pleiotropic genomic regions associated with IMF and several carcass and meat quality traits. For this purpose, data from the two rabbit lines divergently selected for IMF were used and the results were compared with those obtained with univariate GWAS. Further enrichment, whole-genome sequencing analyses and visualization of genes on the associated genomic regions were also performed to integrate the results.

## Methods

### Animals and phenotypes

Samples from two rabbit lines divergently selected for intramuscular fat content (IMF) in the Longissimus thoracis et Lumborum (LTL) muscle were used to perform this study. This divergent selection experiment was carried out at the Universitat Politècnica de València and led to the creation of two lines: one with high-IMF content (H) and one with low-IMF content (L). Each line comprised 8–10 sires and 40–60 dams per generation; they were contemporarily reared under the same environmental and management conditions and were fed the same diet. The selection criterion was the IMF content measured at 9 weeks of age in two full sibs (one male and one female) of the first parity of each doe. Further details about the selection procedure can be found in Zomeño et al. and Martínez-Álvaro et al. [14, 15].

A total of 480 rabbits, 240 from the H line and 240 from the L line, from the 9th generation of selection were used in this study. The animals were slaughtered at 9 weeks of age using exsanguination after electrical stunning. A total of ten carcass traits (including weights, fat depots and colour), and five meat quality traits (including IMF, pH and colour) were recorded. The body weight (BW) of the rabbits was recorded before slaughter. After slaughter, the carcasses were chilled for 24 h at 4 °C, and the chilled carcass weight (CCW) and liver weight (LW) were then recorded. Additionally, the reference carcass weight (RCW) was recorded as the weight of the carcass without the head, liver, lungs, thymus, oesophagus, heart, and kidneys [21]. The scapular (SF) and perirenal fat (PF) were excised and weighted, and the dissectible fat
percentage (DF) was calculated as the sum of SF and PF, relative to the RCW. Finally, the carcass colour parameters lightness (CL*), redness (Ca*) and yellowness (Cb*) were measured on the surface of the fourth lumbar vertebra using a CR-400 Chroma Meter (Konica Minolta, Inc.).

To measure the meat quality traits, the LTL muscle was first dissected. The muscle pH was measured at the level of the fifth lumbar vertebra using a Crison pH-meter Basic +20 (Crison Instruments, Barcelona, Spain), and the colour parameters lightness (LL*), redness (La*), and yellowness (Lb*) were measured at the seventh lumbar vertebra transversal section using the same colorimeter mentioned above. Finally, the muscle was minced and freeze dried, and the IMF content was quantified using near-infrared spectroscopy (NIRS; FOSS NIRSystems Inc., Hilleroed, Denmark), applying the equations developed by Zomeño et al. [22], and expressed as g IMF/100 g of muscle. Twenty percent of the samples were chemically analysed to test NIRS results.

Heritabilities, genetic correlations with IMF, and (co) variance components were estimated by fitting the bivariate animal model described by Martínez-Álvaro et al. [15] and using the TM software [23]. The dataset used comprised 2244 phenotypic records and 3676 animals in the pedigree from all generations.

## Genotypes and quality control

The DNA was extracted from samples of the obliquus abdominis muscle using standard procedures, and the genotyping was performed using the 200 K Affymetrix Axiom OrcunSNP array (ThermoFisher Scientific). A quality control analysis was performed by Axiom Analysis Suite v. 4.0.3 and only the SNPs with a call rate ≥ 0.95, a minor allele frequency (MAF) ≥ 0.05, and a known autosomal chromosome position according to OryCun2.0 assembly, were used in the association analyses. Moreover, animals with more than 3% of missing genotypes or that failed the Mendelian inheritance test were excluded; in the latter case, only the offspring failing the test was removed, while the parent was retained. The remaining missing genotypes were imputed using the software FImpute [24], inferring the haplotype using the complete genealogical information of the lines and also the genomic information of the parents (the 8th generation) obtained using the same aforementioned array. Finally, 477 animals (240 from the H and 237 from the L line) and 89,968 SNPs were retained for the association analyses. The vast majority of SNPs were removed because they were fixed in our populations, likely reflecting the fact that the SNP chip

was not developed using rabbits from our specific lines, resulting in limited polymorphism at many loci.

## Statistical analyses

### GWAS

Both univariate and bivariate GWAS were performed using Bayesian multiple-marker regression (BMMR) under the BayesB model. The bivariate GWAS included one carcass or meat quality trait as the first trait, and IMF as the second one. The model was defined as follows:

$$ 
\begin{bmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \end{bmatrix} = \begin{bmatrix} \mathbf{X}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_2 \end{bmatrix} \begin{bmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \end{bmatrix} + \begin{bmatrix} \mathbf{W}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{W}_2 \end{bmatrix} \begin{bmatrix} \mathbf{c}_1 \\ \mathbf{c}_2 \end{bmatrix} \\ 
+ \begin{bmatrix} \mathbf{Z}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{Z}_2 \end{bmatrix} \begin{bmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \end{bmatrix} + \sum_{j=1}^{k} \mathbf{z}_{ji} \mathbf{D}_j \alpha_j + \begin{bmatrix} \mathbf{e}_1 \\ \mathbf{e}_2 \end{bmatrix}, \quad (1)
$$

where the subscript 1 refers to the analysed carcass or meat quality trait and the subscript 2 refers to the IMF trait; $y_1$ and $y_2$ are the phenotype vectors; $b_1$ and $b_2$ are the vectors of fixed effects, which included month (five levels), sex (two levels), and parity order (two levels); and $c_1$ and $c_2$ are the common litter random effects. Bearing in mind the missing heritability, a pedigree-based random polygenic effect represented by $\mathbf{u}_1$ and $\mathbf{u}_2$ was included in the model; $X_1$, $X_2$, $W_1$, $W_2$, $Z_1$ and $Z_2$ are the incidence matrices relating the corresponding effects to the phenotypes. $k$ is the total number of SNPs after quality control; $z_{ji}$ is the vector containing the genotypic covariate for each SNP or locus $j$ (0, 1 or 2 reference alleles) for individual $i$; $\alpha_j$ is the random substitution effect for SNP $j$; $D_j$ is a diagonal matrix with elements $\text{diag}(D_j) = (\delta_{j1}, \delta_{j2})$, where $\delta_j$ is an indicator variable indicating if the marker effect of locus $j$ for each trait is zero or not. Further details regarding $\delta_j$ are explained below. Finally, $\mathbf{e}_1$ and $\mathbf{e}_2$ are the residual effects. The effects were assumed to be uncorrelated between each other.

The model was assumed to be conditionally distributed as a multivariate normal. Ordering the residuals by individual, the common litter random effects were distributed as $c \sim N(0, I_c \otimes C)$, and the polygenic random effects as $u \sim N(0, A \otimes U)$; in which $C$ and $U$ are $2 \times 2$ (co)variances matrices between the two traits for the common litter ($C$) and the polygenic ($U$) effects; $A$ is the full relationship matrix and $I_c$ is the identity matrix of the same order as the number of levels of common litter effects. Finally, ordering by individual, the residuals were distributed as $e \sim N(0, I_i \otimes R)$, where $I_i$ is an identity matrix with the same order as the number of individuals and $R$ is the residual (co)variance matrix structured as follows:
$$ \mathbf{R} = \begin{bmatrix} \sigma_{e1}^2 & \sigma_{e2,1} \\ \sigma_{e1,2} & \sigma_{e2}^2 \end{bmatrix}. \qquad (2) $$

All (co)variance matrices were a priori assumed to follow an inverse Wishart distribution as described in Cheng et al. [20].

The BayesB model assumes that, at any given MCMC iteration, the majority of SNPs have zero effect on the phenotype (i.e., $\delta_j = 0$, with a probability of $\pi$) and are therefore not included in the model [25]. In a bivariate GWAS, each SNP can affect any combination of traits. More specifically, there are four possible combinations, SNP j: (a) has zero effect on both traits ($\delta_{j1}=0, \delta_{j2}=0$); (b) has an effect only on the first trait ($\delta_{j1}=1, \delta_{j2}=0$); (c) has an effect only on the second trait ($\delta_{j1}=0, \delta_{j2}=1$); or (d) has an effect on both traits ($\delta_{j1}=1, \delta_{j2}=1$). After some exploratory analysis, performed following the methodology described in previous studies [16, 26], the joint probability (a) of SNPs with zero effect on both traits was set to 0.9982, which is the product of $\pi1$ and $\pi2$. The joint probability of the remaining combinations (b), (c) and (d) was proportionally distributed, each with a probability of 0.0006 (i.e., the remaining probability of 0.0018 divided by 3). The same model was used for the univariate GWAS of each carcass and meat quality trait, in which the $\pi$ was set to 0.9988, as the one defined in the previous GWAS for IMF performed using these rabbit lines [16]. Finally, the effects of the SNPs that are not zero were assumed to be multinormally distributed as $N(0, I_j \otimes G)$, where $I_j$ is an identity matrix with the same order as the number of SNPs and G is the (co)variance matrix of the random substitution effect structured as:

$$ \mathbf{G} = \begin{bmatrix} \sigma_{\alpha j 1}^2 & \sigma_{\alpha j 2, 1} \\ \sigma_{\alpha j 1, 2} & \sigma_{\alpha j 2}^2 \end{bmatrix}, \qquad (3) $$

which a priori assumed to follow an inverse Wishart distribution as described in Cheng et al. [20]. All effects were assumed to be independent from each other.

The marginal posterior distributions of the model parameters were estimated by Gibbs sampling, based on MCMC, consisting of 470,000 iterations, with a burn-in of 70,000. Only one sample every 40 iterations was saved to avoid the high correlation between consecutive samples. The prior (co)variance components were those previously estimated using the data from all generations, as explained above. Both univariate and bivariate GWAS were performed with the JWAS software [20]. The percentage of genetic variance was calculated by 1-Mb sliding genomic windows, allocated to the 21 autosomes, as this approach simplifies the tuning of the associated genomic regions. Genomic windows were considered to be associated with the analysed trait when they explained

at least 0.75% of the genetic variance [19]. This threshold represents 14 times the expected genetic variance under an infinitesimal model (89,176 sliding genomic windows). The associated overlapping genomic windows were then grouped in only one genomic region to tune the position of the region in the rabbit genome. This sliding genomic window approach allows to better capture genomic regions associated with a trait than when using single marker regression analysis with Bayes factors [27]. As the selection increased the LD block between the IMF causal variants and the surrounding SNPs, the identification of associated genomic regions may gradually decrease, because the estimation of a causal variant effect would be shared among all surrounding SNPs that unlikely exceed a Bayes factor threshold. The posterior probability of association (PPA) for each SNP was then computed as the proportion of iterations in which the SNP was included in the model with a nonzero effect, relative to the total number of iterations. However, this parameter was not used as a threshold criterion given its dependence on $\pi$ and the marker numbers in BayesB models, both at the chromosome- and genome-wide levels [26].

As the selection criterion was IMF, which is correlated with other meat and carcass quality traits, we hypothesised that a bivariate GWAS analysis [19, 28] would detect the regions affected by IMF selection as well as those moderately correlated with carcass and meat quality traits in rabbits. Thus, the genetic correlation between the two analysed traits for each sliding genomic window was estimated as the mean of the marginal posterior distribution of the correlation between the sampled breeding values for the two traits in the genomic region [19].

### Identification of candidate genes and functional mutations

Candidate genes allocated in the associated genomic windows were obtained from the Ensembl Genes 111 database, using the Oryctolagus cuniculus as the reference genome (OryCun 2.0.) [29]. The biological functions and pathways of the genes were retrieved from the Database for Annotation, Visualization and Integrated Discovery (DAVID) v.6.8 [30]. The genes related to retrieved pathways were also investigated using the gene ontology (GO) [31].

Whole-genome sequencing (WGS) from the 20 breeding males of the 8th generation, corresponding to the sires of the animals used in the GWAS, was performed to identify the genomic variants present in the associated genomic regions that segregated differently between the lines, indicating possible functional mutations. The sequencing, pre-processing, and variant calling were performed at the Centro Nacional de Análisis Genómico (CNAG-CRG Barcelona, Spain). Briefly, sequencing was performed using Illumina Technology; and the WGS data
was pre-processed following the methodology described by Wright, et al. [32]. Illumina adapters and low-quality ends were removed using Trimmomatic v0.39 [33] and the BWA-MEM algorithm was then applied to align the reads to the OryCun v2.0 reference genome [34]. Variant calling was performed applying the GATK Best Practices pipeline [35], using HaplotypeCaller and GenotypeGVCF, to finally obtain the VCF files. Variants were annotated using the snpEff software [36] and only those SNVs and INDELs with moderate to high impact were considered, which included: missense and frameshift variants, start lost, stop gained, conservative in frame deletions and insertions [37]. Finally, we identified as relevant (1) the allele variants that were fixed in one line and segregating at an allele frequency lower than 0.5 in the other, as these suggest directional selection or lineage-specific fixation, and (2) the variants whose alleles were segregating at opposite frequencies between lines (greater than 0.65), indicating substantial differentiation. These thresholds were chosen to balance stringency with sensitivity, allowing us to capture biologically meaningful divergence while minimizing noise from minor frequency shifts.

## Results

Descriptive statistics of the studied traits before the correction for systematic effects are shown in Table 1, while the differences between lines can be found in previous studies [15, 38]. Heritabilities, genetic correlations with

IMF, and the common litter random effects are shown in Table 2.

Regarding the GWAS results, there was a general increase in the number of associated SNPs detected when performing the bivariate GWAS than when performing the univariate one for each trait considered (Tables 3, 4, and 5). As mentioned above, in this study we focussed on the genomic regions found to be associated by genomic window analysis, as only one or few SNPs exceeding a SNP threshold in the single marker regression method are likely to be false positives, even more so if the SNP was not found to be associated by genomic window analysis. The associated overlapping genomic windows have been grouped in only one genomic region to tune the position of the region in the rabbit genome. Therein, this tuned region is defined as associated genomic region. For simplicity and better visualization, the results are presented graphically or in tables based on overlapping sliding windows across GWAS analyses, i.e., the associated genomic regions across all traits have been clustered into: (a) weight traits: BW, CCW, RCW and LW; (b) fat depot traits: PF, SF and DF; and (c) meat quality traits: all colour parameter traits and pH. The bivariate GWAS allowed to detect genomic regions associated with both traits included in the analysis (i.e., regions with pleiotropic effects; Table 6), but also genomic regions associated with only one of the traits (Tables 3, 4, and 5).

<table><caption>**Table 1** Descriptive statistics of carcass and meat quality traits in the 9th generation of divergent selection experiment</caption><thead><tr><th rowspan="2">Trait</th><th colspan="5">High-IMF line</th><th colspan="5">Low-IMF line</th></tr><tr><th>NumA<sup>a</sup></th><th>Mean</th><th>SD<sup>b</sup></th><th>Range<sup>c</sup></th><th>CV%</th><th>NumA<sup>a</sup></th><th>Mean</th><th>SD<sup>b</sup></th><th>Range<sup>c</sup></th><th>CV%</th></tr></thead><tbody><tr><td>IMF (g/100 g)</td><td>240</td><td>1.32</td><td>0.21</td><td>[0.87, 1.94]</td><td>15.85</td><td>237</td><td>0.80</td><td>0.06</td><td>[0.63, 1.08]</td><td>8.01</td></tr><tr><td>BW (g)</td><td>240</td><td>1828</td><td>140</td><td>[1502, 2096]</td><td>7.66</td><td>237</td><td>1777</td><td>144</td><td>[1520, 2168]</td><td>8.13</td></tr><tr><td>CCW (g)</td><td>240</td><td>975</td><td>104</td><td>[703, 1219]</td><td>10.68</td><td>237</td><td>933</td><td>116</td><td>[709, 1312]</td><td>12.46</td></tr><tr><td>RCW (g)</td><td>240</td><td>848</td><td>79</td><td>[636, 1004]</td><td>9.28</td><td>235</td><td>808</td><td>89</td><td>[639, 1088]</td><td>10.96</td></tr><tr><td>LW (g)</td><td>237</td><td>57.81</td><td>8.97</td><td>[39.85, 81.07]</td><td>15.52</td><td>235</td><td>51.24</td><td>7.53</td><td>[37.09, 75.83]</td><td>14.69</td></tr><tr><td>SF (g)</td><td>234</td><td>4.60</td><td>1.24</td><td>[1.5, 8.52]</td><td>26.95</td><td>236</td><td>3.60</td><td>1.48</td><td>[0.96, 10.54]</td><td>41.01</td></tr><tr><td>PF (g)</td><td>239</td><td>11.22</td><td>3.53</td><td>[1.75, 21.58]</td><td>31.49</td><td>237</td><td>5.99</td><td>2.89</td><td>[1.18, 16.19]</td><td>48.28</td></tr><tr><td>DF (%)</td><td>234</td><td>1.85</td><td>0.44</td><td>[0.48, 3.04]</td><td>23.93</td><td>234</td><td>1.16</td><td>0.42</td><td>[0.36, 2.58]</td><td>36.22</td></tr><tr><td>CL*</td><td>239</td><td>51.65</td><td>1.89</td><td>[46.08, 56.72]</td><td>3.65</td><td>236</td><td>51.68</td><td>1.92</td><td>[46.79, 58.77]</td><td>3.72</td></tr><tr><td>Ca*</td><td>239</td><td>1.98</td><td>0.68</td><td>[0.65, 4.61]</td><td>34.33</td><td>236</td><td>2.04</td><td>0.79</td><td>[0.69, 6.09]</td><td>38.79</td></tr><tr><td>Cb*</td><td>239</td><td>-1.99</td><td>1.41</td><td>[-5.74, 1.94]</td><td>NA</td><td>236</td><td>-2.18</td><td>1.69</td><td>[-7.46, 2.28]</td><td>NA</td></tr><tr><td>pH</td><td>239</td><td>5.72</td><td>0.10</td><td>[5.38, 5.97]</td><td>1.73</td><td>237</td><td>5.73</td><td>0.10</td><td>[5.37, 6.05]</td><td>1.78</td></tr><tr><td>LL*</td><td>240</td><td>50.63</td><td>1.82</td><td>[45.55, 57.08]</td><td>3.60</td><td>237</td><td>50.72</td><td>2.19</td><td>[44.87, 58.34]</td><td>4.33</td></tr><tr><td>La*</td><td>240</td><td>2.24</td><td>0.92</td><td>(-0.14, 5.11)</td><td>NA</td><td>237</td><td>2.28</td><td>1.06</td><td>[0.23, 7.53]</td><td>46.60</td></tr><tr><td>Lb*</td><td>240</td><td>1.69</td><td>0.66</td><td>(-0.83, 3.72)</td><td>NA</td><td>237</td><td>1.62</td><td>0.72</td><td>(-1.95, 3.91)</td><td>NA</td></tr></tbody></table>

NumA*: number of animals after quality control; SD*: standard deviation; Range*: includes the minimum (left) and the maximum values

BW: body weight, CCW: commercial carcass weight, RCW: reference carcass weight, LW: liver weight, SF: scapular fat, PF: perirenal fat, DF: dissectible fat percentage, CL*: carcass lightness, Ca*: carcass redness, Cb*: carcass yellowness, LL*: muscle lightness, La*: muscle redness, Lb*: muscle yellowness
<table><caption>**Table 2** Heritabilities of carcass and meat quality traits and their genetic correlations with intramuscular fat content</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">Heritability</th><th colspan="3">Genetic correlation</th><th colspan="2">C<sup>2</sup> of common litter effect<sup>c</sup></th></tr><tr><th>Mean</th><th>HPD<sub>95%</sub><sup>a</sup></th><th>Mean</th><th>HPD<sub>95%</sub><sup>a</sup></th><th>P<sub>0</sub><sup>b</sup></th><th>Mean</th><th>HPD<sub>95%</sub><sup>a</sup></th></tr></thead><tbody><tr><td colspan="8"><strong>Carcass-weight traits</strong></td></tr><tr><td>BW</td><td>0.28</td><td>[0.13, 0.43]</td><td>0.06</td><td>[-0.24, 0.34]</td><td>0.67</td><td>0.28</td><td>[0.22, 0.35]</td></tr><tr><td>CCW</td><td>0.26</td><td>[0.13, 0.41]</td><td>0.09</td><td>[-0.20, 0.39]</td><td>0.72</td><td>0.29</td><td>[0.22, 0.36]</td></tr><tr><td>RCW</td><td>0.30</td><td>[0.15, 0.46]</td><td>0.06</td><td>[-0.23, 0.35]</td><td>0.65</td><td>0.27</td><td>[0.21, 0.34]</td></tr><tr><td>LW</td><td>0.38</td><td>[0.15, 0.64]</td><td>0.18</td><td>[-0.16, 0.52]</td><td>0.84</td><td>0.24</td><td>[0.14, 0.33]</td></tr><tr><td colspan="8"><strong>Fat depot traits</strong></td></tr><tr><td>SF</td><td>0.43</td><td>[0.29, 0.57]</td><td>0.28</td><td>[0.04, 0.51]</td><td>0.98</td><td>0.14</td><td>[0.08, 0.20]</td></tr><tr><td>PF</td><td>0.58</td><td>[0.43, 0.75]</td><td>0.31</td><td>[0.09, 0.53]</td><td>0.99</td><td>0.15</td><td>[0.09, 0.21]</td></tr><tr><td>DF</td><td>0.66</td><td>[0.51, 0.80]</td><td>0.34</td><td>[0.12, 0.55]</td><td>1.00</td><td>0.11</td><td>[0.06, 0.17]</td></tr><tr><td colspan="8"><strong>Meat quality traits</strong></td></tr><tr><td>IMF</td><td>0.49</td><td>[0.35, 0.63]</td><td>–</td><td>–</td><td>–</td><td>0.17</td><td>[0.11, 0.23]</td></tr><tr><td>CL*</td><td>0.18</td><td>[0.08, 0.31]</td><td>– 0.11</td><td>[-0.41, 0.21]</td><td>0.75</td><td>0.31</td><td>[0.24, 0.37]</td></tr><tr><td>Ca*</td><td>0.36</td><td>[0.20, 0.51]</td><td>0.00</td><td>[-0.27, 0.28]</td><td>0.50</td><td>0.19</td><td>[0.13, 0.26]</td></tr><tr><td>Cb*</td><td>0.11</td><td>[0.01, 0.22]</td><td>0.11</td><td>[-0.27, 0.47]</td><td>0.73</td><td>0.28</td><td>[0.21, 0.36]</td></tr><tr><td>pH</td><td>0.07</td><td>[0.00, 0.14]</td><td>0.10</td><td>[-0.25, 0.46]</td><td>0.71</td><td>0.32</td><td>[0.26, 0.37]</td></tr><tr><td>LL*</td><td>0.23</td><td>[0.09, 0.37]</td><td>0.00</td><td>[-0.29, 0.32]</td><td>0.51</td><td>0.33</td><td>[0.26, 0.40]</td></tr><tr><td>La*</td><td>0.32</td><td>[0.19, 0.47]</td><td>0.09</td><td>[-0.19, 0.35]</td><td>0.73</td><td>0.13</td><td>[0.08, 0.19]</td></tr><tr><td>Lb*</td><td>0.20</td><td>[0.09, 0.32]</td><td>0.27</td><td>[-0.01, 0.53]</td><td>0.96</td><td>0.21</td><td>[0.15, 0.27]</td></tr></tbody></table>

IMF: intramuscular fat content, BW: body weight, CCW: commercial carcass weight, RCW: reference carcass weight, LW: liver weight, SF: scapular fat, PF: perirenal fat, DF: dissectible carcass fat percentage, CL*: carcass lightness, Ca*: carcass redness, Cb*: carcass yellowness, LL*: muscle lightness, La*: muscle redness, Lb*: muscle yellowness

## Weight traits

The bivariate analyses allowed us to identify three genomic regions associated with all weight traits (BW, CCW, RCW and LW); two located on Oryctolagus cuniculus chromosome (OCC) 1 and one on OCC19 (Table 3). In addition, genomic regions found on OCC 7, 14, 17 and 20 were associated with BW, CCW, and RCW, as expected due to the high correlation between these traits. From the aforementioned regions, those located at 35.42–38.06 Mb on OCC1 and at 9.86–12.05 Mb on OCC17 were not detected in the corresponding univariate analyses, likely due to the low percentage of genetic variance that they explained. On the other hand, the associated genomic regions explaining the highest percentage of genetic variance matched between the univariate and bivariate model analyses. This was the case for the regions located at 28.90–35.23 Mb on OCC1 (from 1.56 to 2.79% of genetic variance explained considering both univariate and bivariate results), and 30.10–32.82 Mb on OCC19 (from 0.93 to 2.06%), associated with the four traits; the region located at 19.04–20.59 Mb on OCC20 (from 2.08 to 3.21%) associated with BW, CCW and RCW traits; and the region located at 30.38–33.74

Mb on OCC12 (2.62% in the univariate GWAS and 3.69% in the bivariate GWAS) associated with LW (Table 3).

## Fat depot traits

The genomic regions associated with fat depot traits presented the highest percentages of explained genetic variance of all traits analysed in the current study (Table 4). The regions located at 33.71–38.20 Mb on OCC1 and at 148.79–151.57 Mb on OCC3 were related to all fat depot traits (PF, SF and DF), and were also the regions explaining the highest percentage of genetic variance (from 0.80 to 4.35% for the region on OCC1 and from 3.15 to 10.90% for the region on OCC3, considering both univariate and bivariate results). The highest percentage of genetic variance explained on OCC3 corresponded to the SF (10.32% in the univariate GWAS and 10.90% in the bivariate GWAS), reaching up to 0.97 of PPA. Another relevant genomic region was that located at 3.94–6.10 Mb on OCC6, which accounted for up to 3.14% of genetic variance for DF.

The results of weight and fat depot traits showed two common regions associated with most traits: the one located at 33.71–38.20 Mb on OCC1 (0.80–4.35%) and at 30.53–33.06 Mb on OCC19 (0.85–2.32%), suggesting that

<footnote>a HPD95%: the highest posterior density interval at 95%</footnote>

<footnote>b P0: probability of the correlation being greater than zero when positive or lower than zero when negative</footnote>

<footnote>c C² of common litter effect: proportion of phenotypic variance explained by the common litter effect</footnote>
<table><caption>**Table 3** Genomic regions associated with weight traits in rabbits</caption><thead><tr><th rowspan="2">OCC</th><th colspan="2">Position in genome</th><th rowspan="2">Trait</th><th rowspan="2">P_Gen_Var<br>UNIVAR</th><th rowspan="2">(PPA) UNIVAR</th><th rowspan="2">P_Gen_Var<br>BIVAR</th><th rowspan="2">(PPA) BIVAR</th></tr><tr><th>Start</th><th>End</th></tr></thead><tbody><tr><td rowspan="10">1</td><td rowspan="4">28901821</td><td rowspan="4">35233047</td><td>CCW</td><td>2.79</td><td>(0.3899)</td><td>2.40</td><td>(0.5195)</td></tr><tr><td>RCW</td><td>1.98</td><td>(0.3574)</td><td>1.93</td><td>(0.4582)</td></tr><tr><td>BW</td><td>2.40</td><td>(0.4551)</td><td>2.13</td><td>(0.3860)</td></tr><tr><td>LW</td><td>2.58</td><td>(0.4140)</td><td>1.56</td><td>(0.3739)</td></tr><tr><td rowspan="4">35421203</td><td rowspan="4">38061945</td><td>BW</td><td>-</td><td>-</td><td>2.56</td><td>(0.4585)</td></tr><tr><td>CCW</td><td>-</td><td>-</td><td>2.50</td><td>(0.6078)</td></tr><tr><td>LW</td><td>-</td><td>-</td><td>3.36</td><td>(0.4716)</td></tr><tr><td>RCW</td><td>-</td><td>-</td><td>2.29</td><td>(0.5413)</td></tr><tr><td rowspan="3">42411932</td><td rowspan="3">44790700</td><td>RCW</td><td>3.42</td><td>(0.4189)</td><td>2.18</td><td>(0.3882)</td></tr><tr><td>CCW</td><td>2.28</td><td>(0.3165)</td><td>1.50</td><td>(0.3334)</td></tr><tr><td>BW</td><td>0.83</td><td>(0.2198)</td><td>-</td><td>-</td></tr><tr><td>2</td><td>139299216</td><td>140642584</td><td>LW</td><td>0.85</td><td>(0.3537)</td><td>-</td><td>-</td></tr><tr><td>5</td><td>9779841</td><td>11429918</td><td>LW</td><td>1.05</td><td>(0.3703)</td><td>0.83</td><td>(0.3405)</td></tr><tr><td>6</td><td>2242936</td><td>3567876</td><td>LW</td><td>0.88</td><td>(0.3266)</td><td>-</td><td>-</td></tr><tr><td rowspan="6">7</td><td rowspan="2">110543262</td><td rowspan="2">112253533</td><td>CCW</td><td>0.96</td><td>(0.2938)</td><td>0.96</td><td>(0.3742)</td></tr><tr><td>RCW</td><td>1.22</td><td>(0.3369)</td><td>1.15</td><td>(0.3667)</td></tr><tr><td rowspan="2">156958947</td><td rowspan="2">158112218</td><td>BW</td><td>-</td><td>-</td><td>1.02</td><td>(0.2871)</td></tr><tr><td>CCW</td><td>-</td><td>-</td><td>0.75</td><td>(0.2780)</td></tr><tr><td rowspan="3">158765428</td><td rowspan="3">160790627</td><td>BW</td><td>0.93</td><td>(0.2659)</td><td>1.53</td><td>(0.3596)</td></tr><tr><td>CCW</td><td>0.88</td><td>(0.26830)</td><td>1.48</td><td>(0.4290)</td></tr><tr><td>RCW</td><td>0.81</td><td>(0.2694)</td><td>1.15</td><td>(0.3726)</td></tr><tr><td>9</td><td>74758537</td><td>77293167</td><td>LW</td><td>1.27</td><td>(0.3906)</td><td>0.96</td><td>(0.3367)</td></tr><tr><td>10</td><td>16545201</td><td>17708416</td><td>BW</td><td>-</td><td>-</td><td>0.83</td><td>(0.2814)</td></tr><tr><td>12</td><td>30386038</td><td>33742686</td><td>LW</td><td>3.69</td><td>(0.4474)</td><td>2.62</td><td>(0.4278)</td></tr><tr><td rowspan="2">13</td><td rowspan="2">29328268</td><td rowspan="2">30685073</td><td>BW</td><td>0.96</td><td>(0.2676)</td><td>-</td><td>-</td></tr><tr><td>RCW</td><td>0.96</td><td>(0.2766)</td><td>-</td><td>-</td></tr><tr><td rowspan="4">14</td><td rowspan="2">70440101</td><td rowspan="2">72889757</td><td>BW</td><td>0.98</td><td>(0.2296)</td><td>2.29</td><td>(0.3851)</td></tr><tr><td>CCW</td><td>1.50</td><td>(0.3432)</td><td>1.91</td><td>(0.5745)</td></tr><tr><td rowspan="2">134303242</td><td rowspan="2">136330626</td><td>RCW</td><td>1.36</td><td>(0.3808)</td><td>1.83</td><td>(0.5166)</td></tr><tr><td>BW</td><td>0.80</td><td>(0.2834)</td><td>1.04</td><td>(0.3514)</td></tr><tr><td rowspan="2">15</td><td rowspan="2">10831534</td><td rowspan="2">12836884</td><td>LW</td><td>1.71</td><td>(0.4325)</td><td>1.29</td><td>(0.3986)</td></tr><tr><td>BW</td><td>1.12</td><td>(0.3081)</td><td>-</td><td>-</td></tr><tr><td rowspan="3">17</td><td rowspan="3">9861967</td><td rowspan="3">12054383</td><td>CCW</td><td>0.76</td><td>(0.2815)</td><td>3.42</td><td>(0.6866)</td></tr><tr><td>RCW</td><td>-</td><td>-</td><td>2.37</td><td>(0.5925)</td></tr><tr><td>BW</td><td>-</td><td>-</td><td>1.29</td><td>(0.4535)</td></tr><tr><td rowspan="5">19</td><td rowspan="2">28333037</td><td rowspan="2">30004320</td><td>BW</td><td>-</td><td>-</td><td>1.04</td><td>(0.3294)</td></tr><tr><td>LW</td><td>1.26</td><td>(0.3997)</td><td>1.51</td><td>(0.4418)</td></tr><tr><td rowspan="3">30099870</td><td rowspan="3">32818846</td><td>BW</td><td>1.52</td><td>(0.3747)</td><td>2.06</td><td>(0.4210)</td></tr><tr><td>CCW</td><td>1.12</td><td>(0.3134)</td><td>1.89</td><td>(0.5585)</td></tr><tr><td>RCW</td><td>0.93</td><td>(0.3010)</td><td>1.44</td><td>(0.4459)</td></tr><tr><td rowspan="3">20</td><td rowspan="3">19038679</td><td rowspan="3">20588204</td><td>BW</td><td>3.21</td><td>(0.4975)</td><td>2.55</td><td>(0.4620)</td></tr><tr><td>CCW</td><td>2.20</td><td>(0.4175)</td><td>2.27</td><td>(0.6231)</td></tr><tr><td>RCW</td><td>2.27</td><td>(0.4465)</td><td>2.08</td><td>(0.5655)</td></tr></tbody></table>

OCC: Oryctolagus cuniculus chromosome. P_Gen_Var: percentage of genetic variance accounted for the associated genomic region; (PPA) = posterior probability of association or inclusion in the model; UNIVAR: univariate model analysis; BIVAR: univariate model analysis between this trait and intramuscular fat content. "-" = non associated genomic region for this model as the genomic region does not exceed the genomic window threshold (0.75). BW: body weight; CCW: commercial carcass weight; RCW: reference carcass weight; LW: liver weight
<table><caption>**Table 4** Genomic regions associated with fat depot traits in rabbits</caption><thead><tr><th rowspan="2">OCC</th><th colspan="2">Position in genome</th><th rowspan="2">Trait</th><th rowspan="2">P_Gen_Var<br>UNIVAR</th><th rowspan="2">(PPA) UNIVAR</th><th rowspan="2">P_Gen_Var<br>BIVAR</th><th rowspan="2">(PPA) BIVAR</th></tr><tr><th>Start</th><th>End</th></tr></thead><tbody><tr><td rowspan="7">1</td><td>27570818</td><td>28755852</td><td>DF</td><td>0.83</td><td>(0.44)</td><td>-</td><td>-</td></tr><tr><td>30428887</td><td>32483725</td><td>PF</td><td>1.65</td><td>(0.6388)</td><td>-</td><td>-</td></tr><tr><td>33710991</td><td>38203204</td><td>PF</td><td>1.24</td><td>(0.4977)</td><td>4.35</td><td>(0.7668)</td></tr><tr><td></td><td></td><td>DF</td><td>0.80</td><td>(0.4684)</td><td>3.14</td><td>(0.7470)</td></tr><tr><td></td><td></td><td>SF</td><td>-</td><td>-</td><td>2.31</td><td>(0.8742)</td></tr><tr><td>174113797</td><td>175711375</td><td>DF</td><td>-</td><td>-</td><td>0.85</td><td>(0.7048)</td></tr><tr><td>PF</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.78</td><td>(0.6847)</td></tr><tr><td rowspan="4">3</td><td>138874133</td><td>140152676</td><td>PF</td><td>0.77</td><td>(0.3838)</td><td>1.00</td><td>(0.4973)</td></tr><tr><td>148788842</td><td>151571092</td><td>SF</td><td>10.32</td><td>(0.9347)</td><td>10.90</td><td>(0.9687)</td></tr><tr><td></td><td></td><td>DF</td><td>4.54</td><td>(0.9361)</td><td>5.13</td><td>(0.9589)</td></tr><tr><td></td><td></td><td>PF</td><td>3.15</td><td>(0.8676)</td><td>3.50</td><td>(0.9247)</td></tr><tr><td>4</td><td>71082447</td><td>72639613</td><td>SF</td><td>-</td><td>-</td><td>0.77</td><td>(0.5960)</td></tr><tr><td rowspan="3">5</td><td>10154337</td><td>11665709</td><td>DF</td><td>0.91</td><td>(0.7625)</td><td>-</td><td>-</td></tr><tr><td>29840263</td><td>31819891</td><td>PF</td><td>-</td><td>-</td><td>1.41</td><td>(0.8708)</td></tr><tr><td></td><td></td><td>DF</td><td>-</td><td>-</td><td>1.06</td><td>(0.8357)</td></tr><tr><td rowspan="3">6</td><td>3936191</td><td>6108224</td><td>DF</td><td>3.84</td><td>(0.9663)</td><td>3.14</td><td>(0.9573)</td></tr><tr><td></td><td></td><td>PF</td><td>1.79</td><td>(0.7989)</td><td>1.15</td><td>(0.6995)</td></tr><tr><td></td><td></td><td>SF</td><td>1.17</td><td>(0.6058)</td><td>-</td><td>-</td></tr><tr><td rowspan="2">9</td><td>15473709</td><td>17491458</td><td>SF</td><td>1.24</td><td>(0.7700)</td><td>1.63</td><td>(0.7409)</td></tr><tr><td></td><td></td><td>DF</td><td>1.13</td><td>(0.7814)</td><td>0.89</td><td>(0.6678)</td></tr><tr><td>12</td><td>42683826</td><td>44688701</td><td>PF</td><td>1.74</td><td>(0.6513)</td><td>-</td><td>-</td></tr><tr><td rowspan="3">13</td><td>100353890</td><td>101950026</td><td>DF</td><td>-</td><td>-</td><td>0.86</td><td>(0.6097)</td></tr><tr><td>139448688</td><td>141261644</td><td>DF</td><td>1.46</td><td>(0.7405)</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td>SF</td><td>1.70</td><td>(0.7965)</td><td>-</td><td>-</td></tr><tr><td rowspan="2">14</td><td>42451260</td><td>43654926</td><td>SF</td><td>-</td><td>-</td><td>0.80</td><td>(0.4723)</td></tr><tr><td>134489752</td><td>135928030</td><td>SF</td><td>-</td><td>-</td><td>0.86</td><td>(0.4756)</td></tr><tr><td>17</td><td>10061927</td><td>11798813</td><td>PF</td><td>-</td><td>-</td><td>0.87</td><td>(0.7564)</td></tr><tr><td rowspan="2">19</td><td>30530951</td><td>33065887</td><td>PF</td><td>1.41</td><td>(0.7623)</td><td>2.32</td><td>(0.8118)</td></tr><tr><td></td><td></td><td>DF</td><td>0.85</td><td>(0.6803)</td><td>1.40</td><td>(0.7668)</td></tr><tr><td></td><td>54592191</td><td>56054881</td><td>SF</td><td>-</td><td>-</td><td>0.86</td><td>(0.4305)</td></tr></tbody></table>

OCC: Oryctolagus cuniculus chromosome. P_Gen_Var: percentage of genetic variance accounted for the associated genomic region; (PPA) = posterior probability of association or inclusion in the model; UNIVAR: univariate model analysis; BIVAR: bivariate model analysis between this trait and intramuscular fat content. "-" = non associated genomic region for this model as the genomic region does not exceed the genomic window threshold (0.75). SF: scapular fat; PF: perirenal fat; DF: dissectible carcass fat percentage

these genomic regions have relevant effects for both sets of traits.

## Meat quality traits

The genomic regions associated with meat quality and carcass colour traits are shown in Table 5. These associated genomic regions were scattered across the whole rabbit genome, with the only exception of OCC21. Several genomic regions were clearly associated when the bivariate GWAS analysis was performed, most likely due to their relationship and association with IMF, whereas no associations were found in the univariate GWAS. Among these regions, those located at 35.42–37.88 Mb

on OCC1, 26.74–28.90 Mb on OCC7 and 67.28–70.26 Mb on OCC9 stand out because they were associated with most of the studied traits. The region on OCC1 presented association with pH (1.65% of genetic variance explained), La* (1.57%), Lb* (0.76%), Ca* (1.46%), Cb* (1.17%), and CL* (1.04%), and it overlapped with the associated region found for weight and fat depot traits. On the other hand, the region on OCC7 was associated with pH (8.48%), Ca* (1.13%), Cb* (0.79%) and La* (0.80%), while that in OCC9 was associated with Lb* (2.57%), Ca* (1.81%), La* (1.84%), pH (1.44%) and Cb* (1.02%). The regions located at 120.34–122.10 Mb on OCC1 (0.90–0.95%), 70.69–72.82 Mb on OCC4
<table><caption>**Table 5** Genomic regions associated with meat quality traits in rabbits</caption><thead><tr><th rowspan="2">OCC</th><th colspan="2">Position in genome</th><th rowspan="2">Trait</th><th rowspan="2">P_Gen_Var<br>UNIVAR</th><th rowspan="2">(PPA) UNIVAR</th><th rowspan="2">P_Gen_Var<br>BIVAR</th><th rowspan="2">(PPA) BIVAR</th></tr><tr><th>Start</th><th>End</th></tr></thead><tbody><tr><td rowspan="10">1</td><td rowspan="10">35421203</td><td rowspan="10">37880078</td><td>La*</td><td>-</td><td>-</td><td>1.57</td><td>(0.3732)</td></tr><tr><td>pH</td><td>-</td><td>-</td><td>1.65</td><td>(0.5932)</td></tr><tr><td>Ca*</td><td>-</td><td>-</td><td>1.46</td><td>(0.4990)</td></tr><tr><td>Cb*</td><td>-</td><td>-</td><td>1.17</td><td>(0.3757)</td></tr><tr><td>CL*</td><td>-</td><td>-</td><td>1.04</td><td>(0.3486)</td></tr><tr><td>Lb*</td><td>-</td><td>-</td><td>0.76</td><td>(0.3041)</td></tr><tr><td>La*</td><td>0.84</td><td>(0.2501)</td><td>0.78</td><td>(0.2734)</td></tr><tr><td>La*</td><td>0.75</td><td>(0.2193)</td><td>-</td><td>-</td></tr><tr><td>pH</td><td>-</td><td>-</td><td>0.95</td><td>(0.5754)</td></tr><tr><td>Cb*</td><td>-</td><td>-</td><td>0.90</td><td>(0.3920)</td></tr><tr><td rowspan="2">2</td><td>19814073</td><td>20883718</td><td>LL*</td><td>0.90</td><td>(0.2862)</td><td>-</td><td>-</td></tr><tr><td>152094935</td><td>153292895</td><td>Ca*</td><td>0.84</td><td>(0.2522)</td><td>0.76</td><td>(0.2584)</td></tr><tr><td rowspan="10">3</td><td>34588926</td><td>36743147</td><td>Lb*</td><td>1.97</td><td>(0.2419)</td><td>1.67</td><td>(0.2756)</td></tr><tr><td>38625731</td><td>41531359</td><td>Lb*</td><td>1.57</td><td>(0.241)</td><td>1.27</td><td>(0.2699)</td></tr><tr><td>123774002</td><td>125777015</td><td>Lb*</td><td>-</td><td>-</td><td>0.95</td><td>(0.2456)</td></tr><tr><td>133620124</td><td>135530513</td><td>Lb*</td><td>-</td><td>-</td><td>0.99</td><td>(0.2108)</td></tr><tr><td>138874133</td><td>140306258</td><td>CL*</td><td>-</td><td>-</td><td>1.46</td><td>(0.2849)</td></tr><tr><td>Lb*</td><td>1.09</td><td>(0.2461)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>LL*</td><td>0.87</td><td>(0.2621)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ca*</td><td>0.82</td><td>(0.2675)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ca*</td><td>-</td><td>-</td><td>2.98</td><td>(0.5346)</td><td>-</td><td>-</td></tr><tr><td rowspan="2">4</td><td>62983442</td><td>64506834</td><td>pH</td><td>-</td><td>-</td><td>0.93</td><td>(0.4124)</td></tr><tr><td>66806370</td><td>68324187</td><td>Lb*</td><td>-</td><td>-</td><td>0.83</td><td>(0.3446)</td></tr><tr><td rowspan="10">5</td><td>70693711</td><td>72818170</td><td>LL*</td><td>0.89</td><td>(0.2443)</td><td>-</td><td>-</td></tr><tr><td>74338850</td><td>75611556</td><td>Lb*</td><td>-</td><td>-</td><td>0.80</td><td>(0.2807)</td></tr><tr><td>7298347</td><td>8860898</td><td>Ca*</td><td>2.97</td><td>(0.3919)</td><td>2.24</td><td>(0.4636)</td></tr><tr><td>15837999</td><td>19050703</td><td>CL*</td><td>1.33</td><td>(0.2612)</td><td>0.87</td><td>(0.233)</td></tr><tr><td>20180124</td><td>22111503</td><td>Lb*</td><td>1.71</td><td>(0.2535)</td><td>1.87</td><td>(0.3275)</td></tr><tr><td>5452</td><td>1979376</td><td>Lb*</td><td>1.09</td><td>(0.2181)</td><td>1.18</td><td>(0.2521)</td></tr><tr><td>2518414</td><td>4805150</td><td>Lb*</td><td>0.76</td><td>(0.2166)</td><td>-</td><td>-</td></tr><tr><td>8256520</td><td>9284111</td><td>pH</td><td>-</td><td>-</td><td>8.48</td><td>(0.6734)</td></tr><tr><td>26743882</td><td>28900797</td><td>Ca*</td><td>-</td><td>-</td><td>1.13</td><td>(0.482)</td></tr><tr><td rowspan="10">7</td><td rowspan="10">45455815</td><td rowspan="10">47590522</td><td>Cb*</td><td>-</td><td>-</td><td>0.79</td><td>(0.3531)</td></tr><tr><td>La*</td><td>-</td><td>-</td><td>0.80</td><td>(0.3402)</td></tr><tr><td>La*</td><td>1.65</td><td>(0.3593)</td><td>1.23</td><td>(0.3227)</td></tr><tr><td>Ca*</td><td>1.11</td><td>(0.354)</td><td>-</td><td>-</td></tr><tr><td>CL*</td><td>0.84</td><td>(0.2273)</td><td>-</td><td>-</td></tr><tr><td>Ca*</td><td>-</td><td>-</td><td>0.75</td><td>(0.2735)</td></tr><tr><td>69048438</td><td>70596609</td><td>La*</td><td>0.86</td><td>(0.2602)</td><td>-</td><td>-</td></tr><tr><td>72316729</td><td>73316729</td><td>Cb*</td><td>0.81</td><td>(0.1563)</td><td>-</td><td>-</td></tr><tr><td>96996061</td><td>98552759</td><td>CL*</td><td>1.12</td><td>(0.3036)</td><td>0.97</td><td>(0.2968)</td></tr><tr><td>167220424</td><td>168333021</td><td>Ca*</td><td>1.04</td><td>(0.2376)</td><td>1.15</td><td>(0.2807)</td></tr><tr><td>1663933</td><td>3568211</td><td>LL*</td><td>0.85</td><td>(0.2433)</td><td>-</td><td>-</td></tr><tr><td rowspan="2">8</td><td>68365260</td><td>70147889</td><td>La*</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>99133029</td><td>100490608</td><td>LL*</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<table><caption>**Table 5** (continued)</caption><thead><tr><th rowspan="2">OCC</th><th colspan="2">Position in genome</th><th rowspan="2">Trait</th><th rowspan="2">P_Gen_Var<br>UNIVAR</th><th rowspan="2">(PPA) UNIVAR</th><th rowspan="2">P_Gen_Var<br>BIVAR</th><th rowspan="2">(PPA) BIVAR</th></tr><tr><th>Start</th><th>End</th></tr></thead><tbody><tr><td rowspan="12">9</td><td rowspan="12">15140254<br>67289669</td><td rowspan="12">16993305<br>70258747</td><td>LL*</td><td>1.02</td><td>(0.2416)</td><td>–</td><td>–</td></tr><tr><td>Lb*</td><td>–</td><td>–</td><td>2.57</td><td>(0.3983)</td></tr><tr><td>Ca*</td><td>–</td><td>–</td><td>1.81</td><td>(0.3949)</td></tr><tr><td>La*</td><td>–</td><td>–</td><td>1.84</td><td>(0.3640)</td></tr><tr><td>pH</td><td>–</td><td>–</td><td>1.44</td><td>(0.3659)</td></tr><tr><td>Cb*</td><td>–</td><td>–</td><td>1.02</td><td>(0.3564)</td></tr><tr><td>Ca*</td><td>–</td><td>–</td><td>0.86</td><td>(0.2473)</td></tr><tr><td>CL*</td><td>2.06</td><td>(0.2898)</td><td>1.45</td><td>(0.2803)</td></tr><tr><td>CL*</td><td>0.93</td><td>(0.2435)</td><td>–</td><td>–</td></tr><tr><td>La*</td><td>2.12</td><td>(0.3259)</td><td>2.53</td><td>(0.3529)</td></tr><tr><td>LL*</td><td>2.75</td><td>(0.3939)</td><td>–</td><td>–</td></tr><tr><td>CL*</td><td>0.94</td><td>(0.2373)</td><td>–</td><td>–</td></tr><tr><td>LL*</td><td>1.30</td><td>(0.3082)</td><td>–</td><td>–</td></tr><tr><td>CL*</td><td>–</td><td>–</td><td>1.05</td><td>(0.2634)</td></tr><tr><td rowspan="6">10<br>11<br>12</td><td rowspan="6">17680826<br>83081198<br>47446494<br>52181759<br>111962950</td><td rowspan="6">91658263<br>19827740<br>85377090<br>48928659<br>54177785<br>113947302</td><td>LL*</td><td>1.43</td><td>(0.2941)</td><td>–</td><td>–</td></tr><tr><td>LL*</td><td>2.92</td><td>(0.3614)</td><td>1.02</td><td>(0.547)</td></tr><tr><td>Cb*</td><td>–</td><td>–</td><td>0.88</td><td>(0.3942)</td></tr><tr><td>Ca*</td><td>1.44</td><td>(0.3146)</td><td>1.41</td><td>(0.3205)</td></tr><tr><td>La*</td><td>0.79</td><td>(0.2233)</td><td>–</td><td>–</td></tr><tr><td>CL*</td><td>–</td><td>–</td><td>0.77</td><td>(0.2507)</td></tr><tr><td rowspan="6">13<br>14<br>15</td><td rowspan="6">66873358<br>72861641<br>93880669<br>122014443<br>128853115</td><td rowspan="6">68737843<br>75624854<br>95590867<br>123999114<br>130283007</td><td>LL*</td><td>1.05</td><td>(0.2714)</td><td>–</td><td>–</td></tr><tr><td>LL*</td><td>0.77</td><td>(0.2669)</td><td>–</td><td>–</td></tr><tr><td>Cb*</td><td>–</td><td>–</td><td>1.21</td><td>(0.3797)</td></tr><tr><td>Ca*</td><td>0.76</td><td>(0.2232)</td><td>–</td><td>–</td></tr><tr><td>La*</td><td>1.94</td><td>(0.3850)</td><td>1.73</td><td>(0.4783)</td></tr><tr><td>CL*</td><td>–</td><td>(0.2944)</td><td>1.53</td><td>(0.3182)</td></tr><tr><td rowspan="6">16<br>17</td><td rowspan="6">42146089<br>4496821<br>6726456<br>19416243<br>22789228</td><td rowspan="6">44087078<br>5673569<br>9405305<br>20467908<br>24896386</td><td>LL*</td><td>1.20</td><td>(0.2944)</td><td>1.95</td><td>(0.3260)</td></tr><tr><td>LL*</td><td>–</td><td>–</td><td>1.02</td><td>(0.2994)</td></tr><tr><td>Cb*</td><td>–</td><td>–</td><td>0.99</td><td>(0.2948)</td></tr><tr><td>Ca*</td><td>–</td><td>–</td><td>1.87</td><td>(0.3216)</td></tr><tr><td>La*</td><td>–</td><td>–</td><td>0.90</td><td>(0.2394)</td></tr><tr><td>CL*</td><td>–</td><td>–</td><td>1.09</td><td>(0.2843)</td></tr><tr><td rowspan="6">18<br>19</td><td rowspan="6">34695954<br>60905831<br>62484401</td><td rowspan="6">36361523<br>62116156<br>63886847</td><td>LL*</td><td>1.08</td><td>(0.3294)</td><td>0.80</td><td>(0.287)</td></tr><tr><td>LL*</td><td>0.76</td><td>(0.1676)</td><td>–</td><td>–</td></tr><tr><td>Cb*</td><td>0.84</td><td>(0.1377)</td><td>–</td><td>–</td></tr><tr><td>CL*</td><td>–</td><td>–</td><td>1.18</td><td>(0.2739)</td></tr><tr><td>CL*</td><td>1.44</td><td>(0.2859)</td><td>–</td><td>–</td></tr><tr><td>CL*</td><td>0.83</td><td>(0.2019)</td><td>–</td><td>–</td></tr><tr><td>20</td><td>22699320</td><td>26339396</td><td>LL*</td><td>0.93</td><td>(0.2283)</td><td>0.78</td><td>(0.1708)</td></tr></tbody></table>

OCC: Oryctolagus cuniculus chromosome. P_Gen_Var: percentage of genetic variance accounted for the associated genomic region; (PPA) = posterior probability of association or inclusion in the model; UNIVAR: univariate model analysis; BIVAR: bivariate model analysis between this trait and intramuscular fat content. "–": non associated genomic region for this model as the genomic region does not exceed the genomic window threshold (0.75). CL*: carcass lightness; Ca*: carcass redness; Cb*: carcass yellowness; LL*: muscle lightness; La*: muscle redness, Lb*: muscle yellowness
<table><caption>**Table 6** Genomic correlations between regions associated with intramuscular fat (IMF) and carcass-weight, fat depot and meat quality traits in rabbits</caption><thead><tr><th>OCC</th><th colspan="6">Genomic regions associated only with IMF</th></tr><tr><th></th><th>IMF-assGR (Mb)</th><th>IMF-P_Gen_Var</th><th>IMF-(PPA)</th><th>Traits (bivariate analysis)</th><th>Highest GenCor</th><th>Correlation SD</th></tr></thead><tbody><tr><td rowspan="8">1</td><td rowspan="8">120.28-122.19</td><td>1.16</td><td>(0.7555)</td><td>IMF & Ca*</td><td>0.01</td><td>0.5200</td></tr><tr><td>1.20</td><td>(0.7717)</td><td>IMF & Lb*</td><td>0.10</td><td>0.4844</td></tr><tr><td>1.05</td><td>(0.7256)</td><td>IMF & La*</td><td>0.13</td><td>0.5319</td></tr><tr><td>0.92</td><td>(0.6839)</td><td>IMF & BW</td><td>0.15</td><td>0.5285</td></tr><tr><td>0.90</td><td>(0.6691)</td><td>IMF & RCW</td><td>0.06</td><td>0.5062</td></tr><tr><td>0.81</td><td>(0.637)</td><td>IMF & CCW</td><td>0.07</td><td>0.4783</td></tr><tr><td>0.76</td><td>(0.7145)</td><td>IMF & DF</td><td>0.26</td><td>0.4627</td></tr><tr><td>0.76</td><td>(0.808)</td><td>IMF & SF</td><td>0.29</td><td>0.4894</td></tr><tr><td>142.26-143.41</td><td>0.79</td><td>(0.7803)</td><td>IMF & DF</td><td>-0.003</td><td>0.3877</td></tr><tr><td rowspan="2">9</td><td rowspan="2">68.31-69.37</td><td>0.79</td><td>(0.5768)</td><td>IMF & CCW</td><td>0.10</td><td>0.4678</td></tr><tr><td>0.77</td><td>(0.5622)</td><td>IMF & RCW</td><td>0.07</td><td>0.4569</td></tr><tr><td rowspan="4">13</td><td rowspan="4">93.68-95.61</td><td>1.08</td><td>(0.8501)</td><td>IMF & PF</td><td>0.08</td><td>0.4036</td></tr><tr><td>0.80</td><td>(0.6536)</td><td>IMF & RCW</td><td>0.19</td><td>0.5206</td></tr><tr><td>0.76</td><td>(0.6408)</td><td>IMF & CCW</td><td>0.18</td><td>0.5346</td></tr><tr><td>0.78</td><td>(0.7200)</td><td>IMF & Lb*</td><td>0.10</td><td>0.4371</td></tr><tr><td>16</td><td>83.08-84.47</td><td>1.03</td><td>(0.5627)</td><td>IMF & LW</td><td>-0.04</td><td>0.4259</td></tr><tr><td>17</td><td>4.54-5.80</td><td>0.98</td><td>(0.5519)</td><td>IMF & LW</td><td>-0.04</td><td>0.4129</td></tr></tbody></table>

<table><caption>**Genomic regions associated with both IMF and specific trait**</caption><thead><tr><th>OCC</th><th>IMF-assGR (Mb)</th><th>IMF-P_Gen_Var</th><th>IMF-(PPA)</th><th>Trait</th><th>Position in the genome (Mb)</th><th>Highest GenCor</th><th>Correlation SD</th></tr></thead><tbody><tr><td rowspan="17">1</td><td>35.42-37.70</td><td>2.00</td><td>(0.6844)</td><td>BW</td><td>35.42-38.06</td><td>0.29</td><td>0.5699</td></tr><tr><td>35.42-37.63</td><td>1.78</td><td>(0.6484)</td><td>CCW</td><td></td><td>0.31</td><td>0.5561</td></tr><tr><td>35.45-37.54</td><td>1.59</td><td>(0.6190)</td><td>LW</td><td></td><td>0.40</td><td>0.5359</td></tr><tr><td>35.45-37.67</td><td>1.77</td><td>(0.6530)</td><td>RCW</td><td></td><td>0.29</td><td>0.5596</td></tr><tr><td>35.48-37.71</td><td>1.82</td><td>(0.7911)</td><td>PF</td><td>33.71-38.20</td><td>0.60</td><td>0.5034</td></tr><tr><td>35.54-37.90</td><td>1.65</td><td>(0.7675)</td><td>DF</td><td></td><td>0.58</td><td>0.5049</td></tr><tr><td>35.42-37.62</td><td>2.00</td><td>(0.8923)</td><td>SF</td><td></td><td>0.40</td><td>0.5406</td></tr><tr><td>35.42-37.80</td><td>2.38</td><td>(0.7309)</td><td>La*</td><td>35.42-37.88</td><td>-0.17</td><td>0.5646</td></tr><tr><td>35.42-37.78</td><td>2.23</td><td>(0.7389)</td><td>pH</td><td></td><td>0.10</td><td>0.5825</td></tr><tr><td>35.42-37.80</td><td>2.44</td><td>(0.7500)</td><td>Ca*</td><td></td><td>-0.17</td><td>0.5622</td></tr><tr><td>35.42-37.80</td><td>2.19</td><td>(0.7241)</td><td>Cb*</td><td></td><td>0.02</td><td>0.4827</td></tr><tr><td>35.54-37.71</td><td>1.75</td><td>(0.6837)</td><td>CL*</td><td></td><td>-0.09</td><td>0.5508</td></tr><tr><td>35.42-37.78</td><td>2.66</td><td>(0.7585)</td><td>Lb*</td><td></td><td>0.02</td><td>0.3730</td></tr><tr><td>120.30-122.18</td><td>1.10</td><td>(0.7554)</td><td>pH</td><td>120.34-122.10</td><td>-0.08</td><td>0.5727</td></tr><tr><td>120.30-122.16</td><td>1.08</td><td>(0.7296)</td><td>Cb*</td><td></td><td>0.19</td><td>0.5805</td></tr><tr><td>3</td><td>148.94-150.89</td><td>1.21</td><td>(0.9291)</td><td>SF</td><td>148.79-151.57</td><td>0.92</td><td>0.2746</td></tr><tr><td>7</td><td>27.07-28.44</td><td>0.75</td><td>(0.7599)</td><td>pH</td><td>26.74-28.90</td><td>0.63</td><td>0.5144</td></tr><tr><td rowspan="3">9</td><td>68.31-69.36</td><td>0.77</td><td>(0.6397)</td><td>Lb*</td><td>67.29-70.26</td><td>0.34</td><td>0.5044</td></tr><tr><td>67.93-69.37</td><td>0.81</td><td>(0.6128)</td><td>pH</td><td></td><td>-0.15</td><td>0.4935</td></tr><tr><td>67.93-69.37</td><td>0.79</td><td>(0.5951)</td><td>Cb*</td><td></td><td>0.17</td><td>0.4952</td></tr><tr><td>13</td><td>93.88-95.57</td><td>0.84</td><td>(0.7000)</td><td>Cb*</td><td>93.88-95.59</td><td>0.21</td><td>0.5546</td></tr></tbody></table>

OCC: Oryctolagus cuniculus chromosome. P_Gen_Var: percentage of genetic variance accounted for the associated genomic region; (PPA) = posterior probability of association or inclusion in the model; UNIVAR: univariate model analysis; BIVAR: bivariate model analysis between this trait and intramuscular fat content; Highest GenCor: Highest genetic correlation estimate amongst all associated consecutive genomic windows that belong to the genomic region. Correlation SD: correlation standard deviation of the highest genetic correlation estimate BW: body weight; CCW: commercial carcass weight; RCW: reference carcass weight; LW: liver weight; SF: scapular fat; PF: perirenal fat; DF: dissectible carcass fat percentage; CL*: carcass lightness, Ca*: carcass redness, Cb*: carcass yellowness; LL*: muscle lightness, La*: muscle redness, Lb*: muscle yellowness
(0.83–2.98%), 82.76–84.47 Mb on OCC16 (1.02–1.95%), and 45.36–59.70 Mb on OCC17 (0.99–1.87%) were other associated genomic regions identified only with bivariate model analysis.

## Pleiotropic genomics regions associated with intramuscular fat

Table 6 shows the genomic regions found in all bivariate analyses that were associated only with IMF and those associated with both traits (i.e., pleiotropic regions).

Associated genomic regions with pleiotropic effects were found on OCC1, OCC3, OCC7, OCC9 and OCC13. The most important region was the one located at 35.42–37.90 Mb on OCC1, which explained up to 2.66% of the genetic variance of IMF and, as mentioned before, it was also associated with all analysed traits with the only exception of LL*. Other regions to highlight were those located at 148.94–150.89 Mb on OCC3 and 27.07–28.44 Mb on OCC7. The first one was the most relevant region for SF and explained 1.21% of the genetic variance of IMF, while the second one was remarkably important for pH and explained 0.75% of the genetic variance of IMF. These pleiotropic regions were supported by the high genomic correlations between the traits for each specific region (Table 6). However, it is worth mentioning that the magnitude of the genomic correlation of the associated regions depends on the trait analysed and is not necessarily related to the association. For example, in the bivariate analysis for the colour trait Cb*, the region 120–122 Mb on OCC1 presented association with both traits and had a genomic correlation of 0.19, whereas in the bivariate analysis for DF this same region had a higher genomic correlation (0.26) but was only associated with IMF and not with DF.

## Gene search and identification of functional mutations

The associated genomic regions found in all analyses comprised 1687 genes: 500 genes for weight traits, 385 genes for fat depot traits, and 1005 genes for meat quality traits (see Additional file 1: Table S1).

The WGS performed on the sires allowed to identify a large number of allele variants present in the associated genomic regions for each set of traits (see Additional file 2: Table S2). From those, we kept (1) the allele variants that were fixed in one line and segregating at an allele frequency lower than 0.5 in the other, and (2) the variants whose alleles were segregating at opposite frequencies (greater than 0.65) between lines. We finally detected 167 variants (located in 69 genes) in the genomic regions associated to weight traits (see Additional file 3: Table S3); 79 variants (located in 28 genes) in the genomic regions associated to fat traits (see Additional file 4: Table S4); and 134 variants (located in 46

genes) in the genomic regions associated to meat quality traits (see Additional file 5: Table S5).

When analysing the main pleiotropic regions highlighted in the previous section, there were 13 variants detected: 7 on the pleiotropic region of OCC1 (associated to almost all traits) and located in genes *SH3GL2*, *CNTLN*, and *BNC2*; and 6 on the pleiotropic region of OCC3 (associated to all fat traits and explaining a large percentage of genetic variance), all located in the *ENSOCUG00000000157* gene. Unfortunately, we did not detect any relevant variant on the pleiotropic region of OCC7 (explaining a large percentage of the genetic variance of pH) under the assumed parameters.

## Discussion

Understanding how selection for IMF affects the correlated traits is important for making accurate decisions in livestock breeding programmes. Causal variants can affect more than one trait, so the study of their pleiotropic contributions could be relevant in this sense. In fact, bivariate GWAS of two highly correlated traits can increase the power to detect associations on both studied traits [19, 20, 28]. The detection of variants and genomic regions associated with a particular trait also depends on the model and its prior parameters. For instance, novel genomic regions associated with IMF were discovered in these lines when a maternal genetic effect was used in the model [39]. On the other hand, the priors can influence the results of the various “Bayesian alphabet” GWAS methods, especially when the number of unknown parameters (p) exceeds the sample size (n), which is typical in genomic analyses [25, 26, 40]. However, with largely polygenic traits, the identification of genomic regions or variants using regression methods based on whole genome can be difficult even when $n \approx p$. A GWAS for body weight at 35 days performed in broilers using an extremely large sample (137,000 genotyped individuals) could only identify 25 genomic regions that together explained ~30% of the genetic variance [41]. In practice, most common Bayesian methods for GWAS, such as B, C, C-$\pi$ and R, are able to separate the “noise” from those genomic variants with relative relevance to a trait, thus improving the detection of variants and genomic regions associated [20]. Therefore, in this study we assumed that the genomic regions detected could be considered as associated but probably explain a slightly lower variance than those found. Furthermore, given that pleiotropy is being increasingly considered in genomic evaluations, bivariate GWAS would be a relevant genetic analysis to explore in livestock breeding programmes.

The region on OCC1 (~ from 35 to 38 Mb) was found to be associated to all traits analysed in this study (except LL*), demonstrating its importance for carcass,
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Manhattan plots of genome-wide association studies for intramuscular fat bivariate (inner part), scapular fat univariate (mid part) and scapular fat bivariate (outer part). Each blue dot stands for a genomic window exceeding a 0.75% genetic variance threshold, and each orange dot stands for a window exceeding a 1.50% genetic variance threshold</figcaption></figure>

fat depot and meat quality traits. This region has also been reported on these lines in a previous selection sig-nature study for IMF [17], and in GWA studies for fatty acid composition [39, 42]. However, this genomic region was not detected neither in the univariate GWAS for SF, weight, and meat quality traits (see for example Fig. 1 for SF and Fig. 2 for pH) nor in the univariate GWAS pre-viously performed for IMF [16]. This represents a good example of how the analysis of related traits under a bivariate model can aid in the detection of associated genomic regions for both traits. The functional annota-tion of the genes located in this region showed that they

were most likely involved in the glucose, growth and lipid metabolisms, whereas the variants found with the WGS pointed to the *SH3GL2*, *CNTLN*, and *BNC2* genes as candidate genes affecting weight, fat depots and meat quality traits. The *SH3GL2* gene is involved in synaptic vesicle endocytosis and lipid binding and modification [43] and, together with *CNTLN*, has been associated with an obesity index in pigs [44]. On the other hand, the expression levels of *BCN2* have been positively correlated with pH in pigs [45]. This region also harbours the *PLIN2* gene, which has a relevant function in lipids metabolism and has been related to IMF content, growth and lean
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Manhattan plots of genome-wide association studies for intramuscular fat bivariate (inner part), pH univariate (mid part) and pH bivariate (outer part). Each blue dot stands for a genomic window exceeding a 0.75% genetic variance threshold, and each orange dot stands for a window exceeding a 1.50% genetic variance threshold</figcaption></figure>

weight in pigs [46, 47]. In these lines, there was one vari-
ant detected on this gene, which had the reference allele
fixed in the L line, but its frequency in the H line was of
0.6 (i.e., frequency of 0.4 of the alternative allele). Even
though this variant did not exceed the threshold assumed
(i.e., fixed in one line and a frequency lower than 0.5 in
the other), given its relevant function and the evidence in
the literature, it still represents a promising candidate to
evaluate in future functional analyses.

The pleiotropic region associated to SF and IMF on
OCC3 is another interesting example of the advantage of
using a pleiotropic GWAS approach. Both univariate and

bivariate models were able to detect its association to SF,
most likely due to its large effect, explaining over 10% of
the genetic variance in both cases (Fig. 1). However, this
region was not detected in the univariate analysis for IMF
[16] but only in the bivariate GWAS, explaining 1.21%
of its genetic variance, which was consistent with the
high correlation found between these traits in this win-
dow (0.92). This region was also associated with the fatty
acid composition in these lines, especially C14:0, C16:0,
C18:0, SFA, C18:1n9, MUFA, C18:2n6, C20:4n6, PUFA
and PUFA/SFA, explaining between 0.68 and 6.05% of
their genetic variance [42]. A total of 6 allele variants
were detected in this region, all of which were located in the *ENSOCUG00000000157* novel gene, which is also known as *ST3GAL1* in other species like humans, poultry, pigs, and mice. This gene is involved in lipid metabolism and its expression has also been associated to IMF deposition in broilers [48]. Further analysis of this region could help determine its potential for marker-assisted selection aimed at reducing fat depot accumulation with minimal impact on intramuscular fat (IMF) content. However, weighting-SNP-based genomic selection methods would be a better option than using marker assisted selection in rabbits, especially for economically relevant traits, e.g. growth rate or litter size [49–51].

Under the assumption of pleiotropic effects, the genomic region on OCC7 (26.7–28.9 Mb) was associated with most of the meat quality traits and showed an extreme association value for pH (Fig. 2). This region was not detected in the corresponding univariate GWAS. Regardless of the large number of genes retrieved, no allele variants were detected in the WGS under the assumed parameters. Nonetheless, it still represents an interesting candidate region to evaluate on posterior functional analyses. Finally, another interesting pleiotropic region for meat quality traits was located at 120–121 Mb on OCC1, which was associated with pH and Cb* and it was only detected in the bivariate analyses. This region has consistently appeared in the previous univariate GWAS analysis for IMF performed on these lines (explaining up to 2.03% of the genetic variance) [16], as well as in the selection signature analysis [17] and fatty acid composition analyses [42]. In this region, 4 allele variants were detected, which were located on the novel gene *ENSOCUG00000036784* and on *MAML2*. No evidence was found regarding their involvement in meat quality traits like those discussed in this paper. However, the *MAML2* is a transcription factor that could be affecting several pathways, and has been related to other relevant production traits, like feed efficiency in pigs [52], thus representing an interesting candidate gene for further studies.

## Conclusions

Animal breeding programmes are based on multiple trait selection. In this study, bivariate models assuming pleiotropic effects allowed us to identify relevant genomic regions simultaneously associated with IMF and weight, fat depot and meat quality traits. Additionally, the whole genome sequencing performed allowed us to identify candidate genes with putative functional mutations affecting the mentioned traits. Altogether, our results can provide relevant insights into the correlations and relationships between the traits at the

genomic level, which could improve the genomic selection potential. The main pleiotropic region was found on OCC1 between 35.4 Mb and 38.2 Mb, which was associated to almost all analysed traits, with the only exception of the LL* colour parameter. The variants found with the WGS pointed to *PLIN2*, *SH3GL2*, *CNTLN*, and *BNC2* as the main candidate genes in this region affecting weight, fat depots and meat quality traits. Other relevant pleiotropic regions were found on OCC3 between 148.79 and 151.57 Mb, showing the highest association value with SF, and on OCC7 between 26.7 and 28.9 Mb, with the highest association with pH. The allele variants on OCC3 suggested the novel gene *ENSOCUG00000000157* as the main candidate gene affecting fat deposition, whereas on OCC7 there were no allele variants segregating differently between the lines in the identified regions.

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12711-025-00971-5.

Additional file 1: Table S1. Genes retrieved from the genomic regions associated with weight, fat and meat quality traits identified with the bivariate approach. Complete list of genes found in the associated genomic regions to weight, fat and meat quality traits identified in the bivariate genome-wide association analyses

Additional file 2: Table S2. Allele variants found in the genomic regions associated to weight, fat, and meat quality traits. Complete list of allele variants found in the associated genomic regions to weight, fat and meat quality traits identified in the bivariate genome-wide association analyses

Additional file 3: Table S3. Relevant allelic variants located in the genomic regions associated to weight traits. The allelic variants considered relevant were those that were fixed in one line and segregating at an allele frequency lower than 0.5 in the other, and whose alleles were segregating at opposite frequencies between lines

Additional file 4: Table S4. Relevant allelic variants located in the genomic regions associated to fat depot traits. The allelic variants considered relevant were those that were fixed in one line and segregating at an allele frequency lower than 0.5 in the other, and whose alleles were segregating at opposite frequencies between lines

Additional file 5: Table S5. Relevant allelic variants located in the genomic regions associated to meat quality traits. The allelic variants considered relevant were those that were fixed in one line and segregating at an allele frequency lower than 0.5 in the other, and whose alleles were segregating at opposite frequencies between lines

## Acknowledgements

The authors thank Federico Pardo and Marina Morini for their technical assistance.

## Author contributions

AB, NIE, and PH conceived and designed the experimental procedures and supervised the study. AB and PH secured substantial funding. BSSM and PH collected the samples and performed the phenotypic data recording. AZG, BSSM, and NIE organized, edited and performed the genomic analyses. NIE, AB, PH, AZG and BSSM helped with the interpretation of the results. BSSM and AZG wrote the original manuscript. All co-authors provided manuscript editing and feedback. All authors read and approved the final manuscript.
## Funding

This study was funded by the Spanish Ministry of Science and Innovation, project number PID2020-115558GB-C21, and by the Conselleria for Innovation, Universities, Science and Digital Society, project number CIAICO/2022/016.

## Availability of data and materials

The datasets used and analysed during the current study are available from the corresponding author on reasonable request. Genomic data of the ninth generation of intramuscular fat lines were deposited in Figshare Repository (https://doi.org/10.6084/m9.fishshare.9934058.v1).

## Declarations

### Ethics approval and consent to participate

Animal manipulations and the experimental procedures were approved by the Ethical Committee of the Polytechnic University of Valencia (reference number 2017/VSC/PEA/00212), according to Council Directives 98/58/EC and 2010/63/EC (European Economic Community, 1998).

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹Institute for Animal Science and Technology, Universitat Politècnica de València, 46022 Valencia, Spain.

Received: 9 August 2024 Accepted: 7 May 2025

Published online: 11 July 2025

## References

1. Biswas AK, Mandal P. Meat quality analysis: advanced evaluation methods, techniques, and technologies. New York: Academic Press; 2019.

2. Matarneh SK, Silva SL, Gerrard DE. New insights in muscle biology that alter meat quality. Annu Rev Anim Biosci. 2021;9:355-77.

3. Warris PD. Meat Quality. Meat Sci an Introd text. 2nd Ed. Wallingford, UK: CABi; 2010. p. 77-96.

4. Wood JD, Enser M, Fisher AV, Nute GR, Sheard PR, Richardson RI, et al. Fat deposition, fatty acid composition and meat quality: a review. Meat Sci. 2008;78:343-58.

5. Wood JD, Nute GR, Richardson RI, Whittington FM, Southwood O, Plastow G, et al. Effects of breed, diet and muscle on fat deposition and eating quality in pigs. Meat Sci. 2004;67:651-67.

6. DeVol DL, McKeith FK, Bechtel PJ, Novakofski J, Shanks RD, Carr TR. Variation in composition and palatability traits and relationships between muscle characteristics and palatability in a random sample of pork carcasses. J Anim Sci. 1988;66:385-95.

7. Sapp RL, Bertrand JK, Pringle TD, Wilson DE. Effects of selection for ultrasound intramuscular fat percentage in Angus bulls on carcass traits of progeny. J Anim Sci United States. 2002;80:2017-22.

8. Ochsner KP, Macneil MD, Lewis RM, Spangler ML. Economic selection index development for beefmaster cattle I: terminal breeding objective. J Anim Sci. 2017;95:1063-70.

9. Gao Y, Zhang R, Hu X, Li N. Application of genomic technologies to the improvement of meat quality of farm animals. Meat Sci. 2007;77:36-45.

10. Schwab CR, Baas TJ, Stalder KJ, Nettleton D. Results from six generations of selection for intramuscular fat in Duroc swine using real-time ultrasound. I. Direct and correlated phenotypic responses to selection. J Anim Sci. 2009;87:2774-80.

11. Liu L, Cui H, Xing S, Zhao G, Wen J. Effect of divergent selection for intramuscular fat content on muscle lipid metabolism in chickens. Animals (Basel). 2019;10:4.

12. Zhao GP, Chen JL, Zheng MQ, Wen J, Zhang Y. Correlated responses to selection for increased intramuscular fat in a Chinese quality chicken line. Poult Sci. 2007;86:2309-14.

13. Pannier L, Gardner GE, O'Reilly RA, Pethick DW. Factors affecting lamb eating quality and the potential for their integration into an MSA sheepmeat grading model. Meat Sci. 2018;144:43-52.

14. Zomeño C, Hernández P, Blasco A. Divergent selection for intramuscular fat content in rabbits. I. Direct response to selection. J Anim Sci. 2013;91:4526-31.

15. Martínez-Álvaro M, Hernández P, Blasco A. Divergent selection on intramuscular fat in rabbits: responses to selection and genetic parameters. J Anim Sci. 2016;94:4993-5003.

16. Sosa-Madrid BS, Hernández P, Blasco A, Haley CS, Fontanesi L, Santacreu MA, et al. Genomic regions influencing intramuscular fat in divergently selected rabbit lines. Anim Genet. 2020;51:58-69.

17. Sosa-Madrid BS, Varona L, Blasco A, Hernández P, Casto-Rebollo C, Ibáñez-Escriche N. The effect of divergent selection for intramuscular fat on the domestic rabbit genome. Animal. 2020;14:2225-35.

18. Hughes KA, Leips J. Pleiotropy, constraint, and modularity in the evolution of life histories: insights from genomic analyses. Ann N Y Acad Sci. 2017;1389:76-91.

19. Cheng J, Fernando R, Cheng H, Kachman SD, Lim K, Harding JCS, et al. Genome-wide association study of disease resilience traits from a natural polymicrobial disease challenge model in pigs identifies the importance of the major histocompatibility complex region. G3 (Bethesda). 2022;12:441.

20. Cheng H, Kizilkaya K, Zeng J, Garrick D, Fernando R. Genomic prediction from multiple-trait Bayesian regression methods using mixture priors. Genetics. 2018;209:89-103.

21. Blasco A, Ouhayoun J. Harmonization of criteria and terminology in rabbit meat research. Revised proposal. World Rabbit Sci. 1993;4:93-9.

22. Zomeño C, Hernández P, Blasco A. Use of near infrared spectroscopy for intramuscular fat selection in rabbits. World Rabbit Sci. 2011;19:203-8.

23. Legarra A, Varona L, López de Maturana E. TM: Threshold Model (available https://alegarra.github.io/manualtm.pdf). 2011.

24. Sargolzaei M, Chesnais JP, Schenkel FS. A new approach for efficient genotype imputation using information from relatives. BMC Genomics. 2014;15:478.

25. Wolc A, Dekkers JCM. Application of Bayesian genomic prediction methods to genome-wide association analyses. Genet Sel Evol. 2022;54:31.

26. Sosa-Madrid BS, Santacreu MA, Blasco A, Fontanesi L, Pena RN, Ibáñez-Escriche N. A genomewide association study in divergently selected lines in rabbits reveals novel genomic regions associated with litter size traits. J Anim Breed Genet. 2020;137:123-38.

27. Beissinger TM, Rosa GJM, Kaepper SM, Gianola D, de Leon N. Defining window-boundaries for genomic analyses using smoothing spline techniques. Genet Sel Evol. 2015;47:30.

28. Buss CE, Afonso J, de Oliveira PSN, Petrini J, Tizioto PC, Cesar ASM, et al. Bivariate GWAS reveals pleiotropic regions among feed efficiency and beef quality-related traits in Nelore cattle. Mamm Genome. 2023;34:90-103.

29. Cunningham F, Achuthan P, Akanni W, Allen J, Amode MR, Armean IM, et al. Ensembl 2019. Nucleic Acids Res. 2019;47:D745-51.

30. Jiao X, Sherman BT, Huang DW, Stephens R, Baseler MW, Lane HC, et al. DAVID-WS: a stateful web service to facilitate gene/protein list analysis. Bioinformatics. 2012;28:1805-6.

31. Carbon S, Douglass E, Good BM, Unni DR, Harris NL, Mungall CJ, et al. The gene ontology resource: enriching a GOld mine. Nucleic Acids Res. 2021;49:D325-34.

32. Wright MN, Gola D, Ziegler A. Preprocessing and quality control for whole-genome sequences from the Illumina HiSeq X Platform. In: Elston RC, editor. Stat Hum Genet Methods Protoc. New York: Springer; 2017. p. 629-47.

33. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114-20.

34. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754-60.

35. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297-303.

36. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80-92.
37. Casto-Rebollo C, Argente MJ, García ML, Pena R, Ibáñez-Escriche N. Identification of functional mutations associated with environmental variance of litter size in rabbits. *Genet Sel Evol.* 2020;52:22.

38. Zubiri-Gaitán A, Blasco A, Ccalta R, Satué K, Hernández P. Intramuscular fat selection in rabbits modifies the fatty acid composition of muscle and liver tissues. *Animals* (Basel). 2022;12:893.

39. ElNagar AG, Heddi I, Sosa-Madrid BS, Blasco A, Hernández P, Ibáñez-Escriche N. Genome-wide association study of maternal genetic effects on intramuscular fat and fatty acid composition in rabbits. *Animals* (Basel). 2023;13:3071.

40. Gianola D. Priors in whole-genome regression: the Bayesian alphabet returns. *Genetics*. 2013;194:573-96.

41. Dadousis C, Somavilla A, Ilska JJ, Johnsson M, Batista L, Mellanby RJ, et al. A genome-wide association analysis for body weight at 35 days measured on 137,343 broiler chickens. *Genet Sel Evol.* 2021;53:70.

42. Laghouaouta H, Sosa-Madrid BS, Zubiri-Gaitán A, Hernández P, Blasco A. Novel genomic regions associated with intramuscular fatty acid composition in rabbits. *Animals* (Basel). 2020;10:2090.

43. Huttner WB, Schmidt A. Lipids, lipid modification and lipid-protein interaction in membrane budding and fission—insights from the roles of endophilin A1 and synaptophysin in synaptic vesicle endocytosis. *Curr Opin Neurobiol.* 2000;10:543-51.

44. Kogelman LJA, Pant SD, Fredholm M, Kadarmideen HN. Systems genetics of obesity in an F2 pig model by genome-wide association, genetic network, and pathway analyses. *Front Genet.* 2014;5:214.

45. Gondret F, Le Floc'h N, Batonon-Alavo DI, Perruchot MH, Mercier Y, Lebret B. Flash dietary methionine supply over growth requirements in pigs: Multi-facetted effects on skeletal muscle metabolism. *Animal*. 2021;15:100268.

46. Gandolfi G, Mazzoni M, Zambonelli P, Lalatta-Costerbosa G, Tronca A, Russo V, et al. Perilipin 1 and perilipin 2 protein localization and gene expression study in skeletal muscles of European cross-breed pigs with different intramuscular fat contents. *Meat Sci.* 2011;88:631-7.

47. Gol S, Ros-Freixedes R, Zambonelli P, Tor M, Pena RN, Braglia S, et al. Relationship between perilipin genes polymorphisms and growth, carcass and meat quality traits in pigs. *J Anim Breed Genet.* 2016;133:24-30.

48. Zhang M, Lin W, Wu Q, Li Y, Chen D, Liu L, et al. Effects of dietary vitamin E on intramuscular fat deposition and transcriptome profile of the pectoral muscle of broilers. *J Poult Sci.* 2023;60:2023006.

49. Mancin E, Sosa-Madrid BS, Blasco A, Ibáñez-Escriche N. Genotype imputation to improve the cost-efficiency of genomic selection in rabbits. *Animals* (Basel). 2021;11:803.

50. Migdał L, Migdał A, Pałka S, Kmiecik M, Otwinowska-Mindur A, Semik-Gurgul E, et al. Polymorphism within IGFBP genes affects the acidity, colour, and shear force of rabbit meat. *Animals* (Basel). 2023;13:3743.

51. Ning C, Xie K, Huang J, Di Y, Wang Y, Yang A, et al. Marker density and statistical model designs to increase accuracy of genomic selection for wool traits in Angora rabbits. *Front Genet.* 2022;13:968712.

52. Ramayo-Caldas Y, Mármol-Sánchez E, Ballester M, Sánchez JP, González-Prendes R, Amills M, et al. Integrating genome-wide co-association and gene expression to identify putative regulators and predictors of feed efficiency in pigs. *Genet Sel Evol.* 2019;51:48.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.