Integrative Analysis of lncRNA-miRNA-mRNA Regulatory Network Reveals Candidate lncRNAs Responsible for Meat Quality in Goats

Changsheng He Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Yong Wang Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Jiangjiang Zhu Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Yan Xiong Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Yanyan Li Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Juan Chen Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu Yaqiu Lin (  linyq1999@163.com ) Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization of Education Ministry, Southwest Minzu University, Chengdu


Background
Fat is generally divided into subcutaneous fat, intermuscular fat and intramuscular fat. The amount of intermuscular fat is small and its composition is similar to subcutaneous fat, so it is generally considered that fat is mainly divided into subcutaneous fat and intramuscular fat [1]. The distribution of fat in meat animals and the fat content of each part have an important impact on meat production and quality.
Subcutaneous fat determines the carcass lean rate to a certain extent. The thicker the fat means the lower the lean meat rate. Intramuscular fat has an important in uence on the avor and tenderness of goat meat. Within a certain range, the quality of meat gradually improves with the increase of intramuscular fat content. One of the current research directions for meat is to study how to regulate the content of intramuscular fat and subcutaneous fat to improve meat quality of goats.
Adipogenesis is a complex process regulated by various transcription factors, non-coding RNA and signal pathways [2]. RNA sequencing (RNA-seq) is a revolutionary tool to identify differentially expressed genes (DEGs) regulating various biological processes. It enables us to discover new genes and therefore to describe unannotated transcriptional activity by identifying numerous noncoding transcripts [3]. Long non-coding RNAs (lncRNAs) are non-coding RNAs with a length greater than 200 nucleotides and low protein coding ability. MicroRNAs (miRNAs) are a type of endogenous non-coding RNA with a length of approximately 22 nucleotides. Previous studies have shown that lncRNAs and mRNAs containing the same miRNA binding site can regulate mutual expression levels by competitively binding miRNAs. The recent explosion in knowledge demonstrating the importance of miRNAs and lncRNAs in the regulation of multiple major biological processes, which impacts the development, differentiation, and metabolism have brought these neglected molecular players [4][5][6]. Some studies have shown that the adipogenic differentiation ability of intramuscular adipocytes is signi cantly lower than that of subcutaneous adipocytes. The expression of related genes is low, indicating that there are speci c regulatory mechanisms for intramuscular fat and subcutaneous fat in animals [7,8]. In mammals, the differentiation of preadipocytes has been well studied, especially in bovine and porcine [9][10][11][12][13]. Recently, goat meat is gradually welcomed by consumers because of its high protein content, low fat and cholesterol content. Moreover, the consumers demand for goat meat and quality requirements continue to rise [14,15]. However, there are few studies on the molecular mechanism of the difference in the deposition of different fat tissues in goats.
Herein, we provided a comprehensive transcriptome pro le on intramuscular and subcutaneous adipocytes in before and after differentiation of Jianzhou Daer goat. The results revealed the expression patterns of mRNAs and lncRNAs, which were important in the developmental stages of two different adipocytes. An integrated analysis of differentially expressed lncRNAs (DELs) and mRNAs (DEGs) was performed, and the lncRNA and its target miRNA/mRNA that could potentially regulate the differentiation of intramuscular and subcutaneous preadipocytes were screened through bioinformatics analysis. At the same time, a lncRNA-miRNA-mRNA interaction network was constructed based on miRNAs known to play a role in adipogenesis. The lncRNA that may be combined with miRNA and mRNA was screened out and veri ed by RT-qPCR technique. The results from the present investigation provided a theoretical basis for in-depth analysis of the regulation mechanism of goat fat cell differentiation and improvement of meat quality in goats.
Finally, 12519 assumed non-coding transcripts were retained by CNCI, CPC2 and PFAM softwares (Fig.   1A), including 16.9% lncRNAs, 76.1% intronic lncRNAs, and 7% anti-sense lncRNAs. In addition, 43767 mRNAs were identi ed. The expression level of each RNA was standardized to fragments per kilobase of exon model per million mapped reads (FPKM), and it was found that a small part of them was not expressed or was expressed at a relatively low level. The expression levels of a large number of RNAs were mainly between 0.1 and 2 of log10 (FPKM + 1), and a small number of genes had very high expression level (log10 |FPKM + 1| > 2) (Fig. 1B). The extremely signi cant differences of genes were helpful to explore the molecular mechanism of fat deposition in different tissues. . Therefore, in order to observe the characteristics of lncRNAs and the difference between lncRNAs and mRNAs, we compared the length of lncRNAs and mRNAs, including the number of exons, and the open reading frame (ORF). The results showed that most of lncRNAs tended to be shorter in length, and they contained less exons than mRNAs (Fig. 1C). Also, the length of ORFs in the lncRNAs was shorter than those of the mRNAs (Fig. 1D & 1E).

Differentially expressed lncRNAs and mRNAs in goat different adipocytes
The correlation coe cient can represent the degree of similarity between samples. We found that the data between SPA and SA, and between IMPA and IMA were highly correlated in expression ( Fig. 2A ). We found 143 lncRNAs and 33258 mRNAs differentially expressed between IMA and IMPA (1451 were upregulated, 1950 were down-regulated), 47 lncRNAs and 1659 mRNAs differentially expressed between SA and SPA (851 were up-regulated, 855 were down-regulated) (Fig. 2B). To further explore the differences in lncRNAs expression between subcutaneous and intramuscular adipocytes at different developmental stages, we performed a clustered heat map on the differentially expressed genes. The results of the cluster analysis showed that distinct lncRNAs expression patterns are associated with the differentiation of both subcutaneous and intramuscular adipocytes in goats (Fig. 2C).

Function prediction of lncRNAs and corresponding genes during intramuscular and subcutaneous preadipocyte differentiation
LncRNAs not only can regulate the expression of neighboring protein-coding genes through a cis mechanism [17,18], but also regulate the expression of genes located on other chromosomes via a trans mechanism [19,20]. In this study, bioinformatic analysis was used to predict the potential target genes of DELs via trans.The absolute value of Pearson's correlation coe cient was set as greater than 0.95. There were 39012 lncRNA-mRNA relationship pairs. Among them, LNC_004374, LNC_007272, LNC_010037, LNC_004066, LNC_010143 and other lncRNAs can target genes related to adipocyte differentiation or lipid metabolism, CD36 [21], FABP3 [22], FGF11 [23], FOXO6 [24], SMAD1 [25], TGFB226, FGFR2 [27]. The results are shown in Table 3.
We drew Venn diagrams (Fig. 3A) to visualize the number of differential lncRNAs that are common versus unique to the two adipocyte differential differentiation processes. Through the use of Goseq, we compare the GO classi cations of DELs (adjusted P < 0.05) [28]. The DELs target genes are enriched in biological process (BP), cellular component (CC) and molecular function (MF)( Table 4). Unique DELs target genes of goat intramuscular adipocytes before and after differentiation are mainly enriched in binding (Fig. 3B), and the unique DELs target genes of subcutaneous adipocytes are mainly enriched in glucose metabolism and chemokine (Fig. 3C), while the shared DELs target genes are mainly enriched in binding (Fig. 3D).
The KEGG annotation results of the differentially expressed target genes of lncRNAs are classi ed according to the pathway types of the KEGG database. As shown in the Fig. 3E, the unique DELs target genes of intramuscular adipocytes before and after differentiation are most signi cantly enriched in hypertrophic cardiomyopathy (HCM), dilated cardiomyopathy, focal adhesion, arrhythmogenic right ventricular cardiomyopathy (ARVC), ECM-receptor interaction, proteoglycans in cancer and cardiac muscle contraction (P < 0.05). In addition, the analyzed genes also was enriched in fat acid metabolism and TGF-beta signaling pathway, which are related to adipocyte formation. The unique DELs target genes of subcutaneous adipocytes before and after differentiation are most signi cantly enriched in biosynthesis of amino acids, glycolysis / Gluconeogenesis, carbon metabolism, valine, leucine and isoleucine degradation, histidine metabolism, fatty acid degradation, melanoma, fatty acid metabolism(P < 0.05), among which fatty acid degradation, fatty acid metabolism, TGF-beta signaling pathway, PI3K-Akt signaling pathway have a signi cant impact on the adipogenesis (Fig. 3F). The shared DELs target genes are most signi cantly enriched in biosynthesis of amino acids, carbon metabolism, glycolysis / gluconeogenesis and alzheimer's disease(P < 0.05) (Fig. 3G).

Construction of Regulation Networks for adipocytes in goats
In the previous study, we obtained 146 differentially expressed miRNAs (DEMs) unique to intramuscular adipocytes, 159 (DEMs) unique to subcutaneous adipocytes, and 83 (DEMs) shared by two adipocytes (Fig. 4A). After targeted binding prediction, 306571 mRNAs and 32,029 lncRNAs were found to have potential targeted binding relationships with DEMs. Combined with the sequenced intramuscular adipocyte-speci c 2724 DEGs and 135 DELs, 6572 miRNA-mRNA relationship pairs and 1240 miRNA-lncRNA relationship pairs were obtained (Fig. 4B). The 1125 DEGs and 39 DELs unique to subcutaneous adipocytes were combined to obtain 3287 miRNA-mRNA relationship pairs and 307 miRNA-lncRNA relationship pairs (Fig. 4B). At the same time, the intersection with a total of 412 DEGs and 8 DELs yielded 1681 miRNA-mRNA relationship pairs and 52 miRNA-lncRNA relationship pairs. The lncRNA-miRNA-mRNA interaction network was constructed by Cytoscape software, and the top 20 miRNAs were selected after calculating the degree value of each factor (Table 5). Among them, we found miR-20, miR-194, miR-335, miR-363, miR-200, miR-199, and miR-302 related to fat formation, thereby, we constructed a ceRNA network view. As shown in the Fig. 5A, the unique lncRNA-miRNA-mRNA network of intramuscular adipocytes includes 30 lncRNAs, 6 miRNAs and 426 mRNAs, and the lncRNA-miRNA-mRNA network unique to subcutaneous adipocytes includes 8 lncRNAs, 4 miRNAs and 182 mRNAs (Fig. 5B). Whereas, two types of adipocyte-common lncRNA-miRNA-mRNA networks include 1 lncRNA, 4 miRNAs and 96 mRNAs (Fig. 5C).
Veri cation of lncRNA expression pro les using qRT-PCR To verify the reliability of RNA-seq results, 7 candidate lncRNAs were randomly selected from the DELs obtained from the screening, and UXT was used as an internal reference for qRT-PCR analysis. The results showed that XR_001917557.1, and LNC_004191 were signi cantly down-regulated in expression in intramuscular adipose tissue. Whereas, the XR_001918647.1, XR_001917728.1 were signi cantly upregulated in expression in intramuscular adipose tissue. XR_001295810.1, XR_001917637.1 and XR_001297263.2 were signi cantly down-regulated expression in subcutaneous adipose tissue, and the LNC_004191 was signi cantly up-regulated expression in subcutaneous adipose tissue, These results were consistent with the trend of RNA-seq, indicating the credibility of the RNA-seq results (Fig. 6).

Discussion
The distribution and deposition of adipose tissue in different parts of the body are the key factors affecting carcass quality and meat avor. Subcutaneous fat mainly affects carcass quality.
Intramuscular fat (IMF) is the material basis of marbling, and an important factor affecting meat avor. A large number of studies have shown that IMF is directly involved in the formation of meat tenderness, juiciness and avor [29,30]. Goat is an indispensable animal in China's agricultural production, and the molecular regulation mechanism of its lipid deposition has not been fully elucidated yet.
LncRNA is a kind of noncoding RNA longer than 200 nt, which has attracted substantial attention in the last few years. Studies have shown that lncRNAs regulate metabolic tissue development and function, including adipogenesis, hepatic lipid metabolism, islet function, and energy balance [31][32][33][34][35]. Despite the fact that many studies have indicated the importance of lncRNAs in different tissues, little is known about their biological function in goat fat deposition, especially in the differentiation of goat intramuscular and subcutaneous preadipocytes. To the best of our knowledge, our study is the rst to screen for lncRNAs and mRNAs regulating goat preadipocyte differentiation by sequencing and annotating the transcriptome of intramuscular and subcutaneous preadipocytes. A total of 1,118,110,544 reads were successfully mapped to the goat reference genome assembly. We identi ed 12,519 lncRNAs.
The average sequence length of lncRNAs was shorter than that of mRNAs, and the number of exons was less than that of mRNAs, with the ORF length being shorter than that of mRNAs. Our results indicated that the predicted lncRNAs were shorter with fewer exons than mRNAs, which are in agreement with the those reproted in previous studies [36][37][38]. The Pearson correlation (R 2 ) of each sample is greater than 0.8, which indicated that our experiment was reliable and the sample selection was reasonable.
LncRNA functions by regulating mRNA. At present, the mechanism of interaction between lncRNA and mRNA is not clear. We predict the biological function of lncRNAs through its co-expression with protein coding genes. Consequently, we found that many target genes of DELs were also differentially expressed in goat intramuscular and subcutaneous preadipocytes. This suggested that lncRNAs may function through complementary target genes, which can play critical roles in the differentiation of goat intramuscular and subcutaneous preadipocytes. For example, SMAD1 is a target gene of the differentially expressed lncRNAs LNC_009792, LNC_007731, LNC_000706, LNC_008467, LNC_006192 and LNC_004878, and it has been reported to regulate the differentiation of preadipocytes [39]. These ndings suggest that these lncRNAs might be involved in the differentiation of intramuscular preadipocytes by affecting the expression of SMAD1. However, there is no DEGs in the differentiation process of subcutaneous preadipocytes in our selected fat development related genes, so we speculate that the network regulating intramuscular adipogenesis is more complex than subcutaneous fat. The higher number of DELs, DEMs and DEGs in IMF than that of subcutaneous fat also supported this, which is consistent with the fact that intramuscular preadipocytes have stronger ability to deposit fat than that of subcutaneous preadipocytes [40].
To explore the similarities and differences of different adipocytes, DELs target genes (IMPA vs IMA and SPA vs SA ) were subjected to GO and KEGG pathway enrichment analyses. We found that few common term was found between the IMPA vs IMA and SPA vs SA comparisons. Several pathways involved in preadipocyte differentiation were previously identi ed, including the TGF-β signaling pathway(IMF and subcutaneous fat) [41], PI3K/AKT signaling pathway (subcutaneous fat) [42], and arachidonic acid metabolism(subcutaneous fat) [43]. For example, the LncRNA GAS5 inhibits lipogenesis in 3T3-L1 cells through the miR-21a-5p / PTEN signaling pathway [44]. FDNCR1 affects porcine lipogenesis by competitively binding miR-204 to regulate the TGF-β pathway [45]. However, for some pathways identi ed here, their involvement in the goat preadipocyte differentiation process is being reported for the rst time.
Interestingly, in the pathway analysis, we found that the components of two pathways, fatty acid metabolism(IMF and subcutaneous fat) and fatty acid degradation (subcutaneous fat), which have been reported to be involved in lipid metabolism, and were enriched in the entire process of differentiation of intramuscular and subcutaneous preadipocytes [46,47]. The common enrichment pathways during differentiation of both adipocytes involve amino acid metabolism, gluconeogenesis and carbohydrate metabolism, suggesting that IMF and subcutaneous fat are largely different in differentiation pathways and lipid metabolism pathways, while there are similarities in communication with other metabolic pathways. In addition to the speci c pathways of the two adipocytes, there are also differences in the common components of the two pathways. Just as IMF is enriched in 19 target genes of TGF-β signaling pathway, 8 target genes of fatty acid metabolism, and subcutaneous fat is enriched in 8 and 7 targat genes. Here, we hypothesize that there are differences in the pathways regulating intramuscular and subcutaneous adipose differentiation, and that there are differences in the downstream target genes in the same pathways, which indirectly demonstrates that gene expression is tissue-speci c in goats.
To date, many genes have been reported to regulate the differentiation of preadipocytes. However, few studies have been conducted on the roles of lncRNAs in intramuscular and subcutaneous preadipocytes differentiation. The molecular and cellular mechanisms regulating goat preadipocytes differentiation are thus still poorly understood. Here, we constructed a miRNA-lncRNA-mRNA interaction network, and calculated the degree (the number of times each factor interacts with other factors) of each factor through centerscape. We selected the top 20 miRNAs in the dgree, and found that miR-20 [48], miR-194 [49], miR-335 [50], miR-363 [51], miR-200 [52], miR-199 [53], and miR-302 [54] were related to fat development, and visualized them in Cytoscape. We identi ed a number of highly connected lncRNAs and mRNAs in the three modules, including the two kinds of adipocytes are unique and common. In conclusion, we rst generated the expression pro les of lncRNAs from intramuscular and subcutaneous adipocytes of Jianzhou Daer goat (IMA vs IMPA and SA vs SPA) based on RNA-Seq technique. We found that the number of lncRNAs regulating IMF differentiation was more than that of subcutaneous adipocytes. Our results suggest that those lncRNAs might play important roles in meat quality. Collectively, this study takes the rst step toward understanding the molecular mechanisms underlying variations in goat meat quality. Also, our results provided a theoretical basis for molecular breeding to improve the meat quality in goats.

Materials And Methods
Experimental animals and sample collection Total RNA extraction and sequencing A total of 20 cell samples were successfully collected. Total RNA was extracted using Trizol reagent (Takara, Dalian, China). RNA quality was determined using NanoPhotometer spectrophotometer and Agilent 2100 bioanalyzer, which analyzes the integrity of the RNA. The lncRNA library is constructed using a chain-speci c library. The method for synthesizing the rst strand of cDNA by reverse transcription is the same as the normal method of NEB library construction. The difference is that when the second strand is synthesized, dTTP in dNTPs is replaced by dUTP. After that, cDNA end repair, A-tailing, ligation of sequencing adapters, and length screening were also performed, and then the second strand of cDNA

Identi cation of lncRNAs and their expression analysis
To ensure the quality of the obtained lncRNAs, three criteria were used to identify the desired lncRNAs in the transcriptome assemblies: 1) transcripts with length > 200 bp and exon number ≥ 2 were selected; 2) Cuffcompare was used to calculate the read coverage of every transcript, and transcripts with an FPKM value(Cuffquant) of more than 0.05 were removed; and 3) the coding potential of the transcripts was predicted using the coding potential calculator (CPC < 1) [69], Coding Noncoding Index (CNCI) [70], and Pfam [71] protein domain families to further remove coding genes. Subsequent steps were performed based on the intersection of the results obtained using the twenty databases.
The FPKMs values of lncRNAs and mRNAs in IMPA, IMA, SPA and SA libraries were calculated by StringTie, balltown and Cuffdiff. The transcripts in the 20 libraries were analyzed for differential expression, and the P value was used to screen the DELs and DEGs between two different adipocyte samples. When P < 0.05, the lncRNAs and mRNAs are considered to be differentially expressed.

Gene functional annotation
The trans function describes the co-expression relationship between lncRNAs and mRNAs. Pearson correlation coe cient (R > 0.95 or R < 0.95) was calculated by custom scripts. David was used to cluster the target genes among 20 samples for functional enrichment analysis of lncRNA target genes [72]. Pvalue < 0.05 was set as signi cant threshold. Go (gene ontology) can enrich and analyze the target genes of DELs. Goseq is used to enrich and analyze the target genes of DELs. P < 0.05 is considered to be the signi cant enrichment of GEGs. KEGG is a database (http://www.genome.jp/kegg/) for understanding the advanced functions and utilities of biological systems such as cells, organisms, and ecosystems. We used KOBAS software to detect the pathway enrichment analysis of DELs target genes in the KEGG pathway.
Download of miRNA sequencing data and construction of lncRNA-miRNA-mRNA networks The raw data were processed using our previous sequencing to remove adaptor dimers, garbage, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeats. Then we used miRBase to de ne known miRNAs. Finally, unmapped sequences were predicted using RNAfold software. The edgeR software [73] was used to calculate DEMs, in two adipocytes. The screening criterion for signi cant differences is P-value < 0.05. miRWalk, The miRanda, RNAhybrid, and Targetscan were used for screening. Then, the intersections of co-expressed lncRNAs, miRNAs, and mRNAs in the 4 programs were used to construct a lncRNA-miRNA-mRNA network.
DELs and DEGs were combined with the target lncRNAs and mRNAs of DEMs respectively to obtain the miRNA-mRNA relationship pair and miRNA-lnRNA relationship pair. Finally, the lncRNA-miRNA-mRNA network was visualized using Cytoscape v3.7 software.

Validation of gene expression of by RT-qPCR technique
Primers were designed using Primer-BLAST on the NCBI website (Table 6). The rst cDNA strand was synthesized using the RevertAid First Strand cDNA Synthesis Kit (Thermo, Waltham, MA, USA), in accordance with the user manual. Ubiquitously expressed transcript gene (UXT) was used as a housekeeping gene [74]. The qPCR reaction procedure was composed of four steps, including predegeneration (95℃, 3 min), degeneration (95℃, 10 s), annealing (60℃, 10 s) and extension (72℃, 15 s), of which degeneration, annealing and extension were running for 40 cycles. Quanti cation of selected gene expression was performed using the comparative threshold cycle (2 −ΔΔCT ) method [75].
Declarations and Use Committee of Southwest Minzu University (No. 18032). This study was carried out in compliance with the ARRIVE guidelines. All goats were humanely sacri ced and all efforts were made to minimize the suffering.

Consent for publication
Not applicable.
Identi cation and characterization of lncRNAs in goat, Capra hircus. (A) Venn diagrams of coding potential analysis by using stringent criteria. Three tools (CPC2, CNCI and PFAM) were employed to analyze the coding potential of lncRNAs. Those simultaneously shared by three analytical tools were   enriched for differentially expressed lncRNAs in both intramuscular and subcutaneous adipocytes before and after differentiation. The abscissa represents the richness factor, and the ordinate represents the enriched pathway terms. The Q-value represents the corrected P value.   Validation of lncRNAs in intramuscular and subcutaneous adipocytes using RT-qPCR technique. Data were analyzed by the 2-ΔΔCt method where the UXT was used as a reference gene. Each column represents the mean ±SE. White bars represent read from RNA-Seq, and blue bars represent the results from qRT-PCR analysis.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.