Transcriptome analysis of messenger RNA and long noncoding RNA related to different developmental stages of tail adipose tissues of sunite sheep

Abstract The tail fat of sheep is the most typical deposited fat, and it can be widely used in human daily life, such as diet, cosmetics, and industrial raw materials. To understand the potential regulatory mechanism of different growth stages of tail fat in Sunite sheep, we performed high‐throughput RNA sequencing to characterize the long noncoding RNA (lncRNA) and messenger RNA (mRNA) expression profiles of the sheep tail fat at the age of 6, 18, and 30 months. A total of 223 differentially expressed genes (DEGs) and 148 differentially expressed lncRNAs were found in the tail fat of 6‐, 18‐, and 30‐month‐old sheep. Based on functional analysis, we found that fat‐related DEGs were mainly expressed at 6 months of age and gradually decreased at 18 and 30 months of age. The target gene prediction analysis shows that most of the lncRNAs target more than 20 mRNAs as their transregulators. Further, we obtained several fat‐related differentially expressed target genes; these target genes interact with different differentially expressed lncRNAs at various ages and play an important role in the development of tail fat. Based on the DEGs and differentially expressed lncRNAs, we established three co‐expression networks for each comparison group. Finally, we concluded that the development of the sheep tail fat is more active during the early stage of growth and gradually decreases with the increase in age. The mutual regulation of lncRNAs and mRNAs may play a key role in this complex biological process.

representative of local elite breeds in Inner Mongolia, and it is the primary meat resource with the benefits of great numbers, delicate flesh, and good flavor, which is very popular among consumers. SS is mainly raised in the Xilingol grassland of Inner Mongolia. These sheep are accustomed to voluntary movement and typically free-feed (or naturally graze) throughout their lives. SS is a meat breed, and the important phenotype of SS is the fat tail. With the growth of SS, the tail fat increases continuously and reaches the weight of about 3-4.5 kg at 30 months of age (30 M). Tail fat can be used by humans as an important source of dietary fat (Kashan et al., 2005;Moradi et al., 2012) and provides the energy needed by the human body. As a by-product of mutton, it can also be used as a raw material for daily-use products, such as soap, cosmetics, and medicinal materials.
Adipose tissue plays a vital role in maintaining the balance of homeostatic metabolic processes in domestic animals. During severe conditions, such as food scarcity resulting from migration, drought, and winter, the tail fat can provide energy (Ben Sassi-Zaidy et al., 2014). According to a previous study (Kashan et al., 2005), fattailed sheep have a low percentage of intramuscular fat and provide good quality lean meat. In contrast, short-tailed sheep have higher intramuscular fat storage. Thus, the mechanism of tail fat deposition is worth studying. Many studies have employed RNA sequencing (RNAseq) to explore differentially expressed genes (DEGs) in the adipose tissues in fat-tailed sheep recently. To gain a better understanding of fat deposition, Li et al. (2018aLi et al. ( , 2018b performed RNA-seq of perirenal, subcutaneous, and tail fat tissues from Guangling Large-Tailed and Small-Tailed Han (STH) sheep to determine their transcriptome profiles. The result showed that a total of 4131 DEGs were identified in tail fat tissue, and 49 genes were shown to be involved in the peroxisome proliferator-activated receptor (PPAR) signaling pathway, which is the key pathway to balance fat metabolism (Corrales et al., 2018). Wang et al., (2014) used transcriptome sequencing to compare the transcriptome profiles of tail fat tissue between Kazak and Tibetan sheep. This study identified 646 DEGs between the two breeds, and the top two genes with the largest fold change (NELL1 and FMO3), which may be relevant to fat metabolism in adipose tissues.
Further, the lncRNAs and mRNAs associated with tail fat deposition and development in Lanzhou fat-tailed sheep (long fat-tailed sheep), STH sheep (thin-tailed sheep), and Tibetan sheep (short thin-tailed sheep) were analyzed; 407 DEGs and 68 differentially expressed (DE) lncRNAs were identified (Ma et al., 2018). It was shown that the DEGs and target genes of DE lncRNAs were enriched in fatty acid metabolism and fatty acid elongation-related pathways through gene ontology (GO) analysis and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis, which contribute to fat deposition. Network contribution based on DE mRNA and lncRNAs shows that some DE lncRNAs (TCONS_00372767, TCONS_00171926, TCONS_00054953, and TCONS_00373007) may play an important role in tail fat deposition processes. Bakhtiarizadeh and Salami (2019) have performed the transcriptome analysis in fat-tailed (Lori-Bakhtiari) and thin-tailed (Zel) Iranian sheep breeds and identified 7 DE lncRNAs and 311 DEGs between the two breeds. Further, the target prediction analysis shows that the novel lncRNAs can regulate the expression of genes involved in lipid metabolism through cis-or transregulation. In addition, the animal quantitative trait loci database suggested 1 intronic and 6 intergenic lncRNAs as candidates of sheep fat-tail development.
Transcriptome analyses were performed in specific sheep tissues to reveal the potential regulatory roles of lncRNAs, such as in the skeletal muscle (Chao et al., 2018;Li et al. 2018aLi et al. , 2018bWei et al., 2018), pituitary Yang et al., 2020;Zheng et al., 2019), testis Zhang et al., 2017), and ovaries Miao, Luo, Zhao, & Qin, 2016a, 2016b. However, related research on mRNA and lncRNA in SS tail fat is still lacking, including the regulation mechanism of fat deposition and related molecular pathways of tail fat development. To better understand the potential role of mRNAs and lncRNAs in fat-tailed sheep, we explored the transcriptomic differences in SS's tail fat at three different growth stages, 6 months of age (6 M), 18 months of age (18 M), and 30 M. This facilitated the characterization of the mRNA and lncRNA expression profiles in the fat tail of SS and elucidate the molecular mechanism of fat deposition. Our findings may lay a foundation for further studies in fattailed sheep. In particular, our study provides some information on the mechanism of fat development in fat-tailed sheep during different growth processes, which is of great significance for the development and utilization of by-products of meat breeds of sheep.

| Animal and tail fat tissue collection
Nine castrated Sunite rams were selected from three different growth stages, 6 M (n = 3), 18 M (n = 3), and 30 M (n = 3), respectively. All sheep were raised under the same conditions, including food, water source, and environment. After slaughtering, adipose tissue was sampled from the tail fat (top 1/3) and cut into small pieces of 2 mm × 2 mm × 2 mm (Miao et al., 2015). These small pieces were immediately placed into cryotube (sterile without enzyme), frozen in liquid nitrogen, and transferred to −80°C until RNA extraction.

| RNA extraction and RNA-seq
Total RNA from the nine adipose tissue samples was extracted using

| Transcripts assembly
The FastQC (http://www.bioin forma tics.babra ham.ac.uk/proje cts/ fastq c/) software was used to verify the sequence quality, and adaptor contamination, low-quality bases, and undetermined bases in the raw data were removed by the Cutadapt software (Martin, 2011).
The clean reads were mapped into the genome of sheep (Ovis aries v3.1) using Bowtie2 (Langmead & Salzberg, 2012) and Tophat2 (Kim et al., 2013), and mapped reads were assembled using the StringTie software . To reconstruct a comprehensive transcriptome, all transcriptomes from sheep samples were merged using Perl scripts. After the generation of the final transcriptome, the expression levels of all transcripts were estimated using the StringTie and R package Ballgown (Frazee et al., 2015).

| lncRNA identification and different expression analysis
First, transcripts that overlapped with known mRNAs and transcripts smaller than 200 bp were excluded. Subsequently, the Coding Potential Calculator (CPC) (Kong et al., 2007) and Coding-Non-Coding Index (CNCI) software tools (Sun et al., 2013) along with Pfam database (Punta et al., 2012) were utilized to predict transcripts with coding potential. Transcripts that scored CPC < −1 and CNCI < 0 were discarded. The remaining transcripts with class code

| Target gene prediction and functional analysis of lncRNAs
In order to explore the functions of lncRNAs, the DE lncRNAs were analyzed for target prediction. In this study, coding genes 100,000 bp upstream and downstream of the target gene were considered as the cis-target genes. The targets in trans were defined by calculating the expressed correlation with lncRNAs. Then, we performed GO and KEGG analysis of the DE lncRNA targets and mRNAs, respectively, using the in-house scripts. The significance was expressed as FDR <0.05.

| Construction of the co-expression network
To gain a better understanding of interactions between the DEGs and DE lncRNAs, the Pearson correlation coefficient (COR) of mRNA-lncRNA co-expression network was calculated. Finally, the mRNA-lncRNA co-expression network was constructed using Cytoscape (version 3.7.2) with an absolute value of COR ≥0.7.

| RNA-Seq analysis
The results of the RNA-Seq reads mapping are shown in Table 1. To identify the potential function of lncRNAs in tail fat tissues, the nine cDNA libraries were sequenced using the Illumina Hiseq 4000 platform. A total of 139.82 G raw data were generated from the nine adipose tissues. After filtering out low-quality reads, 131.61 G valid data were obtained, and the average valid ratio (reads) was 94%.
In detail, the valid reads obtained were as follows: (1)

| Summary of lncRNA and mRNA expression
To

| Different expression analysis
We compared the expression profiles between any two stages (

| Functional analysis of DEGs
The to the development of the sheep tail fat. We found five highly expressed DEGs, namely EHHADH, LPIN1, ACACA, THRSP, and GPAT4, which were related to these functions. The previous study shows that LPIN1 deficiency will lead to a significant decrease in adipose tissue and abnormal expression of adipogenic genes. Conversely, increased expression of LPIN1 in skeletal muscle or adipose tissue will promote obesity in mice (He et al., 2009). EHHADH is associated with the expression of genes involved in the tricarboxylic acid cycle, mitochondrial and peroxisome fatty acid oxidation, and is indispensable for the production of medium-chain dicarboxylic acids in mice during fasting (Houten et al., 2012). ACACA is considered to be a key regulator of fat production and a limiting factor in the synthesis of long-chain fatty acids. Acetyl-CoA can be converted to malonyl-CoA (Pena et al., 2019), which may play a key role in energy metabolism and homeostasis in sheep tail fat cells. THRSP is involved in the process of adipogenesis in rodents, and it may be a potential marker gene for bovine intramuscular fat. Studies have shown that THRSP is mainly expressed in adipocyte nuclei, intramuscular adipocytes, and related cells and expressed in mature adipocytes rather than in the early stages of adipogenesis (Schering et al., 2017). In our study, the THRSP gene with a higher expression level ( (Cooper et al., 2015), and ACSM1, ACSM3, and ACAT1 were related to fat deposition and fatty acid metabolism (Berton et al., 2016;Guo et al., 2014;Huang, Guo, et al., 2017). ECHS1 was shown to be associated with the fatty acid beta-oxidation (Peng et al., 2016).
Studies have shown that TKT expression affects fatty acid oxidation and mitochondrial function (Tian et al., 2020). On the other hand, a total of 8 KEGG pathways were significantly enriched in three different stages. They were mainly focused on metabolism processes, including carbon metabolism, mineral absorption, glutathione metabolism, butanoate metabolism, and some related amino acid metabolism. Based on the KEGG pathway analysis, those highly expressed DEGs were related to butanoate metabolism, fatty acid metabolism, glycerolipid metabolism, and PPAR signaling pathway, which may contribute to the fat deposition in sheep tail fat.

| Hierarchical clustering analysis of fatrelated DEGs
Based on the above analysis and further screening, we obtained several DEGs that may be related to fat-tail development. We performed a hierarchical clustering analysis to show the expression patterns of these DEGs (Figure 8). It is not difficult to find that most of the genes are active in the early months of age, especially during the 6 M, and the expression of DEGs decreased gradually with the increase in age.
Therefore, we indicate that the vitality of fat development weakens with the increase of age, that is to say, the development of tail fat will be more active at the age of 6 M, but will gradually decrease at the age of 18 M and 30 M. There is a significant difference between 6 M and 30 M of age. Further, the expression pattern at 18 M, as the middle month, plays the role of transitioning from high metabolic activity to low metabolic activity. However, our findings are only possible in theory, and the mechanism needs to be further identified.

| Target gene prediction and functional analysis
In order to explore how lncRNA participates in regulation, we predicted DE lncRNAs based on cis-and transregulation in three different stages of the fat-tail development. In our study, 148 DE lncRNAs (68 upregulated and 80 downregulated) were obtained and the target genes prediction analysis was performed in these DE lncRNAs. with obesity, the concentration and activity of CA3 in rat adipose tissue decreased (Wang et al., 2006  . CYP1A1 is only expressed at 6 M in our study, and it was reported to be expressed in brown adipose tissue (Galván et al., 2005). The study showed that using a specific anti-LBP antibody to inhibit LBP activity can improve the adipogenic status of fully differentiated adipocytes, which makes LBP is a novel adipokine that might display an essential role in inflammation and obesity-associated adipose tissue

| Co-expression network construction
We constructed three co-expression networks based on DEGs

| CON CLUS IONS
In our study, we used transcriptome analysis to explore the un-