Identification of Long Non-coding RNAs Expressed During Early Adipogenesis

Corresponding Author: Rongkun Shen Department of Biology, The College at Brockport, State University of New York, Brockport, NY 14420, USA Tell: 585-395-2808 Fax: 585-395-2741 Email: rshen@brockport.edu Abstract: Long non-coding RNAs (lncRNAs) are transcripts with a length of more than 200 nucleotides that lack protein-coding capacity. LncRNAs play important roles in the regulation of genes during many biological processes, development and disease progression, such as obesity. Obesity is an increasing health concern around the world. Although multiple studies linking lncRNAs and fully differentiated adipocytes have been published, the systematic analysis of those lncRNAs involved in early preadipocyte differentiation, when fate determination decisions are made, has not been reported. In order to fill this gap, we conducted strandspecific RNA-Seq on mouse 3T3-L1 preadipocyte cells and compared the expression profiles before and after early differentiation for 2 days. We identified 82 lncRNAs that significantly changed their expression after early differentiation, 98% of which were newly discovered in association with early adipogenesis. The most remarkable lncRNAs were U90926, Wincr1, Kcnq1ot1, Malat1 and Hotairm1. Expression patterns for these identified genes were also highly correlated with markers of adipogenesis, including PPAR, CEBP, FABP4 and FASN. Further analysis implies that many of those altered lncRNAs might coordinately inhibit Wnt/-catenin signaling pathway to promote preadipocyte differentiation. This work strongly suggests that lncRNAs are critical in the proliferation and differentiation of adipose tissue.


Introduction
Long non-coding RNAs (lncRNAs) are transcripts possessing longer than 200 nucleotides that do not code for proteins. Most lncRNAs are transcribed by RNA polymerase II in the same way as coding genes and undergo posttranscriptional processing, such as splicing, 5'-capping and polyadenylation at 3'-end (Quinn and Chang, 2016). Poly(A) tails provide the premise for purifying those transcripts and converting mRNAs into cDNAs for sequencing library construction. LncRNAs, together with other non-coding RNAs, e.g., microRNAs, used to be regarded as expression noise. In the recent decades of the genome era, however, many lncRNAs were discovered and found to play important roles in gene regulation, cell differentiation, development and disease progression (Zarkou et al., 2018). LncRNAs may function in multiple diverse ways. Some lncRNAs may serve to imprint specific genomic loci. For example, KCNQ1 overlapping transcript 1 (Kcnq1ot1), is a lncRNA located on the opposite strand of potassium voltage-gated channel, subfamily Q, member 1 (KCNQ1). Kcnq1ot1 is involved in imprinting a cluster of genes, including KCNQ1, on the sense-strand (Mancini-Dinardo et al., 2006). Other lncRNAs are involved in epigenetic regulation, such as HOX Transcript Antisense RNA (HOTAIR). HOTAIR is transcribed from HOXC locus, but represses the HOXD locus (~40kb) by binding to the Polycomb Repressive Complex 2 (PRC2) and altering chromatin modifications (Tsai et al., 2010). Finally, aberrant lncRNA expression has been closely associated with disease. For example, metastasis associated lung adenocarcinoma transcript 1 (Malat1) was found to be highly up-regulated during metastasis of early-stage nonsmall cell lung cancer and it is used as a biomarker for prognosis of lung cancer (Ji et al., 2013).
Recently, the role of lncRNAs in adipogenesis and the development of obesity is of great interest (Sun et al., 2013;Wei et al., 2016). While obesity has been a worldwide health issue for several decades, still about 40% of adults in the United States were classified as obese (Hales et al., 2017). Obesity rates in the U.S. dramatically increase yearly based on a study from 1986 to 2010 (Sturm and Hattori, 2013). Moreover, the obesity were associated with increased mortality rate (Flegal et al., 2013). Zeng et al. (2018) reviewed the role of some individual lncRNAs in regulating adipogenesis including HOTAIR and U90926. HOTAIR promotes preadipocyte differentiation by enhancing the expression of AdipoQ, FABP4 and PPAR genes while the upregulation of U90926 inhibits adipogenesis by decreasing the expression of those same genes. Knoll et al. (2015) reviewed the studies of lncRNAs as regulators of the endocrine system including adipose tissue. It contained lncRAP from Sun et al. (2013) and several other lncRNAs, such as Sra1 and Blnc1. Sra1 was turned out to be a coding gene in the later annotation (Frankish et al., 2019). Wang et al. (2019) reported that lnc-OAD (lncRNA associated with osteoblast and adipocyte differentiation, from 1700018A04Rik gene) was upregulated during the differentiation of mouse 3T3-L1 preadipocytes (Wang et al., 2019). Its knockdown inhibited the cell differentiation through promoting -catenin expression and suppressed the expression of PPAR and CEBP genes, i.e., the key marker genes of adipogenesis. Sun et al. (2013) systematically studied lncRNAs during adipogenesis using RNA-Seq and microarray; they identified 175 lncRNAs that are regulated during adipogenesis between preadipocytes and mature adipocyte after 4-day differentiation (Sun et al., 2013). Existing lncRNA studies focusing on adipogenesis utilize samples during their transition from preadipocytes to mature adipocytes usually after differentiation induction for more than 4 days, as these cells have the greater capacity for lipid accumulation (Sun et al., 2013). We are interested in the decision-making signals that initiate preadipocytes to proliferate and/or differentiate, since expansion of adipose tissue is a hallmark of obesity (Rosen and Spiegelman, 2014).
On a genome-wide scale, many genes are working together forming gene regulatory networks and their control is critical for biological processes (Imani and Braga-Neto, 2019a;2019b;Imani et al., 2019). Those networks usually contain the coding genes, which is already of great complexity. With the lncRNAs added to the network, the study becomes even more complicated. For example, how the lncRNAs function within a known pathway and regulate it needs to be explored.
This study focuses on early initiation of adipogenesis in the 3T3-L1 mouse cell line. These cells are committed to the adipocyte lineage, yet can be cultured in a pre-determined state. When adipogenesis is initiated with insulin, dexamethasone and isobutylmethylxanthine, lipid droplets appear at about the 5th day and most experiments are performed around Day 10-14 (Green and Meuth, 1974). Interestingly, Day 2 differentiating 3T3-L1 preadipocytes form a transient primary cilium that may inhibit Wnt signaling, allowing adipogenesis to continue (Ross et al., 2000). We induced differentiation of 3T3-L1 preadipocytes and conducted strand-specific RNA-Seq to compare lncRNA expression profiles between Day 0 and Day 2 (after 2-day differentiation induction). We identified 82 lncRNAs that were significantly changed in their expression levels after differentiation induction. Of these, 98% are newly associated with early adipogenesis and had not been previously reported. This work strongly suggests that lncRNAs play a significant role in the expansion and differentiation of adipose tissue.

Preadipocyte Differentiation
3T3-L1 cell cultures with <10 passages were grown to 100% confluence. Day_0 samples were harvested at this point in triplicates. Media in Day_2 cultures were replaced with differentiation media consisting of DMEM, 10% fetal bovine serum, 10 g/ml insulin, 400 ng/ml dexamethasone and 115 g/ml isobutylmethylxanthene. Day_2 samples were incubated for 2 days and harvested in triplicates.

RNA Extraction and RNA-Seq Sample Preparation
Total RNAs from 3T3-L1 preadipocyte cells, harvested at Day_0 and Day_2, were extracted with a RNeasy Mini Kit (Qiagen) using a QIAshredder (Qiagen). Total RNAs were stored in -80C freezer. RNA samples were checked for quality with an Agilent 2100 Bioanalyzer at The University of Rochester Functional Genomics Center and all RNA samples had a RIN number of >= 9.3.
Total RNAs were sent to the Center for Genome Research and Biocomputing (CGRB) at Oregon State University for library construction and sequencing. Libraries were prepared at CGRB using a WaferGen PrepXTM mRNA Strand-Specific Library Preparation Protocol (WaferGen Biosystems). In short, mRNAs were purified using the WaferGen Robotic Poly(A) Selection approach and then further fragmented by RNase III. The adapters were ligated one at a time on the mRNA fragment, first the 3' adaptor and then the 5' adaptor. The 3' adaptor was ligated with an ATP-independent RNA ligase followed by blocking the 3' end with a primer to make it double-stranded. The 5' adaptor was ligated next with an ATP-dependent ligase. This assures that minimal adapter dimers form and the adapters are ligated in a directional manner. The 1st sequence read should be from the 5' to 3' direction. That is, all the libraries are strand-specific and the directions of reads are identical to the original direction of transcripts or mRNAs. All libraries were barcoded and sequenced in one lane on an Illumina HiSeq3000 platform with 151bp single-end reads.

Bioinformatic Analysis
The sequencing reads from the Illumina HiSeq3000 were first analyzed with FastQC (Andrews, 2010) for quality control. For all the reads, based on the quality scores, the leading and tailing low quality stretches were trimmed, keeping a 100bp stretch from positions 7 to 106. Those trimmed reads were mapped to the mouse 10 mm assembly using STAR v2.5.2a (Dobin et al., 2013) with the following parameters: --outSAMtype BAM SortedByCoordinate --outSAMstrandField intronMotif --outFilterMismatchNmax 3 --alignIntronMax 300000 --alignSJDBoverhangMin 1 --alignEndsType EndToEnd --chimSegmentMin 2. Reads for each gene were counted using HTSeq-count (Anders et al., 2015) with the default parameters including strand-specific and union modes and the corresponding annotations. The annotations used in this study are UCSC RefSeq (version 20170804, https://genome.ucsc.edu) for RefSeq genes and GENCODE mouse vM22 long non-coding RNA annotation (Frankish et al., 2019). Gene differential expression analyses were conducted with DESeq2 (Love et al., 2014). We used prefiltering criteria of >=10 reads summing the counts of all samples for each gene. All data were normalized by a medium-to-ratio approach and dispersion was estimated by DESeq2. Data were fitted in a negative binomial generalized linear model. Wald statistics were used to determine significantly changed genes. For functional analysis, we extracted promoter regions of significantly altered lncRNAs with ±5 kb of transcription start sites and formatted them into BED files (UCSC Genome Browser). The regions were analyzed using the GREAT web tool (McLean et al., 2010). All the above steps were summarized in the workflow in Fig.  1. Data were visualized by UCSC Genome Browser (Casper et al., 2018) and Integrative Genomics Viewer (IGV) (Robinson et al., 2011).

Global Gene Expression Altered during Early Preadipocyte Differentiation
Unlike previous studies comparing the preadipocytes and differentiated adipocytes (usually more than 4 days after differentiation induction), we focused on early differentiation at Day 2. Morphological changes, such as lipid accumulation, can be readily observed at 4 days after differentiation, but not at Day 2. In order to identify the genes that significantly changed their expression levels, we extracted RNAs from Day_0 and Day_2 (2-day differentiation) 3T3-L1 cells in biological triplicates and made strand-specific sequencing libraries using an oligo-dT probe approach to capture both mRNAs and lncRNAs. After sequencing and demultiplexing, we obtained a total number of reads ranging from 31 to 48 million for the 6 samples. After a careful quality check with FastQC (Andrews et al., 2010), we trimmed all reads and kept the bases from 7 to 106 for each read. Reads with low overall sequencing quality were filtered and the remaining reads were mapped against the mouse GRCm38 mm10 assembly (Fig. 1). Eighty-five percent of the reads in each sample were mapped to the mouse genome. Since the libraries were constructed as strand-specific, it facilitates the assignment of reads to the original transcript based on their directions. For example, in Fig. 2A, the Mkrn2 and Raf1 genes are overlapping at their 3'-ends. We can clearly distinguish the red reads belonging to the Mkrn2 gene and the blue reads belonging to the Raf1 gene, including the reads at splice junctions. Another example is shown in Fig. 2B, where the proteincoding gene Bach2 and lncRNA Bach2os (BTB and CNC homology 2 opposite strand) are overlapping and Bach2os locates inside Bach2. We can easily find the blue reads belonging to Bach2os through this strand-specific RNA-Seq approach.
The number of reads for each gene was determined using an HTSeq-count script and the NCBI RefSeq gene annotation (latest version for mm10) (Anders et al., 2015). After differential expression analysis with DESeq2 (Love et al., 2014), there were 2,471 genes that significantly changed expression during early preadipocyte differentiation from Day 0 to Day 2 ( Fig. 3) It is striking that so many genes altered their expression during early differentiation even though phenotypical changes were slight. This is a large percentage (~10%) among the 25,000 total genes in the mouse genome. When excluding genes that were not or extremely lowly expressed in 3T3-L1 cells, the percentage went up to about 18%. Key marker genes, i.e., PPAR, CEBP, FABP4 and FASN, were differentially expressed, consistent with reprograming preadipocytes (Fig. 3). Red blocks represent the reads on the positive strand, which came from Mkrn2 gene with the same direction pointing to the right. Blue blocks represent the reads on the negative strand coming from Raf1 gene. The most distinguishable characteristics is that in the red oval region, only the splicing junctions detected from Mkrn2; in the blue oval, only Raf1 splice junctions detected and more remarkably alternative splicing events detected too; (B) Coding gene Bach2 and lncRNA Bach2os are overlapping genes but have opposite directions. Red blocks represent the reads on the positive strand, which came from Bach2 gene with splicing junctions detected. Blue blocks represent the reads on the negative strand coming from Bach2os lncRNA with splice junctions, which might be novel alternative splicing. Red oval is for forward strand transcripts and reads, blue oval for backward strand transcripts and reads.

Overall Expression Profile of lncRNAs
The overall expression of lncRNAs before preadipocyte differentiation and during early differentiation, was specifically examined using the GENCODE Mouse Release M22 (GRCm38.p6) long-non-coding RNA gene annotation from our RNA-Seq data. The Version M22 lncRNA annotation contains 13,450 lncRNAs, which is a subset of the GENCODE M22 comprehensive gene annotation (Frankish et al., 2019). Most of the lncRNAs contain poly(A) tails as coding mRNAs do and therefore lncRNAs were included in the sequencing libraries during sample preparation and library construction. Violin plots show that before prefiltering, the majority of lncRNAs had zero or extremely low expression (Fig. 4A). We applied prefiltering to remove the unexpressed or very low expressed lncRNAs (Love et al., 2014). The cutoff is set as sum of reads in a row ≥ 10, where each row is a gene or lncRNA. 2,606 (19.4% of total 13,450 lncRNAs) remained after prefiltering. This also indicates that the overall expressional level of lncRNAs are low. The data shape is close to a normal distribution, similar to that seen for coding genes via RNA-Seq (Fig. 4B). Near-normal distribution meets the best expectation for further differentiation analysis (Cambronne et al., 2012). Compared with that of the coding mRNAs, the expression of lncRNAs were remarkably low ( Fig. 4C), which is consistent with previous reports (Cabili et al., 2011). agreed with each other. Therefore, the lncRNA expression profiles after early differentiation (Day 2) were altered compared to non-differentiated cells at Day 0.
Apart from U90926 and Gas5, the remaining 80 lncRNAs have not been reported to be involved in adipogenesis or obesity. Several lncRNAs in the known lncRNA category are noteworthy because they have already been shown to be involved in Wnt-signaling, albeit in other tissues (Table 2). Wincr1, or Wnt signaling-Induced Non-Coding RNA 1, was downregulated by 0.29-fold, FDR 2.0x10-9 after 2 days of differentiation. Mullin et al. (2017) reported that Wincr1 was induced by -catenin activity during skin fibrosis as its name suggests. Wincr1 further regulates the expression of key markers of dermal fibrosis. Kcnq1ot1 is involved in imprinting a cluster of genes on the sense-strand including KCNA1 (Mancini-Dinardo et al., 2006). Our results show that Kcnq1ot1 was downregulated by 0.42-fold, FDR 0.02%. Hotairm1, or Hoxa transcript antisense RNA, myeloid-specific 1, is a lncRNA that is located on the opposite strand of Hoxa gene cluster. It was down-regulated by 0.56-fold, FDR 8%. Zhang et al. (2018) studied Hotairm1 in the pathogenesis of hepatocellular carcinoma and found that it was down-regulated and might have resulted in inhibiting tumor cell proliferation reduction and promoting apoptosis through Wnt signaling pathway. Malat1 is one of the earliest identified lncRNA, which is related to various pathological processes from diabetes to cancers (Wu et al., 2015). Malat1 was down 0.49-fold for 2 days of differentiation. Mir17hg, or Mir17 host gene, is a lncRNA of Pri-miRNA class hosting miR-17 to miR-92 on Chromosome 14. Mir17hg was downregulated by 0.57-fold, FDR 1.1%. While lncRNAs mentioned above were all down-regulated, multiple lncRNAs were up-regulated, including Dubr (Dppa2 upstream binding, 1.95-fold up), Ptprv (protein tyrosine phosphatase, receptor type V, 3.35-fold up), Snhg8 (small nucleolar RNA host gene 8, 2.03-fold up), Snhg1 (small nucleolar RNA host gene 1, 1.57-fold up), Trerf1 (transcriptional regulating factor 1, 1.70-fold up), Brip1os (BRCA1 interacting protein C-terminal helicase 1, opposite strand,1.38-fold up) and Tug1 (taurine upregulated gene 1, 1.22-fold up). Those up-regulated lncRNAs might play important roles during early preadipocyte differentiation, however, the mechanisms of their function remain to be uncovered.

Functional Analysis of lncRNAs Involved in Early Preadipocyte Differentiation
Because the gene ontology information is very limited, we analyzed the promoter regions of lncRNA in our list using GREAT (McLean et al., 2010) in order to find their biological implication. We took ±5 kb region of lncRNA's transcription start site and feed into the GREAT web tool version 3.0.0 with the "two nearest genes" option. Since the number of genes are relatively small, we lowered the threshold of Minimum Regionbased Fold Enrichment to 1.2 and used Raw P-value threshold of 0.05. Among the multiple categories GREAT detected, two of them are most interesting: Molecular Signatures Database (MSigDB) Predicted Promoter Motifs and MSigDB Perturbation (Fig. 6). In the promoter motif analysis, we found that several Transcription Factors (TFs) stood out including RUNX1, AHR, CEBPA and RFX1 (Fig. 6A). Those 4 TFs were also in the list of significantly expressed genes of mRNA-Seq results. RUNX1, or runt-related transcription factor 1, was down-regulated by 0.74-fold, FDR 0.03%. AHR, or aryl-hydrocarbon receptor, was down-regulated by 0.52-fold, FDR 2.110 8 . CEBPA, or CCAAT/enhancer binding protein (C/EBP) alpha, was up-regulated by 1.74-fold, FDR 2.8%. RFX1, or regulatory factor X 1, up-regulated by 1.75-fold, FDR 0.12%. All these lncRNAs might be regulated by those TFs in a coordinated way although the exact mechanism is not known yet. From the previous perturbation studies in MSigDB, we found that expression of some lncRNAs matched the published data (Burton et al., 2004), especially the data from 8-48 h or 8-96 h during differentiation of 3T3-L1 cells into adipocytes, which is consistent with our experimental design (Fig. 6B).

MSigDB Predicted Promote Motif
Genes up-regulated in response to both single dose and fractionated radiation that were common to all three cell line studied.

Progressively up-regulated form 8-48 h during differentiation of 3T3-L1 cells (fibroblast) into adipocytes. Genes down-regulated in NIH3T3 cells (fibroblasts) transformed by activated KRAS vs over-expression of a dominant negative form of CDC25 by rescue
Genes up-regulated in response to hydrogen peroxide in CS-B cells (Cockayne syndrome fibroblast, CS) expressing ERCC6 off a plasmid vector Genes consistently down-regulated in mature mammary luminal cells both in mouse and human species. Genes commonly down-regulated in human alveolar rhabdomyosarcoma (ARMS) and its mouse model overexpressing PAX3-FOXO1 fusion.
Transcripts depleted from pseudopodia of NIH/3T3 cells (fibroblast) in response to haptotactic migratory stimulus by fibronectin, FN1 Genes from 'subtype S1' signature of hepatocellular carcinoma (HCC): aberrant activation of the WNT signaling pathway.

Discussion
Long non-coding RNAs have been reported to be involved in many biological processes, such as adipogenesis, development and disease pathogenesis, especially cancers. We can improve our understanding of how lncRNAs work by specifically identifying the lncRNA involved in biological processes. Here we used cutting-edge strand-specific RNA-Seq technology to study early preadipocyte differentiation. The strandspecific approach provides power and advantage to distinguish overlapping transcripts in opposing directions, such as antisense and opposite strand genes (Fig. 2). We identified 82 lncRNAs that significantly changed their expression after 2-day differentiation of 3T3-L1 preadipocytes ( Fig. 5 and 6). Only 2 (U90926 and Gas5) had previously been reported to be involved in adipogenesis or obesity; they most likely promote adipogenesis through the canonical transcription factors, i.e., PPAR, CEBP, FABP4 and FASN (Chen et al., 2017;. Our data agreed with their findings. Interestingly, both Liu et al. (2018) and  studied Gas5 in adipogenesis in 3T3-L1 cells for 4 days differentiation induction and mesenchymal stem cells for 3 or more days induction, respectively. Both groups concluded that Gas5 suppresses adipogenesis, through other microRNAs, which is controversy to our results Liu et al., 2018). However, this could be explained that both groups studied Gas5 for more than 3 days of differentiation (we did 2 days) and they used overexpression system and knockdown approach, which might result in cell fate during early adipogenesis.
The remaining 80 lncRNAs, however, were relatively unknown with regard to adipocyte differentiation. After further analysis, we found that many known annotated lncRNAs on our list (Table 2) have connections to the Wnt/-catenin signaling pathway. There are many reports linking lncRNA with Wnt signaling pathway (Zarkou et al., 2018). When we were analyzing the promoter regions of the lncRNAs in our list, one item intrigued us, i.e., Wnt signaling pathway in Fig. 6B and LEF1/TCF1/TCF8 (proteins involved in Wnt pathway) in Fig. 6A. The Wnt/-catenin signaling pathway controls many biological processes during animal development and tissue homeostasis. It is tightly regulated as its functions are critical for early tissue development as well as cell proliferation during growth and tissue maintenance. The Wnt pathway is involved in cell differentiation, such as adipogenesis and obesity (Chen and Wang, 2018). Wnt/-catenin signaling pathway was reported to inhibit adipogenesis in 3T3-L1 cells through the key marker genes, i.e., CCAAT/Enhancer Binding Proteins  (CEBP) and Peroxisome Proliferator-Activated Receptor  (PPAR).
Our results show that Gas5 was upregulated during early preadipocyte differentiation (Fig. 5) is consistent with the findings from Song et al. (2019), who found that Wnt/-catenin signaling pathway was activated in colorectal cancer cells. However, lncRNA Gas5 could inhibit the activation of this pathway and repress tumor cell invasion, metastasis and migration. Furthermore, Shi et al. (2019) found that stabilization of lncRNA Gas5 led to restored glucose processing in type 2 diabetic patients and suggested Gas5 to be a therapeutic target for diabetes. Therefore, escalated Gas5 level inhibits Wnt signaling pathway and further promote adipocyte differentiation. Kcnq1ot1 was originally discovered at an imprinting lncRNA (Mancini-Dinardo et al., 2006). Gu et al. (2019) found that lncRNA Kcnq1ot1 delayed bone fracture healing through the Wnt/β-catenin pathway. Gao et al (2018) found that Kcnq1ot1 promotes osteoblast differentiation by activating the Win/b-catenin signaling pathway . Those two papers and our results indicate that Kcnq1ot1 is important for adipogenesis. Hu et al. (2017) and Li et al. (2019) reported that increased Malat1 expression level facilitated nuclear translocation of -catenin and induced its target gene expressions through Wnt/-catenin signaling pathway. On the contrary, decreased Malat1 level reversed nuclear translocation of -catenin, which is consistent with our findings. Hotairm1 was found to be downregulated in hepatocellular carcinoma and inhibited the proliferation and promote the apoptosis by suppressing the Wnt/-catenin signaling pathway and then inhibit the tumor progression. (Mu et al., 2019;Zhang et al., 2018) Yuan et al. found that the overexpression of interferon regulatory factor-1 (IRF-1) downregulated Mir17hg in gastric cancer cells and inhibited the tumor progression by inhibiting Wnt/-catenin signaling pathway. Yuan et al. (2019) Wincr1 was found to be activated by Wnt/-catenin signaling pathway in fibroblast cells. (Mullin et al., 2017) So Wincr1 might be the downstream target of Wnt signaling pathway during adipogenesis. In short, during the early adipogenesis, Wnt/-catenin signaling pathway was inhibited by the down-regulations of Kcnq1ot1, Malat1, Hotairm1 and Mir17hg and the up-regulation of Gas5 (Fig. 7), perhaps together with transcription factors and non-coding RNAs, such as microRNAs Liu et al., 2018).
We also compared with the lncRNAs identified between preadipocytes and mature adipocytes from Rinn laboratory (Sun et al., 2013). To our surprise, the profile of lncRNAs in their data is significantly different from ours (Supplemental Table S2). Malat1 and Kcnq1ot1 were up-regulated by 9.7-and 4.8-fold, respectively, while they were down-regulated in our results. However, several other lncRNAs (Snhg1, Trerf1, Snhg8 and Brip1os) were downregulated in their results, which is opposite to ours. In another study, lnc-OAD (1700018A04Rik) were reported to be required for adipocyte differentiation in 3T3-L1 cells (Wang et al., 2019). Lnc-OAD was rapidly stimulated at day 1 but returned to basal level at day 2, then increased again after day 4. This is consistent with our data at day 2, where expression of 1700018A04Rik was not significantly changed compared with day 1 (data not shown). All these suggest that the expression levels of lncRNAs are dynamic depending on different stages of cell development. It further infers that lncRNAs provide very precise regulation during the adipogenesis process.
To be summarized, in this study we identified a cohort of lncRNAs that play important roles during early preadipocyte differentiation, among which many lncRNAs were discovered for the first time. Besides the key transcription factors, such as PPAR, CEBP, FABP4 and so on, lncRNAs may act as fine-tuning the adipogenesis process. Those lncRNAs may serve as biomarkers for diagnosis or therapeutic targets. We hope all these findings will shed light on obesity research and therapy.

Conclusion
We identified 82 lncRNAs in early preadipocyte differentiation in 3T3-L1 cells and 98% of them are novel. Those lncRNAs promote adipogenesis possibly through inhibiting Wnt/-catenin signaling pathway.