Transcriptomic analysis in the liver of two rabbit lines divergently selected for intramuscular fat content

The intramuscular fat (IMF) content is an important indicator of meat quality, affecting the sensory properties of meat. IMF is a complex trait with polygenic nature. This research aimed to identify differentially expressed (DE) genes and key transcription factors (TFs) associated with IMF deposition in the liver of two rabbit lines divergently selected for IMF (high-IMF: H and low-IMF: L). We used 48 rabbits (24 H and 24 L) belonging to the 9th generation of selection to determine their liver gene expression levels using 3′ RNA sequencing. We found 308 DE genes between H and L lines; 134 upregulated and 174 downregulated in the H line. Among them, ACBD4, ACOT1, ACOT4, AIFM2, CPT1A, CPT1B, CROT, CYP4B1, CYP4A6, HADHB, IGFBP1, IGFBP2, FABP4, GPAT3, MPZL2, MYLIP, MRPL15, NR4A2, PLIN2, PCTP, PRRG4, PAQR9, RAB34, SLC16A11, SLC18A1, SLC2A1, TPMT, and UBC, were related to lipid metabolism, energy metabolism, purine/thiopurine metabolism, and/or molecule transport that could influence the IMF deposition. Of the 308 DE genes, 210 target genes were predicted for 11 DE TFs through the gene regulatory networks. Notably, 3 TFs, namely ETV1, NR4A1, and IKZF3 appear to modulate gene expression of several putative targets (including common targets: CREM, NR4A2, and RRS1 ). Functional analysis of DE genes revealed an overrepresentation of 40 enriched terms, including the PPAR signalling pathway (associated genes: CPT1B, FABP4, and PLIN2 ) and other lipid metabolism bioprocesses. Collectively, our results enhance the understanding of the liver-muscle crosstalk that contributes to IMF deposition and improvement of rabbit meat quality.

The intramuscular fat (IMF) content and fatty acid (FA) composition are important indicators of meat quality in livestock. IMF content influences the sensory properties of meat such as tenderness, juiciness, and flavour (Wood et al., 2008), and it can be modified by genetic selection (Schwab et al., 2009; Martínez-Álvaro et al ., 2016; Liu et al., 2019). A divergent selection experiment for IMF in the Longissimus thoracis et lumborum (LTL) muscle was performed in rabbits at the Universitat Politècnica de València for 10 generations (Zomeño et al., 2013; Martínez-Álvaro et al., 2016; Zubiri-Gaitán et al., 2022). The experiment yielded two rabbit lines: one with high IMF (H) content and the other with low IMF (L) content. Divergent selection was successful, obtaining a direct response to selection of 0.49 g/100 g in the 10th generation, equivalent to 3.8 standard deviations (SD) or to 47% of the mean of the trait. The FA composition of both muscle and liver was modified after selection as well (Martínez-Álvaro et al., 2018a; Zubiri-Gaitán et al., 2022). In addition, a correlated response was observed in blood plasma metabolome (Zubiri-Gaitán et al., 2023). IMF deposition in animals is influenced by multiple factors, including environment, sex, age, nutrition, and genetic factors (Hocquette et al., 2010; Nguyen et al., 2021). Both IMF content and FA composition are complex traits of polygenic nature (Cecchinato et al., 2012; Natacha Pena et al., 2016). Therefore, such traits can be genetically affected by the expression levels of multiple genes in different lipogenic tissues (e.g., fat, muscle, and liver), leading to complex interactions (metabolic crosstalk) between these tissues (Stefan et al., 2008; Alves-Bezerra and Cohen, 2018; Chen et al., 2024). The liver has special interest because it is one of the primary sites for de novo lipogenesis and is the main lipogenic site in growing rabbits (Gondret et al., 1997). The liver is highly metabolically active and, as the central organ for FA metabolism and energy homeostasis, it controls a wide array of functions, including FA uptake, oxidation, and metabolic conversion of non-esterified FAs (Wang et al., 2017; Alves-Bezerra and Cohen, 2018). Similarly, previous experiments with the aforementioned divergent rabbit lines found a greater liver size in the H line compared to the L line, which was attributed to greater lipogenic activity of enzymes FASN, G6PDH, and ME1 in the H line (Martínez-Álvaro et al., 2018b). Furthermore, changes in the liver FA composition were observed in these divergent lines (Zubiri-Gaitán et al., 2022). Numerous studies have focused on comparative analysis of the transcriptome profiles, typically with two groups under scrutiny (e.g., two breeds, two diets, divergent meat quality phenotypes, etc.) (Damon et al., 2012; Hamill et al., 2013; Óvilo et al., 2014). However, an analysis of the liver transcriptome profile of divergent rabbit lines for IMF may help to unravel the biological mechanisms explaining these differences and their link to IMF deposition. To the best of our knowledge, there are no previous studies comparing the liver transcriptome between rabbit lines differing in their IMF genetic background. This knowledge is crucial for understanding liver gene expression regulation and its link with muscle fat content, which is very complex and still not fully elucidated. The control of liver gene expression involves multiple biological processes and pathways, as well as multiple regulatory factors (Schrem et al., 2004; Gatti et al., 2007; Jump et al., 2013). Polyunsaturated fatty acids (PUFAs) are recognized to affect gene expression by regulating the activity or abundance of several families of TFs, including PPARs, LXRs, HNF-4, and SREBPs (Jump, 2002). Besides, n6-methyladenosine (m 6 A) modification in perirenal fat and longissimus lumborum muscle tissues of Hycole rabbit (meat breed) has been detected by Methylated RNA immunoprecipitation sequencing (MeRIP-Seq) (Luo et al., 2023). The analysis revealed that differential methylases (ZC3H13, METTL4, and IGF2BP2) and differential genes (e.g., ABCA1, ADRB1, FABP3, GLI2, PPKAG3, and SOX9 ) between two tissues were modified by m 6 A, both regulating IMF deposition through lipid signaling pathways. Moreover, a recent study of whole transcriptome sequencing of the longissimus dorsi muscles of Hycole and Rex rabbit breeds, identified several candidate genes ( STARD13, HSPB8, and CBLB ) and lncRNAs (e.g., TCONS_00062393 and TCONS_00031457 ) associated to fat deposition at different developmental stages (Wang et al., 2025). Interestingly, the 3’ untranslated regions (UTRs) of the vast majority of mRNA transcripts have been shown as regulatory elements and hubs for the post-transcriptional gene regulation (Matoulkova et al., 2012; Schuster and Hsieh, 2019). 3′ UTRs can establish 3′ UTR-mediated protein-protein interactions (PPIs), which allows them to transmit genetic information encoded in 3′ UTRs to proteins (Mayr, 2019). Ultimately, the 3’ RNA sequencing (3′ RNA-Seq) provides an alternative to whole transcript analysis with similar levels of reproducibility that the standard RNA-Seq method. Moreover, the 3′ mRNA-Seq method is particularly useful for detecting of short transcripts (Ma et al., 2019) and is the approach of choice for accurately determining gene expression at the lowest cost (Moll et al., 2014). Additionally, 3′ RNA-Seq is superior to standard RNA-Seq in cases of sparse data, which relies on enriching for 3′ ends of the transcript and its use is considered as a proxy for expression of the whole gene (McClure et al., 2023). In the current study, we tested the hypothesis that the different IMF deposition in the divergent rabbit lines (H and L) could be associated with changes of gene expression levels in the liver. Our experimental material allows us to explore the genetic link between IMF content and liver gene expression, considering that both lines originated from the same synthetic line, were contemporaneously selected on the same farm and with the same diet and management. Therefore, the aim of this work was to identify differentially expressed (DE) genes in the liver of H and L lines by comparing their hepatic transcriptome obtained by 3′ RNA-Seq. Additionally, we aimed to elucidate the gene regulatory networks and key transcription factors (TFs) that significantly influence the expression levels of DE genes and modulate IMF deposition.

Animals and experimental design Data for this study come from the 9th generation of a 10-generation divergent selection experiment for IMF content measured at 9 weeks of age in the LTL muscle. In each generation, the young rabbits were weaned at 28 days and were then reared jointly within litters until slaughter at 63 days of age. All animals were reared on the same farm, under the same environmental conditions, with the same management, free access to water and fed ad libitum with a cereal-based commercial diet. Details of the divergent selection procedure, including the base population, number of individuals, inbreeding, and selection criterion, are described in Martínez-Álvaro et al. (2016). Further details about ingredient and chemical composition of the diet are in Zubiri-Gaitán et al . (2023). At the 9th generation, the selection response was around 3.1 SD or 41% of the mean of the trait, calculated as the phenotypic difference between lines (Sosa-Madrid et al., 2020b). A selection of 48 animals (24 from H and 24 from L IMF lines) from the 9th generation of selection was made based on the normal distribution of their IMF content within each line. The data were divided into four quartiles, and a balanced number of samples were drawn randomly from each quartile until 24 animals per line of both sexes (12 per sex) were obtained. The animals selected came from 45 litters derived from 45 different mothers and 19 different fathers, with 1 or 2 animals sampled within each litter. Animals were slaughtered after 4 hours of fasting by exsanguination after electrical stunning, and the carcasses were chilled at 4 ºC for 24 h. After chilling, the liver was excised, and liver samples were placed in 2 mL cryogenic tubes with 500 µl of RNA later (QIAGEN). The samples were then stored at ‑80°C until the analysis. The 48 liver samples were processed by the biotechnology company Seqplexing (Valencia, Spain). Total RNA was isolated from the liver samples using the MagMAX mirVana Total RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s recommendations. RNA quality was assessed using an Agilent 2100 Bioanalyzer system (Hercules, CA, USA), and 48 samples with an RNA quality number (RQN) greater than 7 were sequenced. Transcriptome sequencing was performed using the Lexogen QuantSeq 3′ mRNA-Seq Library Prep FWD kit (Lexogen GmbH, Austria) with unique molecular identifiers (UMI). To create RNA libraries incorporating UMIs, the guidelines provided by the manufacturer’s protocol were followed. The quality of the library was evaluated using the QIAxcel Advanced System (Qiagen). Lastly, sequencing was performed at 10× of depth on the Illumina NovaSeq X platform with 150 bp paired-end reads (average of 5.95 million reads/sample), which was ultimately used for the quantitative study of gene expression levels in liver. Bioinformatics and statistical analyses The raw sequences (FASTQ files) were processed using an in-house bioinformatics pipeline developed by the Seqplexing company. The methodology consisted of: (1) quality control with FastQC software v0.12.1 (Andrews, 2010) and trimming with Cutadapt tool v1.2.0 (Martin, 2011), (2) UMI (Unique Molecular Identifier) processing using umi-tools software v0.2.3 (Smith et al., 2017), (3) alignment against the reference genome ”UM_NZW_1.0 (RefSeq GCF_009806435.1)” ( Oryctolagus cuniculus ) using the STAR software v2.7.11a (Dobin et al., 2013), and (4) expression quantification using the HTSeq-count program v2.0.5 (Anders et al., 2015), in order to calculate the counts of each gene. The matrix of transcripts abundance in counts obtained (n = 37667 transcripts by 48 samples) was pre-processed as follows: transcripts that did not have at least 10 counts in 12 or more samples (25%) or were not protein-coding genes were removed. After these filers, we kept 11135 transcripts for further analysis. Analysis of differential expression, regulator-target network, enrichment analysis, regulator scoring and ranking, and annotation of TF-target gene interactions Statistical analyses were performed using R program v4.4.0 (R Core Team, 2024). We applied a comprehensive five-step analytical pipeline to identify DE genes and key TFs associated with IMF deposition, using the RegEnrich R package v1.12.1 (Tao et al., 2022) and the TFLink gateway (Liska et al., 2022). A RegenrichSet object was initially constructed, incorporating the expression data matrix with 11135 transcripts and 48 samples, the metadata information, and a list of TFs identified in the expression data. The TFs annotations for “ Oryctolagus cuniculus ” were corroborated in the AnimalTFDB v4.0 database (Shen et al., 2023). The first step aimed to identify which of the 11135 transcripts were DE between the lines (H and L). For this differential expression analysis (DEA), we used the ”Wald_DESeq2” method (Love et al., 2014) implemented within the RegEnrich approach. The model assumed a negative binomial distribution of the 11135 transcript abundances and included sex (2 levels) and line (2 levels) as fixed effects. Counts within sample were normalized by sequencing depth. Differential expression of each transcript between the lines was evaluated based on p -values estimated using the Wald test and adjusted by the Benjamini–Hochberg (BH) method (Benjamini and Hochberg, 1995) to control the false discovery rate (FDR). Additionally, the logarithmic fold change (log 2 FC) between the lines was also considered for the analysis (Love et al., 2014). All genes with an absolute FC of at least 1.5 (|log 2 FC| ≥ 0.58) and a BH padj < 0.05 were selected as DE genes. We used the normalized gene expression levels of DE genes to construct a complementary heatmap via the ComplexHeatmap v2.14.0 package (Gu et al ., 2016), with the aim of determining the hierarchical clustering for sets of DE genes that exhibit similar expression patterns in all 48 samples. The second step aimed at investigating which of the 11135 transcripts may be putative downstream targets of the list of TFs identified (potential regulators). Here, a regulator-target network inference (Tao et al., 2022) was performed using the gene regulatory network method (Huynh-Thu et al., 2010). This computational method employs the Random Forest (RF) algorithm (Breiman, 2001), to quantify the strength of the relationship (edge weights) between TFs and the DE genes (i.e., of each pair, the nodes). The following parameters were used: RsquaredCut = 0.85, minR = 0.30, edgeThreshold = 0.01, nperm = 10000, nbTrees = 1000, topNetPercent = 10%, qvalueCutoff = 0.05, importanceMeasure = IncNodePurity; among other default parameters. We selected the TFs whose putative targets were identified as DE in step one, and then visualized the regulatory network using Cytoscape v3.10.1 software (Shannon et al., 2003). This step was corroborated using the publicly available annotation for Homo sapiens species in the curated database called TFLink gateway (Liska et al., 2022). Here, we studied the overlap of the prediction of the putative target genes for each DE TFs against the TFLink database of TF-target gene interactions, which includes interactions validated by small-scale experiments, large-scale experiments, or both (Liska et al., 2022). The third and fourth steps are focused on TFs. The third step aimed to identify the key DE TFs together with their DE putative targets that were enriched according to a TF-enrichment analysis (EA) (Tao et al., 2022). In the enrichment task, Fisher’s exact test (FET) was used, and TFs were considered potential regulators of associated genes when the BH padj < 0.05. The fourth step aimed to summarize the importance of the key TFs identified in the previous step by considering information from both DEA and EA. This step is called regulator scoring and ranking (Tao et al., 2022), in which after the EA, the overall ranking scores of regulators were calculated as: \(score\ =\ f(-\log_{10}(P_{E}))\ +\ f(-\log_{10}(P_{D}))\); where\(f\left(x\right)=\ \frac{x-\min(x)}{\max\left(x\right)-\min(x)}\), P D is the vector of p \(-\)values of regulators obtained from DEA, and P E is the vector of p \(-\)values of regulators obtained from EA. Additionally, two thresholds were considered for the TFs to quantify the expression changes between the lines, including a conservative threshold (|log 2 FC| ≥ 0.26 or absolute FC > 1.2) and a more stringent threshold (|log 2 FC| ≥ 0.58 or absolute FC > 1.5). We ranked the TFs by family according to annotation of AnimalTFDB v4.0 database (Shen et al., 2023), and the score calculated above. In the fifth step we carried out a functional analysis with the list of DE genes. This list was submitted to the ClueGO v2.5.10 plug-in in Cytoscape v3.10.1 software (Shannon et al., 2003; Bindea et al., 2009) for gene ontology (GO) and pathways analysis using the Oryctolagus cuniculus database. GO terms (biological processes, molecular function, and pathways) statistically overrepresented were determined using a Fisher test, and this time p -values were corrected by the BH method (Benjamini and Hochberg, 1995). Results with and without the “GO Term Fusion” feature were generated to assess the redundant parent–child terms and extract non-redundant associations of GO-based gene functions. Finally, the Pearson correlation between the normalized expression levels of each DE gene and IMF was calculated using the cor.test function of base R. With this approach, we also tested whether the expression of both the TFs and their targets was related to the trait.

and Discussion Description of the quality control and read alignments using 3′ RNA-Seq The Illumina NovaSeq X platform generated an average of 30.24 million of 150 bp paired-reads per sample, with Q20 > 98%, Q30 > 94% and GC > 38% across all samples, while the mean proportion of unique fragments was 39.87% (Table S1). Furthermore, the unique mapped reads against the reference rabbit genome consisted of an average of 78.53% (ranging from 72.77% to 82.91%). Q20 and Q30 metrics exceeding 90% and 70%, respectively, in all samples indicated reliable data. Additionally, the incorporation of UMIs ensured that the duplicated reads found in the analysis belonged to different molecules, thereby providing greater confidence in the results. Differential expression analysis, gene clustering, and relationships with IMF content A total of 11135 genes were tested to identify DE genes between the H and L IMF samples (Table S2). In total, 308 DE genes were identified between the two divergent lines, based on the previously described threshold. Figure 1a presents a volcano plot displaying these DE genes according to the pairwise comparison (H vs. L), with 134 upregulated (green points) and 174 downregulated genes (red points) in the H line. The remaining genes that did not pass the threshold are shown as grey points. For instance, the DE genes that exhibited the most extreme expression changes (|log 2 FC| ≥ 4.05, Table S2) were LOC100342501 or UBB, LOC127486578 or IFT46, LOC100347496, LOC127486882, TINAG, LOC127489097, NPHS2, and LOC100358583 (upregulated in the H line ) and LOC100343656, MRPL15, LOC100134865 or CDK2, and LOC127489014 or KLHL11 (downregulated in the H line ). Complementary, figure 1b illustrates a heatmap showing the hierarchical clustering of gene expression across all 48 samples. The upper dendrogram separates the samples into two major clusters corresponding to the H and L lines, highlighting the relevant impact of the genetic line on overall gene expression. Conversely, the left dendrogram further classifies the 308 DE genes into two primary clusters corresponding to the genes upregulated (n = 134) and downregulated (n = 174) in the H line, reflecting the distinct expression patterns of these genes between the divergent lines. These findings highlight the prevailing gene expression differences between H and L lines as a consequence of the divergent selection by IMF content. As both rabbit lines were contemporarily reared on the same farm, under identical dietary and management conditions, the observed differences in IMF are likely attributable to changes in allele frequencies and gene expression profiles. The present study contributes to the identification of a panel of candidate genes involved in hepatic lipid metabolism that may be associated with intramuscular fat deposition. Noteworthy putative candidate genes for IMF with upregulated expression in H line that clustered together included ELOVL6 (lipogenesis and purine metabolism), BACH2 and RAB34 (post-transcriptional regulation), MPZL2 (adhesion molecule), PENK (neuropeptide hormone activity), and TPMT (thiopurine metabolism). Conversely, downregulated genes in the H line comprised several involved in FA and lipid metabolism, including ACBD4, CREM, CYP4B1, FABP4, IL1RN, LOC100344509 (or ACOT1 ), PCTP, and PLIN2 ; genes related to FA β-oxidation such as ACAD10, CROT, HADHB, LOC103346593 (or CPT1A ), and CPT1B ; and genes associated with adipogenesis ( GPAT3 and RBM4 ), transport functions ( MYLIP, SLC18A1, SLC2A1, and SLC41A3 ), lipid-related transcription regulation ( ETV1, IKZF3, and NR4A1 ), glucose regulation, insulin metabolism and hormonal status ( IGFBP1 and IGFBP2 ), protein metabolism ( MRPL15 and USP18 ), ubiquitin homeostasis ( UBC ), and ribosome biogenesis/signal transduction ( RRS1 ). The Pearson correlation between IMF content and the abundance of each of the 308 DE genes in liver is shown in Table S3. Correlation coefficients ranged from -0.22 to -0.74 and from 0.20 to 0.72, negative correlations (56.49%) being more common than positive correlations (43.51%). Genes showing the highest correlation with IMF (i.e., |r| ≥ 0.70) included PRRG4, LOC100358950 or RNF5, RFLNA and LOC100344872 or BTN3A3 (all upregulated in the H line ) and LOC127489014 or KLHL11, LOC127482753 or BRI3BP, and LOC100351163 (all downregulated in the H line). We also identified nine genes likely to be involved in hepatic lipid metabolism, energy metabolism, and thiopurine metabolism that correlated with IMF content ( ACAD10, ACBD4, CPT1B, GPAT3, HADHB, MALL, MYLIP, SLC16A11, UBC, and TPMT ), all of which exhibited a |correlation| > 0.50 (Table S3). Given the regulatory role of TFs and co-factors in gene expression and metabolic pathways, we also examined their correlation with the trait. Correlation coefficients ranged from -0.23 to -0.58 for ETV1, IKZF3, NR4A1, NR4A2, RAB11FIP4, and TAL2, and from 0.24 to 0.43 for BACH2, BHLHE22, CREB5, and NFKB2 (Table S3). Several of the DE genes identified in this study have been previously reported in genomics analyses of these rabbit lines. For example, NUSAP1 and TRDMT1, which are upregulated in the H line, and PLIN2 and TYRO3, which are downregulated in the H line, were previously proposed as part of the genes ascribed to selection signatures for IMF content in these divergent rabbit lines (Sosa-Madrid et al ., 2020a). Additionally, PLIN2 plus the upregulated gene ART4 were reported in a GWAS as functional candidate genes influencing IMF content and FA composition (Laghouaouta et al ., 2020; Sosa-Madrid et al., 2020b); while PLIN2 was also detected as a potential candidate gene influencing several weight traits and fat depot traits using a multivariate GWAS (Sosa-Madrid et al ., 2025). Among these five genes, PLIN2 is of particular interest because of its relevant function in lipids metabolism. PLIN2 is known as a member of gatekeepers of intracellular lipolysis and its protein level plays a key role in governing lipid droplet dynamics and their relationship to mitochondria (Xu et al., 2019). The fatty acid binding protein 4 (FABP4), a gene involved in IMF content, was also found to be differentially expressed in these lines, with its expression decreased in the H line. FABP4, also known as A-FABP or aP2, encodes a FA binding protein found in adipocytes. According to Furuhashi et al. (2014), FABP4 is released from adipocytes by a non-classical pathway associated with lipolysis, possibly acting as an adipokine. Furthermore, the ELOVL fatty acid elongase 6 ( ELOVL6 ), an interesting lipogenic gene, showed increased expression in the H line compared to the L line. ELOVL6 is a microsomal enzyme involved in the elongation of saturated and monounsaturated FAs with 12, 14, and 16 carbons, as well as in the biosynthesis of unsaturated FAs (Matsuzaka and Shimano, 2009). Consistent with our results, PLIN2 and FABP4 have been previously reported as candidate genes for IMF deposition in rabbits (Migdał et al., 2018; Ahamba et al., 2024). Likewise, ELOVL6 has been reported as a candidate gene associated with intramuscular FA composition in Iberian × Duroc backcross and Duroc × Luchuan crossbred pigs ( Liu et al., 2021; Valdés-Hernández et al., 2023, 2024). IMF deposition may be highly influenced by the expression of hepatic genes involved in FA and lipid metabolism. Among the DE genes with larger differences in expression, LOC100343656 (or BTN1A1) and LOC100134865 (or CDK2 ) were downregulated in the H line. The butyrophilin subfamily 1 member A1 gene ( BTN1A1 ), which showed the largest expression change, encodes a member of the butyrophilins molecules involved in signaling receptor binding. Notably, BTN1A1 knockout has been shown to alter lipid droplet formation and phospholipid composition in mammary epithelial cells (Han et al., 2020). On the other hand, the cyclin-dependent kinase 2 ( CDK2 ) gene participates in the sequential phosphorylation and activation of the CCAAT/enhancer binding protein-β (C/EBPβ) during adipogenic differentiation (Li et al., 2007). In this study, we also found several DE genes related to FA β-oxidation, which were all downregulated in the H line, namely CPT1A, CPT1B, CROT, and HADHB . CROT, CPT1A, and CPT1B are known to encode different active forms of related carnitine acyltransferases enzymes (Van der Leij et al., 2000). The carnitine octanoyltransferase ( CROT ) is a peroxisomal enzyme catalyzing very long-chain FAs (VLCFAs) conversion to medium-chain FAs (MCFAs) that can be then absorbed by mitochondria during β-oxidation (Sanford et al., 2023). Conversely, the carnitine palmitoyltransferase 1A ( CPT1A ) and carnitine palmitoyltransferase-1B ( CPT1B ) are two of the three isoforms of the carnitine palmitoyltransferase I (CPT1) protein family. These enzymes catalyze the initial and regulated step in the β-oxidation of FAs in the mitochondria. CPT1A is the primary liver isoform and its regulation is complex, involving several layers that include genetic, epigenetic, physiological, and nutritional modulators (Schlaepfer and Joshi, 2020). Instead, CPT1B is primarily known as the skeletal muscle isoform (Maples et al., 2015; Van der Leij et al., 2000). In addition to CPT1B expression, a study performed in humans showed that differential DNA methylation may partially explain the reduced expression of CPT1B in obese individuals (Maples et al., 2015). Moreover, the hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit beta ( HADHB ) catalyses 3-ketoacyl-CoA thiolase (LCKT) activity, participating in the final step of mitochondrial β-oxidation of long-chain fatty acids (LCFAs) (Dagher et al ., 2021). Contrary to what was observed in the muscle of these lines (Martinez-Alvaro et al., 2017; Zubiri-Gaitán et al ., 2023), the present results suggest that the H line, exhibiting lower expression of these key enzymes involved in FA β-oxidation, likely has a reduced capacity in the carnitine shuttle system, leading to decreased FA β-oxidation in the liver of this line. Consistent with these findings, we found decreased expression of insulin like growth factor binding protein 1 ( IGFBP1 ) in the H line. IGFBP1 is known to be associated with glucose regulation and insulin resistance, and has been shown to downregulate hepatic lipogenesis while upregulating lipid β-oxidation (Pan et al., 2021). These regulatory roles align with the previously reported lower hepatic lipogenesis (Martínez-Álvaro et al., 2018a) and the greater abundance of β-oxidation and acyltransferase enzymes found in the L line. Earlier analyses performed in these lines have shown the relevance of β-oxidation in the different fat deposition between the lines (Martinez-Alvaro et al., 2017; Zubiri-Gaitán et al ., 2023). Nonetheless, these studies indicated differences in the oxidation of FA in the muscle suggesting lower β-oxidation in the L line, contrary to what it was observed in the liver. These results reinforce the existence of complex mechanisms and interactions between the liver, muscle, and other lipogenic tissues involved in fat deposition. Finally, we identified other interesting DE genes involved in adipogenesis, such as GPAT3 and RBM4 (downregulated in the H line), and RAB34 (upregulated in the H line). The glycerol-3-phosphate acyltransferase 3 ( GPAT3 ) encodes for a member of the acyltransferase family of proteins (termed GPAT) (Takeuchi and Reue, 2009). The encoded protein is an enzyme that catalyses the conversion of glycerol-3-phosphate to lysophosphatidic acid (LPA) in the synthesis of triacylglycerol. Interestingly, the previous untargeted metabolomics analysis revealed that glycerol-3-phosphate was relevant in discriminating between H and L lines. However, this metabolite showed a 98% probability of being more abundant in the L line (Zubiri-Gaitán et al ., 2023), supporting the complex metabolism involved in fat deposition. Conversely, the ras-related protein rab 34 gene ( RAB34 ) plays important roles during adipocyte differentiation and it is involved in regulating the oligomerization and secretion of adiponectin, a key adipokine with insulin-sensitizing actions, as well as to lipid storage and mobilization from lipid droplets (López-Alcalá et al., 2024). Identification of transcription factors, characterization of their expression changes between the lines, gene regulatory network and regulator enrichment analysis We identified a total of 588 TFs among the 11135 transcripts, which were considered candidate regulators of gene expression in liver (Table S4) and are potentially involved in the regulation of IMF deposition. Of those 588 TFs, 36 were DE between lines based on the threshold of |log 2 FC| ≥ 0.26 (Figure 2), among which 11 ( BACH2, BHLHE22, CREB5, ETV1, IKZF3, NFKB2, NR4A1, NR4A2, RAB11FIP4, TAL2, and ZFP57 ) exhibited a |log 2 FC| ≥ 0.58. Figure 2 shows a summary of these 36 TFs grouped by family (16 families) and ranked by their score as explained in the methodology section. Four families were found with the highest number of DE TFs: C2H2 (n = 11), TF_bZIP (n = 4), ETS (n = 4), and bHLH (n = 3). To be noted, TFs such as BACH2, ETV1, IKZF3, NFKB2, NR4A1, RAB11FIP4, and TAL2 were observed within the top 25 with the highest score (range from 0.97 to 1.46, Figure 2). Using the gene expression of the 11135 identified transcripts, we explored the regulator-target relationships between TFs and putative targets. Results of the gene-regulatory network and regulator EA for the 588 identified TFs, including the putative targets associated with each TF, are presented in Table S5 and Table S6, respectively. Focusing on TF-target gene interactions for the 11 DE TFs mentioned above, and keeping 10% of the best ranked edges, we identified a total of 210 putative target genes that were DE between the lines. Details of the Pearson correlations between the 11 TFs and their predicted target genes are shown in Table S5. The 11 DE TFs with their 210 DE putative targets were represented in the regulator-target network using the weights inferred by GNR approach (Figure 3). The network consisted of 221 nodes and 511 edges. For each TF a variable number of predicted targets was detected (between 24 and 67), with weights (edges) ranging from 0.004 to 0.04 between gene pairs (nodes). Overall, the most interconnected regulators in the network (i.e., exhibiting the highest number of predicted targets) were NR4A1, ETV1, IKZF3, NFKB2, and RAB11FIP4 ( Figure 3 ). Gene regulatory network (GRN) analysis is useful to elucidate the interactions of TFs with their putative targets. In our work, we focused on identifying those key TFs whose expression and that of their putative targets (including co-factors) can be associated with the IMF phenotype. Notably, the most interconnected TFs in the gene networks were ETV1 and NR4A1 (n = 67 each), both downregulated in H line. These TFs had 44 targets in common, which included several candidate genes for lipid and energy metabolism, such as CREM, CROT, CYP4B1, FABP4, HADHB, IGFBP1, IGFBP2, IL1RN, IKZF3, NR4A2, MRPL15, RRS1, PLIN2, LOC103346593 or CPT1A, and SLC41A3 . Among the non-shared target, ETV1 was uniquely associated with genes as ACAD10, CYP4A6, OLAH, and SLC25A30, while NR4A1 uniquely targeted AIFM2, MTHFD1L, LOC100344509 or ACOT1, PAQR9, SLC25A33, SLC25A42, and SLC2A1 . Similarly, IKZF3 TF also had a high number of predicted targets (n = 59), of which 19 are believed to be involved in lipid and energy metabolism ( AGAP2, CREM, CROT, CYP4A6, CYP4B1, ETV1, FABP4, HADHB, IGFBP1, IGFBP2, IL1RN, OLAH, LOC100344509 or ACOT1, LOC103346593 or CPT1A, NR4A1, NR4A2, UBC, SLC41A3, and SLC25A33 ). Subsequently, target overlapping results suggest that CREM, NR4A2, and RRS1 were predicted common targets for the three most interconnected TFs ( ETV1, NR4A1, and IKZF3 ). Furthermore, we identified 43 predicted targets for BACH2, including 4 of the top 12 genes with the highest expression differences between the lines ( LOC100342501 or UBB, LOC100347496, LOC100358583, and MRPL15 ), and other putative targets involved in lipid metabolism and glucose metabolism, such as CPT1B and ELOVL6, PPP1R3B and PRRG4. Three of these TFs, namely ETV1, NR4A1, and IKZF3, showed an interesting pattern of relationship, with NR4A1 being a common target of ETV1 and IKZF3; ETV1 being a target of IKZF3; and IKZF3 being a target of NR4A1 . In addition, in NFKB2 we identified a target known to be TF co-factor (i.e., PCGF2 ). The nuclear receptor subfamily 4 group A member 1 ( NR4A1, also known as Nur77 or NGFIB ) and the BTB domain and CNC homologue ( BACH2 ) have been reported as relevant genes showing strong post-transcriptional signals in porcine skeletal muscle comparing fasted vs. fed Duroc gilts (Mármol-Sánchez et al., 2022). An altered function (hyper- or dys- regulation) of NR4A1 is associated with various metabolic processes, including carbohydrate metabolism, lipid metabolism, and energy balance in key metabolic tissues such as the liver, skeletal muscle, pancreas, and adipose tissue (Zhang et al., 2018). It is worth mentioning that ETV1 has been reported to regulate lipid metabolism, glycolysis and glutamine metabolism (Baena et al., 2013). In addition, Cooper et al . (2015) through structural and biochemical analysis of ETV1 elucidated mechanisms of post-translational regulation in ETS proteins. Another important TF influencing lipid metabolism is the IKAROS family zinc finger 3 ( IKZF3 ). Paniri et al . (2024) have reported that IKZF3 together with long non-coding RNAs cyclin-dependent kinase inhibitor 2B (lncRNA CDKN2B-AS1 ), and microRNA lethal 7-a2 ( let7-a-2 ) are involved in lipid metabolism disorders. Other targets of interest were the polyubiquitin-B ( UBB ; upregulated in the H line) and the myosin regulatory light chain interacting protein gene ( MYLIP ; downregulated in the H line). UBB was the second gene with the largest expression changes between lines and encodes ubiquitin, a protein that serves as an essential core unit of the post-translational modification process known as ubiquitination (Han et al., 2021). Ubiquitin plays a crucial role in regulating lipid metabolism by controlling the turnover of proteins and lipids involved in lipid metabolism (Loix et al., 2024). On the other hand, MYLIP, also known as IDOL, was found to be a common putative target for CREB5, RAB11FIP4, TAL2 and ZFP57 . Interestingly, it is a direct target for regulation by liver X receptors (LXRs), and its expression is responsible to cellular sterol status independent of the sterol-response element-binding proteins (Hong et al., 2010; Weissglas-Volkov et al., 2011). Overall, it is well established that the regulatory mechanisms of lipid deposition in animals include the expression of genes associated with triglyceride biosynthesis pathways, genes encoding transmembrane transport proteins, and FA desaturases and elongases genes (Yi et al., 2023). Our study also suggests that genes involved in signal transduction and post-translational modifications play a pivotal role in regulating FA and lipid metabolism. The RegEnrich approach was useful to identify multiple DE genes between groups, as well as key TFs that play the most important role in changing gene expression profiles using a network-based method (Tao et al., 2022). Overlap between regulator-target prediction by RegEnrich versus TFLink database In order to reinforce the interpretation of the results generated with the RegEnrich approach, we complementarily provide a screening of curated information on TF-target gene regulatory interactions for Homo sapiens organism (Table S7). The analysis revealed that 8 of the 11 TFs mentioned above ( BACH2, BHLHE22, CREB5, ETV1, IKZF3, NFKB2, NR4A1, and NR4A2 ) were annotated in the TFLink gateway. However, for two of them ( CREB5 and NR4A2 ), no matching records were found for their predicted target genes. Table S7 shows the overlap for the remaining 6 TFs and their respective target genes, revealing 54 common targets that were shared by RegEnrich and TFLink approaches. Among these 54 targets, we detected 12 candidate genes of lipid and FA metabolism ( ACAD10, ACBD4, CREM, CROT, CYP4B1, HADHB, IL1RN, NR4A2, PCTP, PLIN2, SLC2A1, and UBC ). On the other hand, similar to what we found using the RegEnrich approach, we observed several TFs as putative targets of other TFs, like NR4A1 for ETV1 and IKZF3, ETV1 for IKZF3, IKZF3 for NR4A1, CREB5 for BACH2, and NR4A2 for IKZF3 and NR4A1 ( Figure 3 ). By exploiting TFLink gateway we can obtain comprehensive and accurate information on TF-target gene interactions derived from the human model, integrating data from several sources such as TF databases, experimental methods and publications (Liska et al., 2022). We utilized this information enabling cross-species (human-rabbit) regulatory analysis and supporting research into gene regulation. This provides indirect inter-species confirmation of our results and those previously identified in humans and, notably, offers the opportunity to discover genetic differences due to the expression of distinctive genes. However, further evidence derived from experimental validation (e.g., using ChIP-seq, CRISPR Knockouts or functional assays) is needed to elucidate the roles of TFs such as BACH2, BHLHE22, CREB5, ETV1, IKZF3, NFKB2, and NR4A1 in muscle and prove their possible involvement in the control of lipid target genes. Among them, NR4A1 and NFKB2 can be investigated due to their known metabolic roles in other tissues, as well as BHLHE22 for investigating proinflammatory immune microenvironment or circadian rhythm interactions. Functional enrichment analysis to assign biological relevance of DE genes To gain deeper insight into the biological function of 308 DE genes in liver associated with IMF deposition, we investigated GO terms and pathways overrepresented. Functional analysis revealed 40 significantly overrepresented terms, including 31 pathways, 8 biological processes and 1 molecular function (Figure 4). A detailed description of the DE genes associated to each of the 40 terms is shown in Table S8. As shown in Figure 4a, the most prominent GO term was the ribosome pathway (KEGG:03010). Other functional annotations of interest found related to lipid metabolism included fatty acid elongation (KEGG:00062), peroxisome proliferator-activated receptor (PPAR) signalling (KEGG:03320), glutathione metabolism (KEGG:00480), and glutathione transferase activity (GO:0004364). Additional lipid metabolism-related terms enriched were arachidonic acid metabolism (KEGG:00590), retinol metabolism (KEGG:00830), biosynthesis of unsaturated fatty acids (KEGG:01040), fatty acid degradation (KEGG:00071), insulin resistance (KEGG:04931), fatty acid transport (GO:0015908), monocarboxylic acid transport (GO:0015718), and regulation of glucose metabolic process (GO:0010906). Seventy of the 308 DE genes were annotated into 12 different functional groups (Figure 4b; Table S8). For instance, genes such as FABP4 and PLIN2 were allocated in three terms (PPAR signalling pathway, fatty acid transport, and monocarboxylic acid transport). ELOVL6 was present in four terms: fatty acid elongation, biosynthesis of unsaturated fatty acids, purine nucleoside bisphosphate metabolic process (GO:0034032) and ribonucleoside bisphosphate metabolic process (GO:0033875). Remarkably, the genes LOC100342501 or UBB and LOC100347496, which exhibited the highest expression differences between IMF-divergent lines, were detected in PPAR signalling pathway, fatty acid elongation, and fatty acid degradation. HADHB was allocated in two terms (fatty acid elongation and fatty acid degradation), while CPT1B was found in four terms: PPAR signalling pathway, fatty acid degradation, insulin resistance and glucagon signalling pathway (KEGG:04922). Other genes of interest allocated in multiple terms were CROT (purine nucleoside bisphosphate metabolic process, ribonucleoside bisphosphate metabolic process, fatty acid transport, and monocarboxylic acid transport), CYP4A6 (arachidonic acid metabolism, retinol metabolism, fatty acid degradation, and PPAR signalling pathway), and CYP4A6 and CYP2C16 (arachidonic acid metabolism and retinol metabolism). As other examples of multiple annotations, we observed the presence of genes such as: CROT, FABP4, and PLIN2 together with SLC16A11 and SLC51B in monocarboxylic acid transport, CPT1B together with CREB5, LOC103346593, PPP1R3B, and SLC2A1 in insulin resistance, as well as PPP1R3B together with GCKR and TIGAR in the regulation of glucose metabolic processing. Finally, three TFs ( CREB5, NR4A1, and NR4A2 ) were detected in the aldosterone synthesis and secretion pathway (KEGG:04925). One of the most interesting findings was the overrepresentation of the peroxisome proliferator activated receptor (PPAR) signalling pathway. PPARs are nuclear receptors functioning as ligand-activated TFs that regulate the FA metabolism and modulate gene expression of target genes depending on the presence of co-repressors or co-activators (Wagner and Wagner, 2020). These receptors bind and are activated by a broad range of FAs and FA derivatives, and they thereby serve as major FA sensors controlling metabolism (Poulsen et al., 2012). Additionally, PPARs may influence the transcription of genes regulating several lipid metabolism pathways such as lipolysis, FA uptake, oxidation and lipogenesis; thereby modulating hepatic triglyceride accumulation (Tailleux et al ., 2012). Regulation of these pathways is complex and involves nuclear receptors, membrane transport proteins and cellular enzymes, and together with the amplification of the FA transporters or FA uptake (e.g., CAV1 and CD36 ), would result in elevated synthesis and accumulation of FAs (Nath and Chan, 2016). On the other hand, the acyl-CoA thioesterase genes (e.g., LOC100344509 or ACOT1 and LOC100344005 or ACOT4, both downregulated in the H line) are of particular interest because they are related to the biosynthesis of unsaturated FAs and energy metabolism. Multiples ACOTs have been identified in peroxisomes, which catalyse the hydrolysis of acyl-CoAs (including short, medium, long and very long-chain), bile acid-CoAs, and methyl branched-CoAs, to free FA and coenzyme A, regulating their respective intracellular levels (Ishizuka et al., 2004). Notably, members of the ACOTs family are targets for TFs of the PPARs family like PPARα (Hunt et al., 2012). Based on these results, we found that the correlated response to selection in liver gene expression was characterized by a downregulation in the H line of genes involved in FA elongation and degradation, transport, and PPAR signalling pathway, as well as by an upregulation in the same line of genes involved in glutathione and purine/thiopurine metabolism. Functional validation studies would be valuable to confirm the roles of the candidate genes identified and to further elucidate the genetic relationship between IMF content and liver gene expression. Based on our findings, it would be particularly interesting to investigate hepatic and muscle-derived secretory proteins, such as the TFs ETV1, NR4A1, IKZF3, and BACH2, along with the genes highlighted above, which may function independently or in combination with others to regulate lipid metabolism, energy homeostasis, insulin levels, and mitochondrial function.

A correlated response to IMF selection in gene expression in the liver was found. In this study, we identified 308 DE genes in rabbit liver. These include candidate genes associated with fat deposition, transporter activity, and energy metabolism. Well-known genes such as ACAD10, ACBD4, ACOT1, ACOT4, CPT1A, CPT1B, CROT, CYP4B1, ELOVL6, FABP4, GPAT3, HADHB, PCTP, PLIN2, and SLC16A11, as well as less-known genes including BRI3BP, BTN1A1, BTN3A3, IGFBP1, MYLIP, NR4A2, PAQR9, PRRG4, RAB34, RBM4, RNF5, SLC16A6, SLC18A1, SLC25A33, SLC25A42, SLC2A1, and USP18 were identified. All these genes are nominated as putative candidates to modulate the variation in IMF content of the divergent lines. Our analysis identified 11 TFs as putative regulators of gene expression in rabbit liver. Notably, six of these TFs ( BACH2, BHLHE22, ETV1, IKZF3, NFKB2, and NR4A1 ) were also annotated as regulators in the human TFLink database, supporting their potential regulatory roles across species. ETV1, NR4A1, and IKZF3 emerged as particularly relevant regulators of lipid-related genes, exhibiting the highest number of predicted targets in the regulatory network. These targets included several shared candidate genes involved in lipid metabolism, suggesting its potential role for IMF deposition. Functional analysis of gene ontology and pathways indicated an overrepresentation of 40 functional annotation terms, including 31 pathways, 8 biological processes and 1 molecular function, among which the PPAR signalling pathway, FA transport, fatty acid elongation, fatty acid degradation, and biosynthesis of unsaturated FAs were detected. Author contributions A.B. and P.H. planned and coordinated the study, conceptualized and designed the divergent selection experiment for IMF content, and contributed to manuscript revision. P.H. recorded phenotypic data, collected animal samples, and contributed to data analysis. J.V.H. proposed the methodology and structure of the manuscript, performed bioinformatics analyses, analysed gene expression data, conducted gene functional analysis, and wrote the manuscript. M.M.A. and A.Z.G. contributed to animal sample collection, worked on the divergent selection experiment, performed formal data analysis, and assisted to the manuscript writing and revision.

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. MM-A thanks the Spanish Ministry of Science and Innovation for a Ramon y Cajal grant (RYC2021-032618-I) funded by MCIN/AEI/ 10.13039/501100011033 and by European Union NextGenerationEU/PRTR. We want to thank for all people that contributed to the generation of the animal material used in the current study. The author thanks Dr. Carmen Ivorra for the feedback on the 3’ RNA-Seq methodology, Dr. E. Mármol‐Sánchez and Dr. B.S. Sosa-Madrid for their comments. Conflict of interest disclosure The authors declare no conflict of interest. Ethics approval statement All experimental procedures used in this study were reviewed and approved by the Ethical Committee for Experimentation with Animals of the Universitat Politècnica de València, Spain, according to Council Directives 98/58/EC and 2010/63/EU (reference number 2017/VSC/PEA/00212). This study was also performed in accordance with the ARRIVE guidelines (https://arriveguidelines.org/). Data availability statement All relevant datasets produced in this investigation are disclosed in the paper, as well as its supplementary information files. Sequencing data generated and/or analysed during this study are available to the corresponding author upon reasonable request.

Ahamba, I. S. et al. (2024) ‘Unraveling the genetic and epigenetic landscape governing intramuscular fat deposition in rabbits: Insights and implications’, Food Chemistry: Molecular Sciences, 9, p. 100222. Alves-Bezerra, M. and Cohen, D. E. (2018) ‘Triglyceride Metabolism in the Liver’, Comprehensive Physiology, 8(1), pp. 1–22. Anders, S., Pyl, P. T. and Huber, W. (2015) ‘HTSeq—a Python framework to work with high-throughput sequencing data’, Bioinformatics, 31(2), pp. 166–169. Andrews, S. (2010) Babraham bioinformatics-FastQC a quality control tool for high throughput sequence data .Baena, E. et al. (2013) ‘Abstract B157: Understanding the role of transcription factor ETV1 in metabolic reprogramming of prostate cells as a route to novel therapeutics.’, Molecular Cancer Therapeutics, 12(11_Supplement), pp. B157–B157. Benjamini, Y. and Hochberg, Y. (1995) ‘Controlling the false discovery rate: a practical and powerful approach to multiple testing’, Journal of the Royal Statistical Society, 57(1), pp. 289–300. Bindea, G. et al. (2009) ‘ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks’, Bioinformatics, 25(8), pp. 1091–1093. Breiman, L. (2001) ‘Random Forests’, Machine Learning, 45(1), pp. 5–32. Cecchinato, A. et al. (2012) ‘Genetic analysis of beef fatty acid composition predicted by near-infrared spectroscopy’, Journal of Animal Science, 90(2), pp. 429–438. Chen, C. et al. (2024) ‘The interplay between the muscle and liver in the regulation of glucolipid metabolism’, Journal of Molecular Cell Biology, 15(12), p. 73. Cooper, C. D. O. et al. (2015) ‘Structures of the Ets protein DNA-binding domains of transcription factors Etv1, Etv4, Etv5, and Fev: Determinants of DNA binding and redox regulation by disulfide bond formation’, Journal of Biological Chemistry, 290(22), pp. 13692–13709. Dagher, R., Massie, R. and Gentil, B. J. (2021) ‘MTP deficiency caused by HADHB mutations: Pathophysiology and clinical manifestations’, Molecular Genetics and Metabolism, 133(1), pp. 1–7. Damon, M. et al. (2012) ‘Comparison of muscle transcriptome between pigs with divergent meat quality phenotypes identifies genes related to muscle metabolism and structure’, PLoS ONE, 7(3), p. e33763. Dobin, A. et al. (2013) ‘STAR: ultrafast universal RNA-seq aligner’, Bioinformatics, 29(1), pp. 15–21. Furuhashi, M. et al. (2014) ‘Fatty acid-binding protein 4 (FABP4): Pathophysiological insights and potent clinical biomarker of metabolic and cardiovascular diseases’, Clinical Medicine Insights: Cardiology, 2014, pp. 23–33. Gatti, D. et al. (2007) ‘Genome-level analysis of genetic regulation of liver gene expression networks’, Hepatology, 46(2), pp. 548–557. Gondret, F., Mourot, J. and Bonneau, M. (1997) ‘Developmental changes in lipogenic enzymes in muscle compared to liver and extramuscular adipose tissues in the rabbit (Oryctolagus cuniculus)’, Comparative Biochemistry and Physiology - B Biochemistry and Molecular Biology, 117(2), pp. 259–265. Gu, Z., Eils, R. and Schlesner, M. (2016) ‘Complex heatmaps reveal patterns and correlations in multidimensional genomic data’, Bioinformatics, 32(18), pp. 2847–2849. Guo, S. et al. (2023) ‘Metabolic crosstalk between skeletal muscle cells and liver through IRF4-FSTL1 in nonalcoholic steatohepatitis’, Nature Communications 2023 14:1, 14(1), pp. 1–16. Hamill, R. M. et al. (2013) ‘Transcriptome analysis of porcine M. semimembranosus divergent in intramuscular fat as a consequence of dietary protein restriction’, BMC Genomics, 14(1), pp. 1–14. Han, S. W., Jung, B. K. and Ryu, K. Y. (2021) ‘Regulation of polyubiquitin genes to meet cellular ubiquitin requirement’, BMB Reports, 54(4), pp. 189–195. Hocquette, J. F. et al. (2010) ‘Intramuscular fat content in meat-producing animals: Development, genetic and nutritional control, and identification of putative markers’, Animal, 4(2), pp. 303–319. Hong, C. et al. (2010) ‘The E3 ubiquitin ligase IDOL induces the degradation of the low density lipoprotein receptor family members VLDLR and ApoER2’, Journal of Biological Chemistry, 285(26), pp. 19720–19726. Hunt, M. C., Siponen, M. I. and Alexson, S. E. H. (2012) ‘The emerging role of acyl-CoA thioesterases and acyltransferases in regulating peroxisomal lipid metabolism’, Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease, 1822(9), pp. 1397–1410. Huynh-Thu, V. A. et al. (2010) ‘Inferring Regulatory Networks from Expression Data Using Tree-Based Methods’, PLOS ONE, 5(9), p. e12776. Ishizuka, M. et al. (2004) ‘Overexpression of human acyl-CoA thioesterase upregulates peroxisome biogenesis’, Experimental Cell Research, 297(1), pp. 127–141. Jump, D. B. (2002) ‘Dietary polyunsaturated fatty acids and regulation of gene transcription’, Current opinion in lipidology, 13(2), pp. 155–164. Jump, D. B., Tripathy, S. and Depner, C. M. (2013) ‘Fatty acid-regulated transcription factors in the liver’, Annual Review of Nutrition, 33(Volume 33, 2013), pp. 249–269. Laghouaouta, H. et al. (2020) ‘Novel genomic regions associated with intramuscular fatty acid composition in rabbits’, Animals, 10(11), p. 2090. Van der Leij, F. R. et al. (2000) ‘Genomics of the Human Carnitine Acyltransferase Genes’, Molecular Genetics and Metabolism, 71(1–2), pp. 139–153. Li, X. et al. (2007) ‘Role of cdk2 in the sequential phosphorylation/activation of C/EBPβ during adipocyte differentiation’, Proceedings of the National Academy of Sciences of the United States of America, 104(28), pp. 11597–11602. Liska, O. et al. (2022) ‘TFLink: an integrated gateway to access transcription factor–target gene interactions for multiple species’, Database, 2022, p. baac083. Liu, L. et al. (2019) ‘Effect of Divergent Selection for Intramuscular Fat Content on Muscle Lipid Metabolism in Chickens’, Animals, 10(1), p. 4. Loix, M. et al. (2024) ‘The ubiquitous role of ubiquitination in lipid metabolism’, Trends in Cell Biology, 34(5), pp. 416–429. López-Alcalá, J. et al. (2024) ‘Localization, traffic and function of Rab34 in adipocyte lipid and endocrine functions’, Journal of Biomedical Science, 31(1), pp. 1–28. Love, M. I., Huber, W. and Anders, S. (2014) ‘Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2’, Genome Biology, 15(12), pp. 1–21. Luo, G. et al. (2023) ‘The Characterization and Differential Analysis of m6A Methylation in Hycole Rabbit Muscle and Adipose Tissue and Prediction of Regulatory Mechanism about Intramuscular Fat’, Animals, 13(3), p. 446. Ma, F. et al. (2019) ‘A comparison between whole transcript and 3’ RNA sequencing methods using Kapa and Lexogen library preparation methods 06 Biological Sciences 0604 Genetics’, BMC Genomics, 20(1), pp. 1–12. Maples, J. M. et al. (2015) ‘Differential epigenetic and transcriptional response of the skeletal muscle carnitine palmitoyltransferase 1B (CPT1B) gene to lipid exposure with obesity’, American Journal of Physiology - Endocrinology and Metabolism, 309(4), pp. E345–E356. Mármol-Sánchez, E. et al. (2022) ‘Modeling microRNA-driven post-transcriptional regulation using exon–intron split analysis in pigs’, Animal Genetics, 53(5), pp. 613–626. Martin, M. (2011) ‘Cutadapt removes adapter sequences from high-throughput sequencing reads’, EMBnet.journal, 17(1), pp. 10–12.Martinez-Alvaro, M. et al. (2017) ‘Muscle lipid metabolism in two rabbit lines divergently selected for intramuscular fat.’, Journal of animal science, 95(6), pp. 2576–2584. Martínez-Álvaro, Marina et al. (2018a) ‘Correlated responses to selection for intramuscular fat in several muscles in rabbits’, Meat Science, 139, pp. 187–191. Martínez-Álvaro, M. et al. (2018b) ‘Liver metabolism traits in two rabbit lines divergently selected for intramuscular fat’, Animal, 12(6), pp. 1217–1223. Martínez-Álvaro, M., Hernández, P. and Blasco, A. (2016) ‘Divergent selection on intramuscular fat in rabbits: Responses to selection and genetic parameters’, Journal of Animal Science, 94(12), pp. 4993–5003. Matoulkova, E. et al. (2012) ‘The role of the 3’ untranslated region in post-transcriptional regulation of protein expression in mammalian cells.’, RNA Biology, 9(5), pp. 563–576. Matsuzaka, T. and Shimano, H. (2009) ‘Elovl6: A new player in fatty acid metabolism and insulin sensitivity’, Journal of Molecular Medicine, 87(4), pp. 379–384. Mayr, C. (2019) ‘What Are 3′ UTRs Doing?’, Cold Spring Harbor Perspectives in Biology, 11(10), p. a034728. McClure, R. S. et al. (2023) ‘3′ RNA-seq is superior to standard RNA-seq in cases of sparse data but inferior at identifying toxicity pathways in a model organism’, Frontiers in Bioinformatics, 3, p. 1234218.Migdał et al. (2018) ‘Single nucleotide polymorphisms within rabbits (Oryctolagus cuniculus) fatty acids binding protein 4 (FABP4) are associated with meat quality traits’, Livestock Science, 210, pp. 21–24. Natacha Pena, R. et al. (2016) ‘Genetic marker discovery in complex traits: A field example on fat content and composition in pigs’, International Journal of Molecular Sciences, 17(12), p. 2100. Nguyen, D. V., Nguyen, O. C. and Malau-Aduli, A. (2021) ‘Main regulatory factors of marbling level in beef cattle’, Veterinary and Animal Science, 14, p. 100219. Nath, A., & Chan, C. (2016). Genetic alterations in fatty acid transport and metabolism genes are associated with metastatic progression and poor prognosis of human cancers. Scientific reports, 6(1), p. 18669.Óvilo, C. et al. (2014) ‘Longissimus dorsi transcriptome analysis of purebred and crossbred Iberian pigs differing in muscle characteristics’, BMC Genomics, 15(1), pp. 1–24. Pan, J. et al. (2021) ‘Insulin-like growth factor binding protein 1 ameliorates lipid accumulation and inflammation in nonalcoholic fatty liver disease’, Journal of Gastroenterology and Hepatology, 36(12), pp. 3438–3447. Paniri, A. et al. (2024) ‘Genetic variations in IKZF3, LET7-a2, and CDKN2B-AS1: Exploring associations with metabolic syndrome susceptibility and clinical manifestations’, Journal of Clinical Laboratory Analysis, 38(1–2), p. e24999. Poulsen, L. la C., Siersbæk, M. and Mandrup, S. (2012) ‘PPARs: Fatty acid sensors controlling metabolism’, Seminars in Cell & Developmental Biology, 23(6), pp. 631–639. Sanford, J. D. et al. (2023) ‘Carnitine o-octanoyltransferase is a p53 target that promotes oxidative metabolism and cell survival following nutrient starvation’, Journal of Biological Chemistry, 299(7), p. 1. Schlaepfer, I. R. and Joshi, M. (2020) ‘CPT1A-mediated Fat Oxidation, Mechanisms, and Therapeutic Potential’, Endocrinology, 161(2), p. bqz046. Schrem, H., Klempnauer, J. and Borlak, J. (2004) ‘Liver-Enriched Transcription Factors in Liver Function and Development. Part II: the C/EBPs and D Site-Binding Protein in Cell Cycle Control, Carcinogenesis, Circadian Gene Regulation, Liver Regeneration, Apoptosis, and Liver-Specific Gene Regulation’, Pharmacological Reviews, 56(2), pp. 291–330. Schuster, S. L. and Hsieh, A. C. (2019) ‘The Untranslated Regions of mRNAs in Cancer’, Trends in Cancer, 5(4), pp. 245–262. Schwab, C. R. et al. (2009) ‘Results from six generations of selection for intramuscular fat in Duroc swine using real-time ultrasound. I. Direct and correlated phenotypic responses to selection’, Journal of Animal Science, 87(9), pp. 2774–2780. Shannon, P. et al. (2003) ‘Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks’, Genome Research, 13(11), pp. 2498–2504. Shen, W. K. et al. (2023) ‘AnimalTFDB 4.0: a comprehensive animal transcription factor database updated with variation and expression annotations’, Nucleic Acids Research, 51(D1), pp. D39–D45. Smith, T., Heger, A. and Sudbery, I. (2017) ‘UMI-tools: Modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy’, Genome Research, 27(3), pp. 491–499. Sosa-Madrid, B. S. et al. (2020a) ‘The effect of divergent selection for intramuscular fat on the domestic rabbit genome’, Animal, 14(11), pp. 2225–2235. Sosa-Madrid, B. S. et al. (2020b) ‘Genomic regions influencing intramuscular fat in divergently selected rabbit lines’, Animal Genetics, 51(1), pp. 58–69.Sosa-Madrid, B. S. et al. (2025). 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, 57(1), p.36.Stefan, N., Kantartzis, K. and Häring, H. U. (2008) ‘Causes and Metabolic Consequences of Fatty Liver’, Endocrine Reviews, 29(7), pp. 939–960. Tao, W., Radstake, T. R. D. J. and Pandit, A. (2022) ‘RegEnrich gene regulator enrichment analysis reveals a key role of the ETS transcription factor family in interferon signaling’, Communications Biology 2022 5:1, 5(1), pp. 1–12. R Core Team. (2024) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing . Vienna, Austria.Tailleux, A., Wouters, K., & Staels, B. (2012). Roles of PPARs in NAFLD: potential therapeutic targets. Biochimica et Biophysica Acta (BBA)-Molecular and Cell Biology of Lipids, 1821(5), p. 809-818.Valdés-Hernández, J. et al. (2023) ‘Global analysis of the association between pig muscle fatty acid composition and gene expression using RNA-Seq’, Scientific Reports 2023 13:1, 13(1), pp. 1–12. Valdés-Hernández, J. et al. (2024) ‘Identification of candidate regulatory genes for intramuscular fatty acid composition in pigs by transcriptome analysis’, Genetics Selection Evolution, 56(1), pp. 1–18. Wagner, N. and Wagner, K. D. (2020) ‘The Role of PPARs in Disease’, Cells, 9(11), p. 2367. Wang, X. et al. (2017) ‘The comprehensive liver transcriptome of two cattle breeds with different intramuscular fat content’, Biochemical and Biophysical Research Communications, 490(3), pp. 1018–1025. Wang, H. et al. (2025). Transcriptome analysis reveals key genes and regulatory networks underlying intramuscular fat deposition in rabbits. BMC genomics, 26(1), p. 785.Weissglas-Volkov, D. et al. (2011) ‘The N342S MYLIP polymorphism is associated with high total cholesterol and increased LDL receptor degradation in humans’, The Journal of Clinical Investigation, 121(8), pp. 3062–3071. Wood, J. D. et al. (2008) ‘Fat deposition, fatty acid composition and meat quality: A review’, Meat Science, 78(4), pp. 343–358. Xu, S. et al. (2019) ‘Perilipin 2 and lipid droplets provide reciprocal stabilization’, Biophysics Reports 2019 5:3, 5(3), pp. 145–160. Yi, W. et al. (2023) ‘Lipo-nutritional quality of pork: The lipid composition, regulation, and molecular mechanisms of fatty acid deposition’, Animal Nutrition, 13, pp. 373–385. Zhang, L. et al. (2018) ‘The Orphan Nuclear Receptor 4A1: A Potential New Therapeutic Target for Metabolic Diseases’, Journal of Diabetes Research, 2018(1), p. 9363461. Zomeño, C., Blasco, A. and Hernández, P. (2013) ‘Divergent selection for intramuscular fat content in rabbits. II. Correlated responses on carcass and meat quality traits.’, Journal of animal science, 91(9), pp. 4532–4539. Zubiri-Gaitán, A. et al. (2022) ‘Intramuscular fat selection in rabbits modifies the fatty acid composition of muscle and liver tissues’, Animals, 12(7), pp. 1–12. Zubiri-Gaitán, A., Blasco, A. and Hernández, P. (2023) ‘Plasma metabolomic profiling in two rabbit lines divergently selected for intramuscular fat content’, Communications Biology, 6(1), p. 893. Figure legends Figure 1 . Liver transcriptomic analysis. (a) Volcano plot showing differentially expressed genes in liver of high and low lines for IMF content. Each dot represents one gene (n = 11135). Dashed lines indicate the values of the thresholds. On the y-axis the log 2 fold change (FC) values with ± 0.58 and limits of −10 and 10 were represented. On the x-axis the significance level -log 10 of the P -value equal to 0.05 (absolute value of 1.3) and limits of 0 and 80 were represented. In the significance legend the dots corresponding to the labels Down, Up and Unchanged indicate respectively: genes with decreased expression (or downregulated), genes with increased expression (or upregulated), and genes that did not pass the threshold (|log 2 FC| ≥ 0.58 and a BH padj ≤ 0.05). (b) Complementary heatmap of differentially expressed genes using their normalized expression levels and then escalated. Legends on the left indicate: Gene expression values range from -4 to 4 and the colour scale: shades of blue represent negative values, shades of red represent positive values, and neutral shades (white) indicate values close to zero. L (n = 24) = low line samples and H (n = 24) = high line samples. Figure 2. Summary of 36 regulators with an absolute FC value of 1.2 and BH adjusted P-value < 0.05. The scores and ranking summarize the importance of regulators by family and delineate those members that were differentially expressed with an absolute FC value of 1.5 and BH adjusted P-value < 0.05 (star symbol within the bars). Figure 3. Gene regulatory network inference with target prediction. Yellow and red circles indicate 11 transcription factors, and 210 putative predicted targets interconnected by the weights calculated in the construction of the network. These 222 genes were corroborated with the results of differential expression analysis and regulator enrichment analysis. Figure 4. Functional analysis of the differentially expressed genes between high (H) and low (L) lines for IMF content. (a) Bar graph representing the 40 functional annotation terms identified by ClueGO plug-in, including 31 pathways, 8 biological processes and 1 molecular function with their respective number of associated genes found. The color of the bars was coded according to the significance of the BH corrected P-value < 0.05. (b) Network built with the functional terms detected by ClueGO plug-in. The 12 different functional groups encoded by color correspond to the classification based on GO term similarities. Supplementary Material File (figure 3.tiff) - Download - 45.89 MB Information & Authors Information Version history Copyright This work is licensed under a Non Exclusive No Reuse License.

Authors Metrics & Citations Metrics Article Usage 415views 172downloads Citations Download citation Jesús Valdés-Hernández, Agostina Zubiri-Gaitán, Marina Martínez-Álvaro, et al. Transcriptomic analysis in the liver of two rabbit lines divergently selected for intramuscular fat content. Authorea. 17 October 2025. DOI: https://doi.org/10.22541/au.176069730.09333094/v1 DOI: https://doi.org/10.22541/au.176069730.09333094/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu.