Analysis of long non-coding RNA expression profiles in high-glucose treated vascular endothelial cells

Background Diabetes mellitus is often associated with microvascular and macrovascular lesions, and hyperglycemia-induced vascular endothelial cell damage is a key factor. Methods We investigated long non-coding RNAs (lncRNAs) and mRNAs that are affected by hyperglycemia-induced damage using human umbilical vein endothelial cells (HUVECs) as a model. HUVECs were cultured under high (25 mmol/L) or normal (5 mmol/L) glucose conditions for 6 d, and then lncRNAs and protein-coding transcripts were profiled by RNA-seq. Result Among 40,379 lncRNAs screened, 214 were upregulated (log2 [fold-change] > 1, FDR < 0.05) and 197 were downregulated (log2 [fold-change] < − 1, FDR < 0.05) in response to high-glucose. Furthermore, among 28,431 protein-coding genes screened, 778 were upregulated and 998 were downregulated. A total of 945 lncRNA/mRNA pairs were identified, including 126 differentially expressed lncRNAs predicted to target 201 mRNAs, among which 26 were cis-regulatory interactions. The corresponding lncRNA-mRNA network was composed of 354 lncRNA nodes, 1167 mRNA nodes and 9735 edges. Dozens of lncRNAs with high degree may play important roles in high-glucose-induced HUVEC damage, including ENST00000600527, NONHSAT037576.2, NONHSAT135706.2, ENST00000602127, NONHSAT200243.1, NONHSAT217282.1, NONHSAT176260.1, NONHSAT199075.1, NONHSAT067063.2, NONHSAT058417.2. Conclusion These observations may provide novel insights into the regulatory molecules and pathways of hyperglycemia-related endothelial dysfunction in diabetes-associated vascular disease.


Background
Diabetes-related microvascular and macrovascular complications are directly correlated with the severity and duration of hyperglycemia [1,2]. Stable and healthy endothelial cells are the basis of normal blood vessels, whereas dysfunction of endothelial cells is a risk indicator of diabetic angiopathy [3,4]. Protein kinase C, hexosamine pathways, polyol pathways, and advanced glycation end-products are believed to be responsible for dysfunction associated with endothelial cell diabetes mellitus. These factors inhibit the production of nitric oxide by promoting the bioavailability of reactive oxygen species, thus altering the structure and physiology of endothelial cells [4][5][6][7][8]. Hyperglycemia is the main contributor of endothelial dysfunction [9]. Advanced glycation end products of hyperglycemia compromise the bioavailability of nitric oxide, which essentially drives endothelial dysfunction [10]. The effect of advanced glycation end products on monocytes, macrophages and vascular smooth muscle cells enhances the inflammatory response and oxidative stress in the endothelial system [9,10]. Thus, unlike other cells or tissues, the vascular endothelium is extremely sensitive to blood glucose and is the direct target of hyperglycemia injury [11].
Long non-coding RNAs (lncRNAs) are non-coding RNA molecules that contain over 200 nucleotides. LncRNAs play important regulatory roles in various diseases, such as cancers [12][13][14][15], atherosclerosis [15], neurodegeneration [11], autoimmune disorders [16][17][18][19][20], and chronic vascular disease [21]. Data from both in vitro hyperglycemia induction and in vivo diabetes mellitus experiments have shown that lncRNA MALAT1 is highly expressed, and that inhibition of its expression can effectively improve endothelial cell inflammation and diabetic retinopathy [22,23]. Furthermore, high glucose (HG) treatment of vascular endothelial cells is accompanied by increased expression of lncRNA MIAT, and interference with MIAT can promote the proliferation, migration and abnormal angiogenesis of vascular endothelial cells induced by HG [24]. In human endothelial cells, the knockdown of lncRNA AGAP2-AS1 inhibits cell proliferation, tubule formation and acetylated LDL uptake [25]. Additionally, HG or oxidative stress inhibits the expression of lncRNA MEG3 in endothelial cells and diabetic mice. Knockout of lncRNA MEG3 can aggravate retinal vascular dysfunction and increase microvascular leakage and inflammation [26]. These results suggest that lncRNAs may comprise novel therapeutic targets in hyperglycemia-related endothelial dysfunction or diabetes-induced vascular disease. In order to further understand the role of lncRNA in endothelial cell injury in diabetes mellitus, we carried out high throughput sequencing using human umbilical vein endothelial cells (HUVECs) cultured under normal (CN) or HG conditions.

Cell culture
HUVECs were purchased from the China Center for Type Culture Collection (CCTCC, Wuhan, China), and were maintained in Modified Eagle's Medium (MEM; Hyclone/ Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, Gaithersburg, MD USA), penicillin (100 U/mL; Gibco), and streptomycin (100 μg/mL; Gibco) at 37°C in a humidified 5% CO 2 chamber. Confluent HUVECs were starved overnight before exposure to media with normal or HG. HUVECs assigned to the Control (CN) group were maintained for 6 d in 5 mmol/L glucose and 20 mmol/L mannitol to account for changes that may be triggered by osmolarity differences. The remaining HUVECs of the HG group were maintained for the same duration in glucose-enriched (final glucose concentration 25 mmol/L) media. There were three biological duplicates in each group of HUVECs.

RNA preparation and RNA-seq
TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate and purify total RNA according to the manufacturer's instructions. The quantity and quality of the RNA were determined using a NanoDrop 2000 instrument (Thermo, Fisher Scientific). The RNA integrity was assessed by electrophoresis with denaturing agarose gels. Libraries were constructed according to the standard TruSeq protocol. Sequencing was performed on the Illumina HiSeq 2500 according to the manufacturer's protocol at Ao-Ji Biotech (Shanghai, China).

Prediction and functional analysis of target genes of differentially expressed lncRNAs
Chromosome localizations, sequence complementarity and correlation coefficients between lncRNA and mRNA pairs were analyzed in order to identify the lncRNAs' cis-and trans-target genes. In brief, a distinction was made according to chromosome location: if lncRNAs were located in the range of 10 KB upstream and 20 KB downstream of the coding gene, they were considered likely to be cis-regulatory [26,27]. The trans-regulatory interaction potential between lncRNAs and mRNAs was analyzed by RNAplex software [28], with binding energy of <− 30 as the threshold. Furthermore, the Pearson correlation coefficients (PCCs) between lncRNAs and mRNAs were calculated, with a cutoff of PCC ≥ 0.6. The functions of these candidate coding genes were assessed using gene ontology (GO) function and pathway analysis using ClusterProfiler.
The lncRNA-mRNA-coexpression network The lncRNA and mRNA co-expression network was constructed according to the normalized fragments per kilobase of transcript per million mapped reads of the unit genes. We calculated the PCC between differentially expressed lncRNAs and differentially expressed mRNAs. LncRNA-mRNA pairs with significant correlations (PCC > 0.99 and p < 0.01) were chosen to build the coexpression network and visualized by Cytoscape. The number of directly linked neighbors for each node was calculated and was defined as the nodes degree.

RT-PCR validation
Seven lncRNAs were randomly selected for verification of the RNA-Seq results by quantitative real-time PCR (qRT-PCR), which was performed on a Roche LightCycler 480 machine (Roche Applied Science, Germany) with the SYBR green assay (TaKaRa, Japan). The qRT-PCR amplification reactions were carried out via the following program: 95°C for 10 min, 40 cycles with 95°C for 15 s and 60°C for 20 s. The primer sequences for qRT-PCR are provided in Table 1. With GAPDH as an internal control, the relative expression was computed according to the 2 −ΔΔCt method.

LncRNA-sequencing data analysis
We characterized the lncRNA landscape of expression by performing deep RNA-seq experiments on three CN and three HG-induced HUVECs samples. After SEQTK quality assessment, more than 33 million total original reads for each sample were obtained, and the proportion of bases with quality values greater than 20 (Q20) was > 94%. These results indicate that the quality of the sequencing results was acceptable ( Table 2). After filtering out the adaptor sequence and low quality reads, the percentage of clean reads within the raw reads accounted for 94% of the total sequences in the two groups. Hisat2 software was used to map the clean reads to homo sapiens reference genome. As shown in Table 2, approximately 97% of the trimmed reads were mapped onto the reference genome. In total, we screened 40,380 lncRNAs from the six samples, including 387 novel and 39,993 known lncRNAs, of which 36,550 were shared lncRNAs detected in both the HG and CN HUVEC groups (Fig. 1a, Supplemental Table S1). Most of the identified lncRNAs were transcribed from protein-coding exons (sense and antisense); others were from introns and intergenic regions (Fig. 1b). Furthermore, 24,304 lncRNA transcripts could be found in all chromosomes, with the majority located on chromosome 1 (Fig. 1c).

Identification of differentially expressed lncRNAs
EdgeR was used to filter differentially expressed lncRNAs (DELs) between the HG-induced and CN HUVEC groups. Among the lncRNAs, 214 were significantly upregulated (log2 (fold-change) > 1, FDR < 0.05) and 197 were significantly downregulated (log2 (foldchange) < − 1, FDR < 0.05) in response to HG exposure (Fig. 2). Additionally, several of the DELs had a fold change value equal to positive infinity and negative infinity, meaning that these lncRNAs were completely switched-on or off with HG induction. The top five  Table S2).

qRT-PCR verification of DELs
To verify our findings, the expression profiles of six differentially expressed lncRNAs were randomly selected for qRT-PCR analysis. There were three repeats per group and five repeats per sample in the qPCR. The results show that the expression of the lncRNAs had similar trends as with the sequencing results, indicating that our sequencing results were reliable (Fig. 3).

Regulatory analysis of DELs and DEGs
lncRNAs act via cis-and trans-regulation of target genes for biological function. To evaluate the regulatory pathways associated with the lncRNAs, we assessed the differentially expressed genes (DEGs) in the same HUVEC samples. Of 28,431 protein-coding genes that were screened, 778 were upregulated and 998 downregulated by HG treatment. By comparing the DELs and the DEGs, a total of 945 matched lncRNA-mRNAs pairs for 126 DELs and 201 DEGs were predicted, of which 26 lncRNA/mRNA interactions were cis-regulatory, with either positive or negative correlations of the lncRNAs with their predicted target genes. An additional 715 interactions were trans-regulatory, including 2 that were both cis-and trans-regulatory (Supplemental Table S3).
To further understand the regulatory functions of the differentially expressed lncRNAs, all predicted target genes were annotated according to GO and pathway function entries using ClusterProfiler. Among the GO Enrichment terms (Fig. 4a) the most abundant in the biological process categories were Mitotic cell cycle, Cell cycle, Cell division, Microtubule cytoskeleton organization, DNA replication, Chromosome segregation, Spindle organization, Cytoskeleton organization, Cholesterol biosynthetic process, and Centromere complex assembly. The most abundant GO terms in the cellular component categories were, Molecular function Centromeric region, Chromosome, spindle, Chromosome, Replication fork, Nuclear chromosome, Condensed nuclear chromosome, Microtubule, Microtubule cytoskeleton, Cytoskeleton, and Nucleoplasm. Among the Pathway Enrichment terms (Fig. 4b), the most abundant were beta-Alanine metabolism, Primary immunodeficiency, Carbohydrate digestion and absorption, Arginine and proline metabolism, Histidine metabolism, Fatty acid elongation, Homologous recombination, Colorectal cancer, Mucin type

Discussion
Using microarray analysis, a variety of lncRNAs have been determined to be dysregulated in HG-induced HUVECs [36]. In addition, mechanistic studies indicate a role for specific lncRNAs in endothelial cell dysfunction induced by diabetes or HG [22][23][24][25][26]. Both clinical and experimental studies indicate that impairment of vascular smooth muscle cells by diabetes and HG contribute to the increased incidence of diabetic cardiomyopathy [37]. Furthermore, a large number of studies have shown that lncRNAs are involved in the injury of vascular smooth muscle induced by diabetes mellitus. For example, lncRNA-ES3 can regulate calcification/senescence of vascular smooth muscle cells through an miR-34c-5p/BMF axis that is activated upon HG induction [38]. The lncRNA SENCR shows high abundance of expression in vascular cells [37,39,40], and can regulate fork-head box protein O1 and transient receptor potential cation channel 6, thereby promoting the proliferation and migration of smooth muscle cells. HG exposure of T2DM db/db mice can inhibit the expression and function in smooth muscle cells, and overexpression of SENCR reverses the inhibitory effect of HG on vascular smooth muscle cells [40]. Given the important role of these lncRNAs, we rationalized that it is likely that other uncharacterized lncRNAs participate in the processes of endothelial cell pathogenesis caused by diabetes. Further, RNA-seq-as a highly sensitive approach towards identifying differentially expressed RNAs-could be used for the purpose of identifying the roles of these yet uncharacterized lncRNAs.
In the present study, a total of 214 upregulated and 197 downregulated lncRNAs (|FC| > 2, FDR < 0.05) were identified through RNA-seq using an HG-induced HUVEC model. Quality control assays were performed to ensure the reliability of the experiments, and several lncRNAs were validated using RT-qPCR. Furthermore, DEGs in the same cells were identified, and the lncRNA and mRNA expression profiles were compared. Using a variety of bioinformatics approaches, we revealed pathways and processes associated with dysregulation of lncRNA expression. It is well known that the AGE-RAGE signaling pathway in diabetic complications plays an important role in endothelial cell injury induced by diabetes mellitus. In Enrichment Analysis based on lncRNA-mRNA co-expression networks, up to 45 lncRNAs were found to be associated with this pathway, especially ENST00000600527, NONHSAT217282.1, NON-HSAT067063.2, NONHSAT093248.2, NONHSAT118785.2, ENST00000444438, NONHSAT200447.1, NON-HSAT108582.2, ENST00000508000, ENST00000594119, NONHSAT186088.1. Therefore, our findings are consistent with current understanding of processes by which diabetes promotes endothelial cell injury. On the basis of our findings, the information obtained in this study should be useful for identifying additional pathways and processes that contribute to the pathogenesis of diabetes. Moreover, we identified some known lncRNAs may participate in high-glucose treated HUVEC. DAPK1-IT1, DAPK1 intronic transcript 1, has proved involved in cancer progression [41,42], chemotherapy insensitivity [43], Atherogenesis [44,45], etc. In our study, we constructed the coexpression subnet of DAPK1-IT1 and function was analyzed with KEGG pathway (Fig. 6), found may function with Terpenoid backbone biosynthesis, Steroid biosynthesis, Vitamin digestion and absorption, Protein processing in endoplasmic reticulum, Phagosome, and some of these pathways have been proved to be involved in the occurrence and  [46]. We will conduct further clinical studies on other lincRNAs in the future.

Conclusion
The potential roles of lncRNAs in diabetic complications were investigated by bioinformatics analysis. These findings may help us to understand the possible molecular mechanism of HG-induced HUVECs and may provide a more comprehensive understanding of the lncRNA expression profile that is dysregulated during diabetes.
Additional file 1: Table S1. List of 40,380 lncRNAs screened from HG and CN HUVEC groups.
Additional file 2: Table S2. The top upregulated and downregulated DELs.