Construction of a long non‑coding RNA-mediated competitive endogenous RNA network reveals global patterns and regulatory markers in gestational diabetes.

Gestational diabetes mellitus (GDM) is a common disease affecting pregnant women. Recent studies have suggested that competing endogenous RNAs (ceRNAs), which compete with long non‑coding RNAs (lncRNAs) for microRNA (miRNA or miR) binding and indirectly regulate miRNA targets through competing interactions, play a critical role in disease. In this study, we present a computationally integrated approach with which to construct a lncRNA‑mediated ceRNA network (LCEN) in GDM by integrating RNA interactions and expression data. lncRNAs exhibited specific features and played critical roles in GDM‑associated LCEN. The construction of a global functional score profile revealed that ceRNAs had a high activity in GDM. We extracted several ceRNA modules and demonstrated that these modules had increased close interactions. We further discovered that these ceRNA modules may be utilized as specific and effective circulating biomarkers for GDM. Finally, functional analyses demonstrated that the GDM‑associated ceRNAs participated in the regulation of irisin and the thyroid hormone signaling pathway. It was suggested that there were close associations between the thyroid hormone and GDM. Collectively, ceRNAs may accelerate biomarker discovery and therapeutic development in GDM.

Introduction

Gestational diabetes mellitus (GDM) is defined as glucose intolerance with onset or first recognition during pregnancy and is one of the most common complications during pregnancy (1). Metabolic and immunological changes can occur during pregnancy and frequently present clinically as increased insulin resistance and immune tolerance against the fetus and placenta (2). The risk of adverse pregnancy outcomes is increased in patients with GDM and affects the mother and the child; common complications include ischemic heart disease, hypertension, predisposition to obesity, metabolic syndrome and type 2 diabetes mellitus (T2DM) following pregnancy (3). Risk factors include pre-pregnancy weight gain and obesity, a family history of diabetes, an advanced maternal age, a poor diet and low physical activity (4,5). In addition, studies have reported that low thyroid hormone levels during early pregnancy are associated with an increased risk of GDM and may affect pregnancy outcome and the intellectual level of the infant (6,7). However, the etiology of and mechanisms responsible for GDM are unknown for the majority of patients and thus, the identification of novel signatures or biomarkers that enhance clinical behavior in the treatment of GDM are essential.

Long non-coding RNAs (lncRNAs) do not code for proteins and are pervasive across the genome (8). The dysregulation of lncRNA expression is associated with various human diseases (9-11). A previous study suggested metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) lncRNA expression may be a novel biomarker to predict GDM (12). Shi et al (13) performed a microarray expression profile analysis to reveal that lncRNAs were differentially expressed in the blood of the umbilical cord from patients with GDM and may play a role in macrosomia development. Recently, an increasing number of studies have suggested that lncRNAs participate in competing endogenous RNA (ceRNA) regulatory processes and in communicating with other RNAs, such as microRNAs (miRNAs or miRs), coding genes and circular RNAs (circRNAs) (14,15). lncRNAs compete with miRNA target genes for binding by sharing common miRNA-binding sites, thus attenuating miRNA-associated target repression. UICLM lncRNA has been shown to promote liver metastasis in colorectal cancer by acting as a ceRNA for miR-215, regulating zinc finger E-box-binding homeobox 2 (ZEB2) expression (16). MT1JP lncRNA has been shown to function as a ceRNA in regulating FBXW7 by competitively binding to miR-92a-3p in gastric cancer (17). In addition, studies have utilized high-throughput expression profiles and interaction data to construct a global ceRNA network in diseases (18,19). These data highlight the important role of lncRNA in interacting with ceRNAs and suggest that the integration of expression profiles and network analysis contributes to the identification of problematic lncRNAs and to the elucidation of the mechanisms of disease. However, a lncRNA-mediated ceRNA network in GDM has not yet been constructed and analyzed, at least to the best of our knowledge.

In this study, we used a comprehensive computational approach to construct a lncRNA-mediated ceRNA network (LCEN) in GDM by integrating the expression profiles of mRNAs, miRNAs and lncRNAs from patients with GDM and experimentally verifying the interactions. We observed that the lncRNAs exhibited specific topological features in the LCEN, consistent with a regulatory association with coding mRNAs in GDM. The LCEN in GDM presents specific and highly competing activity profiles. We further identified a core GDM-associated subnetwork and 3 modules to characterize the properties of LCENs. ceRNA expression in these modules may be able to distinguish between patients that are normal glucose-tolerant (NGT) and patients with GDM. Functional analysis revealed that GDM-associated ceRNAs participated in glycolytic and hormone metabolic processes. Additionally, we identified the thyroid hormone pathway to be associated with ceRNAs in GDM. On the whole, these results suggest that GDM-associated LCENs may provide new insight into the mechanisms of GDM and may aid in the discovery of novel molecular biomarkers and GDM therapeutic strategies.

Materials and methods

Collection of high-throughput data for lncRNAs, miRNAs and genes

lncRNA, miRNA and gene expression profiles for GDM were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo). The study collected data for 8 NGT patients and 8 patients with GDM, which were matched by body mass index and age. Paired data were adjusted for mid-pregnancy weight gain and by pregnancy week (accession no. GSE92772; unpublished data). The RNA sequencing of whole blood cells of these samples were produced by Illumina HiSeq 2500.

Experimental validation of miRNA targets and GDM-associated genes

Gene-miRNA association data were obtained from a public database (miRTarBase 7.0) that contained >360,000 miRNA-target interactions (MTIs) (20). Only MTIs supported by strong experimental evidence were extracted for the current study. lncRNA-miRNA association data was obtained (RAID 2.0) containing experimentally and computationally predicted RNA-RNA and RNA-protein interactions. We further extracted lncRNA-miRNA association data supported by strong experimental evidence (21). GDM-associated gene data were downloaded from DisGeNET, a discovery platform containing one of the largest publicly available collections of genes and variants associated with human diseases (22).

Identification of candidate ceRNA interactions

A hypergeo-metric test was performed to identify candidate competing mRNA-lncRNA interactions. We evaluated the significance of shared miRNAs between lncRNAs and mRNAs. The human genome contains a total number of K miRNAs, with L and M representing the number of miRNAs associated with the current lncRNA or mRNA, respectively. The number of common miRNAs shared by lncRNA and mRNA was defined as ‘i’. P-values were calculated using the hypergeometric test and by evaluating the enrichment significance for competing function as follows:

A false discovery rate (FDR) correction was applied. FDR<0.01 was the threshold to select candidate competing mRNA-lncRNA pairs.

Identification of GDM-associated functional ceRNAs by integrating expression and candidate ceRNA interactions

Pearson correlation coefficients (PCC) were used to identify GDM-specific mRNA-lncRNA pairs based on the expression values of competing mRNA-lncRNA pairs. We filtered GDM-associated functional ceRNAs from expression datasets based on following rules: PCC (lncRNA, gene) >0.5 and P<0.05; PCC (lncRNA, miRNA) <−0.5 and P<0.05; and PCC (gene, miRNA) <−0.5 and P<0.05, with PCC representing interactions based on expression values. A total of 77 GDM-associated functional ceRNAs, comprising 17 lncRNAs, 14 miRNAs and 40 coding genes, were retained for further analysis. A GDM-associated ceRNA network was constructed using Cytoscape 3.0 (http://www.cytoscape.org/).

Functional score of GDM-associated functional ceRNAs

The correlation of mRNA and miRNA expression is an effective statistically method for distinguishing between direct and indirect interactions (23). The functional score was applied to determine the strength of competition in GDM-associated functional ceRNAs and was defined as (|PCC (lncRNA, gene)|+|PCC (lncRNA, miRNA)|+|PCC (gene, miRNA)|)/3. A higher functional score indicated a greater competition between lncRNA and gene for miRNA binding.

Dissecting topological features for LCEN in GDM

We performed topological analyses, including degree and topological coefficients for all the nodes in LCEN using Cytoscape 3.0 (http://www.cytoscape.org/).

Identifying core modules from subnetworks of LCEN in GDM

We extracted a subnetwork form LCEN in GDM following the verification GDM-associated genes. All the network modules were identified using ClusterONE with the default parameters based on the subnetwork of LCEN constructed by GDM-associated genes (http://apps.cytoscape.org/apps/ClusterONE). ClusterONE is a tool that clusters a given network based on topology to identify densely connected regions. A total of 14 modules were identified and the 3 modules with the highest numbers of nodes were extracted.

Classification power of the core modules in GDM

For core modules, the expression data was used to classify 16 samples (NGT, n=8 and GDM, n=8) applying a consensus clustering method (24). This process is performed by ConsensusClusterPlus package in R (https://www.r-project.org/). We defined the optimal category number as the smallest increase in the area under the cumulative distribution function (CDF) curve. Chi-square tests were used to evaluate whether the disease and control samples can be classified using this method. A value of P<0.05 was considered to indicate a statistically significant difference. Gene, miRNA and lncRNA expression were used to perform this analysis and the respective expression was integrated.

Functional enrichment analysis

With the online Enrichr tool and by applying default parameters, a functional enrichment analysis was performed for genes in LCEN (25). We identified enriched GO terms (P<0.01) and KEGG pathways (P<0.05).

Results

Construction and global properties of the GDM-associated ceRNA network

A total of 480,023 candidate ceRNA interactions (gene-miRNA-lncRNA) were identified. An integrated pipeline was used to construct a GDM-associated LCEN based on experimentally verified RNA interactions and expression data (Fig. 1A). The network consisted of 71 nodes (17 lncRNAs, 14 miRNAs and 40 coding genes) and 150 edges. A total of 76 gene-lncRNA interactions, 30 miRNA-lncRNA interactions and 44 gene-miRNA interactions were identified in the GDM-associated LCEN. The GDM-associated LCEN exhibited a scale-free distribution (R2=0.817) and was similar to the majority of biological networks (Fig. 1B). It was further indicated that this network had a similarity to the small-world network (26). In addition, it was discovered that the topological coefficient decreased when the degree increased, which suggested that the GDM-associated LCEN described a hierarchical modularity phenomenon (Fig. 1C) (27). This topological feature has been presented in ceRNA networks of various types of cancer (19). Furthermore, we determined the degrees of coding genes, miRNAs and lncRNAs. A node was considered as a hub, with increased ceRNA interactions representing higher node degrees. In the current study, lncRNA nodes exhibited higher mean degrees compared with coding gene and miRNA nodes (Fig. 1D-F). The mean degree of coding for genes, miRNAs and lncRNAs was 3, 5.54 and 9.40, respectively. These results suggested that although lncRNA do not code for proteins, they exhibited more specific topological properties than mRNAs in the GDM-associated LCEN. The highest degree nodes were determined for phosphatase and tensin homolog (PTEN), hsa-miR-19b-3p and deleted in lymphocytic leukemia 1 (DLEU1).

Figure 1

Construction and global characteristics of the LCEN in GDM. (A) A global LCEN in GDM. Coding genes, lncRNAs and miRNAs are colored blue, red and yellow, respectively. (B) Degree distribution of nodes in the LCEN. (C) Topological coefficient of nodes in the LCEN. The degree distributions of (D) genes, (E) miRNAs and (F) lncRNAs; the horizontal axis represents the degree values and vertical axis represents the number of nodes with a certain degree value. LCEN, lncRNA-mediated ceRNA network; GDM, gestational diabetes mellitus; lncRNA, long non-coding RNA; miRNA, microRNA; ceRNA, competitive endogenous RNA.

Variable competing activity profiles in the GDM-associated ceRNA network

To characterize the compactness of each ceRNA interaction in the GDM-associated LCEN, we evaluated the PCC values for the gene-lncRNA, miRNA-lncRNA and miRNA-gene interactions. Density curves for these interactions were determined and observed to be similar (Fig. 2A). The majority PCC values were concentrated between 0.75-0.80 and the results indicated that interactions were compact. The mean PCC value of gene-lncRNA interactions was highest with 0.83 and the miRNA-gene and miRNA-lncRNA interactions had mean PCC values of −0.8 and −0.79, respectively (Fig. 2B). In addition, a functional score was proposed to evaluate the GDM-associated functional ceRNAs. Based on the functional score, we constructed a competing activity profile for the GDM-associated LCEN (Fig. 2C). We discovered that while certain interactions exhibited lower PCC values, the overall ceRNAs exhibited a high functional score. The result indicated that the functional score may be used to evaluate the properties of ceRNA interactions. We further discovered that the gene expression for VEGFA competed with lncRNA expression of MZF1-AS1 and hsa-miR-378a-3p expression (Fig. 2D). The expression heatmap exhibited differences between GDM and NGT (Fig. 2E). The expression levels of VEGFA, MZF1-AS1 and hsa-miR-378a-3p in GDM were increased compared with NGT (Fig. 2F). We further observed that the difference in lncRNA levels between GDM and NGT was not significant, potentially due to decreased lncRNA expression in blood. The results suggested that gene and miRNA associations with ceRNA may potentially be used to explore the role and function of lncRNA, especially in blood sample.

Figure 2

Construction of a functional score profile for the GDM-associated LCEN. (A) Density distribution curves of gene-lncRNA, gene-miRNA and miRNA-lncRNA. (B) Mean correlation values of gene-lncRNA, gene-miRNA and miRNA-lncRNA interactions. (C) Functional score profile for functional ceRNAs in GDM. (D) An exemplary ceRNA. (E) Heatmap of ceRNA expression in GDM and NGT. (F) Mean expression values in GDM and NGT. GDM, gestational diabetes mellitus; LCEN, lncRNA-mediated ceRNA network; lncRNA, long non-coding RNA; miRNA, microRNA; ceRNA, competitive endogenous RNA; NGT, normal glucose-tolerant.

Characterization of a core GDM-associated subnetwork

We extracted a subnetwork based on verified GDM-associated genes to explore how ceRNAs provide insight into GDM pathogenesis. We observed that 57.5% of the genes in the GDM-associated LCEN were verified as GDM-associated genes (Fig. 3A). We further counted the degree of nodes in the subnetwork (Fig. 3B). DLEU1 lncRNA exhibited the highest degree, which indicated that it may play an essential role in the GDM-associated LCEN. The subnetwork consisted of 50 nodes (23 coding genes, 12 miRNAs and 15 lncRNAs) and 103 edges (Fig. 3C). A module analysis of the subnetwork was performed to further investigate the interactions between various RNA transcripts. Three modules were extracted, and the number of nodes and features were analyzed (Fig. 3D). The first module included 8 nodes (5 coding genes, 1 miRNA and 2 lncRNAs); the nodes formed 10 ceRNAs. The second module included 8 nodes (4 coding genes, 1 miRNA and 3 lncRNAs); the nodes formed 10 ceRNAs. The third module included 7 nodes (4 coding genes, 1 miRNA and 2 lncRNAs); the nodes formed 7 ceRNAs. A small quantity of nodes formed multiple ceRNAs indicating these nodes had strong associations. We further considered the compactness of each ceRNA interaction in modules and discovered the ceRNA interactions in modules were more compact. In the third module, all 7 ceRNAs had high functional scores (Fig. 3E) suggesting that the nodes in this module were close.

Figure 3

Core subnetwork and ceRNA modules in GDM. (A) Percentage of known GDM-associated genes in all genes. (B) Degree of nodes in the subnetwork. (C) A subnetwork extracted from LCEN, which only included known GDM-associated genes and interacting RNAs. (D) Three ceRNA modules extracted from the subnetwork. (E) Co-expression pattern of each ceRNA in the third module; circle size reflects correlation strength; red represents a positive and blue a negative correlation; the functional score is presented in green. ceRNA, competitive endogenous RNA; GDM, gestational diabetes mellitus.

GDM-associated ceRNA distinguish patients with NGT and GDM

To evaluate whether ceRNAs are potential classifiers in GDM, we used gene, miRNA and lncRNA expression values from each module to classify the samples using a consensus clustering method. We further used a Chi-square test to validate the classification power of the modules. No significant differences in gene, miRNA or lncRNA expression between GDM and NGT samples were observed in these 3 modules; however, significance was determined for the integration of first and second module (P=0.03; Fig. 4A). The results indicated that ceRNAs may be an effective classifier for patients with GDM. We then analyzed the classification power to evaluate the integrated expression of 9 coding genes [cullin 5 (CUL5), DNA methyltransferase 1 (DNMT1), endothelial PAS domain protein 1 (EPAS1), nuclear factor of activated T cells 3 (NFATC3), protein phosphatase 2 regulatory subunit B′ epsilon (PPP2R5E), protein kinase AMP-activated catalytic subunit alpha 1 (PRKAA1), PTEN, stromal interaction molecule 1 (STIM1) and thyroglobulin (TG)], 5 lncRNAs (HCG27, LINC00954, MBNL1-AS1, PTENP1 and SND1-IT1) and 2 miRNAs (miR-19b-3p and miR-185-5p). The GDM and NGT samples could be classified into several groups by the above-mentioned genes, lncRNAs and miRNAs. According to the CDF and relative change in the area under the curve (Fig. 4C and D), we determined that the final number of groups was 4. We discovered that the 4 sample groups had respective expression patterns (Fig. 4D). We further determined that a majority of samples was classified accurately, particularly in the second group (C2), where all samples were from GDM (Fig. 4E). Collectively, these results suggested that the integration of expression in ceRNA interactions can distinguish between NGT and GDM samples.

Figure 4

ceRNA modules distinguish NGT and GDM samples. (A) P-values from the Fisher-test classifying the power of miRNA, genes and lncRNA. (B) Cumulative distribution function of the consensus index. (C) Relative changes in the area under the CDF curve by group number. (D) Consensus cluster heatmap of samples. (E) Gene expression heatmap, with sublabels referring to group type classified by the consensus cluster method and the sample type referring to disease status of samples. The names of lncRNAs, genes and miRNAs are represented by blue, red and yellow. ceRNA, competitive endogenous RNA; GDM, gestational diabetes mellitus; NGT, normal glucose-tolerant; miRNA, microRNA; lncRNA, long non-coding RNA.

GDM-associated ceRNA is associated with critical biological functions and the thyroid hormone

We performed GO enrichment analysis based on the genes in the GDM-associated LCEN. The enrichment of these genes in different GO terms was observed (Fig. 5A). Genes were associated with some critical biological functions, including the positive regulation of gene expression and the positive regulation of transcription. Certain genes exhibited metabolism-associated functions, such as the enrichment of the regulation of macromolecular, carbohydrate and ATP processes. In addition, we discovered associations to glycolytic processes. Significantly enriched GO functions were associated with hormones, including intracellular steroid production and hormone metabolism. We further performed pathway enrichment analysis and identified insulin-associated pathways to be relevant in the GDM-associated LCEN (Fig. 5B). The thyroid hormone signaling pathway was identified in the analysis and this hormone may be a potential biochemical marker in predicting GDM and effects on insulin sensitivity in GDM (28,29). As the thyroid hormone status has previously been demonstrated to be associated with GDM during pregnancy (29,30), we further analyzed this signaling pathway and identified 4 enriched genes (Fig. 5C). These genes participated in 8 ceRNA interactions and the pathway was identified to be associated with the gluconeogenesis pathway, which is an essential pathway in diabetes (31). We further focused on irisin, which is an exercise-induced myokine that drives brown-fat-type thermogenesis in murine white fat (32). In this analysis, we observed associations with the MAPK signaling pathway and importance arises from irisin potentially suppressing obesity and associated T2DM via the MAPK signaling pathway (33). A ceRNA module identified in the current study, MAPK1/miR-378/HCG27, may be associated with irisin based to literature known MAPK1 and miR-378 associations with irisin (34,35). Additionally, hyper- and hypothyroidism are associated with the upregulation of serum irisin in male rats (36). These results suggested GDM-associated ceRNAs may participate in the regulation of the thyroid hormone and irisin and influence GDM development.

Figure 5

Functional analysis of genes in the GDM-associated LCEN. (A) Enriched GO terms for genes in the GDM-associated LCEN ranked by −log10(P). (B) KEGG pathway enrichment for genes in the GDM-associated LCEN, ranked by −log10(P). (C) Association of the thyroid hormone signaling pathway with GDM-associated ceRNAs. Some dysregulated ceRNAs are showed in figure. GDM, gestational diabetes mellitus; LCEN, lncRNA-mediated ceRNA network; ceRNA, competitive endogenous RNA.

Discussion

Coding genes, miRNAs and lncRNAs can form complex regulatory clusters, including ceRNAs, which may influence pathogenic mechanisms in GDM. lncRNAs regulate mRNA expression through an indirect post-transcriptional mechanism of competing with miRNA for mRNAs binding sites (14). In various diseases, individual ceRNA crosstalk interactions were verified and global ceRNA networks were constructed (37-39). Lin et al (40) constructed a ceRNA network to reveal the regulatory role of lncRNAs in T2DM. However, Lin et al (40) solely relied on disease-associated genes to construct ceRNA network and expression information was not integrated. The lack of expression information reduced specificity and accuracy in the disease biomarker analysis. In this study, a system was created that integrated experimentally-verified RNA interactions, expression data and verified GDM-associated genes. A global ceRNA network was constructed and GDM-associated ceRNAs were identified. A candidate GDM-associated ceRNA list were provided and some lncRNAs, genes and miRNAs in this list have not been verified by strong evidence. In future studies, researchers could extracted these molecules to construct relations between non-coding RNAs and GDM or identify more accurate biomakers and treatment targets for GDM. Collectively, the current study provides a comprehensive resource for studying the regulation of GDM by non-coding molecules.

Various lncRNAs function as molecular biomarkers of different types of diseases. lncRNA-HEIH levels in serum and exosomes are potential biomarkers for HCV-associated hepatocellular carcinoma (41). Plasma lncRNA GAS5 is a novel biomarker for coronary artery disease (42). Survival analyses further demonstrated that ceRNA network modules are potential prognostic biomarkers in GBM (43). In the current analysis, we found that although single lncRNAs could not perfectly distinguish GDM and NGT (Chi-square test, P>0.05), some combinations among lncRNAs, genes and miRNAs could distinguish GDM and NGT (Chi-square test, P=0.03). The result also verified the advantage of considering the ceRNA as specific and effective biomarkers to distinguish between NGT and GDM. We inferred that lncRNAs could interact with other molecules (including miRNAs and genes) to play their roles in some diseases. Blood-based biomarkers are critical for disease prediction. In this study, it was suggested that certain candidate ceRNAs in the blood may serve as circulating biomarkers in GDM and more experiments focused on whole ceRNA motif in cell lines or animal models should be constructed.

The frequency of thyroid dysfunction in diabetic patients is increased compared with the general population and ≤1/3 of patients with type 1 diabetes develop thyroid dysfunction (44). Thyroid hormones are positively associated with insulin resistance in the early development of T2DM (45). However, little attention has been paid to the diagnosis of thyroid diseases in diabetics, as they are diagnosed in only approximately half of the patients (46). In the current study, we discovered various GDM-associated ceRNAs associated with the thyroid hormone signaling pathway. For example, SLC9A1 is a key upstream gene of thyroid hormone signaling pathway and form two dysregulated ceRNAs with miR-185-5p and lncRNA DLEU1 and LINC00954 in GDM. DLEU1 was associated with BMI-adjusted adiponectin, which is related to diabetes (47). In this study, lncRNA DLEU1 had the highest degree in GDM-associated subnetwork (Fig. 3B). We provide novel insight into how lncRNA DLEU1 influences thyroid hormone synthesis by competing miR-185-5p with SLC9A1. The MAPK signaling pathway is an essential part of the thyroid hormone signaling pathway (34). In this study, we found that lncRNA HCG27 competes for miR-378a-3p with MAPK1 in GDM, and MAPK1 and miR-378 have both been shown to be associated with irisin in previous studies (34,35). Irisin is expressed and produced by human muscle and adipose tissue in association with obesity and insulin resistance (48). Thus, we inferred that there were close relations among irisin, thyroid hormone, obesity and GDM. In summary, novel insight was provided into the study of the associations between thyroid hormone dyscrasia and GDM. Further studies are warranted to focus on investigating an increased number of GDM samples to validate the accuracy and stability of the method presented in the current study.

In conclusion, in the present study, we constructed a GDM-associated LCEN and analyzed its features. A functional score was used to evaluate the activity of each ceRNA in GDM and certain ceRNA modules were extracted analyzed for strict interactions. In addition, these ceRNA modules demonstrated to distinguish between GDM and NGT samples. We further observed a close association between thyroid hormone dyscrasia and GDM. The current analysis provides insight into the identification of novel biomarkers for GDM.