Microarray Analysis of Long Noncoding RNAs in Female Diabetic Peripheral Neuropathy Patients

Diabetes mellitus (DM) causes a wide variety of complications, such as neuropathy, nephropathy, retinopathy, cardiomyopathy, vasculopathy, dermatopathy and encephalopathy [1]. Diabetic peripheral neuropathy (DPN) is characterized by the progressive loss of peripheral nerve axons, resulting in pain, paraesthesia, and sensory loss [2]. It affects approximately 50% of people with a considerable morbidity, mortality and diminished quality of life. Despite the development of diagnostic methods and therapeutic modalities, DPN is still not diagnosed or managed properly in most patients [3].

Over the past two decades, the pathogenesis of DPN has remained the subject of research as well as controversy because of the multifactorial underlying causes of DPN [1]. Previous studies supported the hypothesis that DPN develops as a result of chronic hyperglycemia [4]. However, more recent data seem to indicate that the development of DPN is not necessarily strictly glucocentric but also involves other risk factors that cause nerve and microvascular injury, especially in type 2 DM patients [5].

Long non-coding RNAs (LncRNAs) represent a class of transcripts longer than 200 nucleotides without a protein-coding capacity [6]. LncRNAs function as important modulators of the expression of protein-coding genes located nearby (cis-regulation) or at distance (trans-regulation). Moreover, lncRNAs may also participate in post-transcriptional regulation, the organization of protein complexes, cell-to-cell signaling, and the allosteric regulation of proteins [7, 8]. The misexpression of lncRNAs occurs in a large number of human diseases [9], including diabetes and its complications. Human beta cell transcriptome analyses uncovered lncRNAs that were tissue-specific, dynamically regulated, and abnormally expressed in type 2 diabetes [10]. Moreover, MIAT and PVT1 were found to mediate high glucose-induced renal injury [11, 12], MEG3 was involved in diabetic microvascular dysfunction [13], and the aberrant expression of MALAT1 was found both in diabetic retinopathy and cardiomyopathy [14, 15].

Inspired by these findings, we speculated that lncRNAs might be potential regulators of DPN. In the present study, we used high-throughput microarrays to identify differentially expressed lncRNAs, both common and distinct, between type 2 DM and DPN. The results may provide relevant information for the future development of new mechanism-based diagnostics and therapies.

The participants were recruited from hospitals affiliated with Ningbo University from April 2015 to June 2016. Type 2 DM cases were defined as meeting at least one of the following criteria: a fasting plasma glucose ≥ 7 mmol/L; a 2-h plasma glucose after oral glucose tolerance test (OGTT) ≥ 11.1 mmol/L; or the use of glucose-lowering drugs or physician-diagnosed diabetes. DPN was diagnosed by a positive assessment through neurologic examinations and nerve conduction studies as previously described [16]. Briefly, the criteria for DPN were as follows: (1) confirmed type 2 DM patients; (2) decreased sensation and positive neuropathic sensory symptoms (including pricking, burning, stabbing, or aching pain) in the toes, feet or legs; (3) decreased distal sensation, unequivocally decreased or absent ankle reflexes; and (4) abnormal motor and sensory nerve conduction. The DPN patients involved in the present study were all new cases without any anti-neuropathy medication. Patients with any other major acute or chronic complication associated with diabetes or vitamin B12/folic acid deficiencies were excluded from the DPN group. The healthy controls had no family history of diabetes or neurologic disorders, exhibited normal glucose tolerance and were free of any major chronic diseases. Patients with other types of DM, cardiovascular diseases, neurologic disorders, peripheral vascular occlusive disease, autoimmune disease, or any other possible causes of peripheral neuropathy were excluded from all the groups. The protocol of this study was approved by the medical ethics committee of Ningbo University.

Aliquots of 2 mL of blood samples from each participant were collected with informed written consent. A total of 6 healthy controls, 6 DM patients, and 6 DPN patients were recruited for the microarray analysis. In the validation stage, another 26 healthy controls, 26 DM patients, and 26 DPN patients were included. An overview of the clinical and demographic characteristics of the participants can be found in our previous study.

Total RNA was extracted from peripheral blood mononuclear cells (PBMCs) by using a TaKaRa MiniBEST Universal RNA Extraction Kit according to the manufacturer’s instructions. The RNA quality was determined using a NanoDrop ND-2000 spectrophotometer on the basis of an optical density (OD) 260/280 ratio ≥ 1.8 and the OD260/230 ratio ≥ 2.0. The RNA integrity was measured using an Agilent 2100 Bioanalyzer. The intensity of the 18S and 28S rRNA bands was examined on a 1% formaldehyde denaturing gel electrophoresis. RNA samples with an RNA integrity number (RIN) of ≥ 6.0 and 28S/18S > 1.5 were subjected to the microarray analysis.

The microarray analysis, including RNA amplification, probe labeling, hybridization, and data extraction was performed by BGI Inc. (Shenzhen, China). Briefly, aliquots (100 ng) of total RNA were amplified and transcribed into fluorescent Cy5-labeled antisense RNAs (aRNAs) using the OneArray Amino Allyl aRNA Amplification Kit (Phalanx Biotech, San Diego, CA) according to the manufacturer’s instructions. Subsequently, the Cy5-labeled aRNAs were fragmented and hybridized to the Human LncRNA OneArray Plus (Phalanx Biotech, San Diego, CA). The arrays were scanned by an Agilent G2505C Microarray Scanner (Agilent Technologies, Wilmington, DE). The fluorescence intensities of each spot were analyzed with Feature Extraction software (Agilent Technologies, Wilmington, DE). Overall, 18 microarray chips were analyzed in this study.

The microarray data were normalized using the R/Bioconductor Limma package [17]. An empirical Bayes model was used to compare the differentially expressed lncRNAs between the groups, and the criteria were set for FDR < 0.05 and fold change ≥ 2 or ≤ 0.5.

The coding-non-coding gene coexpression network (CNC network) was constructed based on the correlation analysis between the differentially expressed lncRNAs and the mRNAs. LncRNAs and mRNAs with Pearson correlation coefficients not less than 0.98 were selected to draw the network using the Cytoscape software. In the network analysis, we used our previous mRNA microarray data from same RNA samples with the same protocol.

To investigate the potential role of the lncRNAs–coexpressed mRNAs, the functional classification and biological pathway interpretation of these mRNAs were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 [18]. The pathways that were significantly represented were identified based on having at least 3 differentially expressed lncRNAs and a hypergeometric test with P < 0.05.

Six lncRNAs were further validated by qRT-PCR in another 78 independent samples from CN, DM, and DPN groups. Total RNA was extracted as described above, and double-stranded cDNA was synthesized using a PrimeScript RT Reagent Kit (TaKaRa Biotechnology, Dalian, China) according to the manufacturer’s instructions. Subsequently, the qRT-PCR was performed using the LightCycler 480 SYBR Green I Master mix (Roche, Mannheim, Germany). The 2-ΔΔ CT method was used to quantify the relative expression of each lncRNA, using GAPDH as an internal control. All the experiments were conducted in triplicate and were repeated three times. Differences in the lncRNA expressions among the groups were evaluated with a oneway analysis of variance (ANOVA) using SPSS 18.0 software. Significance was considered to be P < 0.05.

The lncRNA expression data from this study are available in Gene Expression Omnibus (GEO, accession number GSE95849). Cluster assays of the 1000 randomly selected lncRNAs indicated that the expression pattern of the DPN group was quite different than those of the CN and DM groups (Fig. 1). Volcano Plots were constructed using the fold-change values and P-values and thus allowed the visualization of the differential expression between the groups (Fig. 2). Compared with the healthy individuals, patients with DM and DPN exhibited 319 and 256 up-regulated lncRNAs, respectively, whereas 127 and 1071 lncRNAs were expressed at a lower level in the DM and DPN patients compared with the healthy controls (Table 1). In the DPN patients, compared with the DM patients, 1003 lncRNAs were differentially expressed, including 64 up-regulated and 949 down-regulated lncRNAs (Table 1). The Venn diagram further indicated that the DM and DPN groups shared 127 up-regulated and 112 down-regulated genes (Fig. 3).

Although an increasing number of studies report the functional significance of lncRNAs, the biological roles of most lncRNAs remain unknown. The pathway analysis of the lncRNAs-coexpressed mRNAs may reveal the potential functions of these lncRNAs. KEGG pathway analysis indicated that, relative to the CN group, the DPN and DM groups shared three pathways (Table 2), including synaptic vesicle cycle (hsa04721), leishmaniasis (hsa05140), and measles (hsa05162). In contrast, for the differentially expressed lncRNAs-coexpressed mRNAs between the DPN and DM groups, the most enriched pathway was the MAPK signaling pathway (hsa04010). Interestingly, our previous mRNA profiling also indicated that down-regulation of the neurotrophin-MAPK signaling pathway might be the major mechanism of DPN pathogenesis.

The lncRNA/mRNA coexpression network may provide potential interplay between mRNAs and lncRNAs, and therefore, the differentially expressed lncRNAs and mRNAs between the DPN and DM groups were used to draw the coexpression network using the Cytoscape program. Since the differentially expressed lncRNAs-co-expressed mRNAs between DPN and DM groups were enriched in the MAPK signaling pathway, we focused on how the lncRNAs regulate genes in the MAPK signaling pathway (Fig. 4). The network indicated that BDNF and TRAF2 correlate with 6 lncRNAs (CCNT2-AS1, RP1-249H1.2, CTD-3239E11.2, RP11-51B23.3, STAM-AS1, and LINC00629), implying that the inter-regulation of lncRNAs and mRNAs in the neurotrophin-MAPK signaling pathway occurs in the pathogenesis of DPN.

To verify the results of the microarray analysis, six previously reported diabetes-associated lncRNAs were further validated in an additional 78 participants by qRT-PCR. These lncRNAs were MALAT1 (metastasis associated lung adenocarcinoma transcript 1), H19, PVT1 (plasmacytoma variant translocation 1), MEG3 (maternally expressed gene 3), MIAT (myocardial infarction-associated transcript), and MIR143HG (the human homolog of mouse E330013P06, also known as cardiac mesoderm enhancer-associated non-coding RNA, CARMN). The results of the qRT-PCR showed similar trends to those observed in the microarray data. In detail, the expression level of MALAT1 increased in both the DM and DPN groups. The expression level of H19, PVT1, and MIR143HG decreased only in the DPN group, while MEG3 decreased only in the DM group. The MIAT expression remained unchanged in all the groups (Fig. 5).

Since lncRNAs are associated with diabetes and its complications, microarray analyses have been used to reveal the potential roles of lncRNAs in diabetes [19], diabetic retinopathy [14], and diabetic nephropathy [20]. However, to the best of our knowledge, this is the first study to use lncRNA expression profiles to elucidate the molecular mechanisms underlying the pathogenesis of DPN. DPN is the most common complication of diabetes, affecting up to 50% of DM patients, and it contributes significantly to pain, loss of sensation, numbness, injury, and lower extremity amputation [2]. Because complex pathways are implicated in the pathophysiology of DPN, there are still no specific treatments and no means of predicting or preventing the onset or profession of DPN [1]. In our recent transcriptome analysis, we identified that the differentially expressed genes (DEGs) between the DPN and DM groups were significantly enriched in the “MAPK signaling pathway,” the “NODlike receptor signaling pathway,” and the “neurotrophin signaling pathway.” A stratification analysis further indicated that the down-regulation of the neurotrophin-MAPK signaling pathway might be the specific mechanism of DPN. In the present study, the KEGG pathway analysis of the differentially expressed lncRNAs-coexpressed mRNAs between the DPN and DM groups was also significantly enriched in the MAPK signaling pathway. Moreover, brain-derived neurotrophic factor (BDNF), an essential member of the neurotrophin family, plays important roles in regulating the growth of neurons and influences synaptic efficiency and plasticity [21]. Therefore, deficiency in BDNF has been reported in the peripheral nerves of DPN patients [22]. The current coexpression network further indicates the potential interaction of 6 lncRNAs (CCNT2-AS1, RP1-249H1.2, CTD-3239E11.2, RP11-51B23.3, STAM-AS1, and LINC00629) with BDNF. However, the lncRNA BDNF-AS (BDNF antisense RNA) remained unchanged in the present study. Altogether, these findings strongly suggested that these lncRNAs may participate in the pathogenesis of DPN through regulation of the neurotrophin-MAPK signaling pathway.

Accumulating evidence highlights a growing list of lncRNAs that control glucose homeostasis and diabetic pathologies and complications. However, information on lncRNAs in DPN is still very limited. Thus, we selected 6 previously reported diabetes-associated lncRNAs (MALAT1, H19, PVT1, MEG3, MIAT, and MIR143HG) for further validation in another dataset.

MALAT1, also known as NEAT2 (noncoding nuclear-enriched abundant transcript 2), is a lncRNA consisting of more than 8700 nt located on chromosome 11q13. It is one of the first lncRNAs that was demonstrated to be associated with non-small cell lung cancer [23]. Subsequently, a growing number of evidence indicates that MALAT1 is also closely related to various human diseases, including diabetes and diabetic complications [24]. The expression of MALAT1 is significantly up-regulated in the retinas of diabetic mice, high-glucose-treated retinal endothelial cells, and in the fibrovascular membranes of diabetic patients [14, 25]. Moreover, MALAT1 also affects other cells and tissue types under high glucose circumstances. For instance, high-glucose-treated human umbilical vein endothelial cells and renal tissues from diabetic animals have increased MALAT1 expression along with increased levels of serum amyloid antigen 3 and inflammatory factors [26]. In the present study, the expression level of MALAT1 increased significantly in both the DM and DPN groups. Considering it is involved in multiple types of physiological processes and is associated with various diseases [24], MALAT1 cannot be used as a specific biomarker for DM or DPN.

H19 is transcribed from a conserved imprinted gene cluster that also contains the nearby IGF2 gene, encoding insulin-like growth factor 2 [27]. H19 is highly expressed during fetal life and is strongly down-regulated after birth, except for a persistent expression in adult skeletal muscle and heart. It is highly expressed in muscle cells and is reduced approximately 5-fold in the muscle of both human type 2 DM patients and insulin-resistant mice [28]. Acting as a molecular sponge, the decreased expression of H19 leads to the increased bioavailability of miR-let-7, causing diminished expression of miR-let-7 targets, which results in impaired insulin signaling and glucose metabolism [29]. Moreover, H19 acts as the precursor of miR-675, and the latter, through the down-regulation of its target gene VDAC1 (voltage-dependent anion channel 1), regulates cardiomyocyte apoptosis in the pathogenesis of diabetic cardiomyopathy [30]. Unlike previous studies, the present microarray analysis was conducted in PBMCs, and the expression level of H19 decreased significantly only in the DPN group. However, whether it is a reliable biomarker, and what the underlying mechanism is, needs further investigation.

PVT1 is located in the well-known cancer-related region 8q24. It was first discovered as an activator of MYC in murine plasmacytoma variant translocations [31]. Follow-up studies have brought PVT1 under the spotlight, suggesting interesting models of functioning, including participating in DNA rearrangements, encoding microRNAs, and interacting with MYC [32]. In addition to cancers, variations of PVT1 are associated with end-stage renal disease attributed to both type 1 and type 2 DM [33, 34]. PVT1 is expressed in different types of cells in the kidney, and its expression is increased by glucose in mesangial cells [33]. PVT1 is a host gene of miR-1207-5p in renal cells, which acts as a mediator of extracellular matrix accumulation in the kidney, implicating the pathological processes underlying the development of diabetic nephropathy [35]. However, in the present study, the expression level of PVT1 remained unchanged in the DM group but was significantly decreased in the DPN group. Although, the mechanisms of PVT1 in cancers are quite clear, more studies are needed to decipher its roles in diabetes.

MEG3 is a non-coding transcript belonging to the imprinted DLK1-MEG3 region located at chromosome 14q32. 3 in humans. A large-scale genome-wide association study identified a SNP (rs941576) in the DLK1-MEG3 region that alters susceptibility to type 1 DM [36]. Epigenetic studies further revealed that the regulation of the DLK1-MEG3 miRNA cluster is altered in human type 2 diabetic islets [37]. MEG3 levels are down-regulated in type 1 and type 2 diabetic mice, and the inactivation of MEG3 in normal mice causes decreased insulin synthesis and secretion [38]. Furthermore, MEG3 is down-regulated in diabetic retinas and in retinal endothelial cells upon high glucose stress, which is mediated by the activation of PI3k/Akt signaling [13]. We also found a down-regulation of MEG3 in the DM group, and this evidence suggests potential roles for MEG3 in the pathological processes of diabetes.

A recent genome-wide lncRNA expression study using RNA-seq identified that E330013P06 (the human homolog is MIR143HG) was up-regulated in macrophages from db/db and diet-induced insulin resistant type-2 diabetic mice, as well as in monocytes from type-2 diabetes patients [39]. Interestingly, a further functional study revealed that it modulates the expression of macrophage inflammatory genes and foam cell formation in mouse macrophages treated with high glucose, thus implicating lncRNAs in vascular complications [39]. In our previous whole-genome mRNA analysis, dozens of inflammatory genes were aberrantly expressed, and at least two direct immune-related functions, the “B cell receptor signaling pathway” and the “toll-like receptor signaling pathway” were enriched in both the DM and DPN groups. However, unlike the previous RNA-seq analysis, the present study identified that the expression of MIR143HG decreased significantly in DPN, suggesting that the versatility of MIR143HG may exist in different diabetic complications.

In summary, the current study provides the first demonstration that the interplay between lncRNAs and mRNA may be involved in the pathogenesis of DPN, especially the neurotrophin-MAPK signaling pathway. Several lncRNAs were differentially expressed only in the DPN patients, and thus more studies are required to decide whether any of these could be used as reliable biomarkers and delineate the underlying mechanisms.

This work was supported by the Applied Research Project on Nonprofit Technology of Zhejiang Province (2015C33309), the National Natural Science Foundation of China (81541039, 81402747), the Natural Science Foundation of Ningbo (2016A610085, 2017A610199), and the K.C. Wong Magna Fund in Ningbo University.

The authors declare there are no conflicts of interest.

L. Luo and L.-D. Ji contributed equally to this work.