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

BACKGROUND: Preeclampsia (PE) is an idiopathic hypertensive disorder of pregnancy and is the leading cause of maternal death, fetal malformation, and premature birth. The purpose of this study is to identify the key molecules and lncRNA-related competitive endogenous (ceRNA) regulatory network in PE.
METHODS: The differentially expressed mRNAs (DEGs), lncRNAs (DELs), and miRNAs (DEMs) were identified between PE and control using the Deseq R package. In addition, we performed Geno ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) on DEGs and DELs-target genes to explore their function. The ceRNA network was established based on lncRNA-miRNA and miRNA-mRNA interactions and visualized using Cytoscape software. LINCO2532, SLCO4A1-AS1, miR23a-5p, and DYNLRB1 were selected for qRT-PCR assay.
RESULTS: Using microarray analysis, 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. Based on lncRNA-miRNA pairs and miRNA-mRNA pairs, the ceRNA network was constructed, which contained 16 lncRNA, 1 miRNA (miR-23a-5p), and 1 mRNA (DYNLRB1). LncRNA (LINCO2532 and SLCO4A1-AS1) and DYNLRB1 were downregulated and the expression of miR23a-5p was upregulated in PE patients compared with healthy controls.
CONCLUSION: In this study, the novel ceRNA network was established in the placentas of 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.

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.1 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.2 In addition, PE may cause long-term cardiovascular and metabolic diseases in the mother.3 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.5 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.6 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.7 The lncRNA-H19 increased in PE patients and promoted invasion and autophagy via activating PI3K/AKT/mTOR pathways in trophoblast cells.8

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).9 Moreover, miR-134 suppressed the invasion of trophoblast cells in the placenta of PE patients by decreasing ITGB1 expression.10 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.14

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 |log2FC| >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 − CtGAPDH) – mean value of the control group (Ct target− CtGAPDH).

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

	Age(year)	Week Gestation (Week)	Systolic BP (mmHg)	Diastolic BP (mmHg)	Proteinuria (g/24h)	Birth Weight (g)
C1	34	38+2	130	84	–	3550
C2	30	38	126	72	–	3200
C3	22	38+6	120	76	–	2900
P1	30	37+5	150	98	392	3200
P2	26	38+1	158	95	1320	2350
P3	38	32+1	168	106	2945	1550

Abbreviations: C, Control; PE, preeclampsia; BP, blood pressure.

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.

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.

Table 2

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

lncRNAs			miRNAs			mRNAs
lncRNA	Log2(FC)	p-value	miRNA	Log2(FC)	p-value	mRNA	Log2(FC)	p-value
lnc-EGFL6-9	10.98	7.41×10−13	chr1_232	5.07	3.19×10−2	SPANXN5	8.33	3.80×10−6
lnc-CA7-1	9.30	1.14×10−11	chr1_230	5.07	3.19×10−2	ATXN3L	7.98	1.71×10−3
lnc-EGFL6-8	9.15	2.14×10−8	chr1_3167	5.07	3.19×10−2	KRT24	7.58	1.80×10−3
lnc-NCEH1-3	9.09	2.00×10−2	hsa-miR-451b	3.59	4.64×10−2	NPIPA7	7.35	2.88×10−3
lnc-EGFL6-4	9.05	4.4×10−9	chrX_40898	3.45	1.45×10−2	SLCO1B3-SLO1B7	7.21	2.48×10−2
lnc-CADPS-1	8.82	2.41×10−2	hsa-miR-371a-5p	2.82	5.58×10−10	RFLNA	7.18	1.32×10−4
lnc-EGFL6-11	7.90	1.97×10−5	hsa-miR-372-5p	2.73	5.62×10−6	TMPRSS11E	7.15	4.80×10−3
lnc-TNFRSF11B-4	7.70	4.89×10−2	hsa-miR-372-3p	2.61	4.93×10−10	KERA	6.99	1.01×10−2
LINC01203	7.12	4.84×10−4	chr2_21689	2.40	1.42×10−2	FOXA1	6.96	7.08×10−4
LINC02405	7.09	9.40×10−4	chr2_21690	2.40	1.42×10−2	DNAJC12	6.87	9.61×10−3
lnc-ERMN-4	−9.63	1.37×10−2	chr7_34302	−5.11	2.20×10−2	SFT2D3	−10.66	6.37×10−3
lnc-MANSC4-10	−9.31	8.44×10−12	chr7_35189	−5.11	2.20×10−2	URGCP-MRPS24	−8.07	1.44×10−6
lnc-OR7C2-1	−8.93	2.23×10−2	chr7_35545	−5.11	2.20×10−2	MUSTN1	−6.69	4.60×10−2
lnc-DCTD-24	−7.21	1.62×10−5	hsa-miR-4753-3p	−4.63	1.69×10−2	SHISA3	−6.14	3.26×10−2
lnc-SUPT6H-3	−6.98	1.02×10−3	hsa-miR-2276-3p	−3.25	3.72×10−2	TMEM72	−6.05	6.89×10−3
lnc-HSFX4-1	−6.97	7.32×10−3	chr11_6446	−2.70	3.48×10−2	PRTN3	−6.02	3.78×10−2
lnc-MFRP-4	−6.91	1.84×10−3	hsa-miR-138-1-3p	−1.96	2.09×10−2	C7orf33	−5.93	4.44×10−2
lnc-CSMD1-18	−6.90	1.17×10−3	hsa-miR-135a-5p	−1.83	5.48×10−3	H1-1	−5.73	1.87×10−2
lnc-FAM153B-3	−6.75	2.38×10−3	chr5_31003	−1.67	2.29×10−2	ASB12	−5.50	2.94×10−2
lnc-RPL13-4	−6.53	3.36×10−3	hsa-miR-138-5p	−1.42	1.49×10−2	MYOD1	−5.43	3.67×10−2

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).

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.

Table 3

Top 10 Enriched GO Terms of the DEGs

Terms	Pathway Description	Count	p-value
GO.BP:0050900	Leukocyte migration	21	2.07×10−6
GO.BP:0031214	Biomineral tissue development	11	5.84×10−6
GO.BP:0030282	Bone mineralization	10	1.36×10−6
GO.BP:0070167	Regulation of biomineral tissue development	9	2.27×10−6
GO.BP:0002062	Chondrocyte differentiation	9	2.85×10−5
GO.BP:0030500	Regulation of bone mineralization	8	3.86×10−6
GO.BP:0035630	Bone mineralization involved in bone maturation	6	2.39×10−9
GO.BP:0043931	Ossification involved in bone maturation	6	8.88×10−8
GO.BP:0070977	Bone maturation	6	1.81×10−7
GO.BP:0048799	Animal organ maturation	6	6.02×10−7
GO.CC:0031012	Extracellular matrix	17	4.28×10−4
GO.CC:0044445	Cytosolic part	13	4.16×10−5
GO.CC:0022626	Cytosolic ribosome	11	5.01×10−7
GO.CC:0044391	Ribosomal subunit	10	3.03×10−4
GO.CC:0022627	Cytosolic Small ribosomal subunit	6	3.16×10−5
GO.CC:0045334	Clathrin-coated endocytic vesicle	6	2.15×10−4
GO.CC:0015935	Small ribosomal subunit	6	4.81×10−4
GO.CC:0005833	Hemoglobin complex	4	1.53×10−5
GO.CC:0042788	Polysomal ribosome	4	7.95×10−4
GO.CC:0031838	Haptoglobin-hemoglobin complex	3	3.81×10−4
GO.MF:0048018	Receptor ligand activity	16	1.42×10−3
GO.MF:0043177	Organic acid binding	12	3.72×10−5
GO.MF:0005125	Cytokine activity	10	1.21×10−3
GO.MF:0005539	Glycosaminoglycan binding	10	1.34×10−3
GO.MF:0031406	Carboxylic acid binding	9	1.76×10−3
GO.MF:0019825	Oxygen binding	5	1.48×10−4
GO.MF:0008009	Chemokine activity	5	6.43×10−4
GO.MF:0042379	Chemokine receptor binding	5	2.48×10−3
GO.MF:0031720	Haptoglobin binding	3	3.18×10−4
GO.MF:0005344	Oxygen carrier activity	3	9.24×10−4

Abbreviations: BP, biological process; CC, cellular component; MF, molecular function.

Table 4

Top 10 Enriched KEGG Pathways of the DEGs

ID	Pathway Description	Count	p-value
hsa04060	Cytokine-cytokine receptor interaction	15	1.37×10−4
hsa05323	Rheumatoid arthritis	9	2.69×10−5
hsa03010	Ribosome	9	1.47×10−3
hsa00980	Metabolism of xenobiotics by cytochrome P450	8	4.57×10−5
hsa05204	Chemical carcinogenesis	8	7.86×10−5
hsa05150	Staphylococcus aureus infection	8	2.19×10−4
hsa04640	Hematopoietic cell lineage	7	1.46×10−3
hsa00982	Drug metabolism cytochrome P450	6	1.37×10−3
hsa00860	Porphyrin and chlorophyll metabolism	5	6.88×10−4
hsa05144	Malaria	5	1.54×10−3

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).

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.

Table 5

Top 10 Enriched GO Terms of the DELs-Target Genes

Terms	Pathway Description	Count	p-value
GO.BP:0031330	Negative regulation of cellular catabolic process	10	2.02×10−3
GO.BP:0002181	Cytoplasmic translation	7	2.64×10−3
GO.BP:0048813	Dendrite morphogenesis	7	6.43×10−3
GO.BP:0032768	Regulation of monooxygenase activity	5	1.55×10−3
GO.BP:0007257	Activation of JUN kinase activity	4	1.90×10−3
GO.BP:0021932	Hindbrain radial glia guided cell migration	3	2.80×10−4
GO.BP:0021535	Cell migration in hindbrain	3	1.01×10−3
GO.BP:0045725	Positive regulation of glycogen biosynthetic process	3	1.48×10−3
GO.BP:0021533	Cell differentiation in hindbrain	3	2.06×10−3
GO.BP:0070875	Positive regulation of glycogen metabolic process	3	2.06×10−3
GO.CC:0044445	Cytosolic part	9	8.26×10−3
GO.CC:0000151	Ubiquitin ligase complex	9	1.47×10−2
GO.CC:0022626	Cytosolic ribosome	6	5.32×10−3
GO.CC:0005844	Polysome	5	3.41×10−3
GO.CC:0042788	Polysomal ribosome	4	8.40×10−4
GO.CC:0022625	Cytosolic large ribosomal subunit	4	1.13×10−2
GO.CC:0031462	Cul2-RING ubiquitin ligase complex	3	1.05×10−3
GO.CC:0016281	Eukaryotic translation initiation factor 4F complex	2	7.98×10−3
GO.CC:0044666	MLL3/4 complex	2	9.66×10−3
GO.CC:0030130	Clathrin coat of trans-Golgi network vesicle	2	1.15×10−2
GO.MF:0048037	Cofactor binding	15	7.12×10−3
GO.MF:0050662	Coenzyme binding	12	1.15×10−3
GO.MF:0035326	Enhancer binding	8	8.26×10−4
GO.MF:0003730	mRNA 3ʹ-UTR binding	6	1.20×10−3
GO.MF:0016835	Carbon-oxygen lyase activity	5	4.80×10−3
GO.MF:0005089	Rho guanyl-nucleotide exchange factor activity	5	6.67×10−3
GO.MF:0008187	Poly-pyrimidine tract binding	4	1.37×10−3
GO.MF:0032813	Tumor necrosis factor receptor superfamily binding	4	4.70×10−3
GO.MF:0070742	C2H2 zinc finger domain binding	3	1.28×10−3
GO.MF:0030371	Translation repressor activity	3	6.48×10−3

Abbreviations: BP, biological process; CC, cellular component; MF, molecular function.

Table 6

Top 10 Enriched KEGG Pathways of the DELs-Target Genes

ID	Pathway Description	Count	p-value
hsa05131	Shigellosis	10	0.002711
hsa03013	RNA transport	6	0.04739
hsa01200	Carbon metabolism	5	0.027777
hsa04520	Adherens junction	4	0.019671
hsa01230	Biosynthesis of amino acids	4	0.023548
hsa00900	Terpenoid backbone biosynthesis	3	0.003792
hsa00603	Glycosphingolipid biosynthesis - globo and isoglobo series	2	0.019523
hsa00730	Thiamine metabolism	2	0.022102
hsa00910	Nitrogen metabolism	2	0.024814
hsa01210	2-Oxocarboxylic acid metabolism	2	0.027653

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.

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.19 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.20 Leng et al established a ceRNA network to reveal the role of lncRNA HCG27 in gestational diabetes mellitus.21 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.22 Besides, Zhang et al also built a ceRNA regulatory network consisting of 3 genes, 21 lncRNAs and 69 miRNAs in early-onset preeclampsia.23 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.24 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.25 It interacts with transforming growth factors and plays a pivotal role in TGF-β signaling transduction.26 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.32 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.33 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.