Transcriptome sequencing identified the ceRNA network associated with recurrent spontaneous abortion

Abstract

Background

Recurrent spontaneous abortion (RSA) is one of the common complication of pregnancy, bringing heavy burden to the patients and their families. The study aimed to explore the lncRNA-miRNA-mRNA network associated with recurrent spontaneous abortion.

Methods

By transcriptome sequencing, we detected differences in lncRNA, miRNA and mRNA expression in villus tissue samples collected from 3 patients with RSA and 3 normal abortion patients. Differentially expressed lncRNAs, miRNAs and genes (DELs, DEMs and DEGs, respectively) were identified, and Geno Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to determine the functions of DELs and DEGs, which were analysed by Fisher’s test. We also observed the regulatory relationships between miRNA-mRNA and lncRNA-miRNA by Cytoscape 3.6.1.

Results

The results showed that 1008 DELs (523 upregulated and 485 downregulated), 475 DEGs (201 upregulated and 274 downregulated) and 37 DEMs (15 upregulated and 22 downregulated) were identified. And we also constructed a novel lncRNA-related ceRNA network containing 31 lncRNAs, 1 miRNA (hsa-miR-210-5p) and 3 genes (NTNG2, GRIA1 and AQP1).

Conclusions

lncRNA-related ceRNA network containing 31 lncRNAs, 1 miRNA (hsa-miR-210-5p) and 3 mRNAs (NTNG2, GRIA1 and AQP1) was constructed. The results may provide a basic theory for elucidating the mechanism underlying RSA.

Background

Recurrent spontaneous abortion (RSA), a complication of pregnancy, is defined as three or more consecutive spontaneous abortions with the same spouse, and gestational age < 20 weeks of spontaneous abortion [1]. The incidence rate of that was about 5%. To date, chromosomal abnormalities of parents and embryos, anatomical factors, thrombosis, immunological factors, endocrine factors and environmental factors have been reported to be correlated with the occurrence of RSA [2,3,4,5]. RSA seriously endangers women's reproductive health and causes great physical and mental pain for patients and their families. In recent years, new risk factors have been gradually recognized, but the aetiology of many patients with RSA remains unknown. At present, the commonly used treatment methods include immunotherapy, endocrine therapy and anticoagulant therapy. However, the effects of these treatments are not satisfactory. Therefore, it is of great significance to study the aetiology of RSA.

Noncoding RNAs play an important role in the almost all pathological or pathological processes, such as embryonic development, cell proliferation, differentiation, apoptosis, infection and immune response, including RSA [6, 7]. Long noncoding RNAs (lncRNAs), highly conserved noncoding RNAs, have also been found to be involved in RSA related studies [8,9,10]. Gu et al. observed [8] that polymorphisms in lncRNA HULC may be related to the susceptibility to RSA in the Southern Chinese population. Xuan et al. also found that the lncRNA MALAT1 rs619586 G mutation reduced the risk of RSA [9]. Che et al. found that lncRNA CCAT2 rs619586 G mutation may have a potential protective effect and reduce the risk of RSA in southern China [10]. The results described above indicated that lncRNAs played a role in RSA. Furthermore, miRNAs have also been found to be indispensable for the pathogenesis of RSA [11, 12]. By assessing the influence of USP25 on trophoblasts, Ding et al. found that USP25 expression was negatively regulated by miR-27a-3p, and this effect contributed to the pathogenesis of RSA by suppressing the migration and invasion of trophoblasts [11]. It was observed that the upregulation of miR-365 expression may promote the occurrence of RSA by reducing the expression of SGK1, suggesting that miR-365 may be used as a prognostic biomarker and therapeutic target for RSA reported by Zhao et al. [12]. As a new model of gene expression regulation, the large regulatory network of ceRNAs is helpful for exploring the gene function and regulatory mechanisms at a deeper level and for more thoroughly and comprehensively understanding many biological phenomena. However, so far, lncRNA-miRNA interactions and lncRNA-miRNA-mRNA networks have not been reported in RSA.

In the study, we constructed a lncRNA-associated ceRNA network to explore the pathogenesis of RSA in 3 patients with RSA and 3 normal abortion patients, providing a theoretical basis for the elucidation and treatment of RSA in the future.

Material and methods

Subjects

The villus tissue samples were collected from 3 patients with RSA and villus tissue samples from 3 normal abortion patients were served as controls. The fresh tissues were stored in liquid nitrogen tanks for subsequent use. The inclusion criteria for the RSA patients were as follows: (1) patients with RSA suffered three or more consecutive spontaneous abortions at a gestational age of < 20 weeks; (2) female patients with RSA who suffered primary abortion and had no previous history of live births; (3) RSA patients underwent routine examinations, including examination of maternal infection, chromosome aberration, endocrine dysfunction, anatomical factors and autoimmune diseases. Patients who did not meet these conditions were excluded. The controls had at least one childbirth and had no history of spontaneous abortion. Moreover, the controls had no pregnancy-related complications. All the subjects have signed an informed consent form. The study was approved by the Second Affiliated Hospital of Hainan Medical College (2018R005-F01).

Transcriptome sequencing data analysis

Using the TRIzol Reagent (Thermo Fisher Science, USA), we extracted total RNA from the villus tissue samples. Subsequently, we measured the RNA concentration and purity. We performed lncRNA, miRNA and mRNA sequencing with the Illumina transcriptome chip. FastQC software and the R package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) were used to evaluate the quality of the original sequencing data. Using the Trimgalore method (http://www.bioinformatics.babraham.ac.uk/projects/trim_Galore/), we filtered raw reads to obtain clean reads for subsequent analysis. Besides, all the data were processed by quantile normalization.

Analysis of differentially expressed lncRNAs, miRNAs and genes

We used the Cuffdiff version 2.2.1 to identify differentially expressed lncRNAs, miRNAs and genes (DELs, DEMs and DEGs) in the villus tissue samples collected from 3 patients with RSA and 3 controls. p < 0:05 and |log₂FC| > 1 were used as the screening criteria. We completed the heatmap analysis of DELs, DEMs and DEGs with the ComplexHeatmap in the R package.

Functional analyses

In the present studie, Geno Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to determine the functions of DELs and DEGs, which were analysed by Fisher’s test using the clusterProfiler version 2.2.1. Biological process (BP), cell composition (CC) and molecular function (MF) were included in GO annotation analysis. KEGG enrichment analysis mainly focused on the related signaling pathways. p < 0.05 was regarded statistically significant.

Constructing the ceRNA network

We also observed the regulatory relationships between miRNA-mRNA and lncRNA-miRNA by Cytoscape version 3.6.1. The miRNA-mRNA-lncRNA network was constructed and visualized. According to the lncRNAs that directly interacted with mRNAs and regulated their activity as miRNA sponge, we analyze the data through the following three steps: (1) miRNAs targeted by DELs and the interaction between DELs and miRNAs were predicted using the LncTar software; (2) mRNAs targeted by DEMs and the interaction between DELs and miRNAs were predicted by the online tool (MiRDB, miRTarBase and Targetscan databases; (3) lncRNAs and miRNAs negatively regulated by mRNAs were integrated to construct a ceRNA network.

Results

Identification of differentially expressed lncRNAs, miRNAs and mRNAs

According to the quality of the total RNA from 3 patients with RSA and 3 normal abortion personnel (Additional file 1: Figure S1), the transcriptome sequencing showed that 1008 DELs (523 upregulated and 485 downregulated), 475 DEGs (201 upregulated and 274 downregulated) and 37 DEMs (15 upregulated and 22 downregulated), which were shown in the heat map (Fig. 1) and the volcano map (Additional file 3: Figure S3, Additional file 4: Figure S4 and Additional file 5: Figure S5). Table 1 illustrated the DEGs and DELs (top 20), and Table 2 listed the top 15 DEMs. The thresholds of the screening data were p < 0.05 and |log₂FC|> 1.

Fig. 1

Heatmap analysis of DEL, DEM and DEGs. DEL: Differentially expressed lncRNAs; DEMs: Differentially expressed miRNAs; DEGs: Differentially expressed genes

Table 1 Top 20 downregulated DELsand DEGs

Table 2 Top 15 upregulated DEMs

GO and pathway analysis of DELs

To further study the transcriptome differences between the two groups, we performed GO and KEGG pathway analyses of DELs. In the Table 3 and Fig. 2a, the results of the top 10 enriched GO pathways of DELs showed that the biological process (BP) changes were in the regulation of body fluid levels, embryonic skeletal system development, postsynapse organization, carbohydrate derivative transport, activation of JUN kinase activity, mammary gland epithelial cell proliferation, oxygen transport, gas transport, regulation of mammary gland epithelial cell proliferation and pericardium development. Additionally, the cell component (CC) changes of DELs were obviously enriched in transcription factor complex, axon part, postsynaptic specialization, histone methyltransferase complex, clathrin-coated pit, MLL1/2 complex, hemoglobin complex, MLL1 complex, haptoglobin-hemoglobin complex and exocyst. Moreover, molecular function (MF) changes were mainly enriched in DNA-binding transcription activator activity, RNA polymerase II-specific, enhancer sequence-specific DNA binding, enhancer binding, RNA polymerase II distal enhancer sequence-specific DNA binding, oxidoreductase activity, acting on NAD(P)H, molecular carrier activity, kinesin binding, laminin binding, haptoglobin binding and oxygen carrier activity. As shown in the Table 4 and Fig. 2b, the top 10 enriched KEGG pathways of DELs were in Alzheimer’s disease, Thermogenesis, Thyroid hormone signaling pathway, Hippo signaling pathway, Hepatocellular carcinoma, Adherens junction, Arrhythmogenic right ventricular cardiomyopathy (ARVC), Vibrio cholerae infection, Glycosphingolipid biosynthesis—lacto and neolacto series and Antifolate resistance.

Table 3 Top 10 enriched GO pathways of DELs

Fig. 2

GO enrichment items and KEGG pathway analysis of DELs. a Showed that the top 10 enriched GO pathways of DELs were sorted by significance in biological process (BP), cellular component (CC) and molecular function (MF), respectively. b Showed the top 10 enriched KEGG pathways of DELs

Table 4 Top 10 enriched KEGG pathways of DELs

GO and pathway analyses of DEGs

To further study the transcriptome differences between the two groups, we performed the GO and KEGG pathway analysis of DEGs. The results of the top 10 GO pathways of DEGs showed that changes in biological processes (BP) were mainly enriched in regulation of metal ion transport, leukocyte cell–cell adhesion, regulation of leukocyte proliferation, antigen processing and presentation, antibiotic catabolic process, gas transport, cellular extravasation, hydrogen peroxide catabolic process, oxygen transport and eosinophil migration. In addition, cell component (CC) changes were mainly concentrated in extracellular matrix, actin cytoskeleton, contractile fiber, contractile fiber part, myofibril, sarcomere, hemoglobin complex, haptoglobin-hemoglobin complex, MHC protein complex and MHC class II protein complex. Molecular function (MF) changes were mainly distributed in actin binding, actin filament binding, organic acid binding, molecular carrier activity, antioxidant activity, oxygen binding, peroxidase activity, oxidoreductase activity, acting on peroxide as acceptor, haptoglobin binding and oxygen carrier activity (Table 5 and Fig. 3a). The top 10 KEGG pathways of DEGs were mainly enriched in Cell adhesion molecules (CAMs), Chemokine signaling pathway, Staphylococcus aureus infection, Viral protein interaction with cytokine and cytokine receptor, Malaria, B cell receptor signaling pathway, Leishmaniasis, Asthma, African trypanosomiasis and Allograft rejection (Table 6 and Fig. 3b).

Table 5 Top 10 enriched GO terms of DEGs

Fig. 3

GO enrichment items and KEGG pathway analysis of DEGs. a Showed that the top 10 enriched GO pathways of DEGs were sorted by significance in biological process (BP), cellular component (CC) and molecular function (MF), respectively. b Showed the top 10 enriched KEGG pathways of DEGs

Table 6 Top 10 enriched KEGG pathways of DEGs

Construction of the lncRNA-miRNA-mRNA ceRNA network

First, we constructed a lncRNA-miRNA network and miRNA-mRNA network. The lncRNA-miRNA network included 4607 negative interactions (1945 downregulated lncRNAs-upregulated miRNAs and 2662 upregulated lncRNAs- downregulated miRNAs), and the miRNA-mRNA network included 15 negative interactions (6 downregulated miRNAs- upregulated mRNAs and 9 upregulated miRNAs-downregulated mRNAs). Then, we constructed the lncRNA-miRNA-mRNA ceRNA network to identify their relationships based on the lncRNA, miRNA, and mRNA expression profiles, and plotted them using Cytoscape version 3.6.1. First, based on the threshold values (r < 0 and p value < 0.05), we evaluated the relationship between downregulated lncRNAs and upregulated miRNAs shown in Fig. 4a, and the relationship between upregulated lncRNAs and downregulated miRNAs was displayed in Additional file 2: Figure S2. Additionally, the results of the miRNA-mRNA relationship showed a significant link between hsa-miR-210-5p and mRNAs (NTNG2, GRIA1 and AQP1), as shown in Fig. 4b. Besides, we constructed the ceRNA network between DELs, DEMs and DEGs by the Pearson correlation coefficient. Finally, the ceRNA network contained 31 lncRNAs (PSD2-AS1, lnc-ACAN-2, lnc-STON1-1, lnc-HPS4-8, lnc-SHC2-1, lnc-LMO7DN-6, lnc-TPTE-12, lnc-ARRDC3-5, lnc-CHPF-4, lnc-OR1J1-2, lnc-GPAT4-1, lnc-ARPC5L-1, LYPLAL1-DT, lnc-PIWIL4-1, lnc-CCR8-2, lnc-RHBDD3-3, lnc-PPP1R3G-9, RAMP2-AS1, LINC01771, lnc-SFRP4-3, lnc-C1QL3-2, lnc-C6orf223-1, lnc-IGFBP3-2, lnc-CUL2-3, lnc-SRGAP2C-5, PRKCQ-AS1, lnc-C11orf95-5, lnc-IGFBP1-1, lnc-CYP3A7-1, lnc-GPC6-7 and lnc-FUBP1-3), 1 miRNA (hsa-miR-210-5p) and 3 mRNAs (NTNG2, GRIA1 and AQP1) as displayed in Fig. 5 and Table 7 (top20), illustrating that these molecules may be involved in the development of RSA.

Fig. 4

The interactions between lncRNA-miRNA and miRNA-genes were determined, respectively. a Showed the relationship between downregulated lncRNAs and upregulated miRNAs. b Listed the interactions between hsa-miR-210-5p and 3 genes

Fig. 5

The ceRNA network was constructed between DELs, DEMs and DEGs

Table 7 Construction of ceRNA network (Top 20)

Discussion

Recurrent spontaneous abortion is one of the common complications of pregnancy. In the past few decades, the disease has caused heavy psychological burden for couples who want to have children and their families. However, due to the high cost of treatment, many families have failed to realize their desire to have children. The present study showed that we found 1008 DELs, 475 DEGs and 37 DEMs in 3 patients with RSA and 3 normal abortion personnel by transcriptome sequencing of villous tissue samples. We also constructed a novel lncRNA-related ceRNA network containing 31 lncRNAs, 1 miRNA (hsa-miR-210-5p) and 3 mRNAs (NTNG2, GRIA1 and AQP1). The results may provide a theoretical basis for elucidating the mechanism of RSA.

NTNG2 (Netrin G2) the position of which on chromosome is 9q34.13 and encodes the protein NTNG2, a membrane anchor protein. It was found to promote the growth of axons and dendrites. Studies on the correlation of gene polymorphisms in schizophrenia revealed that NTNG1 and its paralogues for NTNG2 gene may be related to the pathophysiology of schizophrenia [13,14,15]. Another paper reported by Maroofian et al. illustrated that NTNG2 played a key role in neurotypical development [16]. Therefore, we speculate that NTNG2 and NTNG1 may play a role in neurological disorders. In addition, based on the bioinformatics analysis of the pediatric onset of multiple sclerosis, genes such as NTNG2 were found to be nodes of the network, and the expression of some miRNAs were significantly correlated with brain volume [17]. But to date, there is no report on the NTNG2-associated network in RSA.

GRIA1 (Glutamate Ionotropic Receptor AMPA Type Subunit 1) is located on the chromosome 5q33.2 and the encoded protein is the main excitatory neurotransmitter receptor in the mammalian brain. It is reported to be a ligand-gated ion channel, regulating the secretion of follicle-stimulating hormone and luteinizing hormone by controlling gonadotropin releasing hormone. Recently, Sugimoto et al. discovered that the gene was linked to the ovulation rate in cattle [18]. Cushman et al. found the correlation between GRIA1 SNPs and cattle infertility [19]. In addition, Sheikhha et al. also observed the relationship between GRIA1 variants and ovarian response to human menopausal gonadotropin in the group of Iranian women [20]. The above studies show that GRIA1 plays an important role in diseases related to pregnancy in women. But so far, the role of GRIA1 in RSA has not been reported. In our study, we used Cytoscape software to construct a network to combine noncoding RNAs to explore its function in RSA.

AQP1 (Aquaporin 1), located on chromosome 7p14.3, contains 4 exons. Some reports have revealed that AQP1 plays an important role in acute lung injury caused by endotoxic shock, delaying the occurrence of renal cyst, and acute lung and brain injury [21, 22]. Su et al. used lipopolysaccharide (LPS)-induced murine model of acute lung injury to detect the function of AQP1, suggesting that AQP1 may be involved in the progression of acute lung injury [23]. Also, noncoding RNAs interacting with AQP1, were involved in the development of acute lung injury. Long noncoding RNA CASC2 can reduce the apoptosis of lung epithelial cells and improve acute lung injury by regulating the miR-144-3p/AQP1 axis [24]. Recent studies have shown that AQP1 participated in the occurrence of diseases through the ceRNA network [25, 26]. Tang et al. observed that lncRNA CASC2 acted as miR-144-3p, and directly interacted with AQP1 after LPS induced A549 cells [25]. After lipopolysaccharide (LPS) induced sepsis, Fang et al. found that AQP1 has been reported to competitively bind to lncRNA H19 and regulated the expression of miRNA-874 [26]. But so far, there has been no study on lncRNA-miRNA-AQP1 in RSA.

At present, some researchers have obtained some results about recurrent abortion through transcriptome sequencing [27, 28]. In this study, we firstly performed transcriptome sequencing analysis on the tissues of 3 patients with RSA and 3 patients with normal abortion, and found key molecules by constructing lncRNA-related ceRNA network, which is helpful to explore the pathogenic mechanism of RSA. However, there are some limitations: (1) the sample size was insufficient; (2) the lncRNAs-miRNAs linked with RSA were not verified; (3) the lncRNA-mediated ceRNA network in RSA was not verified. In the future, we will continue to collect a large number of samples for verification, and further analyze the ceRNA network in RSA by transcriptome analyses, and use molecular biology to verify this network, providing a theoretical basis for the elucidation and treatment of RSA.

Conclusion

In summary, a lncRNA-related ceRNA network containing 31 lncRNAs, 1 miRNA (hsa-miR-210-5p) and 3 mRNAs (NTNG2, GRIA1 and AQP1) was constructed. The results may provide the basic theory for elucidating the mechanism underlying RSA.