Biological networks in gestational diabetes mellitus: insights into the mechanism of crosstalk between long non-coding RNA and N6-methyladenine modification.

BACKGROUND: Gestational diabetes mellitus (GDM) is one of the most common complications of pregnancy. The mechanism underlying the crosstalk between long non-coding RNAs (lncRNAs) and N6-methyladenine (m6A) modification in GDM remain unclear.
METHODS: We generated a lncRNA-mediated competitive endogenous RNA (ceRNA) network using comprehensive data from the Gene Expression Omnibus database, published data, and our preliminary findings. m6A-related lncRNAs were identified based on Pearson correlation coefficient (PCC) analysis using our previous profiles. An integrated pipeline was established to constructed a m6A-related subnetwork thereby predicting the potential effects of the m6A-related lncRNAs.
RESULTS: The ceRNA network was composed of 16 lncRNAs, 17 microRNAs, 184 mRNAs, and 338 edges. Analysis with the Kyoto Encyclopedia of Genes and Genomes database demonstrated that genes in the ceRNA network were primarily involved in the development and adverse outcomes of GDM, such as those in the fatty acid-metabolism pathway, the peroxisome proliferator-activated receptor signaling pathway, and thyroid hormone signaling pathway. Four m6A-related lncRNAs were involved in the ceRNA network, including LINC00667, LINC01087, AP000350.6, and CARMN. The m6A-related subnetwork was generated based on these four lncRNAs, their ceRNAs, and their related m6A regulators. Genes in the subnetwork were enriched in certain GDM-associated hormone (thyroid hormone and oxytocin) signaling pathways. LINC00667 was positively correlated with an m6A "reader" (YTHDF3; PCC = 0.95) and exhibited the highest node degree in the ceRNA network. RIP assays showed that YTHDF3 directly bind LINC00667. We further found that MYC possessed the highest node degree in a protein-protein interaction network and competed with LINC00667 for miR-33a-5p. qPCR analysis indicated that LINC00667, YTHDF3 and MYC levels were upregulated in the GDM placentas, while miR-33a-5p was downregulated. In a support-vector machine classifier, an m6A-related module composed of LINC00667, YTHDF3, MYC, and miR-33a-5p showed excellent classifying power for GDM in both the training and the testing dataset, with an accuracy of 76.19 and 71.43%, respectively.
CONCLUSIONS: Our results shed insights into the potential role of m6A-related lncRNAs in GDM and have implications in terms of novel therapeutic targets for GDM.

Background

Gestational diabetes mellitus (GDM), a condition involving metabolic dysfunction during pregnancy, can cause short-term and long-term adverse effects on both the mother and offspring [1]. It has been reported that GDM affects 2–5% of pregnancies worldwide, due to a diverse genetic background and epigenetic modifications that occur in response to nutritional and environmental factors [1, 2]. Currently, the understanding of the precise etiological mechanisms of GDM remains unclear.

Dysregulation of long non-coding lncRNAs (lncRNAs) has been reported to participate in numerous human diseases [3-5], and aberrant lncRNA expression is associated with GDM pathogenesis [6-9]. It has been reported that metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) expression was elevated in GDM placentas, and its over-expression could suppress the proliferation, invasion, and migration of trophoblast cells via the TGF-β/NF-κB signaling pathway [7]. Reduced expression of lncRNA plasmacytoma variant translocation 1 (PVT1) was observed in the placentas of patients with GDM and preeclampsia, and knocking down PVT1 enhanced apoptosis and inhibited the proliferation, migration, and invasion of trophoblast cells through the PI3K/Akt pathway [6]. lncRNAs exhibit extensive functions and regulate the expression of genes at multiple levels. Several lncRNAs can function as competitive endogenous RNAs (ceRNAs) to regulate the expression of downstream genes by competing for their shared microRNAs (miRNAs) [10, 11]. Ye et al. [12] revealed that lncRNA maternally expressed gene 3 (MEG3) was up-regulated in the human umbilical vein endothelial cells (HUVECs) extracted from GDM pregnancies, and influence the fetal endothelial function through targeting AFF1 via sponging miR-370-3p. Overexpression of MEG3 was also obseved in the blood and placental villous tissue in pregnant GDM patients, and MEG3 could serve as a ceRNA of miR-345-3p to regulate the biological behavior and cell cycle of trophoblasts [9].

Chemical modifications occurring in lncRNAs can modify their secondary structure, splicing, degradation, or molecular stability, which can influence the expression of the lncRNAs and their downstream elements [13, 14]. N6-methyladenine RNA modification (m6A) is the most prevalent type of RNA epigenetic process and is regulated by m6A regulators, including methyltransferases (“writers”), signal transducers (“readers”), and demethylases (“erasers”) [13]. Recent evidence indicates that perturbations of m6A modifications dysregulate glucose/lipid metabolism and the immune/inflammatory response, thereby contributing to obesity, diabetes and cardiovascular diseases [15-17]. However, the knowledge regarding GDM is still in its infancy.

In this study, we focused on a lncRNA-mediated ceRNA network, m6A-related lncRNAs, crosstalk between them, and the potential effects on GDM.

Materials and methods

High-throughput expression profiles for GDM

Expression profiles for lncRNAs and mRNAs were generated by our research team [18] and are referred to here collectively as the GSE Shengjing profile. Samples were extracted from the placentas (maternal side) of women with normal glucose tolerance (NGT; n = 3) and GDM (n = 3). The human mRNA-expression datasets GSE2956, GSE19649, and GSE70493 were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). GSE2956 contains data for placental samples from 6 patients (GDM; n = 3 and NGT; n = 3). GSE19649 consists of data from one pooled GDM placental tissue and one pooled healthy placental tissue. The expression data in GSE70493 were obtained for validation and contain 63 maternal side placental samples (GDM; n = 32 and NGT; n = 31).

Screening of differentially expressed RNAs

Using DEseq2 [19], we identified differentially expressed mRNAs (DEMs) and differentially expressed lncRNAs (DELs) in GSE Shengjing profile, based on the criteria of |log2 fold change| > 1 and adjusted P value < 0.05. A total of 172 lncRNAs and 142 mRNAs were differentially expressed in GSE Shengjing profile, respectively. Of these, 86 lncRNAs and 67 mRNAs were upregulated, and 86 lncRNAs and 75 mRNAs were downregulated. All these differentially expressed genes can also be identified using Limma package [20]. We only extracted lncRNAs with matching transcript IDs for analysis. Genes without official symbols were removed, and all symbols were converted to symbols approved by the HUGO Gene Nomenclature Committee. Heatmaps were conducted using the ‘pheatmap’ package of R software. The aberrantly expressed mRNAs in GSE2956 and GSE19649 that met the criterion of |log2 fold-change| > 1 were also pooled as DEMs for subsequent analysis.

The lncRNA-mediated ceRNA network was constructed on the basis that the miRNA was shared by any pair of ceRNA genes. In such cases, the reliability of miRNA was critical to establish a reliable lncRNA-miRNA-mRNA network. We systematically reviewed PubMed and pooled the dysregulated placenta-specific microRNAs (miRNAs) in GDM with ≥2 supporting experiments [21-27]. Those differentially expressed miRNAs (DEMis) with contradictory findings in different studies were further excluded.

Collection of m6A regulators and identification of m6A-Related DELs

We pooled 21 m6A regulators from previous publications, including writers (METTL14, METTL16, METTL3, RBM15, RBM15B, VIRMA [KIAA1429], WTAP, and ZC3H13), readers (HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3), and erasers (FTO and ALKBH5). Pearson correlation coefficient (PCC) analysis between m6A regulators and DELs was conducted to identify the m6A-related DELs. A |PCC| of > 0.9 and P value of < 0.01 were used as the cutoff criteria.

Construction of a lncRNA–miRNA–mRNA network associated with GDM

miRNA–lncRNA interactions were obtained from DIANA-LncBase v3 and contained > 500,000 experimentally supported miRNA targets on non-coding transcripts [28]. DIANA-TarBase [29] was used to retrieve the mRNA–miRNA-association data with experimental evidence. Finally, a lncRNA–miRNA–mRNA-regulatory network was constructed based on ceRNA theory using matching DEL–DEMi and DEMi–DEM pairs. Cytoscape 3.7.2 software (http://www.cytoscape.org/) was used to visualize the ceRNA network. All node degrees were calculated. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed using the ‘enrichplot’ package of R software. A P value of < 0.05 was set as the cutoff. The Search Tool for the Retrieval of Interacting Genes [30] was applied to predict the protein-protein interactions (PPIs) of DEMs in the ceRNA network with a confidence of > 0.4. We visualized the PPI network using Cytoscape and excavated the most tightly interconnected modules using the MCODE plug-in. We also screened for hub genes in the PPI network using CytoHubb with a node degree of > 15.

Identification of the m6A-related subnetwork and module

To filter the core m6A-related DELs, the m6A-related DELs were overlapped with lncRNAs in the ceRNA network. The m6A-Related subnetwork composed of the core m6A-related DELs, their ceRNAs, and the host m6A regulators was visualized by Cytoscape, and the enrichment assay was further analyzed. An m6A-related module containing a lncRNA with the highest node degree in the ceRNA network, the gene with the highest node degree in the PPI network, their shared miRNA target, and the host m6A gene of the lncRNA were extracted.

Construction of a Support-Vector Machine (SVM) classifier

The support-vector machine (SVM) is comprised of a set of supervised-learning methods for performing binary and multi-class classifications, which has been used to classify and diagnose diseases [31-33]. A SVM classifier based on the selected m6A-related module was established to distinguish samples from subjects with GDM and healthy subjects by analyzing m6A-related module genes using the python package ‘scikit-learn 0.24.0’. Independent training/testing procedures using the SVM classifier were performed to predict the disease status associated with samples in GSE70493.

Clinical samples and cell culture

All protocols for the use of human samples were approved by the Ethics Committee of Shengjing Hospital of China Medical University, in accordance with the Helsinki Guidelines. Thirty-two patients with GDM and 32 healthy controls were recruited from Shengjing Hospital of China Medical University. GDM was diagnosed according to the Chinese Current Care Guidelines for GDM [34]. The control group was matched to the GDM group through a 1:1 pattern according to age and pre-pregnancy body mass index (BMI). Patients who met the criteria as follows were excluded: younger than 18 years old, history of diabetes, infective or inflammatory diseases, hypertension, chronic diseases (thyroid dysfunction, cardio-cerebrovascular diseases, renal failure, etc.), multiple pregnancies, or used assistive reproductive technology. After cesarean section, placental tissues were collected and stored at − 80 °C. The maternal information and data on neonatal outcomes were collected from medical records. HTR-8/SVneo cells were cultured in RPMI-1640 (Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (FBS, Hyclone, Logan, UT, USA) at 37 °C in 5% CO2.

Quantitative real-time polymerase chain reaction (qPCR)

We used TRIzol reagent (Invitrogen, Carlsbad, CA, USA) to extract total RNA from tissues or cells. RNA purity and concentration were determined using a NANO 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) at 260 and 280 nm. cDNA was generated from RNA using the BioTeke super RT kit (Bioteke, Beijing, China) according to the manufacturer’s protocol. Subsequently, qPCR was carried out using a SYBR GREEN mastermix (Solarbio, Beijing, China). Data were analyzed using the 2-ΔΔCt method and normalized to the expression levels of β-Actin or U6. All primers are described in Table S1.

RNA immunoprecipitation (RIP) assay

RIP assay was performed using RIP Kit (Millipore, Billerica, MA, USA) according to the manufacturer’s instructions. Briefly, HTR-8/SVneo cells were lysed with RIP lysate buffer. The lysate was incubated overnight in RIP buffer with magnetic beads coupled to YTHDF3 or IgG antibodies (Abcam, Cambridge, UK) at 4 °C. The expression of LINC00667 was analyzed by qPCR.

Statistical analysis

The Shapiro test and Levene’s test were used to determine the normal distribution and homogeneity, respectively. Quantitative variables are presented as mean ± standard deviation, whereas qualitative variables are presented as counts and/or percentages. Associations between quantitative variables were assessed using the Student’s t-test or the Mann–Whitney U-test. Qualitative variables were compared using the Chi-square test. Graphing and statistical analyses were performed using GraphPad Prism 9, R software (version 4.0.3) and Python (version 3.8). P < 0.05 was considered statistically significant. All reactions in vitro were performed in triplicate.

Results

Filtering results for DELs, DEMis, and DEMs

Based on the criteria described above, we identified 56 DELs, 38 DEMis, and 540 DEMs. Of these, 33 lncRNAs, 22 miRNAs, and 283 mRNAs were upregulated, and 23 lncRNAs, 16 miRNAs, and 257 mRNAs were downregulated. The heatmaps of the DELs and DEMs found in GSE Shengjing profile are shown in Fig. 1. The details of the pooled DEMis and the information related to the DEMs are listed in Tables 1 and 2, respectively. The list of all DEMs identified in this study is shown in Table S2.

Fig. 1

Hierarchical clustering of DELs and DEMs in the GSE Shengjing profile (A) The heatmap of the DELs. B The heatmap of the DEMs. The up- and down-regulation are indicated with orange and green color, respectively

Table 1

Detailed information related to the pooled DEMis

Reference	DEMi name	Location	Trend
Gillet 2019 [26]	miR-122-5p	Placental exosomes	+
	miR-132-3p	Placental exosomes	+
	miR-1323	Placental exosomes	+
	miR-136-5p	Placental exosomes	+
	miR-182-3p	Placental exosomes	+
	miR-210-3p	Placental exosomes	+
	miR-29a-3p	Placental exosomes	+
	miR-29b-3p	Placental exosomes	+
	miR-342-3p	Placental exosomes	+
	miR-520 h	Placental exosomes	+
Ding 2018 [21]	miR-138-5p	Placenta	–
	miR-202-5p	Placenta	+
	miR-210-5p	Placenta	–
	miR-3158-5p	Placenta	–
	miR-4732-3p	Placenta	–
Nair 2018 [22]	miR-125a-3p	Placenta, placental-derived exosomes, circulating exosomes and skeletal muscle	+
	miR-99b-5p	Placenta, placental-derived exosomes, circulating exosomes and skeletal muscle	+
	miR-197-3p	Placenta, placental-derived exosomes, circulating exosomes and skeletal muscle	+
	miR-22-3p	Placenta, placental-derived exosomes, circulating exosomes and skeletal muscle	+
	miR-224-5p	Placenta, placental-derived exosomes, circulating exosomes and skeletal muscle	+
	miR-584-5p	Placenta, placental-derived exosomes	+
	miR-186-5p	Placenta, placental-derived exosomes	+
	miR-433-3p	Placenta, placental-derived exosomes	+
	miR-423-3p	Placenta, placental-derived exosomes	+
	miR-208a-3p	Placenta, placental-derived exosomes	–
	miR-335-5p	Placenta, placental-derived exosomes	–
	miR-451a	Placenta, placental-derived exosomes	–
	miR-145-3p	Placenta, placental-derived exosomes	–
	miR-369-3p	Placenta, placental-derived exosomes	–
	miR-483-3p	Placenta, placental-derived exosomes	–
	miR-203a-3b	Placenta, placental-derived exosomes	–
	miR-574-3p	Placenta, placental-derived exosomes	–
	miR-144-3p	Placenta, placental-derived exosomes	–
	miR-6795-5p	Placenta, placental-derived exosomes	–
	miR-550a-3-3p	Placenta, placental-derived exosomes	–
	miR-411-5p	Placenta, placental-derived exosomes	–
	miR-140-3p	Placenta, placental-derived exosomes	–
Li 2018 [27]	MIR96	Placenta	–
Xu 2017 [23]	miR-503	Placenta	+
Muralimanoharan 2016 [24]	miR-143	Placenta	–
Li 2015 [25]	miR-508-3p	Placenta	+
	miR-27a	Placenta	–
	miR-9	Placenta	–
	miR-92a	Placenta	–
	miR-30d	Placenta	–
	miR-362 -5p	Placenta	–
	miR-502-5p	Placenta	–
	miR-33a	Placenta	–

Table 2

Overall information related to the DEGs identified in this study

	Upregulated	Downregulated	Total
GSE19649	108	35	143
GSE2956	133	185	318
GSEShengjing	42	37	79

Construction of the ceRNA network

A lncRNA–miRNA–mRNA-regulatory network for GDM was constructed (Fig. 2A) and was found to consist of 16 lncRNA nodes, 17 miRNA nodes, 184 mRNA nodes, and 338 edges. The node degrees of the DELs in the ceRNA network are shown in Table 3. LINC00667 was identified as the hub lncRNA as it exhibited the highest degree in the ceRNA network. Based on the GO category, biological process, DEMs in the ceRNA network mainly were enriched for response to drug, response to steroid hormone, fibroblast cell proliferation, positive regulation of cytokin production, and response to steroid hormone (Fig. 2B). KEGG pathways were mainly enriched in fatty acid metabolism, the peroxisome proliferator-activated receptor (PPAR) signaling pathway, and the thyroid hormone signaling pathway (Fig. 2C). A PPI network composed of 152 nodes and 475 edges was constructed and visualized (Fig. 3A). The key module (score = 8.75) was identified from the PPI network (Fig. 3B). The hub genes of the PPI network and the key module nodes are listed in Table 4.

Fig. 2

The lncRNA-mediated ceRNA network in GDM and functional analysis ranked according to the −log (P value). A Depiction of the lncRNA-mediated ceRNA network in GDM. The differentially expressed lncRNAs, mRNAs, and miRNAs are indicated with circle, hexagon, and triangle, respectively. The upregulation and downregulation are indicated with pink and blue shading, respectively. B The top 10 GO biological process categories. C The enriched KEGG pathways

Table 3

List of all DELs in the ceRNA network

DEL number	Name	Node degree	Trend
1	LINC00667	5	Upregulated
2	AC016831.1	4	Downregulated
3	CARMN	3	Downregulated
4	AC090114.2	3	Upregulated
5	AL589745.1	3	Upregulated
6	LINC00639	2	Downregulated
7	AP000350.6	2	Upregulated
8	AL136962.1	1	Downregulated
9	AC106028.3	1	Downregulated
10	LINC01483	1	Downregulated
11	LINC00921	1	Downregulated
12	GMDS-DT	1	Upregulated
13	GASAL1	1	Upregulated
14	LINC00174	1	Upregulated
15	AP002360.3	1	Upregulated
16	LINC01087	1	Upregulated

Fig. 3

Construction of a PPI network. A PPI network of DEMs in the ceRNA network (colored according to the node degree). B The key module was identified using the MCODE plug-in

Table 4

Hub genes and key module nodes in the PPI network

Hub gene		Key module
Annotation	Degree	Annotation	Degree
MYC	36	MYC	36
VEGFA	32	VEGFA	32
FN1	30	EP300	30
EP300	30	FN1	30
ACTB	28	ACTB	28
FOS	22	FOS	22
CREB1	16	PTGS2	16
PTGS2	16	CREB1	16
THBS1	15	CDKN1A	11
SCD	15

Identification of the m6A-related DELs in GDM

To identify m6A-related lncRNAs, we conducted Pearson correlation analysis using the expression matrixes of DELs and 21 m6A regulators from the GSE Shengjing profile. A PPC of > 0.9 and a P value of < 0.01 were set as the cutoff criteria. Ten DELs were significantly correlated with four m6A-related genes, as shown in Fig. 4.

Fig. 4

Correlation heatmap of DELs and m6A regulators

Characterization of the m6A-related subnetwork and module

Four m6A-related DELs were involved in the ceRNA network, including LINC00667, LINC01087, AP000350.6, and CARMN. A m6A-related subnetwork was identified containing 4 lncRNAs, 10 miRNAs, 101 mRNAs, and 3 m6A regulators (Fig. 5A). Functional analysis indicated that the major BP terms were aging, response to the steroid hormone, learning and cognition, and cell proliferation (Fig. 5B). KEGG pathways were mainly enriched in the thyroid hormone signaling pathway, Oxytocin signaling pathway, virus infection, and cancer (Fig. 5C).

Fig. 5

The m6A-related subnetwork in GDM and functional analysis ranked according to the −log (P value). A Depiction of the m6A-related subnetwork in GDM. The differentially expressed lncRNAs, mRNAs, miRNAs, and m6A regulators are indicated with yellow, pink, purple, and green shading, respectively. B The top 10 GO biological process categories. (C) The enriched KEGG pathways

We further identified an m6A-related module to explore molecules associated with GDM, which consisted of four nodes (LINC00667, YTHDF3, MYC, and miR-33a-5p) (Fig. 6A). LINC00667 exhibited the highest node degree among all lncRNAs in the ceRNA network. The heatmap (Fig. 4) shows that the expression level of the hub lncRNA (LINC00667) correlated with that of YTHDF3 (PCC = 0.95), which is an m6A-related reader. MYC exhibited the highest node degree in the PPI network. The subnetwork plot (Fig. 5A) shows that LINC00667 can act as ceRNA for miR-33a-5p to regulate MYC.

Fig. 6

Identification of an m6A-related module. A The m6A-related module consisted of LINC00667, YTHDF3, MYC, and miR-33a-5p. B Venn diagram showing the overlap between modeled (predicted) groups and real groups in the training dataset. C Venn diagram showing the overlap between modeled (predicted) groups and real groups in the testing dataset. D-G The RNA expression levels of LINC00667, YTHDF3, MYC, and miR-33a-5p were analyzed using qPCR. ** P < 0.01. H The interaction between YTHDF3 and LINC00667 was analyzed by RIP from HTR8/SVneo cells. ** P < 0.01

Validation of the m6A-related module

The expression matrix of the four genes in the m6A-related module was extracted from dataset GSE92772, and then the SVM classifier was constructed. We found that 32 of 42 samples and 15 of 21 samples were correctly classified in training set and testing set, with an accuracy of 76.19 and 71.43%, respectively (Fig. 6B, C). We further determined the expression levels of LINC00667, YTHDF3, MYC, and miR-33a-5p in placental tissues (tissues of 32 GDM patients and 32 controls) by qPCR. The clinical data are summarized in Table S3. In addition to indicators of impaired glucose metabolism, we observed more weight gain during gestation on patients with GDM compared with control subjects. LINC00667, YTHDF3, MYC levels were significantly higher in the GDM group than in the NGT group (Fig. 6D-F), while miR-33a-5p was significantly decreased (Fig. 6G). A total of 47 m6A sites in LINC00667 were predicted using RMBase [35]. The details of predicted m6A sites in LINC00667 were shown in Supplementary file 1. Moreover, RIP assays showed that YTHDF3, which is a m6A reader, directly bind LINC00667 (Fig. 6H). These indicated that the stability of LINC00667 may be affected through m6A modification via YTHDF3. Existing literature indicates that both LINC00667 and MYC serve as targets of miR-33a-5p [36-41]. This suggests that YTHDF3/LINC00667/miR-33a-5p/MYC axis maybe potential targets for research on the mechanism of GDM.

Discussion

GDM is a condition of pregnancy-related hyperglycemia with increased morbidity and mortality, for both the mother and fetus. The etiology of GDM, which involves the genetic background and epigenetic modifications, remains unclear. We selected the placenta as the subject of this study based on its specific location between the maternal and fetal bloodstreams. The placenta is exposed to intrauterine conditions that adversely affect placental and fetal development, which explains why most GDM-related adverse pregnancy outcomes originate in the placenta [42].

As a recently discovered type of non-coding RNAs, lncRNAs participate in numerous cellular functions and metabolic diseases [43]. However, only a limited number of studies have been conducted to explore the relationship between lncRNAs and GDM. Cao et al. [44] performed a microarray-based expression-profile analysis, which revealed aberrantly expressed lncRNAs in umbilical cord blood exosomes from patients with GDM. MEG3, MEG8, and MALAT1 upregulation, and PVT1 downregulation were involved in GDM development and the adverse effects on offspring, which provides novel biomarkers or therapeutic targets for GDM [6, 7, 9, 45]. Certain lncRNAs, miRNAs, and genes can form ceRNA motifs, which may participate in the molecular mechanism of GDM. In this study, a lncRNA-mediated ceRNA regulatory network was generated and analyzed. Enrichment analysis was conducted to determine the biological functions enriched for among DEMs in the ceRNA network. Notably, several pathways were closely correlated with the development and the adverse outcomes of GDM, where the fatty acid-metabolism pathway, the PPAR signaling pathway, and the thyroid hormone signaling pathway were the top three pathways [46-49].

As the most abundant epigenetic form of mRNA and non-coding RNA (ncRNA) methylation, m6A modification plays vital roles in various diseases [15, 17, 50]. However, it remains unknown how m6A modification acts in a lncRNA-dependent manner in GDM. We identified m6A-related lncRNAs based on a dataset generated in our previous study [18]. The m6A-related lncRNAs were intersected with the DELs in the ceRNA network to obtain four lncRNAs: LINC00667, LINC01087, AP000350.6, and CARMN. Subsequently, a m6A-related subnetwork was generated based on these four lncRNAs, their ceRNAs, and their related m6A regulators. Genes in the subnetwork were enriched in certain hormone signaling pathways, including thyroid hormone and oxytocin (OT). Recent evidence indicates that perturbations of the thyroid hormone signaling pathway and antibodies are associated with GDM development and mal-outcome [49, 51]. Emerging evidence suggests a role of OT in the pathogenesis of insulin resistance, obesity, and dyslipidemia [52, 53], which are hallmarks or risk factors for GDM. A significant weight loss and improvement of insulin sensitivity, pancreatic β-cell responsivity, and lipid metabolism can be observed in rodents, nonhuman primates, and humans after chronic subcutaneous or intranasal OT treatment [54, 55]. LINC00667 exhibited the highest node degree in the ceRNA network. Moreover, the genes in the LINC00667-mediated ceRNA subnetwork were functionally related to GDM by participating in the thyroid hormone signaling pathway. Consequently, LINC00667 may notably contribute to the development of GDM. A total of 47 m6A modification sites in LINC00667 were predicted using RMBase v2.0 [35]. Pearson correlation analysis suggested that LINC00667 was positively correlated with YTHDF3 (PCC = 0.95), an m6A “reader”. Moreover, RIP assays showed that YTHDF3 directly bind LINC00667 in HTR8/SVneo cells. We further found that MYC possessed the highest node degree in the PPI network. Existing literature indicates that both LINC00667 and MYC serve as targets of miR-33a-5p [36-41]. Through qPCR validation in 64 placental samples, we demonstrated that LINC00667, YTHDF3, MYC levels were significantly higher in the GDM group than in the NGT, while miR-33a-5p was significantly decreased. Furthermore, a SVM classifier for GDM was applied based on a module composed of LINC00667, YTHDF3, MYC, and miR-33a-5p. The expression matrix of these four genes was extracted based on data from 63 samples. Through the training and testing, this module showed good classifying power for GDM. Therefore, our findings suggested that this m6A-related module can be regarded as containing pivotal targets for research on the mechanism of GDM.

This study had a few limitations. A comprehensive analysis of the placenta and peripheral blood is warranted to verify the RNA expression, protein expression and m6A-modification status of the hub genes. The diagnostic ability of this m6A-related module may require further validation using a larger sample size. For subsequent research, more experiments in vitro regarding the biological functions of YTHDF3/LINC00667/miR-33a-5p/MYC axis should be incorporated to shed light into the mechanism of GDM.

Conclusion

In conclusion, this study was focused on potential m6A-related lncRNAs and crosstalk with the lncRNA-mediated ceRNA network that affects GDM. Based on bioinformatics analysis, we identified an m6A-related module consisting of LINC00667, YTHDF3, MYC, and miR-33a-5p, which was used successfully to classify GDM and NGT samples. A comprehensive analysis of the placenta and peripheral blood using a larger sample size, combined with validation of intermolecular interactions in vivo and in vitro, is warranted to validate the accuracy and reliability of our findings in the future.

Additional file 1: Table S1. Primer sequences.

Additional file 2: Table S2. List of all differentially expressed mRNAs (DEMs).

Additional file 3: Table S3. Clinical characteristics of patients included in the study.

Additional file 4. Supplementary file 1.