INTRODUCTION
In China, the incidence of gestational diabetes mellitus (GDM) has increased from approximately 5% to more than 16% since the implementation of diagnostic criteria for GDM in December 2011 [1]. Compared with healthy pregnant women, pregnant women with GDM have a higher incidence of pre-eclampsia and a cumulative incidence of type 2 diabetes within 5–16 years after delivery [2]. In addition, the risk of macrosomia, abortion, hypoglycemia, cardiac dysplasia, obesity, and metabolic syndrome in infants of mothers with GDM is also increased [3]. Due to the high prevalence of GDM and its adverse outcomes, it is necessary to further understand the mechanism of GDM.
In addition to common risk factors, some new environmental factors, such as passive smoking [4] and excessive fruit intake in mid-pregnancy[1], significantly influence the predisposition to GDM. Genomic studies have also revealed that GDM is a multigenic disease. At present, many genes, including protein-coding RNA and noncoding RNA, are involved in the pathogenesis of GDM. MicroRNAs (miRNAs) are a class of small noncoding RNAs with 18–22 nucleotides in length [5]. They modulate gene expression at the posttranscriptional level by interacting with mRNA targets [6]. Serum levels of fibroblast growth factor 23 (FGF23) are upregulated in pregnant women with GDM compared with healthy pregnant women and could be a useful marker of cardiovascular disease in GDM [7]. Serum CysC is significantly and independently associated with insulin resistance, which is helpful in identifying the risk of GDM in Chinese pregnant women [8]. In addition, the combined detection of plasma microRNA-16-5p, -17-5p and -20a-5p is expected to provide early diagnosis for pregnant women with GDM [9]. However, these studies are only based on a single gene. To date, studies providing systematic insight into the composition and expression of the serum mRNA and miRNA transcriptome in GDM are still lacking.
The main objective of this study was to utilize an integrative strategy to construct functional miRNA-mRNA regulatory networks by combining the reverse expression relationships between miRNAs and targets and computational predictions. Our study may provide clues to a better understanding of the mechanism of GDM and help in screening noninvasive molecular markers for the detection of GDM for future biomarkers research.
MATERIAL AND METHODS
Study samples and datasets
A total of three microarray or RNA-seq datasets of plasma samples comparing GDM pregnant women and healthy control pregnant women were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The GSE98043 dataset included two GDM pregnant women and two healthy control pregnant women. The platform used was GPL21575 (Agilent-070156 Human MicroRNA v21.0 Microarray 046064) [10]. In the study containing GSE98043, GDM was diagnosed according to the International Association of the Diabetes and Pregnancy Study Groups (IADPSG) recommendations for the diagnosis and classification of hyperglycemia in pregnancy. All pregnant women were screened for GDM at 24–28 weeks of gestation. Fasting plasma glucose (FPG) was first tested: if the value was ≥ 5.1 mmol/L, a diagnosis of GDM was made; if the value was ≤ 4.4 mmol/L, the woman was diagnosed as no GDM; if the value was ≥ 4.4 mmol/L but < 5.1 mmol/L, a 75-g oral glucose tolerance test (OGTT) was performed, and a diagnosis of GDM was made when at least one glucose value was elevated: FPG ≥ 5.1 mmol/L; 1-hour OGTT ≥ 10.0 mmol/L; or 2-hour OGTT ≥ 8.5 mmol/L. The GSE19649 dataset includes plasma and placental mRNA microarray data, and the GDM was diagnosed under the National Diabetes Data Group criteria for GDM [11]. The 50-g glucose load test was performed at 24–28 weeks of gestation; if the 1-hour plasma glucose ≥ 7.8 mmol/L, a 3-hour 75-g OGTT was performed. The cut-off values for fasting, 1-hour, 2-hour, and 3-hour plasma glucose levels were 5.8, 10.6, 9.2 and 8.1 mmol/L, respectively. If a single threshold was met or exceeded, gestational impaired glucose tolerance was diagnosed. If two or more thresholds were met or exceeded, GDM was diagnosed. Since our study was designed to explore noninvasive plasma markers, we used only plasma mRNA data, which included the pooled blood tissues from 5 GDM women and one healthy control pregnant woman. The platform used was GPL7350 (Aalborg University Illumina Human-6 v2.0 Expression BeadChip). Another independent dataset, GSE92772, was used for validation of the expression of hub genes identified from the GSE19649 dataset [12]. The GSE92772 dataset included blood RNA-seq data from eight pregnant GDM women and eight healthy controls matched for age and BMI. GDM was diagnosed according to the IADPSG criteria. Venous blood was collected during pregnancy at week 24–32. The platform used was GPL16791 (Illumina HiSeq 2500).
Screening of molecular markers
The GEO2R online platform was used to identify the differentially expressed genes (DEmRNAs) in the GSE19649 dataset and the differentially expressed miRNAs (DEmiRNAs) in the GSE98043 dataset between serum samples of GDM pregnant women and healthy control pregnant women. GEO2R is an interactive web tool that allows gene expression analysis of microarray datasets using the limma R package [13]. P < 0.05 and fold change > 1.5 or < 0.5 were set as the cutoff criteria.
Functional enrichment analysis
The target genes of microRNAs were predicted by the mirDIP online platform [14], and the target genes with high scores (top 5%) were selected for further analysis. The computationally predicted target genes and DEmRNAs were represented in a Venn diagram to identify real target genes. KEGG pathway annotation was performed for the real target genes and DEmiRNAs, respectively. All pathway analyses were performed using Funrich 3.1.3 software [15].
Interaction network analysis
Funrich 3.1.3 software was used to systematically analyze and construct the interaction network of DEmRNAs. Then, the hub genes, defined by the genes with the largest number of connecting nodes, were identified.
External hub gene verification
Another dataset (GSE92772) was used as a testing set to verify the differential expression of hub genes between GDM women and healthy controls. Continuous variables were compared using paired t-tests. A receiver operating characteristic (ROC) curve analysis was generated, and the area under the curve (AUC) was calculated to evaluate the ability of the hub gene to predict GDM.
Gene set enrichment analysis (GSEA)
To explore the biological function of hub genes, GSEA was performed in patients with high and low hub gene expression levels in the GSE92772 dataset. The high and low expression groups were separated by the median gene expression. Annotated gene sets c2.cp.kegg.v5.2.symbols. The GMT pathway database was used as the reference. P value < 0.05 and enrichment score (ES) > 0.25 were set as the cutoff values.
A schematic overview of the proposed approach in this study is shown in Figure 1.
RESULTS
By comparing the microarray data of serum samples of GDM patients and those of healthy pregnant women, 264 DEmiRNAs were identified in the GSE98043 dataset, of which 137 genes were upregulated and 127 genes were downregulated. In addition, a total of 1217 DEmRNAs were identified in the GSE19649 dataset. The top 10 DEmiRNAs and DEmRNAs are listed in Table 1. Hsa-miR-146a-3p ranked first because of its lowest p value and was selected for further analysis.
Table 1. Top 10 DEmiRNAs and DEmRNAs |
||||
Differentially expressed genes |
Probe ID |
Gene symbol |
Log(fold change) |
p value |
DEmiRNAs |
9861 |
hsa-miR-146a-3p |
–12.76 |
0.0002 |
55274 |
hsa-miR-146a-5p |
10.82 |
0.0004 |
|
30074 |
hsa-miR-6831-3p |
9.97 |
0.0008 |
|
12711 |
hsa-miR-1285-5p |
9.99 |
0.0008 |
|
32345 |
hsa-miR-146b-3p |
–5.54 |
0.0008 |
|
31670 |
hsa-miR-1911-5p |
6.04 |
0.0009 |
|
17343 |
hsa-miR-4726-5p |
8.89 |
0.0012 |
|
52744 |
hsa-miR-4690-5p |
8.23 |
0.0012 |
|
29799 |
hsa-miR-548w |
8.24 |
0.0012 |
|
56049 |
hsa-miR-4756-5p |
–14.99 |
0.0012 |
|
DEmRNAs |
ILMN_4582 |
KRT2B |
–2.41 |
0.000306 |
ILMN_31624 |
TRPC2 |
–2.45 |
0.000751 |
|
ILMN_25067 |
NR1H4 |
–2.46 |
0.000964 |
|
ILMN_28015 |
LGALS13 |
–2.47 |
0.000993 |
|
ILMN_24642 |
DENND1B |
–2.47 |
0.001123 |
|
ILMN_17357 |
ZNF606 |
–2.48 |
0.001284 |
|
ILMN_26171 |
WDR64 |
–2.48 |
0.001289 |
|
ILMN_11089 |
AGRP |
–2.49 |
0.001331 |
|
ILMN_1266 |
DSC1 |
–2.51 |
0.001709 |
|
ILMN_7518 |
DPYSL4 |
–2.51 |
0.001783 |
The target genes of hsa-miR-146a-3p were identified using two independent and complementary types of information: computational target predictions and expression relationships. For computational target predictions, a total of 1366 mRNAs were predicted to be the target genes of hsa-miR-146a-3p in the mirDIP online platform. These 1366 target genes were further mapped to the 1217 DEmRNAs screened from the GSE19649 dataset. Finally, a total of 47 target genes were shared between the computational target predictions and DEmRNAs and were identified as real target genes of hsa-miR-146a-3p (Fig. 2).
Enrichment analysis of the signaling pathways of DEmiRNAs indicated that GDM-related miRNAs were mainly enriched in the glypican pathway (p < 0.001), proteoglycan syndecan-mediated signaling events (p < 0.001), and syndecan-1-mediated signaling events (p < 0.001) (Fig. 3, Tab. 2). In addition, analysis of the signaling pathways of the real target genes of hsa-miR-146a-3p found that the glypican pathway was also one of the pathways regulated by hsa-miR-146a-3p (p < 0.001). In addition, the DEmRNAs TRAF6, MEF2A, TNFAIP3, SKIL, FLT1 and IRS2 were also targeted by hsa-miR-146a-3p in the present study and were enriched in the glypican pathway.
Table 2. Enrichment analysis of signal pathways involved in DEmiRNAs |
|||
Biological pathway |
Fold enrichment |
p value |
Q value |
Glypican pathway |
1.769653 |
6.84E-13 |
1.99E-09 |
Proteoglycan syndecan-mediated signaling events |
1.760422 |
1.05E-12 |
1.99E-09 |
Syndecan-1-mediated signaling events |
1.763494 |
2.21E-12 |
1.99E-09 |
VEGF and VEGFR signaling network |
1.758073 |
2.81E-12 |
1.99E-09 |
IFN-gamma pathway |
1.757313 |
3.46E-12 |
1.99E-09 |
IFN-gamma — interferon-gamma; VEGF — vascular endothelial growth factor; VEGFR — vascular endothelial growth factor receptor |
To explore the hub DEmRNAs involved in GDM, we further performed an interaction network analysis of DEmRNAs screened from the GSE19649 dataset (Fig. 4). The results indicated that TRAF6, with a total of 16 connecting nodes, was the top hub gene involved in the pathophysiological process of GDM. Other hub genes included CASP8 (number of connecting nodes: 6), PTPN6 (number of connecting nodes: 6), and CHD3 (number of connecting nodes: 6). Next, these four hub genes were selected for independent external validation using the GSE19649 dataset. The expression of TRAF6, CASP8 and CHD3 in 8 pairs of GDM blood samples was confirmed to be higher than that in healthy pregnant women blood samples. However, the expression of PTPN6 in the blood samples of GDM patients is like that of healthy pregnant women blood samples (Fig. 5). The AUC for predicting GDM was 0.813, 0.828, 0.813, and 0.703 for CHD3, PTPN6, TRAF6, and CASP8, respectively (Fig. 6, Tab. 3). Furthermore, TRAF6 was involved in the following three core pathways of GDM: glypican pathway, proteoglycan syndecan-mediated signaling events, and syndecan-1-mediated signaling events.
Table 3. ROC analysis of predicting abilities of hub genes |
||
Gene |
AUC |
95% CI |
CHD3 |
0.813 |
0.592–1.000 |
PTPN6 |
0.828 |
0.621–1.000 |
TRAF6 |
0.813 |
0.600–1.000 |
CASP8 |
0.703 |
0.427–0.979 |
AUC — area under the curve; CI — confidence interval |
To determine the potential mechanism of TRAF6 in GDM, GSEA was conducted (Tab. 4). A total of 9 hallmark gene sets of metabolism processes were enriched, including cysteine and methionine metabolism, pyruvate metabolism, sphingolipid metabolism, beta alanine metabolism, propanoate metabolism, galactose metabolism, starch and sucrose metabolism, butanoate metabolism, and arginine and proline metabolism processes.
Table 4. The gene sets that were significantly associated with TRAF6 by Gene set enrichment analysis (GSEA) |
|||
Name |
Size |
ES |
FWER p value |
KEGG_CYSTEINE_AND_METHIONINE_METABOLISM |
23 |
0.465224 |
0 |
KEGG_PYRUVATE_METABOLISM |
26 |
0.348111 |
0 |
KEGG_SPHINGOLIPID_METABOLISM |
24 |
0.339434 |
0 |
KEGG_BETA_ALANINE_METABOLISM |
15 |
0.332718 |
0 |
KEGG_PROPANOATE_METABOLISM |
26 |
0.331809 |
0 |
KEGG_GALACTOSE_METABOLISM |
18 |
0.324297 |
0 |
KEGG_STARCH_AND_SUCROSE_METABOLISM |
25 |
0.299221 |
0 |
KEGG_BUTANOATE_METABOLISM |
18 |
0.29276 |
0 |
KEGG_ARGININE_AND_PROLINE_METABOLISM |
25 |
0.276111 |
0 |
ES — enrichment score; FWER — family-wise rror rate |
DISCUSSION
The results of the present study provide insights into the molecular mechanisms that underlie GDM. miRNAs are involved in many biological processes through posttranscriptional regulation [16]. They regulate gene expression and transcription by complementary binding to the 3’ end of target genes [16]. However, to date, little is known about the regulatory process between miRNAs and mRNAs in GDM. Thus, we performed an integrative analysis to construct functional miRNA-gene regulatory networks by combining the expression relationships between miRNAs and targets and computational predictions. First, enrichment analysis of the involved signaling pathways was carried out based on DEmiRNAs. The metabolic pathways glypican pathway, proteoglycan syndecan-mediated signaling events, and syndecan-1-mediated signaling events, were the core pathways regulated by miRNAs in GDM. Glypicans belong to the family of membrane-bound heparan sulfate proteoglycans that are linked to the cell surface and are involved in the regulation of growth factor activity [17]. Proteoglycan families are involved in a wide range of diseases, including obesity and diabetes [18]. Impaired glucose metabolism in GDM may alter the expression of proteoglycans, which may impair the biological functions of the placenta [19]. A previous study demonstrated that the expression of a member of the proteoglycan family, endocan, was increased in the human placenta from obese women with GDM and in response to proinflammatory stimuli [18]. In addition, the metabolic effect of high glucose and high osmotic pressure of GDM patients may contribute to the increased proteoglycan perlecan expression in diabetic placentas [20]. Leelalertlauw et al. [21] found that serum glypican four levels were increased with increasing degrees of obesity. Furthermore, Ussar et al. [22] demonstrated that circulating glypican-4 levels correlate with BMI and insulin sensitivity in humans. Moreover, glypican-4 interacts with the insulin receptor, enhances insulin receptor signaling, and increases adipocyte differentiation. These results suggest that abnormal expression of proteoglycan family genes in GDM is helpful for understanding the pathogenesis of GDM.
In the present study, hsa-miR-146a-3p, ranked first in the altered DEmiRNAs of GDM because of its lowest p value, was downregulated with a log (fold change) of –12.76 in the blood samples of GDM patients. Changes in miR-146a in diabetic patients and animals were reported in a previous study. Liu et al. [23] demonstrated that treatment of diabetic mice with miR-146a mimics robustly reduced diabetic peripheral neuropathy. Then, the target genes of hsa-miR-146a-3p were identified based on the parallel DEmRNAs in the present study. The results indicated that TRAF6 was not only identified as a target of hsa-miR-146a-3p but also a hub gene involved in the pathophysiological process of GDM. TRAF6 has been demonstrated to be a target gene of miR-146a at the experimental level in previous immune studies [24]. Regarding diabetes, Kamali et al. found that miR-146a antagomir significantly increased TRAF6 mRNA levels in human umbilical vein endothelial cells under hyperglycemic conditions. In addition, previous research revealed that culture of aortic endothelial cells in high glucose medium significantly increased TRAF6 expression in a time-dependent manner. In addition, TRAF6 mediated high glucose-induced endothelial dysfunction via NF-KB and AP-1-dependent signaling [25]. To date, the role of miR-146a and its target gene TRAF6 in the pathogenesis of GDM has never been explored. In the present study, GSEA revealed that several metabolic processes were enriched in the mechanism of TRAF6 in GDM. Moreover, TRAF6 was involved in all three core pathways regarding glypican metabolism in GDM based on pathway analysis of DEmiRNAs. Furthermore, we confirmed the altered expression of TRAF6 in eight pairs of GDM blood samples to be higher than that in healthy pregnant blood samples in another independent patient group. Thus, serum TRAF6 may serve as a potential novel noninvasive biomarker for GDM. However, further research with sample collection during early pregnancy (before GDM was diagnosed) is needed to evaluate their early predictive value.
Other core genes, such as CASP8 and CHD3, were confirmed to be higher in eight pairs of GDM blood samples than in healthy pregnant women blood samples. At present, there is no research on the relationship between CASP8 or CHD3 and GDM. Wu et al. [25] found that maternal type 2 diabetes mellitus causes heart defects in the developing embryo manifested by excessive apoptosis in heart cells, with increasing apoptosis markers of cleaved caspase 3 and 8. Whether CASP8 can be used as a marker of fetal heart defects in pregnant women with GDM needs to be further studied.
Besides the primary objective of the study to construct miRNA-mRNA regulatory networks, this is also preclinical biomarkers research and certain limitations must be considered. First, although the serum markers screened from the GSE19649 and GSE98043 datasets were independently externally validated using the GSE19649 dataset, its utility for GDM diagnosis warrants further evaluation with sample collection during early pregnancy (before GDM was diagnosed), and the analytical methods must be specified in the clinical standard protocol before clinical application. Moreover, the comparative effectiveness and cost-effectiveness of novel markers versus traditional tests have not yet been established. Secondly, the sample size in this study was relatively small. However, the clinical information, such as age and BMI, of GDM and healthy pregnant women were matched in the GSE19649 dataset, which improved the test effectiveness. Thus, our findings may provide an important clue for future biomarkers research.
In summary, by integrated miRNA-mRNA expression profile analysis, we found that functions of GDM-related miRNAs were mainly enriched in the glypican pathway, proteoglycan syndecan-mediated signaling events, and syndecan-1-mediated signaling events. miR-146a-3p/TRAF6 might play a central role in the pathogenesis of GDM through the above pathways.
Article information and declarations
Data availability statement
All data generated or analyzed during this study are included in this published article. The data of this study are available upon request.
Authors’ contributions
C. M. and Y. J. designed experiments; C. M. carried out experiments, analyzed data and experimental results. C. M. and Y. J. wrote the manuscript.
Funding
C. M. was supported by the Startup Fund for Scientific Research of Fujian Medical University (No. 2019QH1137) and The Natural Science Foundation of Fujian Province of China (No. 2020J01338). Y. J. was supported by the Guide Fund for the Development of Local Science and Technology from the Central Government (2020L3019) and Health research project of Department of Finance (Fujian finance refers to [2019] No. 827) (2020Y183).
Acknowledgements
Not applicable.
Conflict of interest
All authors declare no conflict of interest.
Supplementary material
None.