A meta-analysis approach for annotation and identification of lncRNAs controlling perirenal fat deposition in suckling lambs

Abstract

Long non-coding RNAs (lncRNAs) are being studied in farm animals due to their association with traits of economic interest, such as fat deposition. Based on the analysis of perirenal fat transcriptomes, this research explored the relevance of these regulatory elements to fat deposition in suckling lambs. To that end, meta-analysis techniques have been implemented to efficiently characterize and detect differentially expressed transcripts from two different RNA-seq datasets, one including samples of two sheep breeds that differ in fat deposition features, Churra and Assaf (n = 14), and one generated from Assaf suckling lambs with different fat deposition levels (n = 8). The joint analysis of the 22 perirenal fat RNA-seq samples with the FEELnc software allowed the detection of 3953 novel lncRNAs. After the meta-analysis, 251 differentially expressed genes were identified, 21 of which were novel lncRNAs. Additionally, a co-expression analysis revealed that, in suckling lambs, lncRNAs may play a role in controlling angiogenesis and thermogenesis, processes highlighted in relation to high and low fat deposition levels, respectively. Overall, while providing information that could be applied for the improvement of suckling lamb carcass traits, this study offers insights into the biology of perirenal fat deposition regulation in mammals.

Keywords:

• Non-coding RNAs
• adipogenesis
• RNA-Seq
• meta-analysis
• sheep

Introduction

In recent years, deep sequencing of transcriptomes has revealed the importance of new types of non-coding transcripts, such as long non-coding RNAs (lncRNAs), as regulators of gene expression.¹ LncRNAs are defined as RNA molecules longer than 200 nucleotides with low-coding to no coding potential. These RNAs are tissue-specific and have primary sequences that are less conserved across species than mRNAs.² LncRNAs control gene expression at a supragenomic layer that integrates genomic and epigenomic levels of gene control. In the nucleus, these molecules are involved in chromatin organization, direct transcription, and the modification of nascent RNAs, while in the cytoplasm, they modulate mRNA turnover, storage, and protein translation.^3–5

Growing catalogues of lncRNAs are being characterized⁶ for farm animals, driven by increasing examples of association between lncRNAs and phenotypes of interest. In this sense, adipose tissue, due to its importance in energy storage and other endocrine functions that help in the regulation of energy balance, appetite, immunity and reproduction, has been one of the key organs used to characterize lncRNA regulation in animals.⁷

In sheep, several RNA-Seq studies have focused on comparisons between thin-tailed and fat-tailed sheep breeds to characterize the regulatory role of lncRNAs in fat deposition.^8–11 These studies demonstrated that lncRNAs play important roles in the biosynthesis of fatty acids in the tail and subcutaneous fat of adult sheep.^8,9 However, less research has been conducted on how lncRNAs can regulate fat deposition in lambs. Two previous studies from our research group focused on characterizing the perirenal fat transcriptome of suckling lambs at commercial slaughter weights (10–12 kg).^12,13 The choice of this reservoir is based on the fact that it represents a high proportion of the total fat in lamb during preweaning, and it is used as an indicator of carcass quality.¹⁴

According to several studies, a major factor influencing the carcass tissue composition of suckling lambs is the breed.^15–20 Thus, the first transcriptomic study on perirenal fat from our group compared samples from suckling lambs of two different breeds, Churra and Assaf. These two breeds were selected because at the same standard commercial slaughter weight (11 kg of live weight), they differ in their fat carcass content, with Churra carcasses showing a greater cavity and intermuscular leg fat proportion than Assaf carcasses.¹⁷ In the second mentioned study, Alonso-García et al.¹² analyzed perirenal transcriptome samples from a single breed (Assaf), comparing animals that, at the same slaughter weight, had different proportions of perirenal fat in the half-carcass.

Hence, considering the availability of the two datasets that focus on the study of fat deposition biological processes in suckling lambs, we hypothesize that there may be common regulatory mechanisms associated with different levels of fat deposition in suckling lambs, both between and within breeds. As previously shown in adult sheep for general fat deposition traits^8,9,21, some of these regulatory mechanisms may be related to lncRNAs influencing gene expression regulation.

The lncRNA catalogues in farming animals remain incomplete compared to the mouse or human annotations.^6,22 A common challenge to characterize lncRNAs resides in the limited number of samples analyzed in independent RNA-Seq studies. This is why it has been proposed to utilize expression datasets derived from multiple independent but related studies to accurately annotate and functionally characterize lncRNAs. In sheep, lncRNAs have been characterized across a multitissue dataset²³; however, there is only one study using meta-analysis procedures for that purpose.²⁴ Integration of RNA-Seq datasets allows to better determine transcripts and biological processes that could have been overlooked in individual experiments.²⁴ In addition, meta-analysis techniques have the advantage of increasing the available sample size by integrating related datasets and subsequently increasing the power to detect differentially expressed transcripts.^25,26

Based on the complementary experimental designs of the two previously cited RNA-Seq studies from our research group related to suckling lamb fat deposition^12,13, the objectives of this study are (1) to annotate and characterize the lncRNAs present in the perirenal fat transcriptome of Churra and Assaf suckling lambs using both RNA-Seq datasets and (2) to determine, through a meta-analysis, lncRNA regulatory elements involved in the control of fat deposition in suckling lambs that are commonly shared within breeds at the individual level and between breeds.

Materials and methods

Data collection

RNA-sequencing (RNA-Seq) datasets from the perirenal fat depot of suckling lambs were retrieved from two studies previously published by our own research group. Hence, we obtained the raw sequences related to the study of Suárez-Vega et al.¹³ from the Gene Expression Omnibus (GEO) repository (GSE108314 accession number) and the raw sequences analyzed by Alonso-García et al.¹² from the BioStudies repository (E-MTAB-11839 accession number). In the study performed by Suárez-Vega et al.¹³ the perirenal fat transcriptomes of male suckling lambs from the Churra (n = 6) and Assaf (n = 8) breeds were compared, whereas Alonso-García et al.¹² analyzed perirenal fat transcriptomes of Assaf male suckling lambs (n = 8). Therefore, RNA-Seq datasets used in the present study were obtained from these two Spanish sheep breeds, Churra and Spanish Assaf (Fig. 1). The transcriptomic information of these 14 samples, which will be referred to as DataSet 1 from this point onwards, was analyzed and discussed in relation to the percentage of perirenal fat in the carcass (PF), which was greater in Churra than in Assaf.^13,17 On the other hand, the study reported by Alonso-García et al.¹² contrasted the perirenal fat transcriptome of eight male Assaf lambs selected from a larger population (n = 17) based on their divergent phenotypic values for the proportion of kidney knob and channel fat in the half carcass (KKCF) [High-KKCF (n = 4) and Low-KKCF (n = 4)]. This second dataset of eight transcriptome samples will be referred from this point onwards as DataSet 2.

Figure 1. Breeds analyzed in the in the study: (A) Spanish churra sheep; (B) Spanish assaf sheep.

For the two datasets considered, the RNA sources were perirenal fat samples collected at the slaughterhouse when lambs reached the commercial weight to slaughter (between 10 and 12 kg). After being kept in an RNA-stabilization solution (Ambion RNAlater; Life Technologies), the individual perirenal fat samples were subjected to an RNA extraction protocol using the miRNeasy Mini Kit (Qiagen) for both datasets. The sequencing of the samples is described in detail in the corresponding studies.^12,13 In summary, the samples for DataSet 1 were sequenced on an Illumina HiSeq 2000 Sequencer (Fasteris SA) to a minimum depth of 30 million paired-end reads, producing stranded paired-end reads of 75 bp. On the other hand, DataSet 2 cDNA libraries were sequenced on a NovaSeq 6000 system (Illumina) to a minimum depth of 30 million paired-end reads, producing stranded paired-end reads of 150 bp.

The carcasses of all the animals included in the present study were subjected to phenotypic characterization following the methods of Colomer Rocher et al.²⁶ which have been described in detail in the corresponding studies.^12,13 Among the traits related to fat deposition, the percentage of perirenal fat in the carcass (PF%) was measured in two studies, and this has been the phenotype considered in this group as an indicator of fat deposition level (Supplementary Table 1).

Bioinformatic analysis

The step-by-step bioinformatics workflow applied in this study to 22 perirenal fat transcriptome samples for the detection of lncRNAs is depicted in Fig. 2. A detailed description of the main stages of this workflow is provided below.

Figure 2. Bioinformatic workflow applied in this work for lncRNA annotation and meta-analysis. The different colors represent the bioinformatic analyses steps applied to the different RNA-seq datasets considered: blue for DataSet 1 (14 RNA-seq samples)¹³ green for DataSet 2 (8 RNA-seq samples;¹² and yellow for both RNA-seq datasets (22 samples). The annotation files are indicated in pink color.

Data alignment

The raw sequences obtained from the databases were aligned using the STAR software (v.2.7.6a)²⁷ against the latest version of the ovine reference genome available at the NCBI (ARS-UI_Ram_v2.0). We used the two-pass alignment method for more sensitive novel junction discovery.²⁸

For the first pass, the index was obtained by changing the –sjdbOverhang parameter depending on the length of the reads (75 or 150 bp, depending on the dataset). The same index was also used for the second-pass alignment. In the 1st pass, the novel junctions were identified and included in the genome indices. For this 1st-pass alignment, we used the STAR options –outFilterType BySJout to reduce the number of ‘spurious’ junctions and retain only those reads, including junctions that passed through the SJ.out.tab files, and –outSAMtype BAM Unsorted SortedByCoordinate to obtain both unsorted and sorted bam files as outputs. In the 2nd pass, all reads were mapped again by considering the original GTF file with the –sjdbGTFfile parameter plus the novel junctions detected in the 1st pass mapping using the option –sjdbFileChrStartEnd. This last option allows the use of the list of annotated junctions (SJ.out.tab files) obtained in the previous mapping step.

After the mapping step, we used the StringTie (v.2.1.4)²⁹ transcriptome assembler to reconstruct the transcriptome of each individual sample using the mapped BAM files from the 2nd pass mapping, and the assembly was guided (option –rf) using the ovine reference genome ARS-UI_Ram_v2.0. Then, individual transcriptome annotation files (gtf files) from each sample were merged using Cuffmerge³⁰ to obtain a nonredundant reference annotation file (Oar_1.merged.gtf in Fig. 2).

lncRNA detection, annotation and transcript quantification

Using the new reference annotation file (Oar_1.merged.gtf), we identified potential novel lncRNAs using FEELnc (v.0.2).³¹ The detection of lncRNAs using the FEELnc software requires the use of three interdependent modules. The first module (FEELnc_filter) filters those transcripts found in the analyzed datasets that are longer than 200 nucleotides and overlap protein-coding genes (by using -b transcript_biotype = mRNA). Then, the FEELnc_codpot module was used to identify the potential lncRNAs among the reconstructed transcripts. This module applies random forest to compute the coding potential score (CPS) of the transcripts. For this approach, FEELnc_codpot was used to predict ORF coverage, the multi k-mer frequencies and RNA sequencing length. These predictor scores are incorporated into the Random Forest analysis that computes a CPS for each input training transcripts. Then, the software sets an optimal CPS cutoff based on a threshold that maximized the sensitivity (Sn) and specificity (Sp) using 10-fold cross-validation. In our study, using the shuffle option, we fixed the Sn and Sp thresholds to 0.99. The shuffle strategy is based on the idea that lncRNAs are produced from the ‘debris’ of protein-coding genes.^32,33 To achieve this, FEELnc shuffled mRNA sequences from the reference annotation while preserving the given k-mer frequency (–kmer = 1,2,3,6,9,12) of the input sequences. The annotation file of the lncRNAs detected in this step (Oar_2.lncRNA.gtf in Fig. 2) was merged with the reference annotation file using Cuffcompare, resulting in a third annotation file (Oar_3.merged.gtf in Fig. 2). For the third module (FEELnc_classifier), the novel lncRNAs detected in the previous step were classified into different types (genic, intergenic) and subtypes (containing, nested, overlapping, divergent, convergent, same strand), depending on their closest annotated transcripts in the reference annotation file (ARS-UI_Ram_v2.0). This module aims to identify potential cis-interactions between lncRNAs and target genes.

The Coding Potential Assessment Tool (CPAT) (v.3.0.4)³⁴ was used to confirm that lncRNAs detected by FEELnc had no coding potential and indeed are non-coding transcripts. This software distinguishes between coding and non-coding transcripts from a large pool of candidates using a logistic regression model constructed with four sequence features: open reading frame (ORF) size, ORF coverage, Fickett TESTCODE statistic, and hexamer use bias. First, CPAT was trained with available known protein-coding transcripts (759,400 transcripts annotated as ‘protein_coding’) and non-coding sequences (7,880 transcripts annotated as ‘ncRNAs’) from ARS-UI_Ram_v2.0 ovine genome to build a logistic regression model used to obtain the coding potential cutoff. Then, using 10-fold cross-validation, nonparametric two-graph ROC curves calculated with the ROCR R package (v.1.0.11)³⁵ were used to determine an optimum cutoff value threshold that maximizes the discriminatory power using maximum Sn and Sp (0.975; Supplementary Fig. 1). After obtaining the reference coding potential cutoff, CPAT calculates the CPS for each novel lncRNA detected by FEELnc. Those lncRNAs that had a CPS calculated by CPAT higher than the cutoff were considered coding transcripts and filtered out from the next steps.

Finally, the annotation file containing the lncRNAs identified after the FEELnc_codpot step was merged with the reference annotation file from NCBI (ARS-UI_Ram_v2.0) using Cuffcompare to obtain a complete annotation file that contains the transcripts annotated in the ovine reference genome plus the potential lncRNAs detected in this study (Oar_3.merged.gtf in Fig. 2). This final annotation file (Oar_3.merged.gtf) was subsequently used for mapping and quantification of the total transcripts. Thus, following the scheme represented in Fig. 2, the samples were remapped using STAR (this time using the default one-step mapping method), and the quantification of transcript expression was performed with RSEM (v.1.3.3)³⁶ by implementing the following options: –estimate-rspd, to know if there is any sequencing bias in the data enabling RSEM to learn from the data how the reads are distributed across a transcript; –calc-ci, for calculating 95% credibility intervals and posterior mean estimates; and –forward-prob, set to 0, for a strand-specific protocol where all upstream reads are derived from the reverse strand.

Meta-analysis of differentially expressed genes

Previously, for the meta-analysis, a phenotypic evaluation of the PF% values of the samples in the two datasets was performed by using logit transformation (logPF). For this transformed trait, we performed independent analyses of variance (ANOVA) using the aov function from the R stats package (v. 4.3.1)³⁷ to determine whether logPF was significantly different between the following comparisons: (1) Churra vs. Assaf samples from DataSet 1; (2) Assaf Low-KKCF vs. Assaf High-KKCF from Dataset 2; (3) Assaf samples from DataSet 1 vs. Low-KKCF samples from DataSet 2, which included the two groups of samples with lower fat deposition levels (LFD group, n = 12); (4) Churra breed samples from DataSet 1 and the High-KKCF samples from DataSet 2, which included the samples from each dataset with higher fat deposition levels (HFD group, n = 10); and (5) a final comparison contrasting the samples defining the LFD group (Assaf samples from DataSet 1 vs. Low-KKCF from DataSet 2) vs. the samples defining the HFD group (Churra samples from Dataset 1 vs. High-KKCF from Dataset 2).

We used the ssizeRNA R package (v.1.3.2)³⁸ to estimate the statistical power expected for each individual dataset based on its sample size. To that aim, we used the function ssizeRNA_vary. Considering a two-group contrast RNA-seq experiment, this function calculates appropriate sample sizes for a desired power where the mean and dispersion of gene expression vary among genes. This analysis used gene expression means and dispersions obtained from a suckling lambs perirenal fat transcriptome reference dataset (the four High-KKCF lambs from DataSet 2). The power calculations according to sample size was performed using the following parameters: a desired power threshold to be achieved (power) = 0.8; false discovery rate (fdr) = 0.05; proportion of non-differentially expressed genes (pi0) = 0.8; total number of genes (nGenes) = 18,733; mean counts in the control group (m) = estimated using the considered reference dataset, dispersion (disp) = estimated using the edgeR R package (v.4.0.16)³⁹ for the reference dataset; and three Fold Change scenarios (FCh = 2, 3 and 4).

Once the phenotypic traits were evaluated and the power of the individual analyses was calculated, we performed a differential expression analysis (DEA) using the RNA-Seq data for each of the considered datasets individually. The transcript quantification files obtained from RSEM were imported to the R environment using the tximport R package,⁴⁰ and the DEA for each dataset was performed using the R package DESeq2.⁴¹ The model design used for the DEA was the same as that used in the corresponding manuscript for analyzing each dataset. Hence, for DataSet1, we compared samples clustered based on their breed, comparing the fat transcriptome between Churra and Assaf suckling lambs. To maintain the consistency and coherence of the results in the original study, the factor age was not included in the model due to a lack of uniformity in the age distribution between the two breeds included in DataSet 1.¹³ On the other hand, for DataSet2, we compared the transcriptome between High- and Low-KKCF animals, taking into consideration the age of the animals as a categorical variable. Age was not included in the model for DataSet 1, as age and breed were confounding variables. More details on the models used for the DEA are provided in the original manuscript.^12,13

Then, using the normalized counts, the raw p-values, the P adjusted value (P_adj; using Benjamini and Hochberg⁴² multiple test correction) and the log₂FoldChange (log₂FCh) obtained from DESeq2 for each dataset, a meta-analysis was performed using the R package metaRNASeq.⁴³ For this meta-analysis, the results of the DEA were grouped according to the previous definitions of the LFD and HFD sample groups. The metaRNAseq package uses two different approaches for the detection of new differentially expressed (DE) genes: Fisher’s method and the inverse normal combination technique (INCT). The genes used for subsequent steps, named keepDE, were those genes identified as DE based on both methods (P_adj < 0.05), and the log₂FCh sign (+/−) were consistent between both datasets.

After the meta-analysis, Gene Ontology (GO) overrepresentation analyses were performed with the keepDE genes using the clusterProfiler R package⁴⁴ and the org.Hs.eg.db R package as a database for GO terms.⁴⁵ The functional analyses were performed independently for the keepDE genes overexpressed in each of the sample groups, i.e., HFD and LFD, and for each of the GO categories (Biological Process (BP), Molecular Function (MF), and Cellular Component (CC)). GO terms were considered enriched for a q-value cutoff of 0.05 after applying the multiple test correction defined by Benjamini & Hochberg.⁴²

Detection of co-expressed gene modules

To determine potential trans-interactions between lncRNAs and mRNAs, the R package CEMiTool⁴⁶ was used to identify differentially co-expressed modules of genes. This package uses a weighted gene co-expression network analysis (WGCNA) approach to obtain the co-expressed modules and genes with the highest connectivity, known as hub genes, within each module. Previous to CEMiTool analysis, transcripts were filtered based on their expression using the filter_RNA_seq function from the GWENA R package⁴⁷ to retain genes with at least 10 counts in one sample. After filtering, for the co-expression analysis, we used a total of 16,849 remaining genes. Prior to gene filtering by CEMiTool, we applied variance stabilizing transformation (VST; apply_vst = T). Spearman’s method was used to calculate the correlation coefficients, considering a soft thresholding power (β) of 3.

Once the gene co-expressed modules were identified, CEMiTool was used for gene set enrichment analysis (GSEA), from which a p value and a normalized enrichment score (NES) were obtained for each module. NES was calculated using the genes within a module as the gene set and the median z-score values of each phenotype as the rank. This parameter indicates the enrichment score for a module within each group and is normalized by the number of genes included in the module. For each NES, a Benjamini–Hochberg adjusted p-value was calculated, which indicates whether there were significant differences (p-value <0.05) between the expression of the co-expressed genes for each module between the groups analyzed. The co-expressed gene modules detected in this step were filtered to analyze in depth those modules with co-expressed keepDE-lncRNAs detected in the meta-analysis. This information was assessed to gain further insights into the possible relationships between these regulatory elements and the biological processes underlying perirenal fat deposition in suckling lambs.

Functional analysis and QTL annotation

The genes from co-expressed modules containing keepDE-lncRNAs were used for functional analyses. We performed enrichment analysis using the clusterProfiler R package by applying the same parameters described previously. From the enrichment analysis implemented with clusterProfiler and using the the cnetplot function, we obtained the relationships between the enriched terms and the co-expressed genes in these modules. The results of this analysis were combined with the FEELnc classifier information to identify cis-associations between co-expressed genes and new lncRNAs. In addition, the Pearson correlations between the expression of mRNAs and lncRNAs were calculated using the WGCNA R package⁴⁸ to detect possible trans-associations. To elucidate the biological role of the DElncRNAs detected in the co-expressed gene modules, gene-term networks based on all the information described previously were created using Cytoscape.⁴⁹

The DElncRNAs were also subjected to QTL annotation analysis. For this purpose, we used the GALLO R package⁵⁰ to exploit the information annotated in the AnimalQTLdb,⁵¹ considering the Ovis aries Rambouillet v.2.0 reference genome. We performed QTL annotation using an interval of 5 × 10⁵ bp for all the DElncRNAs detected by the meta-analysis.

Results

Alignment statistics

Approximately 44 million paired-end reads were generated per sample (44 347 682.63 ± 5 738 600.09, mean ± SD) for the samples obtained from DataSet1.¹³ On average, 36 million (36 052 414.5 ± 4 728 699.86) mapped uniquely to the ovine genome assembly (ARS-UI_Ram_v.2.0), corresponding to 81.3% of the total reads. For DataSet2,¹² 31 million paired-end reads were generated per sample (31 661 520.25 ± 1 682 935.43), and an average of almost 30 million reads (29 714 380.5 ± 1 585 934.14) mapped uniquely to the same ovine genome assembly. Those uniquely mapped reads accounted for 93.8% of the total reads.

Transcriptomic identification and characterization of novel lncRNAs in the perirenal fat transcriptome of suckling lambs

FEELnc detected 3,953 novel lncRNAs not previously annotated in the ARS-UI_Ram_v.2.0 reference sheep genome, with a CPS lower than 0.45 (Supplementary Fig. 1). Approximately 11% of the novel lncRNAs detected were located in OAR1, OAR2 and OAR3. Similarly, this distribution was also found for mRNAs and lncRNAs already annotated in the reference genome (Fig. 3A). Based on the analyses performed with CPAT, a total of 40 from the 3,953 novel lncRNAs detected had CPS higher than the coding potential cutoff (0.975) obtained after the 10-fold cross-validation. These 40 lncRNAs (1.01% of the total of novel lncRNAs) were eliminated from further analysis, continuing the following analyses with 3913 novel lncRNAs.

Figure 3. Long non-coding RNA (lncRNA) characterization. The colors used to represent the different transcript categories (green for protein-coding genes, orange for annotated lncRNAs in the sheep reference genome and purple for novel lncRNAs) are the same for A, B, C and D. (A) Bar graph showing the number of genes (Y-axis) per chromosome in the 26 sheep autosomes (X-axis). For better visualization of the bar graph, a square root transformation has been used on the Y-axis. (B) Density plot showing the length of the genes. The X-axis represents the transcript length in base pairs (bp) and the Y-axis the number of transcripts. (C) Density plot showing the number of exons (X-axis) per transcript in each of the categories of the transcripts considered (Y-axis). (D) Expression levels measured in the logarithm of fragments per kilobase million (Log10(FPKMs); Y-axis) of the three transcript categories. (E) Classification of novel lncRNAs into classes by genomic location according to the sheep reference genome.

Compared to those lncRNAs annotated in the sheep reference genome (4,442), the average transcript lengths were similar (2,560.5 and 2,406.8 bp for novel and annotated lncRNAs, respectively), whereas mRNAs annotated in the reference genome were larger (3,841.67 bp on average) (Fig. 3B). In terms of the number of exons, both the novel lncRNAs and those annotated in the reference genome had between two and 40 exons, with ∼46% having fewer than five exons (2.73 exons on average for novel lncRNAs and 2.94 exons for annotated lncRNAs on average), while protein-coding genes had 9.93 exons on average (Fig. 3C).

Interestingly, the novel lncRNAs were slightly more highly expressed than the lncRNAs annotated in the reference genome. Considering expressed transcripts, the median expression for novel lncRNAs was 0.91 FPKMs, and the median expression for annotated lncRNAs was 0.13 FPKMs. However, both the novel and annotated lncRNAs still exhibited lower expression than did the protein-coding genes (median normalized counts = 2.04 FPKMs) (Fig. 3D). Of the 3,953 novel lncRNAs detected in this study, 909 had an average expression for all samples higher than 10 normalized counts. FEELnc_classiffier classified the lncRNAs depending on their location with regard to the genes annotated in the reference genome. Of the 3,953 novel lncRNAs detected in this study, 95 were classified as intergenic lncRNAs. The majority of the lncRNAs were classified as ‘genic containing exonic’ (33% of the total lncRNAs classified by the FEELnc_classifier), as shown in Fig. 3E. The second category most represented by the detected lncRNAs was ‘antisense’ (25%) (both intronic and exonic), followed by ‘intergenic’ (17%) (Fig. 3E).

Meta-analysis of DE genes

In relation to the phenotypic assessment performed for the samples analyzed in the meta-analysis, the average PF% values for the Assaf and Churra animals from DataSet 1 were 1.08 and 1.47, respectively, whereas for the two groups of Assaf suckling lambs contrasted within DataSet 2, the mean PF% values were 1.06 and 2.09 for Low-KKCF and High-KKCF groups, respectively. A box plot distribution for these four groups of samples is shown in Supplementary Fig. 2A. The results of the ANOVA performed for the logPF values in the different considered contrasts are given in Supplementary Fig. 2B. Significant differences were found between the Assaf Low-KKCF and Assaf High-KKCF samples from DataSet 2 (p-value = 0.00277) but not between the Churra and Assaf samples from DataSet 1 (p-value = 0.0665). On the other hand, no significant differences were found when the two groups of samples with lower (Assaf DataSet 1 vs. Assaf Low-KKCF DataSet 2) or higher (Churra DataSet 1 vs. Assaf High-KKCF DataSet 2) fatness levels from the two datasets were compared (p-value = 0.906 and 0.052, respectively). Interestingly, the comparison of the samples grouped as LFD (Assaf from DataSet 1 + Assaf Low-KKCF from DataSet 2) vs. those grouped as HFD (Churra from DataSet 1 + Assaf High-KKCF from DataSet 2) showed a highly significant result (p-value = 0.00109) (Supplementary Fig. 2B). Thus, the phenotypic differences for the logPF trait were maximized when evaluating the two datasets together. Based on these results, we considered the meta-analysis approach proposed here to compare the transcriptome samples included in the LFD group with those included in the HFD group to improve the statistical power of our transcriptomic analysis and to gain a deeper understanding of the gene expression pathways and regulatory lncRNAs underlying fat deposition in suckling lambs.

The power analysis estimated that to achieve a power of 80% in the DEA, thereby detecting 80% of DEGs, when the FCh is set as 2 (log₂FCh = 1), at least 9 replicates per group are necessary; for a FCh of 3 (log2FCh = 1.58), a minimum of 5 samples per group are required, and for FCh = 4 (log2FCh = 2), the minimum drops to 4 samples per group (Supplementary Fig. 3). Focusing on the two datasets analyzed herein, for DataSet 1, which has a minimum of 6 replicates per group, the power analysis estimates that approximately the 68, 94 and 99% of the DEGs are detected with log₂FCh 1, 1.5 y 2, respectively. Meanwhile, for DataSet 2, which comprises 4 replicates per group, approximately 38, 81, and 94% of DEGs are detectable with log2FCh values of 1, 1.58, and 2, respectively.

After the DEA, we obtained 715 and 153 DE genes for DataSet1 and DataSet2, respectively. Then, exploiting the meta-analysis approach with the metaRNASeq R package, we obtained 632 DE genes using Fisher’s method and 463 using INCT (Fig. 4A). For each of the methods used in the meta-analysis, we obtained the Integration-driven Discovery Rate (IDR), which measures the proportion of DE genes found in the meta-analysis that were not identified in any of the individual studies alone. The IDRs were 12.82 and 17.93 for Fisher’s method and INCT, respectively. In addition, we calculated the Integration-driven Revision Rate (IRR) parameter, which measures the proportion of DE genes that were identified in any of the individual studies but not in the meta-analysis. The IRR for Fisher’s method and INCT were 35.63 and 55.61, respectively (Fig. 4A).

Figure 4. Summary of the meta-analysis results using fisher’s method and the inverse normal combination technique (INCT). (A) Table showing the number of differentially expressed (DE) genes for each method, the number of DE genes identified in the meta-analysis that were not identified as DE in individual studies (integration driven discoveries, IDD), number of DE genes identified in individual studies but not in meta-analysis (loss) and the Integration-driven discovery Rate (IDR) and the Integration-driven revision Rate (IRR) which are the corresponding proportions of IDD and loss. (B) Venn diagram with the DE genes identified for each of the individual studies, Suárez-Vega et al.¹³ and Alonso-García et al.¹² and for meta-analyses (fisher’s method and INCT). Those DE genes detected in common by both meta-analysis methods are highlighted in red.

Finally, 250 genes were selected as keepDE genes (genes commonly identified as DE in the meta-analysis by the two analytical methods (p_adj <0.05) whose log2FCh sign (+/−) was consistent for both datasets; Fig. 4B). From the list of 250 keepDE genes, 214 were mRNAs, 20 were novel lncRNAs and 16 were lncRNAs previously annotated in the sheep reference genome. Among the keepDE genes, 148 had higher expression levels in the HFD group (keepHFD) than in the LFD group. Among these genes, 13 were lncRNAs (six novel lncRNAs and seven already annotated). Among the 102 keepDE genes with higher expression in the LFD (keepLFD) group than in the HFD group, 23 were lncRNAs (14 novel lncRNAs and nine already annotated). Almost all lncRNAs detected as keepDE (34) had a cis-associated protein-coding gene. These cis-associated genes, along with the protein-coding keepDE genes, were included in the functional enrichment analysis.

In the functional enrichment analyses, for the keepHFD genes, we obtained a total of 29 significantly enriched GO terms. By category, we obtained 16 significantly enriched terms in the GO: Biological Process category (GO:BP), where the most significant terms were related to angiogenesis processes (e.g., ‘epithelial cell development’, ‘cellular response to vascular endothelial growth factor stimulus’, ‘regulation of vascular permeability’, ‘vascular endothelial growth factor signalling pathway’, ‘positive regulation of angiogenesis’ and ‘positive regulation of vasculature development’), and 13 significantly enriched terms against the GO: Cellular Component category (GO:CC), where the most significant terms were ‘platelet alpha granule lumen’, ‘stress fiber’ and ‘contractile actin filament bundle’. No significant terms were identified in the GO: Molecular Function (GO:MF) category (Supplementary Table 2).

For the enrichment analysis of the keepLFD, we obtained only one significant term for the GO:MF category, ‘C-acyltransferase activity’ (Supplementary Table 2).

Detection of co-expressed gene modules and functional analysis

Six co-expressed gene modules were detected using CEMiTool. Briefly, the size of the modules ranged from 607 genes (Module 1) to 39 genes (Module 6). There were 21 keepDE genes for Module 1 and seven for Module 3, and the number of keepDE genes was less than five for the remaining modules. Although lncRNAs were found in all the modules identified, only modules 1, 2 and 4 retained keepDE lncRNAs. In addition, modules 1 and 2 had novel lncRNAs among the hub genes identified by CEMiTool. All the information obtained from the co-expression analysis is summarized in Fig. 5.

Figure 5. Summary of the results from CEMiTool co-expression network analysis. (A) Table showing the number of protein coding genes (N PGenes), lncRNAs (N lncRNAs), differential expressed Protein-coding genes (DEPgenes) and DElncRNAs detected and the hub genes for each module. (B) Gene set enrichment analysis (GSEA) for each module of co-expressed genes. Circles represent the normalized enrichment score (NES) obtained for each module in each group, High fat deposition (HFD) and low fat deposition (LFD). For Module 1 no significant NES differences between groups were found.

Based on the GSEA, significant differences between the HFD and LFD groups were obtained for all modules except for Module 1 (Fig. 5). As part of the study of the regulatory function of lncRNAs, co-expressed gene modules that included keepDE lncRNAs were further analyzed. Modules 2 and 4 were the only modules containing keepDE lncRNAs for which, based on the NES, significant differences in gene expression levels between the HFD and LFD groups were detected (Fig. 5). For Module 2, which had overall higher gene expression in the HFD group (NES = 2.024), we found 113 significantly enriched GO terms (q-value < 0.05), all of which belonged to the GO:BP category. The most significant terms were ‘cell migration involved in sprouting angiogenesis’ and ‘sprouting angiogenesis’ (see details in Supplementary Table 3). This module included one keepDE lncRNA, LOC114110038 (Fig. 6), which is trans-associated (p-value = 0.003, Pearson’s correlation = 0.6) with the GPLD1 gene. This gene is related to GO terms such as ‘response to peptide hormone’, ‘regulation of lipid metabolic process’ and ‘sprouting angiogenesis’ (Fig. 6). Remarkably, Module 2 included the novel lncRNA MSTRG.18268, which is cis-associated with the VEGFA gene. This gene was identified as a keepDE gene in our meta-analysis (Fig. 6).

Figure 6. Functional enrichment analysis for co-expressed genes in Module 2. (A) Dot-plot of the 10 most significantly enriched terms (q-value <0.05). (B) Heatmap of the pearson’s correlation between the DElncRNA in the X-axis and the co-expressed protein-coding genes in Module 2 in the Y-axis based on their expression level. (C) Gene network and the five most enriched terms. The keepDE genes detected in the meta-analysis are highlighted in bold.

On the other hand, for Module 4, which showed overall greater expression in the LFD group (NES = 1.533) (Fig. 5), 19 terms were found to be significantly enriched. Of these significant GO terms, 10 belonged to the GO:BP category, and the most significant term in this category was ‘cellular glucuronidation’. The remaining significantly enriched terms belonged to the GO:MF category, and the most significantly enriched term was ‘glucuronosyltransferase activity’ (Supplementary Table 3). Module 4 included two keepDE lncRNAs (Fig. 5), MSTRG.11585 and LOC105606890, which were trans-associated with genes included within the GO terms ‘positive regulation of cold-induced thermogenesis’ and ‘temperature homeostasis’ (Fig. 7).

Figure 7. Functional enrichment analysis for co-expressed genes in Module 4. (A) Dot-plot of the 10 most significantly enriched terms (q-value <0.05). (B) Heatmap of the Pearson’s correlation between the DElncRNAs in the X-axis co-expressed and the co-expressed protein-coding genes in Module 4 in the Y-axis based on their expression level. (C) Gene network and the five most enriched terms. The keepDE genes detected in the meta-analysis are highlighted in bold.

Finally, the QTL annotation performed with the GALLO R package for the genomic intervals harboring the 37 keepDE lncRNAs identified corresponded with QTLs previously reported in sheep, especially the QTL class of Meat and Carcass (47.06%), followed by the Milk (26.47%), Health (17.65%), and Production (5.88%) QTL classes. Considering the interest in the correspondence between the 37 identified keepDE lncRNAs and sheep-reported QTLs from the Meat and Carcass categories, detailed information on these correspondences is given in Supplementary Table 4. As a remarkable coincidence from this table, the keepDE lncRNA LOC114110038 included Module 2, which is located in the region of sheep chromosome 21 where several QTLs associated with fatty acid composition (i.e., palmitic, arachidonic, linoleic acids, etc.) have been annotated (see details in Supplementary Table 4). In addition, the novel MSTRG.11585 lncRNA, belonging to Module 4, was located in a region of chromosome 11 associated with QTLs involved in ‘internal fat amount’ in sheep (see details in Supplementary Table 4).

Discussion

In recent years, RNA-Seq has emerged as a valuable tool for identifying key genes involved in adipose tissue development and metabolism. Regulatory elements such as lncRNAs have been implicated in fat deposition in sheep, particularly in relation to tail fat, in some specific breeds.^8,10,11,52 However, their role in perirenal fat deposition has not been thoroughly studied. Perirenal fat, constituting approximately 11.4% of total fat, is significant in suckling lambs and serves as an indicator of carcass quality.¹⁸ Therefore, the characterization of lncRNAs in perirenal fat carried out in the present study provides insight into the regulatory mechanisms influencing fat deposition in lambs and offers valuable information on the molecular factors influencing fatness-related quality traits.

In the present study, a total of 3,953 lncRNAs were identified using FEELnc software; from them, approximately 99% (3,913) were confirmed as non-coding transcripts by the CPAT software, highlighting the accuracy of the detection of lncRNAs with FEELnc. These results are similar to those previously reported for subcutaneous fat in sheep (4,464)⁹. However, these findings differ from those in tail fat, where the novel lncRNAs detected were more heterogeneous (<500 in Bakhtiarizadeh and Salami⁸, 5,722 in He et al.⁵² and 9,082 in Ma et al.¹⁰). The distribution of protein-coding genes and both novel and annotated lncRNAs across ovine chromosomes showed a concentration on the first three chromosomes. There was a clear correlation (Pearson correlation = 0.72) between the number of annotated mRNAs and number of novel lncRNAs on each chromosome, suggesting a relationship in their distribution. This distribution pattern was consistent with studies in adipose tissue and other tissues in sheep.^52–54 The lncRNAs annotated in the perirenal fat also shared typical characteristics already known for lncRNAs, including transcript length, exon number, and expression levels. Regarding transcript length and number of exons, known and novel lncRNAs identified in the perirenal fat of suckling lambs had similar transcript lengths (Fig. 3B and 3C), generally shorter than those of mRNAs, and an average of two exons. These results are consistent with findings in other sheep studies.⁸,²⁴,^54–56 In terms of expression, novel lncRNAs showed higher transcript expression levels than previously annotated lncRNAs. In addition, protein-coding genes were consistently more highly expressed than lncRNAs, which agrees with previous results in sheep.^9,10 Overall, the fact that our results are consistent with common features of lncRNAs observed in other sheep research provides support for the robustness of the lncRNAs identified in the perirenal fat transcriptome of suckling lambs.

Statistical analyses evaluating the power of the DEA underscore the critical importance of the number of replicates in achieving sufficient power to detect DEGs in individual datasets. The power analysis results presented here in reference to perirenal fat transcriptome samples, demonstrated that to achieve a detection power of 80% of the DEGs, for a log₂FCh of one, a high number of replicates per group would be required (n = 9). However, for log₂FCh values of 1.5 or above, with four replicates we would be able to detect more than 80% of the DEGs. When increasing the number of samples in a study is not possible, the use of a meta-analysis provides the opportunity to increase the power for the identification of key genes by combining independent studies focused on a common biological question.⁵⁷ This approach improves the reliability and robustness of the results because it integrates the signals of moderate significance from each individual study while controlling for false positives.^58,59 In the meta-analysis performed herein, the IDRs increased by 12.82 and 17.93%, respectively, for both p-value combination approaches, indicating a clear improved detection power reached by our meta-analysis when compared to DEA (Fig. 4A). In addition, the meta-analysis implemented in the present study appeared to discard false-positive DE genes detected in smaller sample pools, with IRRs of 35.63 and 55.61% for Fisher’s method and INCT, respectively.

Both the keepDE genes, which were more highly expressed in the HFD group, and the Module 2 genes were related to angiogenesis processes (Supplementary Tables 2 and 3). The process of forming new blood vessels, known as angiogenesis, is crucial to adipogenesis⁶⁰ and is a rate-limiting step in fat tissue expansion.⁶¹ Angiogenesis-related GO terms were previously found to be enriched in lambs with a greater fat percentage in our previous work comparing High-KKCF and Low-KKCF Assaf suckling lambs.¹² Despite the fact that in a previous study comparing the perirenal fat transcriptome of Churra and Assaf suckling,¹³ this term was not found to be enriched, several other significantly enriched terms related to vessel formation were found to be enriched in Churra lambs. For example, signaling mechanisms generated as a result of cell adhesion interactions related to ‘focal adhesion’ were enriched in the Churra transcriptome samples analyzed by Suárez-Vega et al.¹³ Focal adhesion is a process that is crucially important as a primary driver of angiogenesis.⁶² In addition, closely related to angiogenesis, the term ‘ERK1 and ERK2 cascade’ was also significantly enriched in the Churra samples analyzed by Suárez-Vega et al.¹³

Moreover, the term ‘regulation of lipid metabolic process’ was significantly enriched in both the Churra lambs in the Suárez-Vega et al.¹³ group and the HFD group in the present meta-analysis. The regulation of lipid metabolism involves genes related to angiogenesis, such as GPLD1, which is trans-associated with the LOC114110038 DElncRNA (Fig. 6). LOC114110038 lncRNA was a keepDE more highly expressed in the HFD group than in the LFD group, and its expression was positively correlated with the expression of the GPLD1 gene. Glycosylphosphatidylinositol phospholipase D1 (GPLD1) is an enzyme that cleaves glycosylphosphatidylinositol and has proven to influence glypican 4, a novel adipokine associated with obesity and insulin resistance in humans.⁶³ Previous studies have shown that serum GPLD1 levels are positively correlated with triglyceride serum levels.⁶⁴ In ruminants, the role of GPLD1 is still unclear. Delosière et al.⁶⁵ hypothesized that this protein, along with AHSG, could mediate increased glucose uptake in non-mammary tissues.

As shown in Fig. 6, LOC114110038 DElncRNA was more highly expressed in the LFD group; however, it interacted with co-expressed genes in Module 2, which had a greater NES in the HFD group. This finding could imply a potential role for LOC114110038 in the inhibition of angiogenesis. This function should be confirmed in further studies, as this lncRNA could be key in the regulation of angiogenesis in the development of adipose tissue in suckling lambs. In addition to its trans-association with GPLD1, LOC114110038 is associated with many QTLs in the sheep genome that are related to meat fatty acid content (i.e., docosapentaenoic acid, arachidonic acid, and gadoleic acid, among others). Hence, LOC114110038 seems to be a promising lncRNA related to fat deposition regulation processes and should be validated in further studies.

On the other hand, the more highly expressed keepDE genes in the LFD group than in the HFD group were associated with ‘C-acyltransferase activity’ and ‘monocarboxylic acid transmembrane transporter activity’. The term ‘C-acyltransferase activity’ was previously found to be associated with a lower amount of perirenal fat in the comparison between Assaf and Churra suckling lambs.¹³ On the other hand, the monocarboxylate transporter 1 (MCT1) is a carrier of short-chain fatty acids, ketone bodies, and lactate in several tissues.⁶⁶ It has been demonstrated that partial inactivation of the MCT1 gene in mice leads to specific adaptations preventing the development of diet-induced obesity, hepatic steatosis and insulin resistance.⁶⁶ In addition, when performing the co-expression network analysis, we detected that Module 4 had a greater NES in the LFD group than in the HFD group. The enriched terms associated with Module 4 were closely related to body temperature regulation via thermogenesis. Previous studies on suckling lambs highlighted thermogenesis as an important process associated with fat tissues in suckling lambs.^12,13,67

For the co-expressed genes in Module 4, we detected two DElncRNAs, MSTRG.11585 and LOC105606890. These DElncRNAs were trans-associated with the G0S2, LCN2, and ACHE genes (Fig. 7).

The G0S2 gene was found to be significantly trans-associated with the MSTRG.11585 lncRNA (p-value = 0.003, Pearson correlation coefficient = 0.61). The expression of the G0/G1 switch 2 (G0S2) gene has a positive relationship with adipocyte differentiation⁶⁸ and lipolysis through the inhibition of adipose triglyceride lipase (ATGL).⁶⁹ El-Assaad et al.⁷⁰ demonstrated that G0S2 is a physiological regulator of adiposity and energy metabolism in a mouse model. G0S2-knockout mice presented enhanced ATGL activity, enhanced lipid metabolism and increased thermogenesis. This gene is also a Module 4 hub gene, as shown in Fig. 7, implying that it is a gene with high connectivity within the module and is likely to play an important role in adipose tissue in the LFD group. Closely related to thermogenesis and trans-associated with the MSTRG.11585 lncRNA, is the LCN2 gene (p-value = 0.0005, Pearson correlation coefficient = 0.68). In murine models, overexpression of lipocalin 2 (LCN2) in white adipose tissue led to significant upregulation of UCP1 and genes involved in mitochondrial oxidation, improving thermogenic adaptation.^71,72 The functions of these two genes might be regulated by MSTRG.11585 due to the positive correlation of this DElncRNA with both genes, so MSTRG.11585 could activate the regulatory cascade of thermogenesis in perirenal adipose tissue. Notably, UCP1 is a thermogenic protein that is thought to be associated only with brown fat mitochondria but has also been shown to be present in beige fat.⁷³

The other DElncRNA of Module 4 was LOC105606890, which was trans-associated with the ACHE gene (p-value = 0.002, Pearson correlation coefficient = 0.62). The main function of AChE is the hydrolysis of acetylcholine (ACh). ACh is a parasympathetic neurotransmitter that has the ability to prevent adipose tissue inflammation by inhibiting the synthesis and release of proinflammatory and inflammatory immune cells.⁷⁴ In animal model studies, obesity has been linked to increased AChE activity.⁷⁵ The acetylcholine signaling pathway has an important effect on beige thermogenesis.⁷⁶ Taken together, these findings suggest that the upregulation of AChE function and G0S2/LCN2 by the LOC105606890 and MSTRG.11585 lncRNAs, respectively, in the LFD samples, together with the presence of UCP1 in Module 4, could imply a greater energy expenditure in LFD lambs, potentially marked by activated thermogenic adipocytes.

In conclusion, the present study implemented RNA-Seq technology to explore the role of lncRNAs in perirenal fat deposition in suckling lambs. The findings revealed a total of 3,953 novel lncRNAs in the perirenal fat of suckling lambs. The distribution of these genes across chromosomes, their typical characteristics, and their expression patterns agreed with common features observed in other sheep studies. The meta-analysis approach applied here, which combines data from two independent transcriptomic datasets, demonstrated increased detection power and filtered out false-positive DE genes, enhancing the reliability of the results. Key genes associated with angiogenesis-related processes were identified in HFD lambs by both meta-analysis and the co-expression network analysis. Notably, the lncRNA LOC114110038 has emerged as a potential regulator that inhibits angiogenesis. In contrast, genes more highly expressed in the LFD were linked to specific activities such as ‘C-acyltransferase activity’ and ‘monocarboxylic acid transmembrane transporter activity’. The co-expression network analysis for Module 4, which was more highly expressed in the LFD samples, revealed an association with body temperature regulation via thermogenesis and identified differentially expressed lncRNAs (MSTRG.11585 and LOC105606890) with trans-associations to key genes (G0S2, LCN2, ACHE) involved in thermogenic processes. These results may suggest a potential increase in energy expenditure in LFD lambs.

Overall, this study provides comprehensive insights into the regulatory mechanisms influencing perirenal fat deposition in suckling lambs that might be commonly shared within and between breeds, highlighting the significant role of lncRNAs in these processes.

Disclosure statement

No potential conflicto of interest was reported by the autor(s).

Additional information

Funding

This work has been performed within the framework of the QUALITY-LAMB project funded by the Junta de Castilla and León Regional Government [Ref. LE084P23]. MAG is funded by a predoctoral fellowship from the Junta de Castilla and León Regional Government and the European Social Fund. PASF is beneficiary of a María Zambrano Grant of the University of León funded by the Ministry of Universities (Madrid, Spain) and financed by the European Union-Next Generation EU.