Dysregulated lncRNAs Regulate Human Umbilical Cord Mesenchymal Stem Cells Differentiation into Insulin-Producing Cells by Forming a Regulatory Network with mRNAs


 ObjectiveIn recent years, cell therapy has become a new research direction in the treatment of diabetes. However, the underlying molecular mechanisms of mesenchymal stem cells (MSCs) participate in such treatment has not been clarified. MethodsIn this study, human umbilical cord mesenchymal stem cells (HUC-MSCs) isolated from newborns were progressively induced into insulin-producing cells (IPCs) using small molecules. HUC-MSCs (S0) and four induced stage (S1-S4) samples were prepared. We then performed transcriptome sequencing experiments to obtain the dynamic expression profiles of both mRNAs and long noncoding RNAs (lncRNAs). ResultsWe found that the number of differentially expressed lncRNAs and mRNAs showed a decreasing trend during differentiation. Gene Ontology (GO) analysis showed that the target genes of differentially expressed lncRNAs were associated with translation, cell adhesion, and cell connection. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the NF-KB signaling pathway, MAPK signaling pathway, HIPPO signaling pathway, PI3K-Akt signaling pathway, and p53 signaling pathway were enriched in these differentially expressed lncRNA-targeting genes. We also found that the coexpression of the lncRNA: CTBP1-AS2 with the PROX1, and the lncRNAs AC009014.3 and GS1-72M22.1 with the mRNA JARID2 was related to the development of pancreatic beta cells. Moreover, the coexpression of the lncRNAs :XLOC_ 050969, LINC00883, XLOC_050981, XLOC_050925, MAP3K14- AS1, RP11-148K1.12, and CTD2020K17.3 with p53, regulated insulin secretion by pancreatic beta cells.ConclusionThis research revealed that HUC-MSCs combined with small molecule compounds were successfully induced into IPCs. Differentially expressed lncRNAs may regulate the insulin secretion of pancreatic beta cells by regulating multiple signaling pathways. The lncRNAs: AC009014.3,Gs1-72m21.1 and CTBP1-AS2 may be involved in the development of pancreatic beta cells, and the lncRNAs: XLOC_050969, LINC00883, XLOC_050981, XLOC_050925, MAP3K14-AS1, RP11-148K1.12, and CTD2020K17.3 may be involved in regulating the insulin secretion of pancreatic beta cells, thus providing a lncRNA catalog for future research regarding the mechanism of the transdifferentiation of HUC-MSCs into IPCs. It also provides a new theoretical basis for the transplantation of insulin-producing cells into diabetic patients in the future.


Introduction
Diabetes mellitus is caused by insulin secretion de ciency associated with pancreatic beta cells damage to varying degrees, and the inability to regulate the metabolic balance of blood sugar. Therefore, replacement therapy with pancreatic beta cells has become a therapeutic option for diabetes. The clinical application of islet transplantation is limited by insu cient donors, immune rejection, immunosuppressive drugs side effects and possible beta cells toxicity [1]. Stem cells are an important alternative that could provide innumerable potential islet cell sources for transplants. Umbilical cord mesenchymal stem cells (UC-MSCs ) are obtained from discarded placenta and they are more readily available and have a higher proliferative potential compared to other MSCs [2,3]. UC-MSCs are more primitive MSCs than bone marrow or adipose mesenchymal stem cells(BM-MSCs or AD-MSCs) and do not express major histocompatibility complex (MHC) class II (HLA-DR) antigens [4,5], which makes them good candidates for potential allogeneic therapeutic applications. Previous studies have shown that UC-MSCs could be differentiated into IPCs in vitro and improve the glycemic status of diabetic mice after transplantation in vivo [6,7]; however, the induction e ciency was low and the maintenance time of hypoglycemia was slightly short [6][7][8]. Therefore, the differentiation e ciency and function of IPCs must be improved. Moreover, studies should focus on identifying the underlying molecular mechanisms by which UC-MSCs differentiate into IPCs.
Many small molecule compounds, growth factors, activators and inhibitors can transdifferentiate MSCs into IPCs, which improves the survival of IPCs and the ability to release insulin in vitro [9][10][11][12]. MSCs transdifferentiation e ciency and insulin secretion can be improved by overexpression of Pdx1, Neurog3, MafA, and Pax4 [13][14][15]. A comprehensive understanding of the signaling pathways and temporal transcription factor activation patterns during human and rodent pancreas organogenesis has accelerated the production of IPCs from human pluripotent stem cells (hPSCs) [16][17][18]. Moreover, MSCs differentiation is accurately regulated and coordinated by mechanical and molecular signals from extracellular environments, which involves complex pathways regulated at the transcriptional and posttranscriptional levels [19][20][21]. Growing evidence suggests that lncRNAs are important gene expression regulators through their interactions with DNA, RNA, or proteins [22][23][24]. By acting as guides, scaffolds, and/or decoys, lncRNAs can modulate chromatin organization or epigenetics, regulate transcription or control mRNA progressions [25]. LncRNAs play an important role in many biological processes, including development [26], differentiation [27], and metabolism [28,29]. However, the speci c mechanism by which lncRNAs regulate the transdifferentiation of MSCs into IPCs remains poorly understood.
In this study, a new 4-stage protocol was developed based on previous studies that investigated small molecule compounds [30]. Moreover, cells at each stage were collected for highthroughput sequencing(Illumina PE150)to obtain the expression pro le of lncRNAs and mRNAs, to explore their molecular mechanisms during the differentiation of HUC-MSCs into IPCs, thereby providing potential pathways to improve the e ciency of induction in the future. To our knowledge, this is the rst study to reveal the changes in lncRNA and mRNA expression during the differentiation of HUC-MSCs into IPCs, and to enrich the lncRNA genes pool for the transformation of stem cells into IPCs.

Separation, Cultivation, Identi cation of HUC-MSCs
Human umbilical cord (HUC) specimens were gathered from the obstetric department of the Second A liated Hospital of Nanchang University. Ethic committee of Second A liated Hospital of Nanchang University approved this study, and informed consent were obtained from all donors. HUC tissue was rinsed with phosphate-buffered saline (PBS). HUC tissue was cut into little pieces using ophthalmic scissor and laid in culture ask. UC-MSCs were cultured in ḁ-MEM complete medium (Gibco, USA), containing 10% fetal bovine serum (10099-141, Gibco, USA ), 1% penicillin and streptomycin (Sigma) in a humidi ed atmosphere with 5% CO 2 at 37°C .After 5 days, the uid was changed for the rst time, and the uid was changed every 3 days after adhesion. When the cells were fused to 70%-80%, they were digested with trypsin (0.25%Trypsin, Sigma, USA) and passed in a ratio of 1:3. The third passage was used for subsequent experiments. The inverted microscope (Nikon Corp., Tokyo, Japan) showed a uniform fusionbroblast morphology and adherent phenotype. The third generation of UC-MSCs were suspended in PBS, stained with uorescence-labeled monoclonal antibodies CD34, HLA-DR, CD105, CD73, CD90, CD14,CD44, and CD45 (Biolegend, Beijing), and then incubated in darkness for 30 min at 4°C, including appropriate isotype controls (for each antibody isotype). LSRFortessa Flow Cytometer (BD) was used to detect the expression of the above cell surface markers. Differentiation kit (Cyagon, China) was used to identify the differentiation ability of stem cells. Brie y, UC-MSCs were counted by 2× 10^4/cm 2 and seeded at a six-well plate, and cultured in the adipose and osteoblast differentiation culture media for 14 d, respectively. Adipocytes and osteocytes were stained with oil red O and Alizarin red, respectively. For chondrogenesis differentiation, UC-MSCs were seeded at a six-well plate by 2x10^4/cm 2 in chondroblast differentiation culture media for 21 d. The chondrocytes were stained with Alcian blue.

Induction of HUC-MSCs
UC-MSCs from third passage were used for IPCs differentiation .The differentiation protocol was adopted from the previous study [37], and slightly modi ed (a). After the cells were digested, cells were laid in a sixwell plate by 2x10^5/well in serum--free DMEM complete medium at 37℃ in a humidi ed atmosphere of 5% CO 2 .When the cell fusion reached 70%-80% ,induction was performed as follows: The medium was replaced with stage 1 medium containing DMEM:F12 with 10%FBS, 1% penicillin and streptomycin (Sigma),10mM nicotinamide (NA, Sigma), 20ng/ml EGF (PEPROTECH,USA), 50µM β-mercaptoethanol and 4nM Activin A (Gibco, USA) for 3 days. On fourth day, at stage 2, the media was repalced with the same supplements as stage1 medium, but Activin A and β-mercaptoethanol were removed and cultured for another 5 days. On ninth days, media was changed to stage 3 containing DMEM:F12 (1:1) with 2% FBS, 1% penicillin and streptomycin (Sigma), 20ng/ml bFGF (PEPROTECH,USA), 100xITS, 10nM exendin-4 (Gibco, USA), and 10mM NA and cultured for another 7 days. On sixteenth days, at stage 4, the media was replaced with the same supplements as stage 3, but 2%FBS supplement was removed. The cells were matured in this medium for another 7 days. Change in morphology was observed under light microscope during the differentiation.

Real time quantitative PCR (QPCR)
Total RNA was extracted using Trizol reagent® (Tiangen Biochemical Technology, Shanghai, China) for each stage, and then was reverse-transcribed into cDNA template for QPCR ampli cation using PrimeScript™ RT reagent Kit with gDNA Eraser (RR047B,TaKaRa, Japan), QPCR ampli cation was performed on Master cycler Realplex 2 (Eppendorf, USA), with a SYBR Green Realtime RCR Master Mix (RR82LR, TaKaRa, Japan). The QPCR Primers were synthesized by Invitrogen (Shanghai, China) and are listed in Table1. Table 1 Primer sequences used in the quantitative RT-PCR.

In vitro insulin secretion assay
IPCs from the fourth stage and UC-MSCs were slightly washed twice with Sugar-free KRBH buffer wash.
The cells were then pre-incubated in KRBH culture media containing 5.5 mM and 25 mM glucose at 37°C for 2 h rs, respectively, and supernatants were then collected for insulin quanti cation. Insulin concentration was assessed by an ELISA kit (DINS001*KT, Alpco, Salem, USA). Insulin levels were detected according to the standard curve.

Immuno uorescent
The cells from each stage of the groups described above were xed in 4% paraformaldehyde for 20 min. After washing with PBS, cells were permeabilized with 0.1%Triton-100 for 20 min, then cells were blocked with the goat serum for 30 min. They were then incubated with the corresponding primary antibodies, Pdx1 antibody (rabbit,1:200, #5679S ,CST),C-peptide antibody (rabbit,1:100, #4593S ,CST), Insulin antibody(rabbit,1:200, #ab181547,Abcam) and overnight at 4℃. Next, the cells were washed 3 times with PBST, then further incubated for 1hr in a wet box avoid light with secondary antibodies: Goat anti-rabbit immunoglobulinG heavy and light chains (IgGH&L ,1:250, proteintech). The nucleus was stained with Hochest33342 for 7 minutes after washing with PBST and nally observed under a uorescence microscope.
2.6 Dithizone (DTZ) staining DTZ (Sigma-Aldrich, USA) stock solution was prepared by dissolving 20 mg of DTZ in 2 ml DMSO. IPCs from the fourth stage were added into 1 mL 1× PBS buffer and 10 µL DTZ stock solution, and then incubated at 37°C for 15 min. A phase contrast microscope (Nikon Corp) were used to examined the scarlet red stained clusters.

Sample quali cation and quanti cation
TRIzol® reagent (Invitrogen, USA) was used to extract total RNA at each stage, and genomic DNA was removed using DNase I (Takara, Japan) according to manufacturer's instructions. The RNA samples were rstly quali ed using 1% agarose gel electrophoresis for possible contamination and degradation. NanoPhotometer® spectrophotometer were used to examine the RNA purity and concentration. RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system were used to measure the RNA integrity and quantity.

Library preparation
RNA library for lncRNA-seq was prepared by rRNA depletion and stranded method. Brie y, the rRNA Removal Kit (RS-122-2402,Illumina,USA) was used to deplete the ribosomal RNA from total RNA following manufacturer's instruction. RNA was then fragmented into 250~300 bp fragments, and fragmented RNA and dNTPs (dATP, dTTP, dCTP and dGTP) were used to reverse-transcribed the rst strand cDNA. RNase H was used to degrade RNA, and DNA polymerase I and dNTPs (dATP, dUTP, dCTP and dGTP) was used to synthesis second strand cDNA. Remaining overhangs of double-strand cDNA were converted into blunt ends by exonuclease/ polymerase activities. After adenylation of 3' ends of DNA fragments, sequencing connector were ligated to the cDNA. To optimize the 250-300 bp cDNA fragments, AMPure XP system was used to purify the library fragments. Uracil-N-Glycosylase was used to perform uridine digestion, which was followed by the cDNA ampli cation using PCR.
After library construction, the Qubit® uorometer was used to measure the concentration of library and adjusted to 1ng/uL. the insert size of the acquired library was deployed to examine using Agilent 2100 Bioanalyzer. At last, qPCR was used to again examine the accurate concentration of cDNA library .Once the inserted library is of the same size and concentration, the sample can be sequenced.

Sequencing
After library preparation and pooling of different samples, the samples were sequenced by Illumina. Commonly, the lncRNA-seq use PE150 (paired-end 150nt) sequencing to get 12G raw data.

Reads alignment and differentially expressed gene (DEG) analysis
Clean reads were aligned to the human GRch38 genome [31] (http://support.illumina.com/sequencing/sequencing_software/igenome.html) using Hisat2(http://ccb.jhu.edu/software/hisat2;Kim D et al., 2015). Finally, uniquely mapped reads were used to calculate read number and fragments per kilobase of exon per million fragments mapped (FPKM) for each gene. The expression levels of genes were evaluated by FPKM. The software DESeq2 [32], which is speci cally used to analyze the differential expression of genes (DEGs), was applied to screen the RNAseq data for DEGs. The results were analyzed using the fold change (FC≥2 or ≤0.5) and false discovery rate (FDR≤0.05) to determine whether a gene was differentially expressed.

LncRNAs prediction and direction identi cation
To systematically analyze the lncRNAs expression pattern, a pipeline was used for lncRNAs identi cation similar as previously reported [33], which was constructed based on the cu inks software [34].

Coexpression analysis
To explore the regulatory mode between lncRNAs and their host mRNAs, the Pearson's correlation coe cients (PCCs) were calculated between them and classi ed their relation into three class: positive correlated, negative correlated and non-correlated based on the PCCs value.

K means analysis
To identify genes that exhibit similar gene expression across the different development stages, gene expression data was clustered using the k means clustering algorithm as implemented in R. FPKM expression values of row-scaled genes combined in a single matrix. To determine the appropriate group size k, the k means clustering was repeated six times with k ranging between 3 and 8.

Functional enrichment analysis
To sort out functional categories of DEGs, Gene Ontology (GO) terms and KEGG pathways were identi ed using KOBAS 2.0 server [35], Hypergeometric test and Benjamini-Hochberg FDR controlling procedure were used to de ne the enrichment of each term.

Other statistical analysis
Principal component analysis (PCA) analysis was performed by R package factoextra (https://cloud.rproject.org/package=factoextra) to show the clustering of samples with the rst two components. The pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R was used to perform the clustering based on Euclidean distance.

Statistical analysis
The experimental data were analyzed and processed using GraphPad Prism8 system. The results were expressed as mean ± deviation (X±SD). The t-test method was used for comparison between two groups, and the one-way analysis of variance was used for comparison between multiple groups.

In vitro IPCs differentiation and Dithizone (DTZ) staining
After in vitro induction of UC-MSCs, the cell morphology did not change signi cantly by adding βmercaptoethanol combined with Activin A in the rst stage(S1, Figure 2i B) ,compared with that of UC-MSCs without induction(S0, Figure 2i A), which showed a typical long spindle shape and some cell death. In the second stage (S2, Figure 2i C), the cells were further induced by removing β-mercaptoethanol over time. In the third stage(S3, Figure 2i D), the cells were further shortened from the initial long fusiform shape to a ake shape, by reducing the serum concentration, and then adding bFGF, exendin-4, and NA. Finally, in the fourth stage (S4, Figure 2i E), the cell morphology changed from a ake shape to an isletlike cluster shape by removing fetal bovine serum (FBS). In the fourth stage (IPCs), dithizone staining was scarlet, indicating that UC-MSCs were successfully differentiated into zinc-containing islet beta-like cells (Figure 2i F).

The mRNA expression levels during IPCs differentiation
After collecting the cells of each stage (S0, S1, S2, S3, and S4), RNA extraction and real-time uorescence quantitative PCR were performed. As shown in Figure 4, the marker genes of de nitive endoderm Foxa2 and Sox17 from the rst stage (S1, Figure 2ii A-B) were increased signi cantly compared with that at other stages, and the difference was statistically signi cant (P<0.001). Then, the pancreatic progenitor marker genes PDX1, Ngn3, and Nkx6.1 were increased signi cantly during the third stage (S3, Figure 2ii C-E) compared with other stages, and the difference was statistically signi cant (P<0.05), indicating the successful induction of pancreatic precursor cells. Finally, the speci c marker genes of islet beta-like cells Glut2, MAFA, and insulin were increased signi cantly from the fourth stage (S4, Figure 2ii F-H), and the difference was statistically signi cant (P<0.001), indicating successful induction of islet beta-like cells.
3.4 Immuno uorescent to detect the expression of key proteins during IPCs identi cation.
After collecting cells from different stages (S0, S3, and S4), immuno uorescence experiments were performed. Figure 3i (A-C) showes the changes in the expression of pdx1 protein in the S0, S3, and S4 stages and reveals the nuclear expression of pdx1. The gure shows that pdx1 was almost not expressed in the S0 stage, presented the highest expression in the S3 stage and showed decreased expression in the S4 stage, which was consistent with the gene expression, indicating that it was successfully induced into pancreatic precursor cells. The D-E diagram shows the changes in the expression of C-peptide protein in the S0 and S4 stages. The C-peptide was expressed in the cytoplasm. The picture showed that the Cpeptide was almost not expressed in the S0 stage, however, expressed in the S4 stage, which is consistent with the function of mature beta cells. The F-G diagram shows the changes in the expression of insulin protein at the S0 and S4 stages. Insulin was expressed in the cytoplasm. The gure shows that insulin was almost not expressed in the S0 stage but expressed in the S4 stage, which is consistent with the gene expression and the function of mature beta cells. Figure 3ii shows the supernatant that collected from UC-MSCs (S0) and IPCs (S4) at 22 days and incubated in low-and high-glucose solutions for 2 h. An insulin release assay was performed using a human insulin ELISA kit (Alpco, Salem, USA). The results showed that the amount of insulin released by the cells in the S4 stage in the 5.5 mM glucose solution was higher than that in the S0 stage (P<0.05). The amount of insulin released by the cells in the S4 stage in the 25 mM glucose solution was also higher than that in the S0 stage, and the difference was statistically signi cant (P<0.05). The amount of insulin released in the S4 stage in the high-glucose 25 mM glucose solution was signi cantly higher than that in the S4 stage in the 5.5 mM glucose solution (P<0.05). All of the above results indicated that functional islet beta-like cells were successfully induced.

Identi cation and global view of lncRNA pro les
We generated cDNA libraries of polyadenylated RNA extracted from ten samples and RNA-seq datasets at a sequencing depth of 8.6 million reads per sample and the average number of reads entering the mapping process across all analyzed samples was 84.9 million. The average percentage of uniquely mapped reads was 92.81%, and the average percentage of reads mapped to multiple loci was 3.09% (Table S1). Some results of the quality control of RNA-seq datasets are shown (Fig. S1A-C). We aligned the ltered reads to the reference sequence (human GRch38 genome) by Hisat2, predicted lncRNAs and calculated their expression level using Cu inks. The number of annotated lncRNAs and novel lncRNAs of ve stages were pro led. Among them, 1053 annotated lncRNAs and 48 novel lncRNAs were continuously expressed in 5 stages ( Figure 4A). LncRNAs obtained by screening included lincRNAs (30.4%), antisense lncRNAs (15.8%), and sense-overlapping lncRNAs (53.8%) (Fig. S1D). Most lncRNAs contained approximately 1 or 2 exons and no more than 10 exons; however, most protein-coding transcripts contained more than 10 exons ( Figure 4B). In addition, a length of 0-200 bp represented the dominant portion of the exon length of lncRNAs and protein-coding transcripts ( Figure 4C). The length of the dominant portion of lncRNAs was mainly 1000 bp, while the length of protein-coding transcripts was nearly 2500 bp ( Figure 4D). Overall, lncRNAs have fewer exons than protein-coding transcripts, and the lengths of these lncRNAs were generally shorter than those of mRNAs. To understand the quality control of experimental data, a Pearson correlation analysis was performed for all pairs of RNA-seq samples, and it demonstrated that the correlation coe cient between each repeated sample reached 0.8, indicating that the overall quality of the data was better ( Figure 4E, Fig S2. A).

Differentially expressed lncRNAs and functional enrichment of target genes
To identify the differentially expressed lncRNAs during IPCs differentiation of UC-MSCs, the number of lncRNAs was detected in ve stages. A total of 2615 lncRNAs with differential expression were identi ed, of which 1101 lncRNAs were continuously expressed in ve stages (Figure5A). To further explore the change of dynamic expression of differentially expressed lncRNAs in ve stages during induction, the S0 stage was compared with the S1, S2, S3, and S4 stages, the S4 stage was compared with the S1, S2, and S3 stages, and the S3 stage was compared with the S2, and S1 stages, the S2 stage was compared with the S1 stage. The up-or down-regulated lncRNAs were identi ed, which revealed that the number of differentially expressed lncRNAs among them decreased successively. Moreover, mRNAs showed the same expression trend (Figure5B, Fig. S2B). The overlap of differentially expressed lncRNAs was shown in S1, S2, S3, and S4 stages compared with S0 stage. Among them, 30 lncRNAs were continuously upregulated, 38 lncRNAs were continuously downregulated in the ve stages. The overlap of differentially expressed mRNAs is shown ( Figure 5C, Fig. S2C).
To identify the function of differentially expressed lncRNAs that controlled IPCs differentiation of UC-MSCs, we performed a Gene Ontology (GO) analysis and pathway analysis. Since lncRNAs do not encode proteins, their regulatory effect can only be exerted by regulating the genes coexpressed with them. First, the GO analysis was performed to enrich the signi cant functions of the genes coexpressed with differentially expressed lncRNAs, which were obtained during IPCs differentiation of UC-MSCs. We obtained the top 20 GO functions according to the p value and FDR (p<0.05, FDR<0.05) in the S1, S2, S3 and S4 stages compared with S0 stage. Coexpressed genes of highly expressed lncRNAs in the S1, S2, S3 and S4 stages were mainly concentrated in protein translation-related pathways. In contrast, the coexpressed genes of lncRNAs with low expression were mainly enriched in the cell adhesion and cell connection pathways ( Figure 5D). The top 20 GO functions of differentially expressed mRNAs are shown (Fig. S2D).
We then enriched the signi cantly changed pathways that mediated the functions of the genes coexpressed with differentially expressed lncRNAs based on the KEGG database. The top 5 signi cant pathways in stage S4 compared with stages S1, S2, and S3 were related to the insulin signal transduction pathway of beta cells, and they may play key roles in the function of pancreatic beta cells. Among them, upregulated lncRNAs were mainly enriched in the NF-KB signaling pathway and MAPK signaling pathway, and downregulated lncRNAs were mainly enriched in the HIPPO signaling pathway,PI3K-Akt signaling pathway, p53 signaling pathway and other signaling pathways. In summary, the KEGG analysis revealed that differentially expressed lncRNAs were related to the insulin signal transduction pathway of beta cells (Figure5E).
3.8 The expression pattern of lncRNAs in S1, S2, S3, and S4 stages were demonstrated by K-means analysis.
To characterize the dynamic changes in lncRNA and mRNA expression, we clustered all their expression patterns (97 lncRNAs and 769 mRNAs) by K-means analysis. We identi ed 4 lncRNAs and mRNAs clusters. The upper part of the module showed four clusters identi ed by K-means, and the lower part showed the lncRNA and mRNA expression pro les of four clusters (Figure6A, Fig. S3A). The lncRNA expression patterns in the four clusters are displayed, among which the change trend of lncRNAs in Clusters 1 and 3 was consistent with that of the mRNAs ( Figure 6B, Fig. S3B). The top 5 signaling pathways of coexpressed genes with lncRNAs in four different clusters are displayed (Figure6C). Moreover, we found that mRNA Cluster 1 was mainly enriched in the cell proliferation and cell differentiation pathways, indicating that cell differentiation tended to be active, Cluster 3 was mainly concentrated in the negative regulation of cell proliferation and cell adhesion pathway. Thus, with the development of cell differentiation, the ability to inhibit cell proliferation and the connection between cells became weaker, namely, cell differentiation became increasingly active, which was in line with the process of cell differentiation according to the functional enrichment analysis of the four mRNA Clusters (Fig. S3C). The expression levels of the P53, FIGF, JARID2 and PROX1 genes, which are involved in the development and function of pancreatic beta cells in Clusters 1 and 3 were determined (Fig. S3D). To explore the functions of lncRNAs, the network interactions of lncRNAs and their coexpressed mRNAs in Clusters 1-3 and the enrichment pathways related to cell differentiation are shown( Figure 6D).

Screening of pancreatic development-related lncRNAs and their co-expressed mRNAs
Since the focus of our study was to investigate lncRNAs associated with IPCs differentiation, we screened the mRNAs JARID2, PROX1 and p53, which are related to the development and function of the pancreas in the functional enrichment pathway. Cervanets et al. found that JARID2 plays a role in the late differentiation of embryonic pancreatic beta cells. Pethe et al. found that the expression levels of JARID2 decreased gradually during directed differentiation and spontaneous differentiation of the pancreatic system. Studies have also found that PROX1 is a speci c marker of the endoderm in early pancreatic development. Other studies have found that PROX1 is highly expressed in pancreatic endocrine progenitor cells, but is lacking in mature beta cells. PROX1 activity is essential for the formation of endocrine progenitor cells and the differentiation of α cells [37][38][39][40][41]. Kung et al. demonstrated that p53 regulates insulin secretion and pancreatic beta cell survivals through multiple signaling pathways [42]. The coexpression of the lncRNAs :AC009014.3 and Gs1-72m21.1 with JARID2 and lncRNA CTBP1-AS2 with PROX1 was identi ed through the coexpression analysis (Figure7BDE). Which indicated that the lncRNAs : AC009014.3, Gs1-72m21.1 and CTBP1-AS2 may be involved in the development of pancreatic beta cells. The lncRNAs: XLOC_050969, LINC00883, XLOC_050981, XLOC-050925, MAP3K14-AS1, RP11148K1.12, and CTD2020K17.3 were coexpressed with p53 and may be involved in regulating insulin secretion by pancreatic beta cells (Figure7AC).

Discussion
Long noncoding RNAs (lncRNAs) are the largest component of the mammalian noncoding transcriptome, with a length greater than 200 nt [43]. According to the genomic position of lncRNAs relative to adjacent or overlapping protein-coding genes, they can be divided into sense, antisense, intron, intergenic or enhancer lncRNAs that which mediate short-and long-range interactions between transcriptional enhancers and other regulatory elements in the genome [44]. lncRNAs can control the expression of cis or trans genes by directly interacting with transcription factors or recruiting chromatin modi cation complexes [45][46][47]. The expression levels of lncRNAs are generally lower than those of mRNAs, although a large proportion of the identi ed lncRNAs have high islet speci city, with a dynamic modulated pattern observed in differentiated beta-like cells in vitro, and most detected lncRNAs have increased speci city in functional endocrine differentiated cells inhibited insulin secretion and reduced the insulin content in Endoc cells stimulated by glucose, and the expression levels of many lncRNAs were highly correlated with transcription factors (GLIS3,HNF1A,NKX2.2,PDX1, and MAFB) in beta cells [49]. Therefore, lncRNAs play important roles in the development of pancreatic beta cells and present regulatory functions in pancreatic beta cells.
In this study, differentially expressed lncRNAs and mRNAs were obtained during the differentiation of stem cells into insulin-producing cells. GO analysis of these differentially expressed lncRNAs showed that the functional enrichment pathway of the upregulated genes was mainly enriched in translation-related pathways, such as the initiation, extension and termination of translation, and downregulated genes were mainly concentrated in the cell connection and adhesion and extracellular matrix tissue pathways.
Differentially expressed lncRNAs were enriched in the S4 stage compared with the S1, S2, and S3 stages based on the KEGG analysis. Many target genes were involved in the NF-KB signaling pathway, MAPK signaling pathway, HIPPO signaling pathway, PI3K-Akt signaling pathway, and p53 signaling pathway. In a study of beta cells treated with free fatty acids that mimic pancreatic dysfunction and apoptosis induced by gut obesity, activation of p53 led to the induction of miR34a, a microRNA that sensitizes beta cells to apoptosis and inhibits the insulin exocytosis pathway, thus leading to impaired insulin secretion [50]. It has been shown that in beta cells, elevated glucose levels lead to elevated Ca 2+ levels and enhanced NF-KB activity, which promotes insulin release. Overexpression of IKBa overexpression inhibits intracellular NF-KB activity and leads to impaired glucose-stimulated insulin secretion [51]. Studies have shown that PI3K-Akt is one of the major signaling pathways that maintains the survival and replication of beta cells and the expression and secretion of insulin genes [52][53][54]. The above discussion indicates that these pathways are related to the insulin secretion of islet beta cells.      The expression of LncRNAs at each stage and functional enrichment of target genes of differentially expressed lncRNAs. A Venn diagrams of LncRNAs detected in samples of 5 stages, the p value or corrected p value (padj) was used to determine the level of signi cance for multiple samples; B Bar plots of differential LncRNAs at different stages; C the Venn diagram showed the overlap of differentially upregulated or down-regulated LncRNAs (S0 stage was compared with S1,S2,S3 and S4 stages, respectively); D S0 stage was compared with S1,S2,S3 and S4 stages, respectively, and GO(molecular process) enrichment analysis of the mRNA of differentially up-regulated or down-regulated LncRNAs, the color scale showed the signi cance of these terms scaled by column(-log10 corrected p value); E S4 stage was compared with S1,S2,and S3 stages, respectively, and KEGG enrichment analysis of the mRNA of differentially up-regulated or down-regulated LncRNAs, the color scale showed the signi cance of these terms scaled by column(-log10 corrected p value).

Figure 6
The expression pattern of LncRNAs in S1 stage to S4 stage was demonstrated by K-means analysis. A Differential LncRNAs in S1 stage to S4 stage were clustered by K-means; B LncRNAs expression pro les generated by K-means clustering; C the top 5 GO terms(molecular processes) of co-expressed mRNAs with 4 clusters of differential LncRNAs generated by K-means clustering, and the color scale showed the signi cance of these terms scaled by column(-log10 corrected p value). D Co-expression network of coregulated LncRNAs in gene clusters 1 to 3 and co-expression mRNAs participating in 9 GO terms. Figure 7