Identification of potential gene markers in gestational diabetes mellitus.

This study aims to investigate underlying mechanisms of gestational diabetes mellitus (GDM). In this work, the GSE70493 dataset from GDM and control samples was acquired from Gene Expression Omnibus (GEO) database. Afterward, differentially expressed genes (DEGs) were screened between GDM and control samples. Subsequently, functional enrichment analysis and protein-protein interaction (PPI) network analysis of these DEGs were carried out. Furthermore, significant sub-modules were identified, and the functional analysis was also performed. Finally, we undertook a quantitative real-time polymerase chain reaction (qRT-PCR) with the purpose of confirming several key genes in GDM development. There were totally 528 up-regulated and 684 down-regulated DEGs between GDM and healthy samples. The functional analyses suggested that the above genes were dramatically enriched in type 1 diabetes mellitus (T1DM) process and immune-related pathways. Moreover, PPI analysis revealed that several members of human leukocyte antigen (HLA) superfamily, including down-regulated HLA-DQA1, HLA-DRB1, HLA-DPA1, and HLA-DQB1 served as hub genes. In addition, six significant sub-clusters were extracted and functional analysis suggested that these four genes in sub-module 1 were also associated with immune and T1DM-related pathways. Finally, they were also confirmed by qRT-PCR array. Besides, the four members of HLA superfamily might be implicated with molecular mechanisms of GDM, contributing to a deeper understanding of GDM development.

INTRODUCTION

Gestational diabetes mellitus (GDM) is regarded as one of health problems threatening pregnant women and manifested with glucose intolerance during pregnancy. , Numerous studies have suggested that patients undergoing GDM tended to develop some complicated complications, and there existed an enhanced risk of type 2 diabetes mellitus (T2DM) in pregnancies suffering from GDM. Existing evidence has shown that the incidence and progression of GDM were caused by multiple factors, including age, ethnic group, and body mass index (BMI). , In recent years, although various strategies containing medical nutrition therapy, exercise, and healthy lifestyle have been applied to improve the clinical outcomes of GDM, the results were unsatisfactory. Besides, although increasing investigators have concentrated on exploring the genetic basis for the initiation and development of GDM, there still existed many difficulties such as GDM screening and heritability estimates.

Obviously, previous research has highlighted that there were alterations in glucose tolerance for pregnancies suffering from GDM and these women generally exhibited a higher risk of T2DM. Encouragingly, in recent few years, the advancement of sequencing techniques considerably promotes genetics research of GDM. Wang et al. carried out a bioinformatics analysis based on a gene expression profile from GDM samples. Consequently, they established a promising screening approach for GDM and multiple members of human leukocyte antigen (HLA) superfamily were associated with GDM. Deng et al. pointed out that antigen processing pathway and immune‐associated genes played key roles in GDM progression according to an integrated analysis of gene expression and methylation data which was generated from visceral omental adipose tissue of numerous Chinese pregnancies. Ding et al. preformed RNA sequencing for placentas of subjects with GDM and control individuals. They found that miR‐138‐5p was differentially expressed in samples with GDM, and it can significantly suppress the proliferation and invasion of trophoblasts, suggesting that miR‐138‐5p might be a potential prognostic biomarker for GDM. Although numerous studies have explored the underlying molecular etiology of GDM, a deeper understanding of GDM development and progression is still needed to promote the establishment of effective therapeutic measures for conquering GDM.

In this study, a microarray dataset (GSE70493) from GDM was retrieved and downloaded from GEO database. Then, we extracted the differentially expressed genes (DEGs) between GDM and control followed by functional annotation and enrichment analyses. Moreover, we conducted a PPI network analysis for DEGs to identify and screen significant sub‐clusters. On this basis, the functional analyses of those DEGs in sub‐modules were performed. Finally, a RT‐PCR was undertaken to confirm several potential gene markers in GDM development, which will contribute to the understanding of GDM pathogenesis.

MATERIALS AND METHODS

Data collection and pre‐processing

A gene expression dataset (GSE70493) provided by Binder et al. was downloaded from the National Center for Biotechnology Information GEO database (http://www.ncbi.nlm.nih.gov/geo/). This dataset included a total of 63 samples (32 placentas samples with GDM and 31 healthy controls). These samples were subjected to [HTA‐2_0] Affymetrix Human Transcriptome Array 2.0 [transcript (gene) version]. Subsequently, the raw sequencing data obtained was pre‐processed by affy package in R 3.3.2 software, mainly including background correction and normalization of the gene expression value.

Screening of DEGs

The gene expression matrix was categorized into two groups (GDM group and control group) and used to identify DEGs. The unpaired t test in Limma was employed to calculate the p value of differential gene expression. Then, we employed Benjamini and Hochberg method to adjust all p values. Notably, in accordance with the cutoff of p value < 0.05, the DEGs in the current study were extracted. Finally, the pheatmap package was applied to construct the heatmap of obtained DEGs.

Gene ontology (GO) and pathway enrichment analyses of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID), a freely online resource, offers extensive functional analysis tool suites. In this study, we performed the GO functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for up‐regulated and down‐regulated DEGs with the use of DAVID. Moreover, the p value < 0.05 and count ≥2 were regarded as the thresholds for significant enrichment analyses.

Protein–protein interaction (PPI) network and sub‐network construction

To further investigate whether there existed correlations among proteins of DEGs, this study performed a PPI analysis of DEGs screened using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database. Briefly, all the DEGs were regarded as the input gene set, the species was set as saccharomyces cerevisiae, and the PPI score of 0.9 was considered to be the cutoff for PPI analysis. Afterward, PPI network was built and visualized by adopting Cytoscape software. Furthermore, the node score in PPI network was assessed by degree centrality representing the topological features of the PPI network. Specifically, the high node score indicated that this protein played a critical role in this network and the top 10 nodes were listed.

Generally, the proteins in the same module tend to exert the same or similar roles, that is, this module generates same biological effect. Therefore, we further extracted the significant functional modules using the culsterONE method of cytoscape plugin. Finally, the functional enrichment analyses of each module were performed with DAVID.

The quantitative real‐time PCR (qRT‐PCR)

The critical genes such as several members of HLA superfamily were examined by a qRT‐PCR array by adopting Maxima™ SYBR Green qPCR mix (Thermo Scientific). The specific primers for HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 are exhibited in Table 1. In the meanwhile, the glyceraldehyde 3‐phosphate dehydrogenase (GAPDH) was considered to be internal gene. In this work, the PCR amplification was conducted on Real‐Time PCR System according to following procedures: initial denaturation at 95°C for x min, 30 cycles at 95°C for 15 s, and 60°C for 60 s. Finally, the amplification products were measured with agarose gel electrophoresis.

TABLE 1

Specific primers for qPCR assay

Gene	Primer sequence
HLA‐DQA1	Sense primer: 5′‐ACCTGGGGAGGAAGGAG‐3′
	Antisense primer: 5′‐CGGTAGAGTTGGATCGTTT‐3′
HLA‐DRB1	Sense primer: 5′‐GGGTGGAGAGGGGTCATAG‐3′
	Antisense primer: 5′‐GCTGGAGAACAGGACAGGA‐3′
HLA‐DPA1	Sense primer: 5′‐GCCTGAGACAACGGAGAC‐3′
	Antisense primer: 5′‐GCCAGAACGCAGAGACTT‐3′
HLA‐DQB1	Sense primer: 5′‐CTTATGCCTGCCCAGAAT‐3′
	Antisense primer: 5′‐GGGGATGAAAGGAGATGAC‐3′
GAPDH	Sense primer: 5′‐GTGGATCAGCAAGCAGGAGT‐3′
	Antisense primer: 5′‐AAAGCCATGCCAATCTCATC‐3′

RESULTS

DEGs identification

The gene expression matrix containing 70,523 probes was obtained after the raw microarray data were standardized (Figure 1A). Totally, there were 1212 DEGs containing 528 up‐regulated genes and 684 down‐regulated genes between GDM and control groups. Furthermore, the cluster analysis indicated that these DEGs could distinctly differentiate GDM and healthy control samples, as displayed in Figure 1B.

FIGURE 1

Data standardization and clustering analysis of differentially expressed genes (DEGs). (A) The boxplot after data standardization: the green represents gestational diabetes mellitus (GDM) samples, and the red indicates the control samples. (B) The heatmap of DEGs: the red color represents the high gene expression level while the green color shows the low gene expression level

Functional analyses of DEGs

As shown in the GO functional annotation analysis, these DEGs were enriched in 105 GO‐biological process (BP) terms containing immune response and antigen processing. Meanwhile, they were clustered in 32 GO‐molecular function (MF) terms such as antigen binding and major histocompatibility complex (MHC) class II receptor activity along with 30 GO‐cellular component (CC) terms mainly focusing on MHC class II protein complex and integral component of endoplasmic reticulum membrane (Table 2). In addition, KEGG analysis demonstrated that these DEGs were primarily responsible for 36 pathways including graft‐versus‐host disease, systemic lupus erythematosus, and type I diabetes mellitus (T1DM; Table 2).

TABLE 2

Functional enrichment analyses of differentially expressed genes

Category	Term	Count	p Value
KEGG_pathway (top 10)	hsa05332:Graft‐versus‐host disease	20	1.73E‐16
	hsa05322:Systemic lupus erythematosus	35	7.22E‐15
	hsa05330:Allograft rejection	19	5.95E‐14
	hsa04940:Type I diabetes mellitus	20	6.30E‐14
	hsa04612:Antigen processing and presentation	24	3.59E‐12
	hsa05320:Autoimmune thyroid disease	20	6.55E‐12
	hsa05150:Staphylococcus aureus infection	19	1.41E‐10
	hsa05416:Viral myocarditis	18	3.25E‐09
	hsa05140:Leishmaniasis	19	1.94E‐08
	hsa05168:Herpes simplex infection	31	2.72E‐08
GOTERM_MF (top 10)	GO:0003823~antigen binding	30	2.25E‐16
	GO:0042605~peptide antigen binding	14	3.36E‐11
	GO:0032395~MHC class II receptor activity	11	5.24E‐11
	GO:0004252~serine‐type endopeptidase activity	34	1.76E‐08
	GO:0008009~chemokine activity	13	8.14E‐07
	GO:0031730~CCR5 chemokine receptor binding	5	2.10E‐04
	GO:0048248~CXCR3 chemokine receptor binding	4	7.51E‐04
	GO:0030246~carbohydrate binding	20	8.57E‐04
	GO:0005518~collagen binding	10	0.001
	GO:0046982~protein heterodimerization activity	35	0.002
GOTERM_CC (top 10)	GO:0071556~integral component of lumenal side of endoplasmic reticulum membrane	17	3.59E‐15
	GO:0042613~MHC class II protein complex	14	4.33E‐13
	GO:0012507~ER to Golgi transport vesicle membrane	18	2.00E‐11
	GO:0005576~extracellular region	123	1.02E‐10
	GO:0070062~extracellular exosome	182	1.27E‐09
	GO:0005615~extracellular space	101	1.84E‐08
	GO:0042612~MHC class I protein complex	8	6.66E‐08
	GO:0030669~clathrin‐coated endocytic vesicle membrane	13	7.58E‐08
	GO:0009897~external side of plasma membrane	29	9.90E‐08
	GO:0030658~transport vesicle membrane	12	3.02E‐07
GOTERM_BP (top 10)	GO:0006955~immune response	87	5.38E‐35
	GO:0050776~regulation of immune response	42	3.44E‐19
	GO:0006958~complement activation, classical pathway	29	6.54E‐16
	GO:0006956~complement activation	27	1.49E‐15
	GO:0060333~interferon‐gamma‐mediated signaling pathway	23	1.01E‐13
	GO:0038096~Fc‐gamma receptor signaling pathway involved in phagocytosis	29	5.72E‐13
	GO:0002504~antigen processing and presentation of peptide or polysaccharide antigen via MHC class II	11	2.92E‐10
	GO:0019882~antigen processing and presentation	16	6.33E‐09
	GO:0006898~receptor‐mediated endocytosis	29	6.65E‐09
	GO:0038095~Fc‐epsilon receptor signaling pathway	28	1.05E‐08

Abbreviations: BP, biological process; CC, cellular component; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

PPI network analysis

With the purpose of exploring the interactions among proteins of DEGs, PPI network was established on the basis of STRING database containing 321 nodes and 1437 interactions (Figure 2). The top 20 nodes were regarded as hub genes (Table 3), mainly including B2M (beta‐2‐microglobulin), C3 (complement 3), PIK3R1 (phosphoinositide‐3‐kinase regulatory subunit 1), HLA‐E (MHC, Class I, E), HLA‐A (MHC, Class II, DR beta 1), HLA‐DRA (MHC, Class II, DR alpha), HLA‐DRB1 (MHC, Class II, DR beta 1), HLA‐DRB5 (MHC, Class II, DR beta 5), HLA‐DQB1 (MHC, Class II, DQ beta 1), HLA‐DQA1 (MHC, Class II, DQ alpha 1), HLA‐DQB2 (MHC, Class II, DQ beta 2), and HLA‐DPA1 (MHC, Class II, DP alpha 1).

FIGURE 2

Protein–protein interaction network of differentially expressed genes. The red circular nodes show up‐regulated genes and green circular nodes represent down‐regulated genes. The node size suggests degree value

TABLE 3

Top 20 hub genes in protein–protein network

Name	Degree	Type
B2M	43	Down‐regulation
C3	38	Down‐regulation
PIK3R1	38	Down‐regulation
HLA‐E	37	Down‐regulation
HLA‐A	37	Down‐regulation
HLA‐DRA	34	Down‐regulation
HLA‐DRB1	34	Down‐regulation
FYN	32	Down‐regulation
HLA‐DRB5	31	Down‐regulation
HLA‐DQB1	31	Down‐regulation
HLA‐DQA1	31	Down‐regulation
HLA‐DQB2	31	Down‐regulation
HLA‐DPA1	31	Down‐regulation
HLA‐DQA2	31	Down‐regulation
HLA‐B	30	Down‐regulation
HLA‐C	30	Down‐regulation
OAS2	30	Down‐regulation
HLA‐F	29	Down‐regulation
OAS1	29	Down‐regulation
HIST1H2BJ	25	Up‐regulation

Typically, an individual gene tends to exert its function through interacting with others and genes in a same module play similar biological roles. Therefore, we further extracted the significant modules from PPI network using the cytoscape according to the criteria described in methods. Our results suggested that there were six distinct modules (Figure 3). Furthermore, we performed the functional enrichment analyses of those DEGs in sub‐modules. The results showed that genes like (HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1) in module 1 were obviously enriched in immune‐related pathways and T1DM. Those genes in modules 2 and 3 were primarily involved in oxidative phosphorylation and chemokine signaling pathway, respectively. Additionally, the genes in module 4 played essential roles in ubiquitin‐mediated proteolysis processes and those in module 5 were mainly implicated with spliceosome. Finally, we discovered that genes in module 6 were responsible for systemic lupus erythematosus and alcoholism pathway (Table S1).

FIGURE 3

Sub‐modules of protein–protein interaction network of differentially expressed genes. The red circular nodes represent up‐regulated genes and green circular nodes stand for down‐regulated genes. The node size stands for degree value

The results of qRT‐PCR

To further confirm the underlying roles of four HLA genes, we carried out a qRT‐PCR by adopting the methods presented in methods. The findings showed that these four genes were down‐regulated in GDM samples in comparison with controls, offering the evidence that HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 might be significant therapeutic targets for GDM management (Figure 4).

FIGURE 4

Results of the quantitative real‐time PCR

DISCUSSION

In the current work, we screened 528 up‐regulated DEGs and 684 down‐regulated between GDM and healthy control samples. According to the KEGG results, they were preferentially enriched in T1DM process while GO annotation analysis suggested that these genes were mainly implicated with the MHC class II protein complex binding. Moreover, down‐regulated HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 served as hub genes in PPI network. Apart from that, six significant sub‐clusters were extracted and functional analysis suggested that several genes such as HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 in cluster 1 were significantly enriched in immune and T1DM‐related pathways. Meanwhile, those genes in other clusters were also involved in multiple critical cellular processes such as oxidative phosphorylation and chemokine signaling pathway.

GDM has been an increasingly serious public health issue globally over the past decades. A previous research has suggested that maternal GDM severely influences the living quality of descendants of those mothers suffering from diabetes in terms of attention disorder. Although many researchers argued that the pharmacotherapy such as metformin intervention could effectively retard progression toward type 2 DM, the prevalence of GDM was still continuously rising. Herein, our findings proved that those DEGs identified were predominately associated with immune and T1DM pathways, which were also supported by the existing evidence that the dysregulation of immune responses was responsible for GDM development. Recently, Sifnaios et al claimed that there existed a different T‐cell proportion between patients with GDM and healthy individuals, especially for type 1 and 2 T‐helper cells. More interestingly, Warncke et al found that there were higher numbers of lymphocytes and T‐tregs which were critical for immune system in pregnancies undergoing T1DM. In addition, Negishi et al. also pointed out that patients undergoing GDM tended to present a high risk of developing T1DM and numerous autoimmune disorders. Collectively, our results further provided evidence for previous findings that dysregulated immune responses might participate in the pathogenesis of GDM.

Notably, our findings also found that there were six distinct sub‐modules in PPI network. Moreover, multiple members of HLA gene system were clustered in sub‐module 1 and the functional analyses indicated that genes including down‐regulated HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 were dramatically enriched in immune‐correlated pathways and T1DM process. Accumulating evidence has demonstrated that HLA was located on chromosome 6 and has been identified to be a vital susceptibility gene for developing T1DM. More importantly, increasing investigations showed that many of the pregnancies with GDM had a higher risk of developing T1DM, and the examination of T1D autoantibodies was crucial for women suffering from GDM. Besides, recent studies suggested that HLA class II allele might make contributions to the incidence of T1DM. For example, Rajaei et al. stated that HLA‐DQA1 and HLA‐DQB1 were strongly linked to T1DM, also emphasizing that there existed close correlations between HLA alleles and autoimmune disorders, which indirectly provided support for our result that HLA‐DQA1 and HLA‐DQB1 exerted prominent roles in some immune responses. HLA‐DPA1, another under‐expressed gene, has been reported to signally enhance the risk of developing T1DM. Interestingly, Guo et al performed a meta‐analysis of GDM risk and HLA class II mutations (DQB1 and DRB1), finding that DQB1*02 and DRB1*1302 were related to elevated risks for progression to GDM. Therefore, it can be inferred that four HLA class II variants (HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1) were probably involved in the molecular mechanisms of GDM progression.

In the current work, although we have identified multiple prominent gene markers through the bioinformatics analysis and RT‐PCR verification, there still exist limitations. Firstly, comprehensive bioinformatics studies with larger sample sizes are required to validate our findings. Apart from that, while we found that many critical signaling pathways such as immune responses might participate in the pathogenesis of GDM, the precise regulatory mechanism has not been illuminated. Finally, the clinical information is also required to be integrated into a detailed analysis. Despite the above‐mentioned limitations, the current work offers a solid theoretical foundation for the diagnosis and treatment of GDM. It is well known that GDM, as a prevalent metabolic complication during pregnancy, shows close association with an incremental risk of adverse pregnancy outcomes. Thus, the in‐depth insight into pathogenesis and coping strategies for conquering GDM are urgently required. The present study revealed the close association of dysregulated immune responses with the development of FDM, also strengthening the understanding of the pathogenesis of GDM. Furthermore, we identified HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 as potential biomarkers for diagnosing GDM, which may make contributions to the early prediction of GDM, thereby enhancing the adverse outcomes of both mothers and their offspring.

In summary, we extracted the DEGs between GDM and control samples, finding that pathways of immune responses and T1DM were possibly associated with molecular mechanism of GDM. Moreover, HLA‐DQA1, HLA‐DRB1, HLA‐DPA1, and HLA‐DQB1 exerted vital roles in immune‐related pathways, providing deeper insights into the development of novel therapeutic strategies for managing GDM. However, some cytological evidence will be required with the purpose of confirming our findings.

CONFLICT OF INTEREST

The authors claim that they have no competing interests.

CONSENT FOR PUBLICATION

Not applicable.

Table S1

Click here for additional data file.