Identification of Long Non-Coding RNA-Associated Competing Endogenous RNA Network in the Differentiation of Chicken Preadipocytes

Emerging evidence indicates that long noncoding RNAs (lncRNAs) play important roles in the regulation of cell differentiation by acting as competing endogenous RNA (ceRNA). However, the regulatory mechanisms of lncRNA and the lncRNA-associated ceRNA network involved in adipogenic differentiation of chicken preadipocytes remain elusive. Here, we first constructed the chicken preadipocyte in vitro induction model. Then, we identified differentially expressed lncRNAs (DELs), miRNAs (DEMis), and mRNAs (DEMs) between differentiated and undifferentiated preadipocytes. Furthermore, we constructed the lncRNA associated ceRNA network by gene expression correlation analysis and target prediction of DELs, DEMis, and DEMs. Finally, we determined twelve candidate lncRNA-miRNA-mRNA interactions from the lncRNA associated ceRNA network. Eight out of the twelve interactions were validated by RT-qPCR, indicating their potential role in the regulation of chicken preadipocytes differentiation. Among the eight interactions, TCONS_00026544-gga-miR-128-1-5p-RASD1, TCONS_00055280-gga-miR-135a-5p-JAM3, TCONS_00055280-gga-miR-135a-5p-GPR133, TCONS_00055280-gga-miR-135a-5p-CLDN1, and TCONS_00055280-gga-miR-135a-5p-TMEM123 may promote adipogenic differentiation of chicken preadipocytes while TCONS_00057272-gga-miR-146a-3p-FOXO6, TCONS_00057242-gga-miR-6615-3p-FOXO6, and TCONS_00057242-gga-miR-6615-3p-ENSGALT00000043224 have the opposite effects. Our results not only provide novel insights into ceRNA roles of lncRNAs in chicken preadipocytes differentiation and but also contribute to a better understanding of chicken fat deposition.


Introduction
In the past few decades, the growth rate of broilers has been significantly improved due to continuous selection [1]. However, the high growth rate is accompanied by excessive fat deposition [2]. Excessive fat deposition (mainly abdominal fat) could influence the laying rate, carcass yield, feed conversion ratio, hatching rate, and fertility rate. Therefore, more attention is paid to reducing abdominal fat, which has become a major breeding goal in broilers [3]. In addition, study showed that the adipose tissue of chicken has a relative resistance to insulin, and chicken could be used as an animal model to study human type 2 diabetes [4]. Thus, understanding the biological basis of the development of chicken adipose tissue can benefit the development of poultry industry and also provide new idea for human biomedical research. In mammals, the molecular and cellular mechanisms of adipose tissue development have been well studied. However, the mechanisms underlying chicken fat tissue development are complex and remain elusive.
Long non-coding RNAs (lncRNAs) are a class of non-coding RNA greater than 200 nucleotides in length. MicroRNAs (miRNAs) are a class of non-coding single-stranded RNA of approximately 22 nucleotides in length encoded by endogenous genes. miRNAs play essential roles in various life activities, including metabolism [5], development, and disease [6]. Chicken abdominal fat tissue development is a result of the increase in cell number (hyperplasia) and the increase in cell size (hypertrophy or preadipocytes differentiation). Recent studies have shown that LncRNAs and miRNAs participate in regulating preadipocytes differentiation in mammals [7,8]. Their roles in the differentiation of chicken preadipocytes are still poorly understood.
The competitive endogenous RNA (ceRNA) hypothesis proposes that transcripts can regulate each other at the post-transcription level by competing for shared miRNAs [9]. lncRNA can act as a miRNA sponge via competing endogenous RNA (ceRNA) activity [10], thereby regulating the gene expression of miRNA [11]. Understanding the RNA interaction will lead to significant insights into gene regulatory network in the differentiation of chicken preadipocytes. The role of lncRNA as ceRNA in adipocytes differentiation of mammal is gradually becoming clear. Liu et al. [12] found that lncRNA Gm15290 sponges miR-27b to promote PPARγ-induced fat deposition in mice. Li et al. [13,14] confirmed that lncRNA GAS5 negatively regulates the adipogenic differentiation of MSCs and 3T3-L1 cells by modulating miR-18a and miR-21a-5p as ceRNAs. Li et al. [15] reported that lncRNA ADNCR suppresses adipogenic differentiation by targeting miR-204 in bovine. However, the regulatory mechanisms of lncRNA as ceRNA and the lncRNA-associated ceRNA network involved in adipogenic differentiation of chicken preadipocytes remain elusive.
To identify lncRNAs and lncRNA-associated ceRNA network involved in differentiation of chicken preadipocytes, we systematically analyzed the differentially expressed lncRNAs, miRNAs, and mRNAs using RNA-seq and miRNA-seq. We then constructed the lncRNA-associated ceRNA network following the ceRNA hypothesis. Crucial ceRNA interactions were identified based on the interacted miRNAs function and validated using quantitative reverse transcription polymerase chain reaction (RT-qPCR). Our results not only provide novel insights into ceRNA roles of lncRNAs in chicken preadipocytes differentiation and but also contribute to a better understanding of chicken fat deposition.

Ethics Statement
This experiment was performed in accordance with Chinese guidelines for animal welfare, and the animal protocol was approved by the animal welfare committee of Yangzhou University (permit number SYXK (Su) 2012-0029).

Cell Culture and Induction
Three Haiyang Yellow chickens were anaesthetized with sodium pentobarbital and killed at two weeks old. Primary preadipocytes were isolated from the abdominal fat tissue following methods described by Shang et al. [16]. Complete medium (90% DMEM/F12 medium,10% fetal bovine serum, 100 units/mL penicillin, and 100 µg/mL streptomycin) was used to culture chicken preadipocytes in a cell incubator with 5% CO 2 at 37 • C. when the cell reached 90% confluence, the differentiation medium (supplemented with oleic acid (400 µmol/L)) was used to induce the differentiation of chicken preadipocytes. Undifferentiated (0 days) and differentiated (6 days) cells were collected. Each interval was carried out for three biological replicates. The differentiation of chicken preadipocytes was determined using RT-qPCR and Oil red O staining. After induction for six days, the cells were fixed in 10% formalin for 30 min at room temperature. Then the cells were washed using PBS buffer for three times and stained using Oil red O working solution for eight minutes. Finally, the cells were washed using PBS buffer and observed using a microscope. The expression level of adipocyte marker genes PPARG and FABP4 in undifferentiated and differentiated preadipocytes was determined using RT-qPCR. The HiScript II Q RT SuperMix for qPCR (+gDNA wiper) kit (Vazyme, Nanjing, China) was used to synthesize the first-strand cDNA following the manufacturer's protocol. The AceQ qPCR SYBR Green Master Mix kit (Vazyme, Nanjing, China) was used to perform RT-qPCR following the manufacturer's protocol.

RNA Extraction and Library Construction
Trizol reagent (Tiangen, Beijing, China) was used to extract total RNA from undifferentiated and differentiated chicken preadipocytes (Add 1 mL of Trizol per 2.5 × 10 6 cells). The NEBNext Poly(A) mRNA Magnetic Isolation Module kit (New England Biolabs, Ipswich, MA, USA) and NEBNext Small RNA Library Prep Set for Illumina (Multiplex Compatible) kit (New England Biolabs, Ipswich, MA, USA) were used to construct the lncRNA libraries and miRNA libraries following the manufacturer's protocol, respectively. Six lncRNA libraries (three for each group) and six miRNA libraries (three for each group) were constructed. The 12 constructed libraries were sequenced using the Illumina HiSeqTM 4000 sequencing platform (Illumina Inc., San Diego, CA, USA).

LncRNA Identification
The raw data were filtered by removing low-quality reads (Q-value ≤ 20), reads containing adapters, reads that are all A bases, and reads containing unknown nucleotides ratio greater than 10%. The clean data were aligned with rRNA sequences of chicken using Bowtie2 to remove reads mapped to rRNA [17]. Then the high-quality clean data were mapped to the Gallus gallus reference genome using Tophat2 (version 2.0.3.12) to identify novel transcripts [18] and known lncRNAs. New lncRNAs were predicted based on the novel transcripts using CPC [19], CNCI [20], and SwissProt database [21]. The arguments and the thresholds used to identify lncRNAs were specified in Data S1.

miRNA Identification
We filter the raw data by removing the low-quality reads (Q-value ≤ 20), reads containing adapters, reads shorter than 18 nt, and reads with polyA. Then reads were aligned with small RNAs in the GenBank [22] and Rfam [23] database using blastall software (version 2.2.25). The clean data was obtained by removing reads that are more than 97% identical to rRNA, scRNA, snoRNA, snRNA, and tRNA. The bowtie (1.1.2) software was used to identify known miRNAs by aligning clean data with the miRbase database. Furthermore, the MIREAP_v0.2 software (http://sourceforge.net/projects/mireap/) was used to predict novel miRNAs. The arguments and the thresholds used to identify miRNAs were specified in Data S1.

Identification of Differentially Expressed Genes
Fragments per kilobase of transcript per million mapped reads (FPKM) normalization method was used to quantify the expression levels of lncRNAs and mRNAs. Tags per million (TPM) normalization method was used to quantify the expression levels of miRNAs. The edgeR [24] software was used to identify differentially expressed lncRNAs (DELs), mRNAs (DEMs), and miRNAs (DEMis). Genes with FDR < 0.05 and |log2FC| > 1 were considered as significantly differentially expressed (DEGs). First, we predicted the DEMis targets within DELs and DEMs using MIREAP_v0.2, miRanda [25], and TargetScan [26] software. Secondly, we calculated the expression correlation between DEMis and its targets (lncRNA-miRNA or miRNA-mRNA) using the Spearman rank correlation coefficient (SCC) [27]. Pairs with SCC smaller than −0.7 were selected as candidate lncRNA-miRNA or miRNA-mRNA pairs. Thirdly, the expression correlation between DELs and DEMs was calculated using SCC, and pairs with SCC greater than 0.9 were selected as candidate ceRNA (lncRNA-mRNA) pairs. At last, the hypergeometric cumulative distribution function test in R software was used to identify final ceRNA pairs (p-value < 0.05). The lncRNA-miRNA-mRNA network was visualized using Cytoscape software v3.6.0 (http://www.cytoscape.org/) based on the target prediction, expression correlation, and hypergeometric cumulative distribution test results. The NetworkAnalyzer plug-in in Cytoscape software was used to calculate the connection degree. Nodes with connection degree greater than average degree of the whole ceRNA network were identified as highly-connected genes. The hypergeometric cumulative distribution function is as follow: n is the number of miRNA that is shared by lncRNA and mRNA, U is the number of all miRNAs, M is the number of target miRNAs of lncRNAs, and N is the number of target miRNAs of mRNAs.

Functional Enrichment Analysis
To characterize the underlying function of the network, DEMs involved in the lncRNA-associated ceRNA network were performed to Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using clusterProfiler package [28] in R software based on the GO (http://www.geneontology.org/) and KO (http://www.kegg.jp/kegg/kegg1.html) databases. Gallus gallus was selected as the reference species for enrichment analysis. GO terms and KEGG pathways with p-value < 0.05 were considered as significantly enriched.

Gene
Primer Sequences

Statistical Analysis
The SPSS for Windows software (version 22, SPSS, Inc.) and Excel software (Microsoft Corp.) were used to analyze experimental data. The SPSS was used to test the normal distribution of the RT-qPCR data. The RT-qPCR data conformed to the normal distribution (Data S2). Student t-test was used to analyze the difference in RT-qPCR data. All tests were performed at least in triplicate. The p < 0.05 value was considered to indicate a statistically significant difference.

Induction of Chicken Preadipocytes
The cells were induced to differentiation with differentiation medium supplemented with oleic acid. The oil red o staining and RT-qPCR were used to determine cell differentiation. We found that, after six days of induction, preadipocytes were fully differentiated and filled with large lipid droplets ( Figure 1A,B). The relative expression level of adipocyte marker genes FABP4 and PPARG significantly increased ( Figure 1C). The relative expression level of adipocyte marker genes before and after differentiation. The relative expression level of FABP4 and PPARG genes significantly increased post-differentiation, indicating that the chicken preadipocytes were successfully induced. *denotes p < 0.05, **denotes p < 0.01.

Identification of Differentially Expressed lncRNAs, miRNAs, and mRNAs
A total of 3881 lncRNAs were identified in chicken preadipocytes (Data S3). A correlation heat map was used to show the relationship between samples, which indicates that samples in the same group were highly correlated ( Figure 2A). Differentially expressed genes between undifferentiated and differentiated preadipocytes were identified using edgeR software. Consequently, 235 DELs, 145 DEMis, and 660 DEMs were identified on the basis of RNA-Seq and small RNA-seq data (Table S1). Heatmaps generated from the expression of DELs, DEMis, and DEMs were used to show the expression patterns of these genes between undifferentiated and differentiated preadipocytes ( Figure  2B-D).

Identification of Differentially Expressed lncRNAs, miRNAs, and mRNAs
A total of 3881 lncRNAs were identified in chicken preadipocytes (Data S3). A correlation heat map was used to show the relationship between samples, which indicates that samples in the same group were highly correlated (Figure 2A). Differentially expressed genes between undifferentiated and differentiated preadipocytes were identified using edgeR software. Consequently, 235 DELs, 145 DEMis, and 660 DEMs were identified on the basis of RNA-Seq and small RNA-seq data (Table S1). Heatmaps generated from the expression of DELs, DEMis, and DEMs were used to show the expression patterns of these genes between undifferentiated and differentiated preadipocytes ( Figure 2B-D).

Construction of ceRNA Network Related to Preadipocytes Differentiation
By the microRNA response elements (MREs) prediction, we obtained 2262 interactions between 145 DEMis and 232 DELs, as well as 999 interactions between 139 DEMis and 192 DEMs. Based on SCC, 438 interactions between DEMis and DELs, 297 interactions between DEMis and DEMs, and 2883 interactions between DELs and DEMs were identified (Table S2). Combined the MREs prediction and SCC results, DELs and DEMs that shared common DEMis were selected as candidate ceRNAs and tested using hypergeometric cumulative distribution function. In total, 251 ceRNAs (lncRNA-miRNA-mRNA) (p < 0.05) of 96 lncRNAs, 71 miRNAs, and 84 mRNAs were identified

Construction of ceRNA Network Related to Preadipocytes Differentiation
By the microRNA response elements (MREs) prediction, we obtained 2262 interactions between 145 DEMis and 232 DELs, as well as 999 interactions between 139 DEMis and 192 DEMs. Based on SCC, 438 interactions between DEMis and DELs, 297 interactions between DEMis and DEMs, and 2883 interactions between DELs and DEMs were identified (Table S2). Combined the MREs prediction and SCC results, DELs and DEMs that shared common DEMis were selected as candidate ceRNAs and tested using hypergeometric cumulative distribution function. In total, 251 ceRNAs (lncRNA-miRNA-mRNA) (p < 0.05) of 96 lncRNAs, 71 miRNAs, and 84 mRNAs were identified (Table S3). The constructed lncRNA-miRNA-mRNA network was visualized using the Cytoscape software ( Figure 3). Network analysis showed that the average degree of the whole network was 4.9.
Genes 2019, 10, 795 8 of 16 (Table S3). The constructed lncRNA-miRNA-mRNA network was visualized using the Cytoscape software ( Figure 3). Network analysis showed that the average degree of the whole network was 4.9. Figure 3. Visualization of the lncRNA-miRNA-mRNA network. Green nodes indicate lncRNAs, red nodes indicate miRNAs, and blue nodes indicate mRNAs, respectively. Node size is proportional to the connection degree. Connection degree shows the number of connected nodes with the individual node. The higher the degree of a node, the more important the node in the network. Hub genes are the nodes with higher degree i.e., nodes with more connections (degree > 4.9).

Functional Enrichment Analyses
To further understand the underlying function of lncRNA-associated ceRNA network, GO and KEGG pathway analyses were performed on mRNAs involved in the ceRNA network. GO analyses showed that genes involved in the ceRNA network are associated with response to external stimulus, positive regulation of cell differentiation, and positive regulation of multicellular organismal biological processes (Figure 4). KEGG pathway analyses identified four significant enriched pathways, including pantothenate and CoA biosynthesis, mineral absorption, steroid hormone biosynthesis, and vascular smooth muscle contraction ( Figure 5). Green nodes indicate lncRNAs, red nodes indicate miRNAs, and blue nodes indicate mRNAs, respectively. Node size is proportional to the connection degree. Connection degree shows the number of connected nodes with the individual node. The higher the degree of a node, the more important the node in the network. Hub genes are the nodes with higher degree i.e., nodes with more connections (degree > 4.9).

Functional Enrichment Analyses
To further understand the underlying function of lncRNA-associated ceRNA network, GO and KEGG pathway analyses were performed on mRNAs involved in the ceRNA network. GO analyses showed that genes involved in the ceRNA network are associated with response to external stimulus, positive regulation of cell differentiation, and positive regulation of multicellular organismal biological processes (Figure 4). KEGG pathway analyses identified four significant enriched pathways, including pantothenate and CoA biosynthesis, mineral absorption, steroid hormone biosynthesis, and vascular smooth muscle contraction ( Figure 5).

Identification and RT-qPCR Validation of Crucial Genes
In our study, miRNAs with TPM greater than 4.5 (more than two-thirds of the DEMis expression is lower than 4.5) were identified as highly expressed. Genes with connection degree greater than 4.9 were identified as highly connected. gga-miR-135a and gga-miR-6615 were high-connected in the network (degree > 4.9) and high-expressed in chicken preadipocytes (TPM > 4.5). gga-miR-146a, gga-miR-135a, and gga-miR-128-1 were reported to participate in the regulation of 3T3-L1 differentiation [29][30][31], suggesting their potential roles in chicken preadipocytes differentiation. Thus, ceRNA

Identification and RT-qPCR Validation of Crucial Genes
In our study, miRNAs with TPM greater than 4.5 (more than two-thirds of the DEMis expression is lower than 4.5) were identified as highly expressed. Genes with connection degree greater than 4.9 were identified as highly connected. gga-miR-135a and gga-miR-6615 were high-connected in the network (degree > 4.9) and high-expressed in chicken preadipocytes (TPM > 4.5). gga-miR-146a, gga-miR-135a, and gga-miR-128-1 were reported to participate in the regulation of 3T3-L1 differentiation [29][30][31], suggesting their potential roles in chicken preadipocytes differentiation. Thus, ceRNA

Discussion
In mammals, the molecular and cellular mechanisms underlying adipose tissue development have been well studied. The development of adipose tissue is a result of both increased number of new adipocytes (hyperplasia) and increased deposition of lipid in adipocytes (hypertrophy) [32]. However, adipose tissue development in chickens is still poorly understood. lncRNA was originally thought to be noise of RNA polymerase II transcription because it cannot encode a protein [33]. In recent years, with the development of high-throughput sequencing technology, more and more lncRNAs have been annotated. Increasing evidence has confirmed that lncRNAs are widely involved in the regulation of gene expression [34], cell proliferation [35], cell differentiation [36], cell apoptosis [37], and cancer development [38].
Recent studies suggest a significant number of lncRNAs play vital roles in regulating adipogenesis. For instance, lncRNA Lnc-U90926 attenuates 3T3-L1 adipocyte differentiation via inhibiting the transactivation of PPARγ2 or PPARγ [39]. PU.1 AS lncRNA, a novel AS lncRNA transcripted from the porcine PU.1 gene promotes adipogenesis through the formation of a sense-antisense RNA duplex with PU.1 mRNA [7]. However, the roles of lncRNAs in chicken adipogenesis remain elusive. To investigate the function of lncRNAs in chicken adipogenesis, we firstly isolated and cultured the preadipocytes from chicken abdominal fat tissue. Chicken preadipocytes were induced to differentiation in vitro using differentiation medium supplemented with oleic acid. Oil red O staining and RT-qPCR of adipocyte differentiation marker genes indicated that the chicken preadipocyte in vitro induction model was successfully constructed. We characterised the expression profiles of lncRNAs during the differentiation of chicken preadipocytes using Ribo-zero RNA-seq based on the in vitro induction model. A total of 3874 lncRNAs were identified, in which 235 lncRNAs were differentially expressed in the differentiation of chicken preadipocytes. Of the 235 DELs, 169 lncRNAs were up-regulated after differentiation, suggesting that most lncRNAs might play positive roles in the regulation of chicken preadipocytes differentiation. We also identified 145 DEMis and 660 DEMs between undifferentiated and differentiated preadipocytes. We found that top DEMs, such as FABP4 [40] and KLF5 [41] are involved in adipose tissue development of chicken. In the DEMis, gga-miR-128 [42] and gga-miR-223 [43] participate in regulating chicken preadipocytes differentiation. These differentially expressed genes were considered as potential key regulators in the differentiation of chicken preadipocytes.
Growing evidence suggests that lncRNAs can act as ceRNAs to indirectly regulate mRNAs by competitively binding miRNAs [10,15,44]. lncRNA-associated ceRNA networks in human adipocyte differentiation have been constructed [45,46], and lncRNAs regulating adipocyte differentiation as ceRNA in human [13], bovine [15], and mouse [47] have been identified. Nevertheless, the specific lncRNA-associated ceRNA network involved in the differentiation of chicken preadipocytes is yet largely unknown. In this study, we mainly focused on the construction of lncRNA-associated ceRNA network related to chicken preadipocytes differentiation. We firstly identified 235 DELs, 145 DEMis, and 660 DEMs between differentiated and undifferentiated preadipocytes. By target prediction and expression correlation analysis, 251 ceRNAs (lncRNA-miRNA-mRNA) (p < 0.05) of 96 DELs, 71 DEMis, and 84 DEMs were identified. We then constructed a ceRNA network based on these 251 ceRNA triples. To characterize the underlying function of the ceRNA network, we performed GO and KEGG pathway analysis on the 84 DEMs. GO analysis results revealed that DEMs were significantly enriched in positive regulation of cell differentiation biological process, suggesting the potential role of the constructed ceRNA network in cell differentiation. Pathway analysis also showed that DEMs were enriched in pathways related to preadipocyte differentiation, including pantothenate and CoA biosynthesis and steroid hormone biosynthesis.

Conflicts of Interest:
The authors declare no conflict of interest.