Transcriptome Analysis to Identify the Potential Role of Long non-coding RNAs in Enteric Glial Cells under Hyperglycemia

Background Long non-coding RNAs (lncRNAs) are important mediators in the pathogenesis of diabetic gastrointestinal autonomic neuropathy, which has just been reported to have a relation to enteric glial cells (EGCs). However, the role of lncRNAs in the pathogenesis of diabetic gastrointestinal autonomic neuropathy, especially EGCs-related gastrointestinal dysfunction, has never been reported. Methods RNA sequencing technology (RNA-Seq) was used to screen the differential lncRNAs and mRNAs in EGCs under hyperglycemia (300 mmol L − 1 high glucose). were obtained. GO enrichment analysis and KEGG pathway analysis showed signicant differences. 2910 and 1549 co-expressed mRNAs were respectively expressed in up-regulated and down-regulated DE lncRNA target genes. Several up- or down-regulated lncRNAs were at the key junction points of the regulatory network. Protein-protein interaction networks showed highly connected clusters were TP53, AKT1, Casp9, Casp8, Casp3, TNF, etc, which are known closely related to apoptosis. FLRT3, Fras1, and other related target genes, which revealed the potential function of lncRNAs, may be important targets for differential lncRNAs to regulate the apoptosis of glial cells induced by hyperglycemia.


Introduction
Diabetes Mellitus (DM) is one of the most serious chronic diseases in the world and a continuous rise in DM prevalence is expected. According to World Health Organization statistics, > 422 million adults globally were suffering from DM since 2014 [1] and 8.8% of the global adults have been diagnosed with diabetes [2]. Diabetic gastrointestinal autonomic dysfunction is a common complication of DM, which occurs in 50-55% of diabetic patients. Diabetic gastrointestinal autonomic dysfunction shows various gastrointestinal symptoms, such as abdominal pain or discomfort, abdominal distention, nausea, vomiting, diarrhea, constipation, etc [3][4][5]. Intestinal sensory and motile gastric dysfunction is the main synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After the adenylation of 3' ends of DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were puri ed with the AMPure XP system (Beckman Coulter, Indianapolis, IN, USA). Then 3 µl USER Enzyme (New England Biolabs, Ipswich, MA, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were puri ed (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2100 system [33].

Clustering and sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, California, USA) according to the manufacturer's instructions [34]. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated.

Quality control
Raw data (raw reads) of FASTQ format were rstly processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. All clean reads were aligned to the reference genome and gene model annotation databases using bowtie2 (v2.2.8) and HISAT2 (v2.0.4) [35], respectively. The aligned reads from all samples were assembled by StringTie (v1.3.1) [34]. At the same time, Q20, Q30 and GC content the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Reads mapping to the reference genome
Reference genome and gene model annotation les were downloaded from the genome website directly. The index of the reference genome was built using STAR and paired-end clean reads were aligned to the reference genome using STAR (v2.5.1b) [36]. STAR used the method of Maximal Mappable Pre x (MMP) which can generate a precise mapping result for junction reads.

Quanti cation of gene expression level
HTSeq (v0.6.0) was used to count the reads numbers mapped to each gene [34]. And then FPKM (fragments per kilobase of transcript per million fragments mapped) of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions of base pairs sequenced, considers the effect of sequencing depth and gene length for the reads count at the same time, and is currently the most used method for estimating gene expression levels.
Differential expression analysis Before differential gene expression analysis, for each sequenced library, the read counts were adjusted by the edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the edgeR R package (3.12.1). The P values were adjusted using the Benjamini & Hochberg method [37]. Corrected P-value of 0.05 and absolute foldchange of 2 was set as the threshold for signi cantly differential expression. Heat map clustering analysis of all differentially expressed lncRNAs and mRNAs was calculated using the Heatmaps R package.

GO and KEGG enrichment analysis
Gene Ontology (GO) enrichment analysis was implemented by the cluster pro ler R package, in which gene length bias was corrected (http://geneontology.org/). GO terms with corrected P-value less than 0.05 were considered signi cantly enriched by differential expressed target genes. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used the cluster pro ler R package to test the statistical enrichment of differentially expressed target genes in KEGG pathways. mRNA-lncRNA-miRNA co-expression network analysis A lncRNA-miRNA-mRNA co-expression network was constructed based on the differentially expressed lncRNAs, mRNAs and miRNAs [38], and the interactions in the network were analyzed using Cytoscape software (http://www.cytoscape.org/). Different colors represent different types of RNA.
Prediction and functional enrichment analysis of lncRNA target genes To reveal the potential function of lncRNAs, their target genes were predicted in Cis. For Cis role, it refers to the action of lncRNA on neighboring target genes. In this study, the coding genes from 100 kb upstream and downstream of a lncRNA were searched.

Real-time quantitative PCR (RT-qPCR)
To validate the RNA-Seq results, cDNA was synthesized by reverse transcription using 1μg of the original RNA sample as described previously. RNAs were reverse-transcribed using the PrimeScript™RT reagent Kit (Perfect Real Time) (TAKARA, Dalian, China) according to the manufacturer's instructions. cDNA was analyzed by using TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) (TAKARA, Dalian, China). GAPDH was chosen as the reference gene. The threshold cycle (CT) value was detected by the 7500 Fast Real-Time PCR system (Applied Biosystems, CA, USA) [40]. The expression levels were computed according to the 2 -ΔΔCt method.

Statistical analysis
GraphPad software Prism (5.0) was used for statistical analysis. Differences between the treatment groups were analyzed by one-way ANOVA with Bonferroni's or Dunnett's posttests. A P-value of less than 0.05 was considered signi cantly enriched.

Hyperglycemia induces cytotoxicity and cell apoptosis in EGCs
To investigate the effect of hyperglycemia on EGCs, cells were cultured with different concentrations of D-(+)-Glucose (0, 50, 100, 200, 300, 400 mM) for 24h. Cell viability determined by the MTT assay was used to measure metabolic activity (intracellular NADH concentration) of the cells. As shown in Figure 1A, G300 (Glucose300) treatment induced the viability of cells from 100.00% ± 8.85% to 74.07% ± 8.71%, thus which was selected as representative concentration. The signi cance was determined by the ANOVA test (*P 0.05, **P 0.01, ***P 0.001, ****P 0.0001 compared with G0). Then, cell apoptosis was determined by the Tunnel assay. As shown in Figure 1B, G300 obviously induced shrinkage, and the formation of apoptotic bodies was observed, while the apoptosis of M300 was signi cantly less than that of G300, which meant that apoptosis was caused by hyperglycemia, not by hypertonic. Tunnel assay showed that typical apoptotic cells were labeled with red uorescence (marked by the white arrows).

Sequencing data overview in EGCs under hyperglycemia
A total of 104821486 and 116495776 raw reads were produced using the Illumina PE150 platform in G0 and G300 groups. After discarding adaptor sequences and low-quality sequences, 101754234 and 111958252 corresponding clean reads were obtained, and the percentage of clean reads was 97.07% and 96.10% (Table 1). The raw bases, clean bases, and the percentage of Q20, Q30, GC content of the cleaned reads were calculated and shown in Table 2. The percent of reads to genome regions were shown in Figure 2A and the results revealed that most clean reads in the two periods were located in the exonic region, and about 20% of the total clean reads were in the intergenic region of the genome. The percent of the intronic region in G0 and G300 was 16.23% and 9.03%. The transcript length, exon number, and ORF were calculated and shown in Figure 2B. Most of the lncRNA lengths were less than 2500, while mRNAs were distributed between 0 and > 1000. The exons and ORF lengths of lncRNAs and mRNAs also showed similar patterns with length. The exon lengths of lncRNAs were less than 10, while that of mRNA was between 0 and 30. The ORF length of lncRNAs was less than 500, while most of the exon lengths of mRNAs were evenly distributed between 0-2000. It showed the difference between lncRNAs and mRNAs and it veri es that the predicted novel lncRNAs conformed to the general characteristics.

Differential expression of lncRNAs and mRNAs in EGCs under hyperglycemia
The volcano plots were constructed to show the rapid evaluation of differentially expressed lncRNAs and mRNAs as well as their statistical signi cance ( Figure 3A). Each point in the volcano plot represented a lncRNA or an mRNA, and the abscissa represented the logarithm of the fold change (FC) in the expression of RNAs between the 2 groups. Higher absolute values were correlated with larger differences. Downregulated lncRNAs were represented in green, and up-regulated lncRNAs were shown in red.  Figure 3B). Red indicated high relative expression, and blue indicated low relative expression. The samples were grouped by hierarchical clustering, and a circus plot was generated to show the location and expression values of the identi ed lncRNAs and mRNAs in two groups ( Figure 3C). The outermost ring indicates the genome, the middle one indicates the mRNAs expression, and the innermost circle showing the lncRNAs expression. In all, according to the transcriptional level, there were signi cant differences between the two groups.

GO enrichment analysis of DE mRNAs and DE lncRNA target genes
To investigate the probable function of candidate lncRNAs, the nearby coding genes around 100kb up and downstream from the identi ed lncRNAs, termed as Cis-regulatory function, were searched. Then functional enrichment analysis of the signi cantly regulated genes was performed based on their relative expression and fold changes. As shown in Figure 4A Figure 4A showed GO enrichment analysis of differentially expressed genes (DEGs), such as multicellular organismal development (GO:0007275), single-organism process (GO:0044699), and protein binding (GO:0005515), etc. Figure 4B showed GO enrichment analysis of DE lncRNA target genes, such as metabolic process (GO:0008152), protein binding (GO:0005515), organic substance metabolic process (GO:0071704), etc.

KEGG enrichment analysis of DE mRNAs and DE lncRNA target genes
The top 20 enriched KEGG pathways obtained from the DE mRNAs or DE lncRNA target genes were listed in Figure 5A-B (Padj<0.05). The ECM-receptor interaction, Hypertrophic Cardiomyopathy (HCM) and Transcriptional Mis-regulation in cancer were identi ed in the KEGG pathway analysis as containing the highest rich factor in DE mRNAs ( Figure 5A). The Non-small cell lung cancer, estrogen signaling pathway and endometrial cancer were identi ed in the KEGG pathway analysis as containing the highest rich factor in DE lncRNA target genes ( Figure 5B). Apart from the basic cellular activities, other pathways, including MAPK signaling pathway, T cell receptor signaling pathway and Base excision repair, have been proved to be involved in cell metabolism, synthesis, immune response, apoptosis and other important life activities [22,41]. DE mRNAs and DE lncRNA target mRNAs were shown in the Venn diagrams. 2910 and 1549 co-expressed mRNAs were respectively expressed in up-regulated and down-regulated DE lncRNA target genes ( Figure 5C-D).
Then, 19 KEGG pathways related to apoptosis, such as apoptosis, Jak-STAT signaling pathway, PI3K-Akt signaling pathway, NF-kappa B signaling pathway were analyzed. We listed the expression of DE mRNAs and DE lncRNA target genes in each pathway and marked the co-expressed targets (Table 3). Coexpression of DE mRNAs and DE lncRNA target genes in these pathways can be considered as important targets for apoptosis, such as Bcl2, Tp53, Wnt4, TRADD, Bax, TNF, which are widely believed to be involved in apoptosis regulation.

mRNA-lncRNA-miRNA interaction network analysis
Co-expression is con rmed to be able to predict and evaluate the function of lncRNAs, thus mRNA-lncRNA-miRNA co-expression network with statistically DE mRNA, lncRNA and miRNA was built ( Figure  6A). Figure 6B showed 161 up-regulated lncRNAs, 23 up-regulated mRNAs and 11 down-regulated miRNAs were included in this interaction network. Figure 6C showed 180 down-regulated lncRNAs, 79 down-regulated mRNAs and 13 up-regulated miRNAs were included in this interaction network. Among this network, several up-regulated or down-regulated lncRNAs were at the key junction points of the regulatory network, which can be used as important targets for further research, such as TCONS_00142807, TCONS_00050289, TCONS_00008543, ENSRNOT00000081689 and TCONS_00089281. Several differentially expressed mRNA at network nodes can also be used as key targets for research, such as ENSRNOT00000071420 (Tk1), ENSRNOT00000002119 (RGD1563888), ENSRNOT00000066211 (Tmem243), ENSRNOT00000026392 (Adprhl1), ENSRNOT00000010619 (Timm10), etc.

Protein-Protein interaction network
The target genes of DE lncRNAs related KEGG apoptosis pathway was mapped into PPI. Figure 7A showed the combined bioinformatics analysis of DE mRNA genes and target genes of lncRNAs. Figure   7B showed the PPI between con uent proteins. The redder and larger circle displayed the more degree connection in the network. The bluer and smaller the circle displayed the less degree connection in the network. The thicker the line and the redder the color, the more degrees were connected in the network. As shown in the network, highly connected clusters were TP53, AKT1, Casp9, Casp8, Casp3, TNF, etc.

DE lncRNA target genes prediction
The top 5 up-and down-regulated differentially expressed lncRNAs were selected for predicting the potential target genes ( Table 4). The DE lncRNAs were ENSRNOT00000092009, TCONS_00073996,  TCONS_00113453, TCONS_00053590, TCONS_00102805, ENSRNOT00000086889,   ENSRNOT00000086208, TCONS_00022155, TCONS_00090467, TCONS_00059568. Their potential target genes, including Fras1, Sim1 and Flrt3, have been reported to be involved in the physiological and pathological functions of immunity, barrier protection and apoptosis.

Real-Time quantitative PCR analysis
In order to verify the repeatability and reliability of lncRNAs identi ed by transcriptome sequencing, we randomly selected four lncRNAs from Table 4, namely TCONS_00053590, TCONS_00113453, ENSRNOT00000086889 and ENSRNOT00000086208, for RT-qPCR analysis ( Figure 8A). Several of their potential target genes, which contained Flrt3, Fras1, Mrpl1, Sim1, were also veri ed by RT-qPCR analysis ( Figure 8B). The above mentioned lncRNAs and mRNAs, which were potentially important targets for apoptosis regulation at the intersection of the network graph, have also been veri ed ( Figure 8C-D). GAPDH was used as an internal control and relative quantity of expression (fold change). Each lncRNA or mRNA was calculated with the comparative 2 -ΔΔCT method and the results were consistent with the sequencing data. Values of RT-qPCR shown were mean with Mean±SD. Each RT-qPCR validation was repeated three times. The results showed that these RNAs showed signi cant differential expression, and was consistent with the change of transcriptome sequencing, indicating that the transcriptome sequencing was reliable.

Discussion
Diabetic gastrointestinal autonomic neuropathy is a common complication of diabetes mellitus. The delayed gastric emptying caused by diabetic gastrointestinal autonomic neuropathy will make the caloric intake and absorption of diabetic patients di cult to match with insulin secretion, which will have adverse effects on blood glucose control [42]and may be a potential mechanism for the development of brittle diabetes [43]. The damage of the ENS in diabetic patients has been proved to be closely related to diabetic gastrointestinal autonomic neuropathy [6,7]. Studies have shown that ENS is particularly vulnerable to hyperglycemia [44]. When the blood glucose uctuates obviously (brittle hyperglycemia) or hyperglycemia continues, glucose metabolism in cells changes, which leads to the formation of advanced glycation end products. This eventually leads to cell damage and cell death through oxidative stress, in ammation and other pathways, which is often referred to as glucose neurotoxicity [45][46][47].
As an important part of the ENS, EGCs have been proved to play an important role in the occurrence and development of gastrointestinal dysfunction [5,9,10]. EGCs can protect neurons by releasing reduced glutathione and glial cell-derived neurotrophic factor (GDNF) under hyperglycemia activation [48,49]. However, when EGCs apoptosis through oxidative stress and other pathways, the protective effect of EGCs on neurons is weakened [48]. Lopes et al. observed a reduction of the number of EGCs decreased in intermuscular and submucosal plexus in the diabetic rats [6]. It has been shown that EGCs begin to proliferate in the ileum of diabetic mice 4 weeks after the model establishment. The number of EGCs decreased after 8 weeks compared with 4 weeks [28]. However, few articles clearly reveal the mechanism of EGC apoptosis induced by hyperglycemia. In this study, we found G300 treatment could obviously induce cytotoxicity and cell apoptosis; thus G300 was selected as induced concentration to further perform transcriptome analysis.
Increasing evidence has shown that lncRNAs are key mediators in the pathogenesis of diabetes mellitus and its complications [50,51]. LncRNA-ANRIL participates in the development of diabetes mellitus by promoting β-cell proliferation and inhibiting the expression of CDKN2A. CDKN2A is the key regulator of CREB induced hepatic gluconeogenesis [52,53]. Li et al. con rmed that lncRNA-Sox2OT gene knockout can protect retinal ganglion cells from hyperglycemia-induced injury and play a neuroprotective role in diabetic retinopathy [54]. LncRNA-MALAT1 has been reported to be involved in the pathogenesis of diabetes and its complications [55]. You et al. have found that downregulation of lncRNA-Meg3 affects Insulin synthesis and secretion in mouse pancreatic β cells [56]. However, the role of lncRNAs in EGCs under hyperglycemia has not been reported. In this study, the involvement of lncRNAs in EGCs under hyperglycemia was analyzed using transcriptome analysis. To reveal the potential function of lncRNAs, their target genes were predicted in Cis in this study. Cis acting lncRNAs regulate gene expression in a manner dependent on their own transcriptional sites [57,58]. Cis acting lncRNAs have been proved to inhibit, activate or via other mechanisms to regulate the expression of target genes [59][60][61].
Bioinformatics analysis of lncRNAs target genes in cis suggested that they may participate in some important KEGG pathways, such as T cell receptor signaling pathway, MAPK signaling pathway and Gap junction. Sauer et al. found that T cell receptor signaling can regulate cell function via PI3K/Akt/mTOR pathways [62], which can mediate cell apoptosis as well. It was con rmed by Chyuan et al. that Tumor necrosis factor-related apoptosis-inducing ligand (TRAIL), a ligand which can induce cell apoptosis by transducing apoptosis signals after interacting with its receptor (TRAIL-R), can inhibit T cell activation and in ammatory response through T cell receptor signal pathway [41]. MAPK signaling pathway has been known as one of the key pathways regulating apoptosis [30,63,64]. Du et al. found that inhibition of gap junctional communication can induce cardiomyocyte apoptosis through the mitochondrial pathway [65]. Another study showed that interference of gap junctional communication can induce the formation of ovarian reactive oxygen species and further mediate cell apoptosis [66].
Further, Table 3 showed the expression of DE mRNAs and DE lncRNA target genes in apoptosis related KEGG pathways. 19 KEGG pathways were shown in the table. Co-expressed targets, such as Bcl2, Bax, TP53, Casp9, Casp3 and TNF, are widely believed to be involved in apoptosis regulation. And the target genes of DE lncRNAs related KEGG apoptosis pathway were mapped into PPI (Fig. 7). It has been found that TRADD, the key adaptor molecule of TNF-α signal transduction, plays a key role in the regulation of NF-κB [67]. TRADD is a key signal intermediate connecting TNF-α and NF-κB activation. TRADD adaptor protein binds to the death domain of the receptor, and then forms a multi-protein complex. Finally, it leads to the phosphorylation of IkBα, which makes NF-κB transport to the nucleus and stimulates cytokine transcription [68]. As well, TRADD was found to mediates islet β cell apoptosis in diabetes [69,70]. BIRC family, also known as the inhibitor of apoptosis proteins (IAPs), consists of eight members [71]. BIRC2 and BIRC3 proteins indirectly regulate caspase activation through E3 ligase activity, TNF signal transduction and NF-kB signal transduction [72][73][74][75]. Through the analysis of TNF-α protein network, Jamil et al. found that TNF-mediated proteins BIRC2 and BIRC3 were involved in the pathogenesis of T2DM [76]. RelA/p65 is considered to be a key subunit of NF-κB signal transduction, which is a key target of mitochondrial metabolism and apoptosis [77]. It was known that the regulatory networks of lncRNAs, miRNAs, and mRNAs, also called competing endogenous RNAs (ceRNAs), communicated with each other to regulate gene expression [81]. Studies suggested that ceRNAs can regulate the occurrence and development of diabetes and its complications through protein synthesis, ER stress, RNA binding and protein translation [82]. Li et al. have shown that targeted ceRAN is involved in the regulation of apoptosis induced by hyperglycemia, which provides a new treatment strategy for diabetes and its complications .
In this study, several DE lncRNAs and mRNAs were at the key junction points of the regulatory network.
TCONS_00012166-miR-362-5p-ENSRNOG00000046202 (Metrnl) axis was found to be up-down-up change. Metrnl is a recently discovered adipokine highly expressed in subcutaneous adipose tissue [40], which contains functions such as promoting nerve development [84], regulating the immune system and so on [85,86]. It can regulate triglyceride, cholesterol, low density lipoprotein, high density lipoprotein and other blood lipid components [87]. Chung et al. demonstrated that in human subjects with T2DM, Metrnl levels were elevated and negatively correlated with various metabolic risk factors [88], which has been proved to work by antagonizing insulin resistance through multiple signaling pathways [89,90].
TCONS_00083285-miR-187-3p-ENSRNOG00000018233 (Gas6) axis was found to be down-up-down change. Growth arrest speci c-6 (Gas6) can be linked to phosphatidylserine on the surface of apoptotic cells, thus triggering effector cell action [91]. Nepal et al. found that macrophages can remove apoptotic neutrophils and eliminate in ammation by inducing Gas6 expression [92]. In addition, we selected 5 mRNAs and 5 lncRNAs at the intersection of the network graph and veri ed them by RT-qPCR (Fig. 8C-D).
As shown in the Table 4, we screened out the top 5 up-and down-regulated DE lncRNAs with potential target genes. FLRT3, an axon guidance-related factor, has been found to be related to the regulation of nerve cell growth and morphogenesis [93]. FLRT3 has been proved to regulate the growth of nerve cells after peripheral nerve injury [94,95]. In our results, TCONS_00113453 was signi cantly down-regulated in G300 group. It can promote the repair of nervous system damage and the growth of synapse by regulating potential target gene FLRT3. FRAS1 is related to the regulation of extracellular matrix composition, maintenance of basement membrane integrity, epithelial adhesion and signal transduction [96,97]. Nikolova et al. reported that FRAS1-related extracellular matrix 3 (FREM3)[98] may participate in cell-cell interaction and maintain the structural and functional integrity of nerve tissue. As shown in our results, TCONS_00053590 may up-regulated its potential target gene FRAS1 through down regulation in G300 group, which may play a role in extracellular matrix regulation and neuronal protection. In our study, the expression changes of randomly selected lncRNAs and several potential target mRANs have been veri ed by RT-qPCR ( Fig. 8A-B), showed a strong correlation with that identi ed with RNA-Seq.

Conclusion
In Ethics approval Not applicable.

Consent for publication
Not applicable.

Con icts of interest
The authors declare no con icts of interest.      Con uent proteins of DE mRNA genes and target genes of lncRNAs related KEGG apoptosis pathway.