Crucial lncRNAs associated with adipocyte differentiation from human adipose-derived stem cells based on co-expression and ceRNA network analyses

Background Injection of adipose-derived stem cells (ASCs) is a promising treatment for facial contour deformities. However, its treatment mechanisms remain largely unknown. The study aimed to explain the molecular mechanisms of adipogenic differentiation from ASCs based on the roles of long noncoding RNAs (lncRNAs). Methods Datasets of mRNA–lncRNA (GSE113253) and miRNA (GSE72429) expression profiling were collected from Gene Expression Omnibus database. The differentially expressed genes (DEGs), lncRNAs (DELs) and miRNAs (DEMs) between undifferentiated and adipocyte differentiated human ASCs were identified using the Linear Models for Microarray Data method. DELs related co-expression and competing endogenous RNA (ceRNA) networks were constructed. Protein–protein interaction (PPI) analysis was performed to screen crucial target genes. Results A total of 748 DEGs, 17 DELs and 51 DEMs were identified. A total of 13 DELs and 279 DEGs with Pearson correlation coefficients > 0.9 and p-value < 0.01 were selected to construct the co-expression network. A total of 151 interaction pairs among 112 nodes (10 DEMs; eight DELs; 94 DEGs) were obtained to construct the ceRNA network. By comparing the lncRNAs and mRNAs in two networks, five lncRNAs (SNHG9, LINC02202, UBAC2-AS1, PTCSC3 and myocardial infarction associated transcript (MIAT)) and 32 genes (i.e., such as phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1), protein tyrosine phosphatase receptor type B (PTPRB)) were found to be shared. PPI analysis demonstrated PIK3R1 , forkhead box O1 (FOXO1; a transcription factor) and estrogen receptor 1 (ESR1) were hub genes, which could be regulated by the miRNAs that interacted with the above five lncRNAs, such as LINC02202-miR-136-5p-PIK3R1, LINC02202-miR-381-3p-FOXO1 and MIAT-miR-18a-5p-ESR1. LINC02202 also could directly co-express with PIK3R1. Furthermore, PTPRB was predicted to be modulated by co-expression with LINC01119. Conclusion MIAT, LINC02202 and LINC01119 may be potentially important, new lncRNAs associated with adipogenic differentiation of ASCs. They may be involved in adipogenesis by acting as a ceRNA or co-expressing with their targets.


INTRODUCTION
Autologous adipose tissue grafting has been a widely accepted surgical tool for anti-aging cosmetics (Charles-De-Sá et al., 2015) and reconstructive restoration of various congenital or acquired facial soft tissue deformities (Bashir et al., 2018). However, conventional fat grafting procedure needs to be repeated multiple times to achieve satisfactory results (Bashir et al., 2018), which may be associated with the low graft survival rate and poor revascularization (Ma et al., 2015). To overcome these two limitations, recent scholars propose to combine with additional autologous adipose-derived stem cells (ASCs) which have the ability to differentiate into mature adipocytes to supplement apoptotic cells and secrete angiogenic growth factors to enhance angiogenesis (Bashir et al., 2018;Kotaro et al., 2008;Philips, Marra & Rubin, 2014). The clinical trials also confirm that supplementation of ASCs to adipose grafts is superior to conventional lipoinjection for facial recontouring (Bashir et al., 2018;Kotaro et al., 2008). Nevertheless, the use of autologous ASCs has not been FDA-approved. This may be because there still remains a huge gap in understanding the potential mechanisms of ASCs for adipocyte differentiation.
Increasing evidence has suggested long noncoding RNAs (lncRNAs), a class of noncoding RNAs more than 200 nucleotides, play crucial roles in adipogenesis for ASCs. For example, Nuermaimaiti et al. (2018) demonstrated that knockdown of HOXA11-AS1 inhibited adipocyte differentiation, leading to suppression of adipogenic-related gene transcription, as well as decreased lipid accumulation in ASCs. Huang et al. (2017) observed knockdown of MIR31HG inhibited adipocyte differentiation, whereas overexpression of MIR31HG promoted adipogenesis in vitro and in vivo. MEG3 was also found to be downregulated during adipogenesis of ASCs. Functional analysis showed that knockdown of MEG3 promoted adipogenic differentiation of ASCs . Furthermore, current research shows lncRNAs, on one hand, functions as microRNA (miRNAs) sponges to bind the miRNA response elements and regulate miRNA-mediated gene silencing (i.e., competing endogenous RNA (ceRNA) hypothesis); and, on the other hand, directly influences their neighboring genes expression by chromatin remodeling or transcriptional control (co-expression model) (Huang et al., 2016;Li, Ao & Wu, 2017). These theories have also been reported in ASCs. Li, Ao & Wu (2017) proved downregulated MEG3 may be insufficient to sponge miR-140-5p and lead to its upregulation during adipogenesis in ASCs. Huang et al. (2017) revealed inhibition of MIR31HG reduced the enrichment of active histone markers, histone H3 lysine 4 trimethylation and acetylation in the promoter of fatty acid binding protein 4, resulting in suppression of its expression and adipogenesis. However, the adipogenic differentiation related lncRNAs and its mechanisms of ASCs remains rarely reported.
The present study aimed to identify crucial lncRNAs involved in adipocyte differentiation of ASCs by constructing lncRNA-miRNA-mRNA ceRNA network and lncRNA-mRNA co-expression network using high throughput analysis data. Our findings might offer greater insights into the molecular mechanisms of adipocyte differentiation from ASCs and provide potentially new targets for inducing adipogenesis.

MATERIALS AND METHODS
Collection of microarray data GSE113253 (Rauch et al., 2019) and GSE72429 datasets (Supplemental Information 1) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi. nlm.nih.gov/geo/). GSE113253 dataset applied the high throughput sequencing methodology to simultaneously detect the lncRNA and mRNA expression profiles in two repeats of undifferentiated human ASCs and 10 repeats of adipogenic differentiation cells using an Illumina HiSeq 1500 instrument, which was submitted to GEO on April 17, 2018. GSE72429 dataset analyzed the miRNA expression profile in four undifferentiated human ASCs and two adipogenic differentiation cells using an Agilent-031181 Unrestricted_ Human_miRNA_V16.0_Microarray (miRBase release 16.0 miRNA ID version), which was submitted to GEO on August 27, 2015.

Differential expression analysis
The normalized series matrix files of each dataset were downloaded from GEO. Following re-annotation according to corresponding platform (GPL18460), the expression values of the lncRNAs and mRNAs in GSE113253 were obtained. The differentially expressed genes (DEGs), lncRNAs (DELs) and miRNAs (DEMs) were identified using the Linear Models for Microarray Data method software (version 3.34.0; Ritchie et al., 2015). p-Value was adjusted by using Benjamini-Hochberg method to avoid false positives. The heatmap was constructed to present the expression difference of DEGs, DELs and DEMs in different samples using the pheatmap package (version: 1.0.8; Kolde, 2019) based on Euclidean distance.

Co-expression network between lncRNA and mRNA
The co-expression network was constructed based on the correlation analysis between DELs and DEGs. Pearson correlation coefficients were calculated using the Weighted Gene Correlation Network Analysis (Langfelder & Horvath, 2016) algorithm to assess the correlation. Only the co-expressed pairs with absolute value of Pearson correlation coefficients ≥ 0.9 and p < 0.01 were selected to draw the network using Cytoscape (version 3.4; Shannon et al., 2001Shannon et al., -2008Kohl, Wiese & Warscheid, 2011).

Protein-protein interaction network
Protein-protein interaction (PPI) data of DEGs in the ceRNA network was collected from Search Tool for the Retrieval of Interacting Genes (STRING; version 10.0; Szklarczyk et al., 2019) database (Szklarczyk et al., 2015). Only interactions with combined score >0.4 were selected to construct the PPI network. Several topological features of the nodes (protein) in the PPI network were calculated using the CytoNCA plugin in cytoscape software (Tang, Li & Wang, 2014;Tang et al., 2015) to screen hub genes, including degree, eigenvector, betweenness and closeness centrality. Furthermore, transcription factors were predicted using iRegulon (Janky et al., 2014) in Cytoscape and then integrated to the PPI network.

Function enrichment analysis
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery online tool (version 6.8; http://david.abcc.ncifcrf.gov) (Huang, Sherman & Lempicki, 2009) to reveal the function of DEGs. p < 0.05 was set as the cut-off value.

Construction of co-expression and ceRNA networks
A total of 13 DELs and 279 DEGs with Pearson correlation coefficients > 0.9 and p-value < 0.01 were selected to construct the lncRNA-mRNA co-expression network, which contained 440 positive connections ( Fig. 2; Supplemental Information 3). Based on at least eight database analyses in miRwalk 2.0 and negatively correlated principles, a total of 79 downregulated DEGs were predicted to be regulated by eight upregulated DEMs, while 128 upregulated DEGs were predicted to be regulated by 32 downregulated DEMs. Using the starBase database, 355 miRNAs were predicted to interact with 25 DELs; using the miRcode database, 192 miRNAs were predicted to interact with eight DELs; using the DIANA-LncBase database, 1,343 miRNAs were predicted to interact with 15 DELs. After overlapping the DEMs that interacted with DELs and DEMs that regulated DEGs, 151 interaction pairs among 112 nodes (10 DEMs, four upregulated and six downregulated; eight DELs, four upregulated and four downregulated; 94 DEGs, 46 upregulated and 48 downregulated) were obtained, which were used for constructing the ceRNA network ( Fig. 3; Supplemental Information 4).

PPI network
Protein-protein interaction pairs were predicted for the 94 DEGs in the ceRNA network using the STRING database, which resulted in 80 interaction relationship pairs that were screened between 58 nodes (24 upregulated and 34 downregulated) (Fig. 4). Phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1), FYN proto-oncogene, Src family tyrosine kinase and estrogen receptor 1 (ESR1) were considered as hub genes in the PPI network because they ranked the top 10 in all four topological features (Table 2). In addition, FOXO1, which was included in the PPI network, was predicted as a Note: All the differentially expressed miRNAs and lncRNAs were shown, but only top 25 upregulated and downregulated mRNAs as well as crucial genes were displayed. FC, fold change. p-Value with asterisk indicated their adjusted p-value were also less than 0.05.
differentially expressed transcription factor to regulate the other target genes in the PPI network using IRegulon plug-in (Fig. 4), indicating FOXO1 was also a hub gene. Function analysis showed eight significant KEGG pathways were enriched, including hsa04015:Rap1 signaling pathway (PIK3R1), hsa05200:Pathways in cancer (PIK3R1,

DISCUSSION
In present study, we identified three crucial lncRNAs (MIAT, LINC02202 and LINC01119) for adipogenesis from human ASCs. MIAT may sponge hsa-miR-18a-5p and influence the inhibition of hsa-miR-18a-5p on the expression of ESR1. LINC02202 may function as a ceRNA for hsa-miR-136-5p/hsa-miR-381-3p to respectively regulate the expressions of PIK3R1 and FOXO1; LINC02202 also may directly affect the transcription of PIK3R1. LINC01119 may co-express with PTPRB to impact its transcription. Although all these relationship pairs may be potentially important, LINC01119-PTPRB co-expression axis may be especially verifiable because their expression significance met the criterion of adjusted p-value < 0.05. Although there have studies to show the roles of lncRNA MIAT for stem differentiation, only osteogenic (Jin et al., 2017) and endothelial cell  differentiation were investigated, without evidence to prove its effect on adipogenesis of human ASCs. A recent study revealed MIAT was an estrogen-inducible lncRNA and its expression was positively related to estrogen receptor (Li et al., 2018b). There was accumulating evidence to reveal that exposure of bone marrow stem cells to icariin or flavonoids of Herba Epimedii inhibited adipogenic differentiation, exhibiting decreased adipocyte numbers and downregulated mRNA expression of adipogenic differentiation markers, peroxisome proliferator-activated receptor gamma (PPARγ) and CCAAT/enhancer-binding protein α (C/EBPα) (Li et al., 2018c;Zhang et al., 2015); while treatment of bone marrow stem cells with estrogen receptor antagonist ICI182780 revered the effects of Herba Epimedii ingredient and promoted adipogenesis (Li et al., 2018c;Zhang et al., 2015). The study of Ihunnah et al. (2014) also demonstrated activation of estrogen receptor in ASCs inhibited adipogenesis by decreasing the recruitment of the adipogenic PPARγ onto its target gene promoters, whereas the use of estrogen receptor antagonism ICI 182780 or knockdown of estrogen receptor-α via lentiviral shRNA enhanced adipogenesis by increasing the expression of PPARγ. Thus, it can be hypothesized that MIAT may be lower expressed in adipogenic differentiation cells like ESR1, which was also confirmed in our study. However, the interaction mechanisms between MIAT and estrogen receptor remain unclear. In present study, we predicted that downregulated MIAT may be insufficient to sponge hsa-miR-18a-5p and lead to more hsa-miR-18a-5p to bind with the 3′ untranslated region of ESR1, inducing the lower expression of ESR1. This hypothesis may be indirectly demonstrated by the fact that miR-18a mimic significantly promoted mesenchymal stem cell (MSC) adipogenic differentiation, while the addition of miR-18a inhibitor obtained the negative effects on adipogenic differentiation of MSCs (Li et al., 2018a). The negative regulatory relationship between ESR1 and miR-18a were also validated in human trophoblast cell line by the luciferase assay (Zhu et al., 2015). LINC02202 may be a newly identified lncRNA associated with stem cell differentiation because its role had not been previously mentioned in the literatures. In this study, we predicted upregulated LINC02202 may be involved in ASCs adipogenic differentiation by regulating phosphatidylinositol 3-kinase (PI3K) signaling. It has been reported that PI3K signaling pathway was strongly activated in MSCs under the adipogenesis-inducing hormone cocktail (Kim et al., 2017), and the addition of PI3K specific inhibitor LY294002 severely suppressed lipid accumulation, as well as the expression of adipogenic markers PPARγ and C/EBPα (Yu et al., 2008). PIK3R1 is a critical component of the PI3K signaling pathway and its expression was also demonstrated to be increased after the induction of adipocyte differentiation from preadipocytes 3T3-L1 (Kim et al., 2014a). Thus, theoretically, PIK3R1 may be upregulated in adipogenic differentiation cells compared with undifferentiated human ASCs, which was confirmed in our study. Activated PI3K/AKT signaling may promote adipogenesis through upregulating downstream transcription factors, such as FoxO1 (Yi et al., 2018) which may subsequently enhance the transcription of its target genes, PPAR-γ and C/EBP-α (Ambele et al., 2016;Munekata & Sakamoto, 2009); whereas persistent inhibition of FoxO1 with its antagonist AS1842856 (Zou et al., 2014) or knockdown of FoxO1 (Sun et al., 2017) was also observed to almost completely suppress adipocyte differentiation and lipogenesis. As expected, we also found FoxO1 was significantly high expressed during adipogenic differentiation. In addition to directly affect the transcription of PIK3R1, LINC02202 may function as a ceRNA for miR-136-5p and hsa-miR-381-3p to regulate the expression of PIK3R1a and FoxO1, respectively. Although there was no study to demonstrate these ceRNA interaction axes, the negative correlation between the expression of miR-136 and adipogenic markers C/EBPα and PPARα in subcutaneous adipose tissue of lambs may indirectly illuminate the importance of miR-136 for adipogenic differentiation (Meale et al., 2014). As expected, we also found miR-136-5p was significantly downregulated in adipogenic differentiation cells.
There was only one sequencing study to identify that LINC01119 was downregulated in colorectal cancer cells after hypoxia treatment (Han et al., 2019). Several authors had demonstrated hypoxia exposure was effective to enhance adipocyte differentiation from ASCs (Fink et al., 2004;Valorani et al., 2012;Kim et al., 2014b), which was medicated by the generation of reactive oxygen species (ROS) and activation of PI3K/ Akt/mTOR (Kim et al., 2014b); the addition of ROS scavenger or Akt/mTOR inhibitor prevented adipocyte differentiation (Kim et al., 2014b). Thus, LINC01119 may have anti-adipose differentiation potential and lower expressed in adipogenic differentiation cells compared with undifferentiated human ASCs, which was validated in our study. However, its mechanisms for adipocyte differentiation remain unclear. We predicted LINC01119 may co-express with PTPRB. The study of Kim et al. showed ectopic over-expression of PTPRB inhibited the expression of adipocyte-related genes (such as PPAR-γ) and led to a reduced adipocyte differentiation from preadipocytes. Also, PTPRB was reported to suppress the tyrosine phosphorylation of VEGFR2 during adipocyte differentiation (Kim et al., 2019). Generally, VEGF functions by binding with VEGFR2, while transfection of VEGF to ASCs increased fat cell survival (Zhang et al., 2017). These findings suggest PTPRB may also be downregulated to promote VEGF secretion and activate its mediated pathways, ultimately inducing adipogenic differentiation from ASCs. This hypothesis was in line with our study showing PTPRB was lower expressed in adipocyte differentiation cells and was involved in angiogenesis.
There are some limitations in this study. First, only two datasets were submitted between 5 years until now, and not all were used for this analysis, which may cause some bias in results due to the small sample size and different data platforms. However, we believe the sequencing or microarray technology may be more mature recently and thus the results may be more believable. This was also indirectly reflected by the less overlapped genes if the other datasets were used (only two comparing GSE72429 with GSE25715; Guo & Cao, 2019) and thus, we renounced the use of multiple datasets and only the newly one. Moreover, this work investigated lncRNA co-expression and ceRNA mechanisms, which required the lncRNA and mRNA should be simultaneously analyzed. Thus, some datasets that only independently investigated lncRNA or mRNA were also excluded. Second, the crucial co-expression and ceRNA axes were obtained by database prediction, which may lead to many false positives. Therefore, further in vitro wet experiments (PCR, luciferase assay, knockdown or overexpression) are still indispensable to confirm the interaction between lncRNAs and miRNAs, lncRNA and mRNAs as well as the miRNAs and mRNAs and their roles during adipogenic differentiation of ASCs.