RNA-Seq analysis
The results of the RNA-Seq reads mapping is 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.61G valid data were obtained, and the average valid ratio (reads) was 94%. In detail, the valid reads obtained were: (1) 94,550,880; 98,994,810; and 99,167,036 per fat tail tissue sample from 6M (A1, A2, and A3), (2) 96,734,688; 97,915,792; and 95,384,562 per fat tail tissue sample from 18M (B1, B2, and B3), and (3) 98,530,172; 95,442,058; and 100,597,240 per fat tail tissue sample from 30M (C1, C2, and C3), respectively. The average percentage of Q20 and Q30 base was more than 99% and 98%, respectively, and the percentage of the Guanine-Cytosine (GC) content of each sample on an average was 48%. Above all, we indicate that results of the RNA sequencing were highly reliable, and follow-up analysis can be carried out.
Table 1
Summary of the reads mapped to the tail adipose tissue transcriptomes
Sample
|
Raw Data
|
Valid Data
|
Valid Ratio(reads)
|
Q20%
|
Q30%
|
GC content%
|
|
Read
|
Base
|
Read
|
Base
|
|
|
|
|
A1
|
101449274
|
15.22G
|
94550880
|
14.18G
|
93.20
|
99.74
|
98.48
|
48
|
A2
|
103736482
|
15.56G
|
98994810
|
14.85G
|
95.43
|
99.81
|
98.63
|
48
|
A3
|
105051744
|
15.76G
|
99167036
|
14.88G
|
94.40
|
99.77
|
98.51
|
48
|
B1
|
102382092
|
15.36G
|
96734688
|
14.51G
|
94.48
|
99.76
|
98.47
|
48
|
B2
|
103267850
|
15.49G
|
97915792
|
14.69G
|
94.82
|
99.79
|
98.64
|
49
|
B3
|
102376260
|
15.36G
|
95384562
|
14.31G
|
93.17
|
99.70
|
98.36
|
47
|
C1
|
105462958
|
15.82G
|
98530172
|
14.78G
|
93.43
|
99.70
|
98.34
|
48
|
C2
|
101700360
|
15.26G
|
95442058
|
14.32G
|
93.85
|
99.71
|
98.35
|
49
|
C3
|
106595358
|
15.99G
|
100597240
|
15.09G
|
94.37
|
99.67
|
98.37
|
47
|
A1, A2, and A3 are 6 months of age; B1, B2, B3 are 18 months of age; C1, C2, and C3 are 30 months of age; valid ratio (reads) = (valid reads/raw reads).
Summary of lncRNA and mRNA expression
To understand the expression profile of the lncRNAs in the tail fat tissue of SS, we identified the expression levels of the lncRNAs and compared them with the expression levels of mRNA. First, a total of 20,670 mRNAs and 6,794 lncRNAs were identified. Here, 5,722 lncRNAs were identified as novel, and the remaining 1,702 lncRNAs were identified as known, which is more than the number of lncRNAs present in chicken [30] and cattle [31] adipose tissue. According to the classification rules, we classify novel lncRNAs as 1,395 (21%) class i, 354 (5%) class j, 288 (4%) class o, 3,174 (47%) class u, 511 (7%) class x, and 1,702 (16%) class = as known lncRNA (Fig. 1a). The chromosome distribution of lncRNAs is shown in circos figure (Fig. 1b). We found that most of the lncRNAs were mainly enhanced in chromosomes 1, 2, and 3. Then lncRNAs and mRNAs were compared with exon number, open reading frames (ORF) length, transcript length, and expression levels. The mRNAs and lncRNAs had 9.7 and 1.8 exons on an average; 86% of lncRNAs contained 1–2 exons, and 38% of mRNAs contained more than 9 exons (Fig. 2a). The size of the ORF of lncRNAs and mRNAs are mainly concentrated in the range of 0–200 and 0–600 amino acids, respectively (Fig. 2b). The majority of lncRNAs and mRNAs were > 1000 bp in size, and short-range (≤ 300 bp–600 bp) lncRNAs were more than mRNAs. The average length of lncRNAs and mRNAs was 3184 bp and 1903 bp, respectively. This significant difference might be due to the quantity gap of lncRNA and mRNA under similar distribution patterns (Fig. 3a). The expression levels of lncRNA were higher than the expression levels of mRNAs (Fig. 3b), which suggests that the lncRNAs may play an important role in the development of sheep tail fat tissue.
Different expression analysis
We compared the expression profiles between any two stages (30M vs 18M, 18M vs 6M, and 30M vs 6M) using |log2(Fold Change) | ≥ 1 and FDR < 0.05 to identify DEGs and DE lncRNAs. In the comparison between 30M vs 6M, we found 377 DEGs (167 up-regulated and 210 downregulated genes). In the 30M vs 18M group, 125 DEGs (56 upregulated and 69 downregulated genes) were obtained. In the comparison of 18M vs 6M, 75 DEGs (38 upregulated and 37 downregulated) were found (Fig. 4). Furthermore, 4 DEGs were commonly expressed in the comparison groups of 30M vs 18M and 18M vs 6M, including IFIT5, THBS1, ENSOARG00000004030, and ENSOARG00000018868. Sixty-eight DEGs were commonly expressed in 30M vs 6M and 30M vs 18M, and 35 DEGs were commonly expressed between 30M vs 6M and 18M vs 6M. On the other hand, 151 DE lncRNAs were identified in 30M vs 6M, 30M vs 18M, and 18M vs 6M. Among them, 78 DE lncRNAs including 38 upregulated (36 novel, 2 known) and 40 downregulated (39 novel, 1 known), 71 DE lncRNAs including 30 upregulated (30 novel, 0 known) and 41 downregulated (38 novel, 3 known), 61 DE lncRNAs including 34 upregulated (33 novel, 1 known) and 27 downregulated (25 novel, 2 known), respectively (Fig. 5). Fifteen DE lncRNAs were commonly expressed in the comparison groups of 30M vs 18M and 18M vs 6M, 25 DE lncRNAs were commonly expressed in 30M vs 6M and 18M vs 6M, and 19 DE lncRNAs were commonly expressed in 30M vs 6M and 30M vs 18M.
Functional analysis of DEGs
The top 15 GO terms and KEGG Pathway analysis performed in DEGs of the fat tail tissue at three different growth stages are shown in the scatterplot (Fig. 6). GO terms were determined by three functions, including cellular component, biological process, and molecular function. We found that more than half of the GO terms were enriched in the biological process in the comparison groups, and the cellular component obtained the most number of genes in the three comparison groups. In the three stages, the DEGs were significantly enriched in 80 GO terms, and several fat related functions were obtained, including fatty acid beta-oxidation, triglyceride biosynthetic process, triglyceride homeostasis, lipid homeostasis, lipid biosynthetic process, regulation of fat cell differentiation, suggesting that these functions might contribute 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 [32]. 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 [33]. 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 [34], 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 [35]. In our study, the expression of the THRSP gene was higher in the tail adipose tissue during 6M and 18M, and significantly lower at 30M. We can speculate that in the early fat tissue of sheep’s tail, fat hypertrophy is mainly manifested by the increase in the number of fat cells, and as the age increases, fat hypertrophy is reflected by the increase in the volume of fat cells. In addition, there were also some highly expressed genes related to fat metabolism, such as GPAT4, ACSM1, ACSM3, ACAT1, TKT, and ECHS1. GPAT4 was reported to be responsible for maintaining triacylglycerol stores [36], and ACSM1, ACSM3, and ACAT1 were related to fat deposition and fatty acid metabolism [37, 38, 39]. ECHS1 was shown to be associated with the fatty acid beta-oxidation [40]. Studies have shown that TKT expression affects fatty acid oxidation and mitochondrial function [41]. 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, PPAR signaling pathway, which may contribute to the fat deposition in sheep tail fat.
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 (Fig, 7). It is not difficult to find that most of the genes are active in the early months of age, especially during the 6M, 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 6M, but will gradually decrease at the age of 18M and 30M. There is a significant difference between 6M and 30M of age. Further, the expression pattern at 18M, 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 trans-regulation 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. A total of 186, 113, and 150 GO terms were significantly enriched in target genes of 30M vs 6M, 30M vs 18M, and 18M vs 6M (FDR < 0.05), respectively. The top 15 GO terms and KEGG pathway of target genes of DE lncRNAs in the three comparison groups are shown in the scatterplot (Fig. 8). There were 5 common GO terms enriched in the three comparison groups, namely plasma membrane, extracellular exosome, membrane, extracellular, and cytoplasm. The target genes of DE lncRNAs in 30M vs 6M were significantly enriched in 4 KEGG pathways, including calcium signaling pathway, cell adhesion molecules, oxytocin signaling pathway, and tight junction. Among these DE lncRNA, only one cis-regulated target gene was obtained: MSTRG.13384.1 targets CLDN4. We found that most of the lncRNAs were targets to more than 20 mRNAs as their trans-regulators, MSTRG.20969.1 targets to 53 mRNAs, as the largest number in 30M vs 6M. The most commonly enriched top 5 target genes were SLC7A6 (38 DE lncRNA), CDS2 (32 DE lncRNA), CA3 (31 DE lncRNA), SLC6A2 (31 DE lncRNA), PRTG (30 DE lncRNA). These target genes were mainly enriched in cellular components, such as membrane and integral component of membrane. Previous studies have indicated that with obesity, the concentration and activity of CA3 in rat adipose tissue decreased [42]. The complement and coagulation cascades are the only KEGG pathway that is significantly enriched in 30M vs 18M. In this comparison group, MSTRG.12899.1 and ENSOART00000028120 were connected to 38 mRNAs as the largest number in 30M vs 18M, respectively. The most commonly enriched top 5 target genes were SNORA23 (29 DE lncRNA), ERICH6B (27 DE lncRNA), ENSOARG00000018868 (22 DE lncRNA), FBP2 (22 DE lncRNA), and ENSOARG00000014791 (16 DE lncRNA), among which ENSOARG00000018868 is related to lipid binding. In 18M vs 6M, MSTRG.14210.1 targets to 48 mRNAs as the largest number. The most commonly enriched top 5 targets genes were HECW1 (23 DE lncRNA), CRHR2 (19 DE lncRNA), FRK (19 DE lncRNA), IFIT5 (19 DE lncRNA), and PTPRZ1 (19 DE lncRNA). Based on the GO analysis, these target genes were mainly enriched in molecular function, including ATP binding and Hippo signaling pathway, and Steroid hormone biosynthesis were obtained in KEGG analysis.
Further, we obtained several fat related DE target genes in these three comparison groups. Among these targets DEGs, some of them were regarded as common DE target genes in two comparison groups, including TRIB3, ACSM1, ACSM3, TKT, SPTB, and ASGR in 30M vs 6M and 30M vs 18M, and CYP1A1 and LBP in 30M vs 6M and 18M vs 6M. Based on the GO analysis, these DEGs were mainly enriched in fatty acid ligase activity, fatty acid biosynthetic process, glyceraldehyde-3-phosphate biosynthetic process, negative regulation of fat cell differentiation, and some lipid related functions, such as lipid binding and lipid homeostasis. It has been reported that TRIB3 was might inhibit subcutaneous fat deposition in Large White pig, and lncRNA XLOC_064871 trans-regulates TRIB3, so XLOC_064871 might play an important role in adipocyte differentiation and fatty acid metabolism in pig [43]. CYP1A1 is only expressed at 6M in our study, and it was reported to be expressed in brown adipose tissue [44]. 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 dysfunction [45]. In addition, we found that the same target gene was affected by different amounts and types of DE lncRNAs at different ages. For example, in 30M vs 6M and 30M vs 18M, 13 and 11 DE lncRNAs were connected to ACSM1, and there were four common lncRNA targets to ACSM1 between two comparison groups; 13 and 15 lncRNAs were connected to ACSM3, 22 and 19 lncRNAs were connected to TKT between two comparison groups, however MSTRG.3410.1 is the only one lncRNA that acts as a target to ACSM3 and TKT, which suggest that MSTRG.3410.1 may be related to the fat deposition. It could indicate that different lncRNAs with different regulation patterns may impact the target gene expression pattern and play its role in different growth stages of sheep tail fat.
Co-expression network construction
We constructed three co-expression networks based on DEGs and DE lncRNA in sheep fat tail using Cytoscape (version 3.7.2) (Fig. 9). A total of 538, 158, and 184 pairs of co-expression pairs were obtained in 30M vs 6M, 30M vs 18M, and 18M vs 6M, respectively. In the comparison of 30M vs 6M, 78 DE lncRNAs connected to the 26 mRNAs, and 538 pairs of co-expression pairs were obtained (403 positively and 135 negatively correlated). There were 20 lncRNAs connected to more than 10 mRNAs. MSTRG.20969.1 and MSTRG.12518.1 were connected to mRNA 21 and 20, respectively. In the 30M vs 18M, 52 DE lncRNAs connected to the 9 mRNAs, and 149 pairs (118 positively correlated and 31 negatively correlated) were obtained. ENSOART00000028120 and MSTRG.19382.4 were co-expressed with 7 mRNAs. In 18M vs 6M, 60 DE lncRNAs connected to the 13 mRNAs, and 184 pairs (139 positively correlated and 45 negatively correlated) of co-expression pairs were obtained. MSTRG.15348.1, MSTRG.14210.1, and MSTRG.14211.1 were co-expressed with 8 mRNAs. In these two comparison groups, there were only 7 (30M vs 18M), and 8 (18M vs 6M) mRNAs connected to single lncRNA, at most. This indicate that these co-expression pairs might play a crucial role, and lncRNA may regulate the development of sheep tail fat mainly through positive correlation with multiple mRNAs.