Identification of Key Molecules and lncRNA-miRNA-mRNA ceRNA Network in Preeclampsia

Introduction

Preeclampsia (PE) is an idiopathic hypertensive disorder of pregnancy, which refers to the occurrence of hypertension and proteinuria occurring in patients after the 20th week of gestation.¹ The incidence of PE is reported to be 4%-5% and is the leading cause of maternal death, fetal malformation, fetal intrauterine growth restriction, and premature birth.² In addition, PE may cause long-term cardiovascular and metabolic diseases in the mother.³ Nevertheless, the etiology and pathogenesis of PE have not been completely elucidated. Therefore, exploring key genes for PE is essential to revealing the underlying molecular mechanisms and investigating possible diagnostic and therapeutic targets.

Long non-coding RNAs (lncRNAs), a subclass of non-coding RNAs (ncRNAs) with a length of more than 200 nucleotides, play essential roles in regulating various biological functions.^3,4 LncRNAs may act as competing endogenous RNAs (ceRNAs) to sponge miRNAs, thus influencing the regulatory function of miRNAs and imposing posttranscriptional regulation of gene expression.⁵ Recently, the abnormal expression of lncRNAs has been reported to contribute to the pathogenesis of PE. For example, Xu et al have demonstrated that lncRNA TUG1 was significantly reduced in the placental tissue of PE patients, and knockdown of TUG1 affected cell proliferation, apoptosis, and migration.⁶ A previous study by Li et al uncovered that lncRNA CCAT1 was overexpressed in PE patients and promoted the progression of PE via suppressing CDK4 expression.⁷ The lncRNA-H19 increased in PE patients and promoted invasion and autophagy via activating PI3K/AKT/mTOR pathways in trophoblast cells.⁸

miRNAs are a class of endogenous single-stranded ncRNAs with a length of about 22 nucleotides. Accumulating evidence has shown that miRNAs are involved in PE. Qian and Liu indicated that miR-30b may facilitate the development of PE via MAPK pathway by targeting matrix-remodeling associated 5 (MXRA5).⁹ Moreover, miR-134 suppressed the invasion of trophoblast cells in the placenta of PE patients by decreasing ITGB1 expression.¹⁰ Besides, miR-200a, miR-141, miR-144 could affect the proliferation and invasion of trophoblast cells by targeting VEGF, PTEN, and other angiogenic factors, thus participating in the maintenance of pregnancy and the occurrence and development of PE.^11–13 These results demonstrated that lncRNAs and miRNAs played a pivotal role in the development and occurrence of PE. However, previous studies mainly focused on the effect of a single lncRNA, miRNA, and gene on PE development, and there were few reports on the lncRNA-miRNA-mRNA regulatory network in PE.

In the present study, we aimed to identify the key genes and construct a pivotal lncRNA-related ceRNA regulatory network in PE. We analyzed lncRNA, miRNA and mRNA expression profiles in the PE and healthy controls based on microarray technology. Then, GO and KEGG enrichment analyses were conducted to examine the biological function of differentially expressed genes (DEGs) and lncRNAs (DELs). This finding may provide a new target for the diagnosis and treatment of PE.

Materials and Methods

Participants

Placental tissues from PE patients (n=3) and healthy controls (n=3) were collected from the Second Affiliated Hospital of Hainan Medical College. Placental samples were taken from a representative block of the central portion of tissue below one-third of the placenta, near the maternal side. All samples were taken from the same region. This region was enriched for cytotrophoblast, syncytiotrophoblast, and villous interstitium. Samples were stored at −80°C until use. The exclusion criteria included multiple gestations, gestational diabetes mellitus, chronic hypertension, thyroid dysfunctions, and kidney disease. According to the American College of Obstetricians and Gynecologists (ACOG), PE is defined as systolic blood pressure of ≥ 140 mmHg and diastolic blood pressure of ≥ 90 mmHg with proteinuria of ≥ 300 mg/day (or a protein/ creatinine ratio of ≥ 0.3 mg/dl or proteinuria of ≥ 1+) or without proteinuria but with severe clinical features after 20 weeks of gestation in a woman with previously normal blood pressure.¹⁴

The research protocol complied with the Declaration of Helsinki and was approved by the Ethics Committee of the Second Affiliated Hospital of Hainan Medical College (NO.2018R005-F01). All samples were collected with patients’ informed consent.

RNA Extraction and Amplification

Total RNA of placental tissue was extracted by TRIzol Reagent (Thermo Fisher Scientific, USA). RNA concentration was measured using the Qubit RNA Assay Kit in the Qubit 3.0 Spectrophotometer (Thermo Fisher Scientific, USA). RNA purity was checked using the NanoDrop 2000 (Thermo Fisher Scientific, USA). The RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent 2100 Bioanalyzer system (Agilent Technologies, USA). RNA was amplified and labeled by a Low Input Quick Amp WT Labeling Kit (Agilent Technologies, USA) following the manufacturer’s instructions. Labeled cRNA was purified by the RNeasy mini kit.

Microarray Analysis

To prepare RNA samples, 3μg RNA per sample was used. First of all, ribosomal RNA was eliminated using the Epicentre Ribo-zero rRNA Removal Kit (Epicentre, USA). Next, ethanol precipitation was done to clean up rRNA-free residue. Successively, sequencing libraries were created using NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB) as per the supplier’s instructions. Concisely, fragmentation was done by divalent cations under raised temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First-strand cDNA was created by random hexamer primer as well as M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA was successively created by DNA Polymerase I plus RNase H. In the reaction buffer, dNTPs with dTTP were substituted by dUTP. Outstanding extensions were transformed into blunted ends by exonuclease/polymerase actions. Following the adenylation of 3ʹ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to get ready for hybridization. With the intention of the selection of cDNA fragments of favorably around 150~200 bp in length, the library fragments were cleaned with the AMPure XP system (Beckman Coulter). At that stage, 3μL USER Enzyme (NEB) was used with appropriate adaptor-ligated cDNA at 37°C for 15 minutes trailed by 5 minutes at 95°C. Next, PCR was done with Phusion High-Fidelity DNA. The clustering of the index-coded samples was accomplished on a cBot Cluster Generation System by TruSeq PE Cluster Kit v3-cBot-HS (Illumia) as per the supplier’s protocol. After cluster generation, the libraries were sequenced on an Illumina Hiseq. 2500 platform.

Quality Control of Raw Sequencing

For lncRNA, raw data was obtained by removing rRNA through an Epicentre ribo-0rRNA removal kit (Epicentre, USA) and removing free rRNA residue through ethanol precipitation. Subsequently, sequencing libraries were generated by the rRNA-depleted RNA by NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations. For miRNAs, raw data of fastq format was first processed through custom Perl and Python scripts. In this step, clean data was obtained by removing reads containing ploy-N, with 5ʹ adapter contaminants, without 3ʹ adapter or the insert tag, reads containing ploy-A, ploy-T, ploy-G, or ploy-C, and low-quality reads from raw data. At the same time, Q20, Q30, and GC-content of all the raw data were calculated. Taken together, all the downstream analyses were based on clean data.

Differentially Expressed Analysis

The Differentially expressed mRNA (DEGs), lncRNAs (DELs), and miRNAs (DEMs) were identified between PE and control using the Deseq R package. A threshold value of |log₂FC| >1 with p value <0.05 was determined. The ggplot2 in the R package was used to draw heat maps of the DERNAs.

Functional Enrichment Analysis

We performed Geno ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) on DEGs and DELs-target genes to explore their functions using Cluster Profiler in R package. GO analysis included biological processes (BP), cellular components (CC), and molecular functions (MF). KEGG analysis mainly focused on enrichment analysis of the signaling pathways of DEGs and DELs-target genes. P < 0.05 as the threshold for statistical significance. The ggplot2 in R package was used to draw the GO and KEGG analysis results.

CeRNA Regulatory Network Construction

The LncTar software was used to predict the lncRNA-miRNA interactions and then combined with the identified DEMs to build the DELs-DEMs network. Similarly, we utilized miRanda and RNAhybrid to predict the target genes of the selected DEMs in this study. And a miRNA-mRNA network was constructed by combining DEGs. By integrating the lncRNA-miRNA, miRNA-mRNA pairs, the lncRNA-miRNA-mRNA ceRNA network was established and visualized by Cytoscape software.

qRT-PCR for the Expression Examination

Based on the results of lncRNA-miRNA-mRNA ceRNA network, LINCO2532, SLCO4A1-AS1, miR23a-5p, and DYNLRB1 were selected for expression verification. Placental tissues from 27 PE patients (30.04 ± 4.37 years) and 30 healthy controls (29.97 ± 4.26 years) were collected from the Second Affiliated Hospital of Hainan Medical College for qRT-PCR assay. Total RNA of placental tissue was extracted by TRIzol Reagent (Thermo Fisher Scientific, USA). The reverse transcription (RT) reactions for the expression of LINCO2532, SLCO4A1-AS1, DYNLRB1 and GAPDH were performed by Takara-PrimeScript™ RT Master Mix (Perfect Real Time). SYBR Green-based qPCR profiling was performed using Takara-TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus). The primers sequences were designed as following, forward primer: LINC02532-F: GCAGATCAAGATGTGGCAACC, LINC02532-R:TGCCCTCAGCTCAGAGTTCT, SLCO4A1-AS1-F: AGACATTGGGGCTGAACGTG, SLCO4A1-AS1-R:GTTCAAGTCAGCGACACTCTG, DYNLRB1-F: GACCCCCAGAACGATCTCAC, DYNLRB1-R: AGGAAATAGTCTTTATCTGGTGCAA, GAPDH –F: GGAGCGAGATCCCTCCAAAAT, and reverse primer: GAPDH –R: GGCTGTTGTCATACTTCTCATGG for GAPDH. The levels of miRNAs were normalized using GAPDH as reference RNA. The relative expression quantity (RQ) of mRNA was calculated as RQ = 2^−ΔΔCt. ΔΔCt = mean value of the study group (Ct _target − Ct_GAPDH) – mean value of the control group (Ct _target− Ct_GAPDH).

The reverse transcription (RT) reactions of U6 and miR23a-5p were performed by the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (cat. no. KR211; Tiangen) with Oligo (dT)-Universal Tag primer. SYBR Green-based qPCR profiling was performed using miRcute Plus miRNA qPCR Detection Kit (cat. no. FP411; Tiangen) in ABI 7500 Real-Time PCR System (Applied Biosystems) according to the manufacturer’s instructions. The miR23a-5p specific forward primer sequences (GCGCAGTTCCTGGGGATG) were designed based on the miRNA sequences obtained from the miRBase database. The levels of miRNAs were normalized using U6 as reference RNA. The relative expression quantity (RQ) of miRNA was calculated as RQ = 2^−ΔΔCt. ΔΔCt = mean value of the study group (Ct _miR23a-5p − Ct _U6) – mean value of the control group (Ct _miR23a-5p − Ct _U6).

Results

Subjects

Three patients with PE and three healthy controls were enrolled in the study. The average age of PE and control was 31.33 ± 4.99 years and 28.67 ± 4.99 years, respectively. The clinical features of subjects are shown in Table 1.

Table 1 Clinical Characteristics of PE and Control for Microarray

Differential Expression Analysis

Based on the given threshold with p < 0.05 and |logFC| > 1), we screened 726 DELs (456 upregulated and 370 downregulated), 49 DEMs (37 upregulated and 12 downregulated), and 318 DEGs (230 upregulated and 88 downregulated) between PE patients and control. The heat map indicated that these DELs (Figure 1A), DEMs (Figure 1B), and DEGs (Figure 1C) distinguished the PE from the control samples. In addition, the top 10 DELs, DEMs, and DEGs are exhibited in Table 2.

Table 2 Top 10 Upregulated and Downregulated Differentially Expressed lncRNAs, miRNAs and mRNAs

Figure 1 The expressions of differentially expressed lncRNAs, miRNAs and mRNAs. Heat maps of (A) lncRNA, (B) miRNA, and (C) mRNA (|log2FC| >1 and p<0.05). Red and green indicate up- and downregulation, respectively. Abbreviations: P, preeclampsia; C, control.

Functional Enrichment Analysis

To better understand the major biological function of DEGs and DELs-target genes in PE, we conducted GO enrichment and KEGG pathway analysis. The enriched GO terms of BP, CC, and MF are exhibited in Figure 2A and Table 3. Our results revealed that DEGs were mainly enriched in leukocyte migration, cytosolic ribosome, and cytokine activity. The KEGG analysis enriched these DEGs in pathways of cytokine-cytokine receptor interaction, and metabolism of xenobiotics by cytochrome P450 (Figure 2B and Table 4).

Table 3 Top 10 Enriched GO Terms of the DEGs

Table 4 Top 10 Enriched KEGG Pathways of the DEGs

Figure 2 GO term and KEGG pathway analyses of DEGs. (A) Top 10 GO terms enrichment of DEGs. (B) Top 10 KEGG pathway enrichment of DEGs. Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed gene.

Besides, the GO analysis revealed that DELs-target genes were mainly enriched in negative regulation of the cellular catabolic process, cytosolic part, and cofactor binding (Figure 3A, Table 5). The KEGG pathway analysis indicated that DELs-target genes were significantly enriched in RNA transport (Figure 3B, Table 6).

Table 5 Top 10 Enriched GO Terms of the DELs-Target Genes

Table 6 Top 10 Enriched KEGG Pathways of the DELs-Target Genes

Figure 3 GO term and KEGG pathway analyses of DELs-target genes. (A) Top 10 GO terms enrichment of DELs-target genes. (B) Top 10 KEGG pathway enrichment of DELs-target genes. Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DELs: differentially expressed lncRNAs.

Construction of ceRNA Network

A total of 906 lncRNA-miRNA interactions were identified by lncTar software. Combining with DELs and DEMs, 138 lncRNA-miRNA pairs were selected, including 79 DELs, and 15 DEMs (Figure 4A). The 15 DEMs identified in the above step were predicted their target genes using miRanda and RNAhybrid databases. Subsequently, these target genes were integrated with DEGs. Finally, target genes of 1 DEMs (hsa-miR-23a-5p, Figure 4B) were determined.

Figure 4 lncRNA-miRNA and miRNA-mRNA interaction network. (A) lncRNA-miRNA interaction network. (B) miRNA-mRNA interaction network. In this figure, lncRNA, miRNA, and mRNA were indicated by ellipse, diamond, and triangle, respectively.

Based on “ceRNA hypothesis”, we established lncRNA-miRNA-mRNA regulatory network. The ceRNA network contains 16 lncRNA (SLCO4A1-AS1, lnc-TMEM41B-5, LINC02532, lnc-ZNF487-4, lnc-ZNF791-1, lnc-MFHAS1-1, lnc-MICALL2-2, lnc-RNF208-2, lnc-FAM71E1-2, lnc-NUB1-2, lnc-KIN-9, C8orf37-AS1, lnc-OR4F21-4, lnc-SCRN1-4, lnc-KLHDC7B-3, lnc-KCNE2-2), 1 miRNA (miR-23a-5p), and 1 mRNA (DYNLRB1). The ceRNA network is presented in Figure 5.

Figure 5 Competing endogenous RNA interaction network of lncRNA-miRNA-mRNA. Diamond represents miRNAs, ellipse represents lncRNAs, and triangle represents mRNAs.

qRT-PCR for the Expression Examination

Based on Figure 5, LINCO2532, SLCO4A1-AS1, miR23a-5p, and DYNLRB1 were selected for expression verification, as shown in Figure 6. We found that lncRNA (LINCO2532 and SLCO4A1-AS1) and DYNLRB1 were downregulated in PE patients compared with healthy controls (p = 0.005, p = 0.013, and p < 0.001, respectively). The expression of miR23a-5p (p < 0.001) was upregulated in PE patients compared with healthy controls.

Figure 6 The expression of LINCO2532 (A), SLCO4A1-AS1 (B), miR23a-5p (C), and DYNLRB1 (D) in PE patients and healthy controls. Statistical significance of expression level with *for p < 0.05, **for p < 0.01 and ***for p < 0.001.

Discussion

It is reported that about 3–5% of pregnant women worldwide suffer from PE. Moreover, PE can cause damage to the kidney and liver function of pregnant women and severely can cause the death of pregnant women and fetuses. In recent years, the role of lncRNA and miRNA in the development of PE has received extensive attention.^15,16 However, there were few reports on lncRNA-mediated regulatory networks in PE, as well as reports of new transcripts. To identify the new transcripts and biological functions of lncRNA and miRNA in PE, we performed RNA-sequencing of placental samples from PE patients and healthy controls and elucidated the underlying pathogenesis of PE.

In the present study, we identified 726 DELs, 49 DEMs and 318 DEGs based on RNA-sequencing data from placental tissue. To further explore the regulatory role of DELs in PE, the interaction network of lncRNA, miRNA and mRNA was constructed. This network showed that multiple lncRNAs are associated with a single miRNA (miR-23a-5p), and then, it regulates a single mRNA (DYNLRB1). And miR-23a-5p and DYNLRB1 were determined to be involved in various diseases, including renal cell carcinoma, diabetic nephropathy.^17,18 These results suggested that miR-23a-5p and DYNLRB1 may play an important role in the process of PE.

The results of qRT-PCR displayed that lncRNA (LINCO2532 and SLCO4A1-AS1) and DYNLRB1 were downregulated and the expression of miR23a-5p was upregulated in PE patients compared with healthy controls. Next, the role of miR-23a-5p and DYNLRB1 in PE need to be further investigated in vivo and in vitro.

ceRNA regulatory network suggests that RNA with miRNA response element (MRE) can competitively bind miRNAs and achieve mutual regulation. Recently, the role of ceRNA regulatory network in various diseases has been widely studied. The study of Zhou et al indicated that an intrahepatic cholangiocarcinoma ceRNA network was built with 30 miRNAs, 16 lncRNAs, and 80 genes, which lncRNA ARHGEF26-AS1 and MIAT were hub genes.¹⁹ Sun et al identified 3 lncRNAs, 15 miRNAs, and 138 mRNAs via constructing lncRNA-miRNA-mRNA network as potential markers in the pathogenesis of inflammatory bowel disease.²⁰ Leng et al established a ceRNA network to reveal the role of lncRNA HCG27 in gestational diabetes mellitus.²¹ At present, some studies have reported the role of ceRNA regulatory network in PE. For example, Liu et al established a ceRNA network related to PE, showing that LINC01492, hsa_circ_0088196 and LIF promote PE development through suppressing miR-100-5p.²² Besides, Zhang et al also built a ceRNA regulatory network consisting of 3 genes, 21 lncRNAs and 69 miRNAs in early-onset preeclampsia.²³ Moreover, Hu and his colleagues revealed the global ceRNA network in PE and has_circ_0036877 may serve as a potential novel blood biomarker for early PE.²⁴ However, the etiology and pathogenesis of PE are unclear. To further elucidate the pathogenesis of PE, a new regulatory network of ceRNA was constructed based on transcriptome sequencing data from placenta samples.

In this study, one hub mRNA (dynein light chain roadblock-type 1, DYNLRB1) in the ceRNA network was identified. DYNLRB1 is a member of the roadblock dynein light chain family and is expressed in the placenta.²⁵ It interacts with transforming growth factors and plays a pivotal role in TGF-β signaling transduction.²⁶ Moreover, increasing evidence has highlighted the importance of TGF-β signaling in PE, placenta, and umbilical cord development.^27–29 These results suggested that DYNLRB1 may act in PE through TGF-β signaling.

Moreover, we found DYNLRB1 was targeted by miR-23a-5p in the ceRNA network. miR-23a-5p has been reported to be involved in cell apoptosis and growth.^30,31 It is reported that miR-23a was increased in the plasma of pregnant women affected by diabetes mellitus.³² Moreover, Li et al uncovered that miR-23a was upregulated in placental samples of PE patients and it promoted trophoblast apoptosis by inhibiting XIAP, ultimately facilitated the development of PE.³³ Besides, some researchers reported that the expression of miR-23a-3p showed an upward tendency in women experiencing recurrent miscarriage, preterm birth, and cervical shortening.^34,35

The above molecules and their predicted interactions might unveil the mechanisms of PE and provide novel biomarkers for the identification of PE. Despite these interesting results, there are some limitations in this study. One is that the sample size is small. An additional validation cohort should be included in further studies to analyze the expression of these identified lncRNAs, miRNAs and mRNA. Second, this is a preliminary screening study and further experiments are needed to confirm the interactions in the identified ceRNA axes in PE.

Conclusions

In conclusion, a novel lncRNA-associated ceRNA network was constructed in PE patients. It elucidated the regulatory mechanism of PE, and identified novel PE biomarkers, which have important guiding significance for clinical treatment and further scientific research of PE.