Analysis of long non-coding RNA expression profiles in vascular endothelial cells with hyperglycemia-induced damage

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. Results  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, 1,167 mRNA nodes and 9,735 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. recombination, Colorectal cancer, Mucin type O-Glycan biosynthesis, Arrhythmogenic right ventricular cardiomyopathy (ARVC), Aldosterone-regulated sodium reabsorption, Cardiac muscle contraction, Hypertrophic cardiomyopathy (HCM), Endocrine and other 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 our findings with current of processes by cell injury. of our findings,


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) 4 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 ug/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.

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).

Identification of lncRNAs by using a computational approach
Quality control of the RNA-Seq reads was conducted using FastQC (v0.11.3). Reads were trimmed using the software SEQTK for known Illumina TruSeq adapter sequences, poor reads and ribosome RNA reads. Trimmed reads were aligned to the rat genome (Rn6) using Hisat2 (version: 2.0.4). (27) Transcripts were assembled using Stringtie (v1.3.0) and then were compiled together by gffcompare (v0.9.8), (28,29) Transcripts with class codes "i," "u," and "x," were considered to be potential novel long transcripts. Pfam, (30) CPC (31) and CNCI (32) were used to compute the coding potential of each novel transcript. Transcripts with a Pfam score <0, CNCI <0 and CPC non-significant were considered to lack coding potential. Transcripts were matched with annotation databases, including NONCODE (v5) (http://www.noncode.org) (33) and Ensembl (34) The matched transcripts were considered to be known lncRNAs, and others were considered to be novel lncRNAs. All lncRNAs were quantified using Stringtie. According to the positional association between lncRNA and mRNA in the genome, lncRNA was classified into six types: Bidirectional, exonic_antisense, exonic_sense, intergenic, intronic_antisense and intronic_sense. (35)

Prediction and functional analysis of target genes of differentially expressed lncRNAs
Chromosome localizations, sequence complementarity and correlation coefficients between lncRNA and RNA 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 cisregulatory. (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 co-expression 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 realtime 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, 6 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 the Rattus norvegicus 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 ( Figure 1A, Supplemental Table S1). Most of the identified lncRNAs were transcribed from proteincoding exons (sense and antisense); others were from introns and intergenic regions ( Figure 1B). Furthermore, 24,304 lncRNA transcripts could be found in all chromosomes, with the majority located on chromosome 1 ( Figure 1C).

Identification of differentially expressed lncRNAs
EdgeR was used to filter differentially expressed lncRNAs (DELs) between the HG-induced and CN  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. The results show that the expression of the lncRNAs had similar trends as with the sequencing results, indicating that our sequencing results were reliable (Figure 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  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 ( Figure 4A

lncRNA-mRNA co-expression network
To visualize the co-expression network, pairs of lncRNAs and mRNAs that had PCC > 0.99 and p < 0.01 were assessed using Cytoscape software. As shown in Figure 5 (Fig. 5). In addition, some well characterized lncRNAs, such as SNHG1 (degree = 23) and GACAT2 (degree = 24), were highly represented within the network.

Discussion
Using microarray analysis, a variety of lncRNAs have been determined to be dysregulated in HGinduced 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

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.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. Table S1.xlsx Table S2.xlsx Table S3.xlsx