Clinical characteristics of pregnant women and infants
Anthropometric and metabolic parameters of the 30 pregnant women and newborns after screening are shown in Table 2, apart from significant differences in OGTT indicators (t-test, P value < 0.05) between the two groups, other characteristics such as maternal age, gestational weight gain, maternal height, gestational week and fetal birth weight did not show significant differences.
Distribution and characteristics of DMPs and DMRs in umbilical cord blood of pregnant women affected by gestational diabetes mellitus
The occurrence of DMPs and DMRs may regulate their neighbouring genes. Identifying DMRs can more accurately find the DNA methylation regulatory genes. After analysing the microarray sequencing data of 30 samples, we identified 23331 hyper-methylated DMPs (hyper-DMPs) containing 9327 genes (Additional file 2: Table. S1) and 29501 hypo-methylated DMPs (hypo-DMPs) containing 10908 genes (Additional file 2: Table. S2). 42 genes contained in the hyper-methylated DMRs (hyper-DMRs) and 209 genes contained in the hypo-methylated DMRs (hypo-DMRs) were also identified (Additional file 2: Table. S11). DMPs showed significant changes on all chromosomes (Fig. 1A). GDM cord blood DNA promoter regions (1stExon, TSS200, TSS1500, 5’utr) and IGR regions contained more hypo-DMPs than the control group (Fig. 1B). According to the distance of DNA methylation sites from CpG islands, the methylation sites can be divided into four categories: “cgi”, “shore”, “shelf” and “opensea” [19], In these four regions, even if DMPs exist in the promoter region, they do not necessarily exist in the CpG islands, and even if CpG islands are not necessarily in the gene promoter region, showing the irregular distribution of methylation sites on the gene (Fig. 1C). DMRs are associated with DMPs. Because there are more hypo-DMPs than hyper-DMPs, and the hypo-DMPs are dominant in the promoter region, there are more hypo-DMRs than hyper-DMRs (Fig. 1D). This is because the promoter regions of most genes in the human body have CpG islands, so the change of DMPs in the promoter region plays a decisive role in the calculation of DMRs. The number of hypo-DMRs was also higher than that of hyper-DMRs on each chromosome (Fig. 1E). Due to different regions, different sequencing technologies and other reasons, sometimes the data obtained cannot be reproduced, and some genes are occasionally absent. To solve these problems, we introduce external datasets for multiple validations.
Dataset Integration Analysis Identified Genes Whose Expression Was Altered Due To Changes In Dna Methylation
The exact process was used for DNA methylation microarray data to obtain public data more efficiently and reliably. We identified 41678 hyper-DMPs and 13747 hypo-DMPs from the GSE153219 dataset, which contained 12518 and 6533 genes, respectively (Additional file 2: Table. S3, S4). 176 hyper-DMPs containing 110 genes (Additional file 2: Table.S5) and 66 hypo-DMPs containing 42 genes (Additional file 2: Table. S6) were identified in the GSE141065 dataset. 606 hyper-DMPs containing 524 genes (Additional file 2: Table.S7) and 9 hypo-DMPs containing 8 genes (Additional file 2: Table. S8) were identified in the GSE88929 dataset. 2025 genes in the hyper-DMPs (Additional file 2: Table.S9) and 1646 in the hypo-DMPs (Additional file 2: Table. S10) were identified in the GSE141065 dataset. DMRs have identified fewer genes, and we identified 72 genes containing in the hyper-DMRs and 34 genes containing in the hypo-DMRs from the GSE153219 dataset (Additional file 2: Table. S12), 28 genes in the hyper-DMRs and 19 in the hypo-DMRs were identified in the GSE141065 dataset (Additional file 2: Table. S13), 136 genes in the hyper-DMRs and 15 in the hypo-DMRs were placed in the GSE88929 dataset (Additional file 2: Table. S14). 60 genes in the hyper-DMRs (Additional file 2: Table. S15) and 71 in the hypo-DMRs (Additional file 2: Table. S16) were identified in the GSE203169 dataset. Subsequently, these data were integrated and analysed by upset package in R, and it was found that the genes contained in the DMPs and DMRs in the dataset were sometimes absent in other datasets (Fig. 2A, B, C, D). To make the data more reliable, the genes contained in DMPs and DMRs obtained from the public database were used as the control group, the genes contained in DMPs and DMRs received in the Test dataset were used as the experimental group for integration analysis, and the significant genes that were repeated at least twice were identified. Finally, we identified 8797 genes by DMPs (Additional file 2: Table. S17) and 23 genes by DMRs (Additional file 2: Table. S18). The identified genes were then integrated with the 136 DEGs (Additional file 2: Table. S19) placed in the GSE203346 dataset to identify 52 genes (Table 3) and 0 genes (Fig. 2E, F), respectively. In this way, we identified genes whose mRNA expression was substantially altered in the presence of GDM because of the affecting cord blood DNA methylation.
Insert Table 3
Functional Enrichment Analysis
To verify which genetic pathways and functions are affected by GDM, we used GO enrichment analysis to identify biologically significant gene annotations. A total of 7 GO terms were identified (p. adjust < 0.05) (Fig. 3A, Additional file 2: Table. S20). Mainly includes "regulation of chondrocyte differentiation" (PTHLH/ACVRL1/ADAMTS12/ ZNF219/SHOX2, p. adjust = 8.59e-05), "Regulation of cartilage development" (PTHLH/ACVRL1/ADAMTS12/ZNF219/SHOX2, p. adjust = 2.15e-04), " negative regulation of cellular response to vascular endothelial growth factor stimulus " (DCN/IL12B/ADAMTS12, p. adjust = 1.97e-03), These GO terms are mainly related to the development and differentiation of chondrocytes and connective tissue and the regulation of cytokines. Subsequently, KEGG pathway analysis was performed, and a total of 57 pathways were found (Additional file 2: Table. S21), but there were no significant pathways with a threshold value of p. adjust < 0.05, so we used p-value < 0.05 as the threshold value found 4 pathways (Fig. 3B). "Arrhythmogenic right ventricular cardiomyopathy" (CDH2/SGCD, p-value = 0.01), "Toll-like receptor signalling pathway" (IL12B/TLR3, p-value = 0.03), "Non-homologous end-joining" (DNTT, p-value = 0.03), “Cytokine-cytokine receptor interaction” (IL12B/ACVRL1/IL27, p-value = 0.03), These pathways are associated with congenital heart disease, innate immunity and DNA repair.
Wgcna Identified Gdm-associated Modules And Functional Enrichment Analysis Of Modules
Subsequently, we tried to use STRING (https://cn.string-db.org/) to explore the protein interaction network to explore the interactions among these genes further and find the hub genes in these genes. It was found that there were few associations between these genes, and then we changed our thinking and used WGCNA to find the modules related to GDM. WGCNA is a commonly used bioinformatics method to identify gene sets of interest using the information of thousands or nearly ten thousand genes with the most changes or all genes and to perform significant association analysis with phenotype. We used Fig. 8 as the soft threshold (Additional file 1: Fig. S1A); the results showed that the blue module was significantly associated with GDM (p. adjust = 0.02) (Additional file 1: Fig. S1B). There were 58 genes in the blue module (Table 4), most of which were histone-related protein genes. GO enrichment analysis of these genes identified a total of 91 GO Terms (p. adjust < 0.05) (Fig. 4A, Additional file 2: Table. S23), The top GO terms were, "Defense response to Gram-negative bacterium" (RNASE3 LTF/ELANE/DEFA4/DEFA1/CTSG/AQP1, p. adjust = 1.08e-10), "Antimicrobial humoral response" (RNASE3/PRTN3/LTF/ELANE/DEFA4/DEFA1/CTSG, p. adjust = 6.14e-10), "Antibacterial humoral response" (RNASE3/LTF/ELANE/DEFA4/DEFA1/CTSG, p. adjust = 6.52e-10), " Defense response to fungus " (LTF/ELANE/DEFA4/DEFA1/CTSG, p. adjust = 4.22e-08). Most GO terms were related to antibacterial and immune function in vivo. KEGG pathway enrichment analysis revealed 29 pathways (Additional file 2: Table. S24); there were also no significant pathways with a threshold value of p. adjust < 0.05, so we used p-value < 0.05 as the threshold value found 9 pathways (Fig. 4B), The top KEGG pathways were “Transcriptional misregulation in cancer“ (ELANE/DEFA4/DEFA1, p-value = 0.002), "Staphylococcus aureus infection" (DEFA4/DEFA1, p-value = 0.007), " Systemic lupus erythematosus” (ELANE/CTSG, p-value = 0.014), Some of these pathways are also involved in viral infection, bacterial infection and immune regulation.
Insert Table 4
Interaction Between Hub Genes And Cd4t Cells
We get two gene sets from the integration analysis of DMPs and transcriptome to WGCNA analysis. We found that the functions of these gene sets affect immune-related functions in vivo, and previous studies have also found immune dysregulation in GDM pregnant women. To explore whether GDM pregnant women affect the immune regulation of their offspring, Here, we performed immune infiltration analysis on the self-measured methylation EPIC BeadChip data and found that the proportion of immune cells in the 30 samples was generally consistent (Fig. 5A, Additional file 2: Table. S25). However, the ratio of CD4T in cord blood was significantly lower in the GDM group (Fig. 5B). We further identified two genes, HIST1H1B and PRTN3, from the integrated analysis of the two gene sets. 6 DMPs associated with these two genes in the Infinium Methylation EPIC BeadChip (Table 5), respectively cg24587427, cg00700007, cg27535410, cg27234823, cg05257528, and cg21462596. To explore the relationship between these sites and CD4T cells, we linked six sites in each sample to CD4T cells (Fig. 5C) and found that cg24587427 and cg00700007 were more closely associated with CD4T cells. These two loci are located on chromosome 6 and are present within the promoter region. Although the four loci belonging to PRTN3 were not closely associated with CD4T cells, they appeared significant in both gene sets and were believed to be important in other aspects.
Insert Table 5
Constructing Cerna Network
To further explore the ncRNAs associated with the two hub genes, we identified 69 DELs from the GSE203346 dataset (Additional file 2: Table. S27) and 891 DELs (Additional file 2: Table. S28) obtained from the GSE65737 dataset. Subsequently, “miRcode” (http://mircode.org/) was used to predict the miRNA of lncRNA interaction, and two hub genes related miRNAs were found. Then the ceRNA network was constructed (Fig. 6). The identified miRNA-16, 17, 20 can be used as a biomarker for predicting maternal GDM [20], and lncRNA TTTY10 was also found to be significantly higher in GDM than in the standard group [21]. These miRNAs and lncRNAs interacted with HIST1H1B and PRTN3, further proving these two genes’ importance.
Validation Of Gene Expressions
Finally, we used qPCR to validate HIST1H1B, PRTN3, DNMT1, DNMT3a, and DNMT3b (Fig. 7). The results showed that GDM significantly reduced (p-value < 0.05) the expression of these genes. In general, low expression of DNMT1, DNMT3a and DNMT3b accelerates gene expression. However, DNMT1 and DNMT3b were significantly reduced in obese pregnant women [22]. Compared with GDM, the expression level in obese pregnant women decreased relatively less. These results suggest that GDM is associated with obesity.