High-throughput deep screening and identification of four peripheral leucocyte microRNAs as novel potential combination biomarkers for preeclampsia.

OBJECTIVE: To identify the specific microRNA (miRNA) biomarkers of preeclampsia (PE), the miRNA profiles analysis were performed. STUDY
DESIGN: The blood samples were obtained from five PE patients and five normal healthy pregnant women. The small RNA profiles were analyzed to identify miRNA expression levels and find out miRNAs that may associate with PE. The quantitative reverse transcriptase-PCR (qRT-PCR) assay was used to validate differentially expressed peripheral leucocyte miRNAs in a new cohort. RESULT: The data analysis showed that 10 peripheral leucocyte miRNAs were significantly differently expressed in severe PE patients. Four differently expressed miRNAs were successfully validated using qRT-PCR method.
CONCLUSION: We successfully constructed a model with high accuracy to predict PE. A combination of four peripheral leucocyte miRNAs has great potential to serve as diagnostic biomarkers of PE.

Introduction

Preeclampsia (PE), a disease characterized by high blood pressure and proteinuria, is one of the causes of maternal and fetal morbidity and mortality and occurs in around 8% of pregnancies.[1] Part of the PE may progress quickly and cause eclampsia, placental abruption, HELLP syndrome, neonatal asphyxia, disseminated intravascular coagulation, stroke, heart failure or death despite intensive treatment. Previous studies have shown that placental hypoxia enhanced platelet aggregation and endothelial dysfunction in vascular aspect and immunological dysfunctions were responsible for the occurring of PE.[2] Even though several factors have been investigated, few effective results were revealed.[3, 4, 5, 6] However, the pathogenesis of PE has not been elucidated, although it has been suggested that PE encompasses a combination of genetic, immune and environment aspects. Neither its useful biomarkers nor its risk factor has been found to reach an agreement.

MicroRNAs (miRNAs) are small non-coding RNAs of 18 to 25 nt long that affect the stability and translational efficiency of target mRNAs.[7, 8] Abnormal expression levels of miRNAs have been associated with various diseases, including PE.[9, 10] The specific expressed miRNA might be used as the biomarker of diseases.

Circulating miRNAs have proven as powerful biomarkers for human disease.[11, 12] Many researchers have found a handful of miRNAs expressed as aberrant in PE samples.[13, 14] Anton et al.[15] reported that miR-210 may be a novel predictive serum biomarker for PE. On the other hand, Luque et al.[16] recently declared maternal serum miRNA assessment at the first trimester of pregnancy showed no predictive value for early PE.

The miRNA expression in leukocytes has been confirmed a stable state[17] that can be several diseases' predicted markers.[18, 19] To better address the role of miRNAs in PE, we extracted leukocytes from peripheral blood of PE patients and performed the miRNA sequencing of PE patients and normal healthy pregnant women using the second-generation sequencing technology. We found that a collection of miRNAs was expressed aberrantly in PE samples compared with the healthy pregnant women. Different expressed miRNA were validated by quantitative reverse transcriptase–PCR (qRT-PCR) in a large sample set. The receiver operating characteristic (ROC) analysis illustrated that a combination of four miRNA expression profiles might be used as biological markers for PE.

Methods

Ethics statement and samples

This study has approved by the Clinical Research Ethics Committee of Shanxi Medical College. PE was defined according to the criteria by American College of Obstetricians and Gynecologists practice bulletin on diagnosing and managing PE and eclampsia.[20] A total of 575 cases were sampled from the Second Hospital of Shanxi Medical University between 2011 and 2014. We filtered the cases of drug treatment, severe anemia, kidney failure, severe hypoalbuminemia, disagreement for experiment and failure in RNA extraction, and 57 cases were excluded from the experiment. Written informed consents were obtained from all participates. Five PE patients and five normal pregnant women were included in the first discovery stage (Supplementary Tables S1 and S2). Forty-seven samples were included in the second validation stage (Supplementary Table S3), which were separated into four groups: the control (N) group (13 normal cases of pregnant women); the PE group (13 cases of PE women aged between 21 and 29 years); PE with advanced age (PA) group (13 cases of PE women aged >30 years); and PE women with complications (PC) group (8 cases of PE women with chronic hypertension and gestational diabetes). Considering the ethics and convenience of pregnancy test, the peripheral blood were collected from all participant.

RNA extraction and sequencing

The procedure of RNA extraction was as follows: (a) an equal volume phosphate-buffered saline (1 ×) was added to fresh whole blood with anticoagulant and mixed; (b) the mixture slowly transferred to another centrifuge tube, which contained the lymphocyte separated liquids and mixture was allowed to settle over lymphocyte separated liquid (so as to not mix the two liquids and keep the interface clear), centrifuged at 3000 × g for 30 min at 4 °C; (c) leucocyte was separated with a pipette carefully, washed with phosphate-buffered saline (1 ×), collected using centrifuge and completely removed and discarded the supernate; (d) TRIzol reagent was added with 20 volumes of lymphocyte, washed the lymphocyte until the cell block was broken, kept the entire solution clear but not the viscous state; and (e) the solution was poured into dry ice or stored at −80 °C. RNA quantity and integrity were evaluated using Agilent 2100 BioAnalyzer (Agilent Technologies, Santa Clara, California, USA). Small RNA libraries were constructed using the method described in previous study.[21, 22] Briefly, for each library, 50 μg of the total RNA was size-fragmented on a 15% tris-borate-EDTA (TBE) urea polyacrylamide gel (Invitrogen, Waltham, Massachusetts, USA) and 15 to 30 base pair (bp) fraction was excised, using 10 bp ladder (Invitrogen) as marker. RNA was eluted from the polyacrylamide gel slice in 600 μl of 0.3 M NaCl overnight at 4 °C. The resulting gel slurry was passed through a Spin-X cellulose acetate filter column (Corning, Corning, New York, USA) and precipitated in two 300-μl aliquots by the addition of 750 μl of ethanol and 3 μl of glycogen (5 mg ml−1; Invitrogen). After washing with 75% ethanol, the pellets were allowed to air dry at 25 °C and dissolved in diethylpyrocarbonate (DEPC) water. The RNA was dephosphorylated by alkaline phosphatase and recovered by ethanol precipitation. The small RNA was ligated with 5′ adapter (5′-GUUCAGAGUUCUACAGUCCGACGAUC-3′) using T4 RNA ligase (Promega, Madison, Wisconsin, USA) in the presence of RNase Out (Invitrogen) overnight at 20 °C. The ligation reaction was stopped by the addition of 10 μl of 2 × Gel Loading Buffer II (Ambion, Waltham, Massachusetts, USA). The ligated RNA was size-fractionated on a 15% TBE urea polyacrylamide gel (Invitrogen), and a 40 to 70 bp fraction was excised. The RNA was eluted from the gel and precipitated as described above followed by resuspension in DEPC-treated water. The precipitated RNA was subsequently ligated to the 3′ RNA adapter (5′-pUCGUAUGCCGUCUUCUGC UUGidT-3′ p, phosphate; idT, inverted deoxythymidine) using T4 RNA ligase (Promega) in the presence of RNase Out (Invitrogen) overnight at 25 °C. The ligation reaction was stopped by the addition of 10 μl of 2 × Gel Loading Buffer II (Ambion). Ligated RNA was size-fractionated on a 10% TBE urea polyacrylamide gel (Invitrogen), and the 60 to 100 bp fraction was excised. The RNA was eluted from the polyacrylamide gel and precipitated from the gel as described above and resuspended in 5.0 μl of DEPC water. The RNA was converted to single-stranded cDNA using Superscript II reverse transcriptase (Invitrogen) and Illumina's small RNA RT-Primer (5′-CAAGCA GAAGACGGCATACGA-3′) following the manufacturer's instructions. The resulting cDNA was PCR-amplified with Hotstart Phusion DNA Polymerase (NEB, Ipswich, Massachusetts, USA) in 15 cycles using Illumina's small RNA primer set (5′-CAAGCAGAAGACGGCATA CGA-3′ 5′-AATGATACGGCGACCACCGA-3′).

PCR products were purified on a 12% TBE urea polyacrylamide gel (Invitrogen), and the 100-bp fraction was excised. The DNA was eluted into 100 μl of 1 × NEBuffer 2 at room temperature for 2 h. The resulting gel slurry was passed through a Spin-X filter (Corning) and precipitated by the addition of 325 μl of ethanol, 10 μl of 3 M sodium acetate and 3 μl of glycogen (5 mg ml−1; Invitrogen). After washing with 75% ethanol, the pellet was allowed to air dry at 25 °C and dissolved in 10 ml of resuspension buffer (10 mM Tris-HCl, pH 8.5) by incubation at 4 °C for 10 min. The purified PCR products were quantified on the Agilent DNA 1000 chip and diluted to 10 nM for sequencing on the Illumina Hiseq 2000 platform (San Diego, California, USA) according to the manufacturer's protocol.

Sequencing data analysis

After removing low-quality reads and adapter reads, the remaining reads with at least 18 nucleotides (nt) were aligned against human reference genome (version hg19) using SOAP (version 2.0, Microsoft Corporation, Redmond, Washington, USA).[22] Perfect match sequence tags were then annotated as coding exons, repeats, rRNA, tRNA, snRNA and snoRNA by aligning sequence tags to the databases of UCSC RefGene, Repbase, NCBI Refseq and NCBI Genbank (http://www.ncbi.nlm.nih.gov/). Expression levels of known miRNAs were calculated after mapping the tags to the precursors/mature miRNA sequence in miRBase (version 18).[23] NOISeq[24] was used in the detection of differentially expressed miRNAs between the PE and healthy group under the criteria: probability >0.70 and log2 ratio >1.5. The target genes for differentially expressed miRNAs were predicted by DIANA-microT.[25] To obtain dependable prediction, the DIANA-microT threshold was set to 0.9. Gene Ontology (GO) functional analysis[26] and KEGG pathway[27] analysis were carried out to identify significantly enriched GO terms and pathways associated with the target genes using miRPath,[28] with a cutoff of corrected P-value<0.05.

Real-time PCR validation

miRNA-451a, miRNA-15a-3P, miRNA-31-3P and miRNA-122-5P were selected to conduct real-time PCR validation in a four-group data set as mentioned above. The known miRNA has-U6 was used as the endogenous control during the process. RT reactions were carried out by miRcute miRNA First-Strand cDNA Synthesis Kit (TIANGEN, Beijing, China) according to the manufacturer's protocol. Amplification were conducted under the following condition: 94 °C for 20 s, then 94 °C for 20 s, and finally 60 °C for 34 s with 40 cycles. All reactions were run in triplicate. The relative expression level was calculated using the comparative CT method and normalized by the expression level of control miRNA has-U6.

The differential expression analysis of miRNAs between difference groups were analyzed by the method of analysis of variance tests using the software SPSS (version 18.0, SPSS Inc., Chicago, Illinois, USA). The significance P-value cutoff was set as 0.05.

Disease prediction model construction and ROC analysis

Logistic regression analysis of the four miRNA combination regression was carried out using glm package of R (http://www.r-project.org/). Relative expression level profiles of four miRNAs from qRT-PCR experiments were used as independent variable, and the PE status was used as dependent variable. Leave-one-out cross-validation accuracy was calculated using cv.binary from DAAG package of R (http://www.stats.uwo.ca/DAAG/). ROC analysis was carried out using pROC package of R.[29]

Results

Sequencing data description and small RNA annotations

Ten small RNA libraries were constructed and sequenced using second-generation sequencing technology. Primary sequencing data of each sample contained 12.98 M raw reads on an average. After filtering out about 3% low-quality and adaptor reads, 12.76 M final reads per sample were reserved. The length of most reads was located in 20 to 24 nt. In our study, >70% small RNAs were identified as miRNA.

Differentially expressed miRNAs

Eight known miRNAs were significantly differentially expressed. Among these eight miRNAs, four were significantly upregulated in PE with up to 10-fold changed level, and the expression of the other four miRNAs were downregulated in patients compared with the normal group (Table 1).

Table 1

The most significantly differentially expressed miRNA in preeclampsia patients

miRNA	Log2 ratio (cases/controls)	Probability
hsa-miR-122-5p	5.057686705	0.815618937
hsa-miR-451a	1.722145196	0.799611952
hsa-miR-299-5p	2.053540151	0.753608847
hsa-miR-1299	2.952837121	0.735428793
hsa-miR-15a-3p	−2.060379369	0.755665503
hsa-miR-31-3p	−1.721995698	0.70613116
hsa-miR-4785	−1.575520542	0.726329065
hsa-miR-4752	−8.593514168	0.758963912

Abbreviation: miRNA, microRNA.

Based on the prediction of MircoT, 125 genes were identified as target genes for these eight miRNAs (Supplementary Figure S1). GO terms of sensory perception of chemical stimulus (GO: 0007606, P=0.044) and G-protein coupled receptor activity (GO:0004930, P=0.005) were significantly enriched in these target genes of differentially expressed of miRNAs. The overrepresented pathways of these target genes included mucin type O-Glycan biosynthesis (ko 00512, P=4.85 × 10−4) and endocytosis (ko 04144, P=7.83 × 10−4) (Table 2).

Table 2

Pathway enrichment analyses of target genes of differentially expressed miRNAs between healthy pregnant women and preeclampsia

KEGG pathway	P-value	No. of genes	No. of miRNAs
Mucin type O-Glycan biosynthesis	0.000485253	4	3
Endocytosis	0.000783152	19	7
Phosphatidylinositol signaling system	0.001235167	10	4
Ubiquitin-mediated proteolysis	0.001338773	14	6
Citrate cycle (TCA cycle)	0.002431726	5	2
Carbohydrate digestion and absorption	0.004777441	5	2
B-cell receptor signaling pathway	0.004777441	9	5
Chronic myeloid leukemia	0.005586001	8	5
GABAergic synapse	0.01768863	8	4
mTOR signaling pathway	0.02970234	7	4
Bacterial invasion of epithelial cells	0.02970234	8	5
T-cell receptor signaling pathway	0.02970234	10	6
Inositol phosphate metabolism	0.04573376	7	4

Abbreviations: miRNA, microRNA; mTOR, mammalian target of rapamycin; TCA, tricarboxylic acid.

qRT-PCR validations

Of these eight differential expressed miRNAs we found in the sequencing stage, hsa-miR-31 and hsa-miR-15 were reported to be deregulated in other PE studies either by sequencing method or PCR validations,[30, 31] miR-451 were differentially regulated in primary trophoblasts exposed to hypoxia[32] and miR-122 were reported to be associated with gene GIT1 (G-protein-coupled receptor kinase interactor 1), a gene which is related to the endothelial cell repair.[33] Therefore, we selected has miR-31, miR15, miR 451 and miR122 to perform qRT-PCR validation in other 47 subjects. Analysis of variance test showed that there were significant different means of the expression levels of each miRNA between groups. T-test showed that the expression levels of miR-15a-3p and miR-31-3p were significantly downregulated in the three patients groups compared with the healthy group. The expression levels of miRNA-451a and miRNA-122-5p were upregulated in PE group and downregulated in the PA and PC groups (Figure 1 and Table 3). We also found that expression of miR-122-5p was downregulated in the PA and PC groups compared with the PE group. Meanwhile, the expression of miR-122-5p in the PC group was significantly increased compared with the PA group.

Figure 1

Box-plot diagrams of relative microRNA (miRNA) expression levels in peripheral blood. Relative miRNAs expression values from quantitative reverse transcriptase–PCR are shown in y axis.

Table 3

Relative expression (2−ΔΔCт) of four selected miRNAs in different groups (mean±s.d.)

miRNA ID	N group	PE group	PA group	PC group	P-value (ANOVA)
miRNA-451a	1.00±0.49	2.47±1.19a	2.04±0.76	2.57±1.54b	0.0361
miRNA-15a-3P	1.00±0.79	0.14±0.11b	0.12±0.16b	0.09±0.11b	0.0087
miRNA-31-3P	1.00±0.63	0.17±0.25b	0.04±0.06b	0.19±0.19b	<0.0001
miRNA-122-5P	1.00±0.46	2.15±0.60b	0.05±0.10b,c,d	0.49±0.35a,c	<0.0001

Abbreviations: ANOVA, analysis of variance; miRNA, microRNA; N, normal control; PA, preeclampsia with advanced age; PC, preeclampsia women with complications; PE, preeclampsia.

P-value<0.05 when compared with the N group.

P-value<0.01 when compared with the N group.

P-value<0.05 when compared with the PE group.

P-value<0.05 when compared with the PC group.

PE prediction model construction and ROC analysis

As four miRNAs could significantly separate PE patients from healthy women, these miRNAs can be used as biomarkers for PE. We constructed a logistic regression model to evaluate these four miRNAs and compare the model with single marker using ROC analysis. The results showed that the combination of four miRNAs achieved the highest area under the ROC curve (AUC=0.99, Figure 2) and was better than other four single markers. To prevent over fitting of the logistic regression model, the leave-one-out cross validation were performed, and we got an estimated accuracy of 0.94.

Figure 2

Receiver operating characteristic (ROC) curves to assess the utility of microRNAs (miRNAs) to differentiate preeclampsia and healthy women. ROC curves for miR-122 (yellow, area under the curve (AUC)=0.59), miR-451 (black, AUC=0.90), miR-15a (green, AUC=0.96), miR-31 (blue, AUC=0.95) and all four miRNAs (red, AUC=0.99).

Discussion

The activity of neutrophil may have a significant impact on the pathophysiology of vascular endothelial cells for PE.[34, 35] The supernatants of placental villus cultured from PE patients can significantly increase the expression of CD62L and CD1lb.[36] All these studies indicated that the biomarkers might be located in leukocyte. Thus our study mainly focused on finding circulating miRNAs from leukocyte that can distinguish PE from healthy pregnant women.

GIT1, a targeted gene of hsa-miR-122-5p, controls the activity of endothelial nitric-oxide synthase in endothelial cells, and down regulation of GIT1 can induce damages to endothelial cells.[33] Previous study[37] showed that the changes initiated by endothelial cell injury results in the clinical syndrome of PE. In our study, the upregulation of miR-122-5p in PE may induce the vascular endothelial damage, immune dysfunction and abnormal liver function.

The overactivity of miR-451 might also be responsible for the damage of the patients' circulatory system and nervous system.[38] Previous studies[39, 40] showed that miR-451 has a crucial role in promoting erythroid maturation under conditions of stress. Overexpression of miR-451 could significantly inhibit growth and induce cell apoptosis.[41, 42] The downregulated miR-451 has been reported in other studies.[43, 44] In our study, the upregulated miR-451 might be responsible for placenta hypoxia in the PE group. Meanwhile, the abnormal expression of miR-451a in PA and PC may be complicated. Further study are needed to explore the mechanism of miR-451a in PA and PC.

miR-15a negatively regulate BCL2 at a posttranscriptional level.[45] In our study, downregulation of miR-15a-3p might lead to an increased expression of BCL2 in leukocytes and induce normal leukocyte apoptosis by increasing oxidative stress and eventually lead to vascular endothelial injury.

miR-31 is blood vascular endothelial cell-specific posttranscriptional regulatory mechanism that has negative function in the development and function of lymphatic vascular system.[46] In our study, downregulation of miR-31-3p may enhance immune response and induce endothelial damage and pathogenesis of PE.

Screening by maternal serum pregnancy-associated plasma protein A, free β-hCG, activin A and inhibin A concentrations at 22+0 to 24+6 weeks' gestation detect 75% and 92% of patients who subsequently develop pre-eclampsia for false positive rates of 5% and 10%, respectively.[47] A multivariate predictive model combining 14 metabolites as metabolomic biomarkers in PE obtained an AUC of 0.94.[48] To distinguish PE patients from healthy women, we constructed the prediction model of PE using four miRNAs. Three of them achieved a very high AUC range from 0.9 to 0.96. In our study, the prediction model achieved 0.99 AUC and estimated accuracy of 0.94 when using the combination of these four miRNAs to predict the status of samples. Therefore, the constructed model is fit for predication of PE using these four miRNAs.

In summary, a set of miRNAs were found that significantly differentially expressed in PE samples based on the results of data analysis of next-generation sequencing and qRT-PCR validation. We successfully constructed a logistic regression model to distinguish PE patients from healthy women with extremely high accuracy. We provide information regarding the valuable miRNAs for diagnosing the systematic subtypes of PE and advance our understanding of the pathogenesis of PE.