Microarray Expression Profile Analysis of Long Non-Coding RNAs in Umbilical Cord Plasma Reveals their Potential Role in Gestational Diabetes-Induced Macrosomia

Foetal macrosomia is the predominant adverse outcome of gestational diabetes mellitus (GDM), which is usually defined as a birth weight (BW) above the 90th percentile or BW above 4000g [1,2]. Macrosomia is closely associated with higher rates of cesarean section and sequential use of vacuum and obstetric forceps for vaginal delivery for the mother, and clavicle fracture, intracranial hemorrhage and respiratory distress syndrome for the newborn [3,4]. Furthermore, macrosomic babies may have an increased susceptibility to obesity and type-2 diabetes and/or cardiovascular diseases in later life [5,6]. Recently, the increasing prevalence of GDM in pregnant women in many parts of the world could be associated with a parallel increase in macrosomic births. The continuously increasing incidence of gestational diabetes-induced macrosomia during the past three decades in developed and developing countries suggested that it is critical to understand the mechanisms underlying gestational diabetes-induced macrosomia, which will help to reduce the incidence of obstetrical complications and type 2 diabetes [3,7].

The placenta is an important endocrine organ during human pregnancy, which produces numerous proteins that may promote early embryonic growth [8]. Furthermore, release of these critical proteins from placental tissues into the umbilical circulation presumably contributes to their increased concentration in cord blood and thereby, fetal tissues. These proteins are reported to be involved in accelerating fetal growth [9,10,11]. Since the umbilical cord vein blood is primarily produced from placental sources, the study of the critical gene products in the umbilical cord vein blood might be help to understand the mechanisms of fetal macrosomia.

Recently, genetic studies have focused on non-coding RNAs. Over the past decade, many studies have reported that non-coding RNAs have important regulatory potentials, both transcriptionally and post transcriptionally [12]. Long non-coding RNAs (lncRNAs) are functional RNAs longer than 200 nucleotides in length. LncRNAs are associated with more than 150 diseases [13], such as preeclampsia [14]. Therefore, LncRNAs are important biological molecules in the mechanisms of disease and are useful for exploring biomarkers for disease diagnosis and treatment. However, the expression of lncRNAs in umbilical cord vein blood and their biological functions in gestational diabetes-induced macrosomia remain unknown.

In this study, we examined the umbilical cord vein blood lncRNA expression profiles of 30 cases of gestational diabetes-induced macrosomia compared with 30 matched non-diabetic macrosomia samples. We then performed internal validations using individual quantitative real-time polymerase chain reaction (qRT-PCR) assays to find a class of functional lncRNAs. Our results demonstrated that lncRNA expression profiles may provide new molecular biomarkers or a new basis for the diagnosis and treatment of gestational diabetes-induced macrosomia.

This study was conducted in accordance with the declaration of Helsinki and performed with approval from the Ethics Committee of the Nanjing Maternity and Child Health Care Hospital affiliated to Nanjing Medical University. Written informed consent was obtained from all subjects before enrolment.

Women with GDM, diagnosed by 75g oral glucose tolerance test (OGTT) during their 24th to 28th week of gestation at the Obstetrics Department in Nanjing Maternity and Child Health Care Hospital in 2013, were enrolled in this study. According to the IADPSG (The International Association of Diabetes and Pregnancy Study Groups) and ADA (American Diabetes Association) recommendation [15,16], GDM was diagnosed if any one glucose level was equal to or greater than fasting 5.1 mmol/L, one hour post load 10.0 mmol/L or two hours post load 8.5 mmol/L. Pregnant women with GDM were treated with a strict diet control to maintain the fasting plasma glucose level at less than 5.3mmol/L and 2h postprandial level at less than 6.7mmol/L. Insulin treatment was administrated when dietary control could not achieve this goal.

Clinical data and outcomes of GDM women and their paired controls were collected until delivery. The babies of GDM group were diagnosed as macrosomia when their birthweight was equal to or exceeded 4000 g [17]. The baby birthweight was less than 4000 g in the control group with normal OGTT. They were chosen on the basis of matched gestational age and gender, and to keep the maternal baseline characteristics similar. The exclusion criteria for both groups included hypertensive disorders, history of type 2 diabetes, smoking, chemical dependency, assisted reproductive technology treatment, multiple pregnancies, fetal congenital anomalies and any other confounding pathologies (including intrahepatic cholestasis of pregnancy, hyperthyroidism and hypothyroidism). Thirty control cases and 30 GDM cases were enrolled in this study ultimately. Maternal and neonatal characteristics are presented in Table 1.

Blood samples of the umbilical vein of the control and GDM groups were collected at the time of cesarean or vaginal delivery. Plasma was obtained by centrifugation of EDTA-blood at 1600 × g for 10 min, and supernatant was centrifuged at 12000 × g for 5 min at 4°C. The control or GDM group was divided into three subgroups randomly, such that every subgroup was a mixture of 10 mL plasma (1 mL from every pregnant woman). In total, 2 mL clean plasma from each subgroup was filtered through a 0.22-µm filter (MILLEXGV; Millipore), and then mixed with 2 mL TRIZOL reagent and 0.4 mL chloroform. The mixture was centrifuged at 12000 × g for 5 min at 4°C. The aqueous layer was transferred to a new tube, 1.5 mL 70% ethanol was added, and then applied to an RNeasy minicolumn (Qiagen), according to the manufacturer's recommendations. RNA quantity and quality were measured using a NanoDrop ND-1000. Total RNA was eluted with 20 µL RNase-free water after DNase digestion and then stored in liquid nitrogen at -20°C [18].

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology), with minor modifications. Briefly, mRNA was purified from total RNA after removal of rRNA (mRNA-ONLY™ Eukaryotic mRNA Isolation Kit, Epicentre). Each sample was then amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3' bias, utilizing a random priming method (Arraystar Flash RNA Labeling Kit, Arraystar). The labeled cRNAs were purified by an RNeasy Mini Kit (Qiagen). The concentration and specific activity of the labeled cRNAs (pmol Cy3/µg cRNA) were measured using a NanoDrop ND-1000. 1 µg of each labeled cRNA was fragmented by adding 5 µL 10 × Blocking Agent and 1 µL of 25 × Fragmentation Buffer, then heating at 60°C for 30 min; finally 25 µl 2 × GE Hybridization buffer was added to dilute the labeled cRNA. 50 µL of hybridization solution was dispensed into the gasket slide and assembled into the LncRNA expression microarray slide. The slides were incubated for 17 hours at 65°C in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed and scanned using an Agilent DNA Microarray Scanner (part number G2505C).

Arraystar Human LncRNA Microarray V3.0 is designed for the global profiling of human LncRNAs and protein-coding transcripts. About 30,586 LncRNAs and 26,109 coding transcripts can be detected by our third-generation LncRNA microarray. 33,045 LncRNAs and 30,215 coding transcripts collected from the most authoritative databases such as, UCSC, NCBI RefSeq, RNAdb, Ensembl, H-invDB and many related literatures, could be detected using the microarray. The microarray work was performed by Kang Chen Biotech, Shanghai, P.R. China [19]. The criteria were as follows: fold-change cut-off: 2.0; fold-change: positive value indicates upregulation in GDM group and negative value indicates downregulation in the control group. Log fold-change means a log₂ value of the absolute fold-change. The fold-change and p-value were calculated from the normalized expression.

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies). After quantile normalization of the raw data, LncRNAs and mRNAs that in at least three out of six samples had flags in Present or Marginal (“All Targets Value”) were chosen for further data analysis. Differentially expressed LncRNAs and mRNAs with statistical significance between the two groups were identified through P-value/false discovery rate (FDR) and fold change filtering. Hierarchical Clustering and combined analysis were performed using in-house scripts.

The total RNA was isolated using the TRIZOL reagent (Invitrogen, Carlsbad, CA, USA) and then reverse transcribed using PrimeScript® RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian China), according to the manufacture's instructions. The expressions of five upregulated lncRNAs and three downregulated lncRNAs in 60 pregnant women enrolled in this study were measured by qRT-PCR and GAPDH was used as an internal control. The primers are listed in Table 2. The relative expression of LncRNAs was determined by the 2^-∆∆^Ct method with GAPDH expression to normalize the data. The expression level of lncRNAs differentially expressed between normal control and GDM women complicated with macrosomia were analyzed using Student's t-test in SPSS (Version 17.0 SPSS Inc). P < 0.05 was considered statistically significant [20]. The clustering analysis of differentially expressed lncRNAs was performed using the software Cluster 3.0, and the heatmap of the clustered results was generated using the software TreeView.

The lncRNA expression profiling data showed that 8814 lncRNAs were expressed in umbilical cord venous plasma from women with gestational diabetes-induced macrosomia (Fig. 1). The lncRNA expression levels between 30 women with gestational diabetes-induced macrosomia and their paired control were compared. 349 lncRNAs were significantly upregulated and 892 lncRNAs were significantly downregulated (fold-change ≥ 2.0) in the GDM group. Hierarchical clustering analysis was used to analyze the relationships between differentially expressed lncRNAs and gestational diabetes-induced macrosomia. XLOC_003497 (Fold change = 42.2945021) was the most significantly upregulated lncRNA and XLOC_006112 (Fold change = -68.7441146) was the most significantly downregulated lncRNA (Table 3).

The mRNA expression profiling data showed that 11668 mRNAs were expressed in umbilical cord venous plasma from women with gestational diabetes-induced macrosomia (Fig. 2). The mRNA expression levels between two groups were compared. 252 mRNAs were significantly upregulated in the GDM group and 654 mRNAs were significantly downregulated (fold-change ≥ 2.0). Hierarchical clustering analysis was used to analyze the relationships between differentially expressed mRNAs and gestational diabetes-induced macrosomia. PSMC1 (Fold change = 53.7099895) was the most significantly upregulated mRNA and NCOR1 (Fold change = 181.5956457) was the most significantly downregulated mRNA (Table 4).

Gene Ontology (GO) analysis was performed to determine the gene and gene product attributes in biological process, cellular components and molecular functions for the transcripts derived from genes encoding or associated with the differentially expressed lncRNAs. Fisher's exact test was performed to determine if there was more overlap between the differentially expressed list and the GO annotation list than would be expected by chance. The p-value ≤ 0.05 denoted the significance of GO terms enrichment in the DE genes. The most enriched GOs targeted by the downregulated transcripts were biological regulation (Fig. 3A), cytoplasm (Fig. 3B) and binding (Fig. 3C). The most enriched GOs targeted by the upregulated transcripts were positive regulation of cellular process (Fig. 3D), organelle part (Fig. 3F) and protein kinase binding (Fig. 3E).

Pathway analysis indicated that six pathways corresponded to upregulated transcripts (Fig. 4B). The most enriched network was “β cell receptor signaling pathway” comprising four target genes. Pathway analysis indicated that nine pathways corresponded to downregulated transcripts (Fig. 4A). The most enriched network was “mineral absorption” comprising seven target genes. Among these pathways, the gene category “mineral absorption” has been reported to be involved in fetal growth [21,22]. The gene category “chemokine signaling pathway” has been reported to be involved in insulin resistance and metabolic phenotype variation [23]. The gene category “Ras signaling pathway” (Fig. 4C) has been reported to be involved in insulin resistance and transplacental transportation of glucose, thus affecting fetal intrauterine growth [24].

Five upregulated lncRNAs and three downregulated lncRNAs (fold changes ≥ 10), whose potential target genes could be closely related with growth, glucose and lipid metabolism by bioinformatics analysis, were chosen for validation. QPCR was performed to identify the expression levels of the differentially expressed lncRNAs in GDM women with macrosomia, using all 60 plasma samples from women enrolled in this study. The results suggested that the qPCR results were consistent with the microarray data (Fig. 5).

GDM, first diagnosed in pregnancy, complicates about 5% to 10% of pregnancies [25]. Macrosomia is the predominant adverse outcome in cases of GDM. It is the main factor linked to reported cases of other complications in GDM [4]. Different from non-diabetic macrosomia, gestational diabetes-induced macrosomia is characterized by asymmetric overgrowth and size increases of insulin-sensitive tissues and organs, especially the adipose tissue, skeletal muscle and liver [26,27]. Previous studies suggested that GDM women, characterized by hyperglycemia, expose the fetus to higher concentrations of glucose than normal women, forcing the fetus to increase its own insulin production, thus leading to macrosomia [28]. However, there is increasing evidence that, besides the hyperglycemic factor, the releasing of certain critical proteins or nucleic acids from the placenta into the umbilical circulation contributes to their increased concentration in cord blood and thereby fetal tissues, and these factors are reported to be involved in accelerating fetal growth [9,10,11,29].

LncRNAs are an emerging class of RNA molecules with a length of more than 200 bases. Studies have revealed that LncRNAs play an essential role in gene silencing, as scaffolds for histone modification, transregulatory elements affecting multiple gene transcription, enhancers to modify gene expression, participators in genomic reprogramming involved in stem cell differentiation and as nucleators to generate the dynamic assembly of the nuclear structure [30]. Recently, some studies have analyzed global changes in LncRNA expression in the pathogenesis of cancer [31], development [32], adipogenesis [33] and preeclampsia [15].

To date, there have been no reports describing LncRNA expression in maternal and fetal circulation system. In our study, we identified the global expression changes in lncRNAs in umbilical vein plasma from GDM-induced macrosomia and paired controls. Since GDM-induced macrosomia is closely related with various factors, such as history of diabetes, BMI before pregnancy, age and genetic reasons, we excluded these confounding factors in the subject enrollment process. Our aim was to study the placenta-derived factors that were released into the umbilical cord vein blood and accelerated fetal growth. Using abundant and varied probes accounting for 33,045 lncRNAs and 30,215 coding transcripts in a microarray, we identified 8814 expressed lncRNAs, among which 259 were upregulated and 892 were downregulated. All raw data discussed in this paper have been deposited in NCBI Gene Expression Omnibus and are accessible with the GEO Series accession number GSE65737 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65737). Most of these lncRNAs have not been functionally characterized. Five upregulated lncRNAs and three downregulated lncRNAs were further validated using qRT-PCR, and the results were consistent with those of the microarray. Moreover, GO analysis and pathway analysis were performed and a coexpression network was constructed to study preliminarily the biological functions and potential mechanisms of these differentially expressed lncRNAs in umbilical vein plasma of GDM.

LncRNA ENST00000552367 is an 885 bp lncRNA transcript from the Pleckstrin homology-like domain, family A (PHLDA1) gene, located on Chromosome 12. PHLDA1, also known as T cell death-associated gene 51 (TADG51), plays an important role in the anti-apoptotic effects of insulin-like growth factor-1, which is an important regulator in cell proliferation [34] and gestational diabetes-induced macrosomia [35,36]. Recently, Basseri [37] suggested that PHLDA1 (TDAG51) expression was inversely correlated with fatty liver in multiple mouse models of hepatic steatosis. TDAG51(-/-) mice fed a chow diet exhibited greater body and adipose deposition, had reduced energy expenditure, displayed mature-onset insulin resistance, and was predisposed to hepatic steatosis. Taken together, TDAG51 is involved in energy homeostasis, especially in regulating lipogenesis in liver and adipose tissue, and may constitute a novel therapeutic target for the treatment of obesity and insulin resistance. The fold change of upregulated LncRNA ENST00000552367 was 14.1039315. Coincidentally, PHLDA1 (TDAG51), the target gene of LncRNA ENST00000552367 was also identified in the significantly downregulated mRNAs list in this study (fold change = 3.7422605, p value = 1.65074E-06). This evidence could shed light on the potential roles of LncRNA ENST00000552367 and PHLDA1 (TDAG51) in gestational diabetes-induced macrosomia. LncRNA uc021tzj.1 is a 14,121 bp lncRNA transcript from the HOXB-AS3/HOXB3/HOXB5/HOXB6 gene, located on Chromosome 17. HOXB-AS3, HOXB3, HOXB5 and HOXB6 are members of homeobox transcription factors family, which are important in differentiation phenotypes of pancreatic beta-cells, insulin resistance [38], and for regulating adipose tissue formation and loss [39]. Furthermore, HOXB3, the target gene of LncRNA uc021tzj.1, was also identified among the differentially expressed mRNAs list (fold change = 2.4270548, p value = 2.50138E-06). LncRNA RP11-316M1.12 is a 999 bp lncRNA transcript from the ceramide synthase 2 (CERS2) gene, located on Chromosome 1. Raichur [37] demonstrated that insufficiency of CERS2 led to hepatic steatosis and insulin resistance. CERS2 could become a therapeutic target for treating metabolic diseases associated with obesity. LncRNA XLOC_009474 is a 11,954 bp lncRNA transcript from the metastasis associated lung adenocarcinoma transcript 1 (MALAT1) gene, located on Chromosome 1. Liu [40] showed that MALAT1, also known as noncoding nuclear-enriched abundant transcript 2 (NEAT2), could regulate cell proliferation and migration. MALAT1 upregulation represents a critical pathogenic mechanism for diabetes-induced microvascular dysfunction. Cooper [41] suggested that MALAT1 associates with SRp40 to regulate splicing of PPARγ2 (peroxisome proliferator-activated receptor γ2) during adipogenesis in 3T3-L1 Cells.

To understand the functions of the lncRNAs and the mechanism by which they exert their roles, pathway analysis were performed to study the differentially expressed lncRNAs. Nine pathways corresponded to downregulated transcripts and six pathways corresponded to upregulated transcripts. Among these pathways, the gene category “mineral absorption” has been reported to be involved in fetal growth [21,22]. The gene category “chemokine signaling pathway” has been reported to be involved in insulin resistance and metabolic phenotype variation [23]. The gene category “Ras signaling pathway” has been reported to be involved in the insulin resistance and transplacental transportation of glucose, thus affecting fetal intrauterine growth [24].

This study has highlighted the critical role of the materno-feto-placental communication involved in the fetal weight gain in macrosomic babies. The expression profiles of umbilical cord vein plasma, and GO analysis, signal pathway analysis and lncRNA classification analysis, will help us to understand the potential mechanisms of the development of gestational diabetes-induced macrosomia. Furthermore, this study focused, for the first time, on the molecular mechanisms of gestational diabetes-induced macrosomia at the lncRNA level. In addition, we showed that lncRNA ENST00000552367, which is associated with PHLDA1, plays an important role in the anti-apoptotic effects of insulin-like growth factor-1, which is an important regulator in cell proliferation and gestational diabetes-induced macrosomia. It is suggested that lncRNAs may exert their functions by regulating the transcription of their related protein-coding genes in gestational diabetes-induced macrosomia. However, further work, including functional and mechanistic experiments, is needed to elucidate the potential roles of umbilical cord vein lncRNAs in gestational diabetes-induced macrosomia.

This work was supported by grants from the National Natural Science Foundation of China [grant numbers 81100436], the State Key Laboratory of Reproductive Medicine Fund [grant number SKLRM-KF-201109], the Nanjing Medical Technology Development Project [grant numbers YKK14126] and Nanjing Medical Young Investigators Project [grant number QRX11211].

All authors have no conflicts of interest to declare.

Z. Shi and C. Zhao contributed equally to this work.