Comprehensive Analysis of circRNAs, miRNAs, and mRNAs Expression Profiles and ceRNA Networks in Decidua of Unexplained Recurrent Spontaneous Abortion

The diagnosis and treatment of unexplained recurrent spontaneous abortion (URSA) are subject to debate, because the exact underlying mechanisms remain unclear. To address this issue, we elucidated the expression profiles of dysregulated circRNAs, miRNAs, and mRNAs and constructed circRNA-associated competitive endogenous RNA (ceRNA) networks by comparing the decidua of URSA with that of normal early pregnancy (NEP) using RNA-sequencing. In total, 550 mRNAs, 88 miRNAs, and 139 circRNAs were differentially expressed (DE) in decidua of URSA. Functional annotation revealed that DE mRNAs as well as potential target genes of DE miRNAs and DE circRNAs are mainly involved in immunologic function, such as antigen processing and presentation, allograft rejection, and T cell receptor signaling pathway. In addition, the top hub genes, including CCL4, DDX58, CXCL10, CXCL9, MX1, CD44, RPS2, SOCS3, RPS3A, and CXCL11, were identified. The mRNAs involved in ceRNA network were enriched in complement and coagulation cascades and protein processing in the endoplasmic reticulum. We found that circRNAs in the ceRNA network, which acted as decoys for hsa-miR-204-5p, were positively correlated with MFGE8 expression. Collectively, the results demonstrated that circRNAs, miRNAs, and mRNAs were aberrantly expressed in the decidua of patients with URSA and played a potential role in the development of URSA. Thus, the establishment of the ceRNA network may profoundly affect the diagnosis and therapy of URSA in the future.


INTRODUCTION
Recurrent spontaneous abortion (RSA), which refers to three or more spontaneous abortions with a single spouse, occurs in 1-5% of women of childbearing age and exerts a negative impact on human reproductive health (Medicine, 2008). Although the etiology of RSA has not yet been completely clarified owing to its complicated molecular and cellular regulation, known causes include chromosomal translocations, uterine abnormalities, thrombophilias, endocrine defects, autoimmune diseases, and infectious agents (Srinivas et al., 2006;Cohn et al., 2010;Daher et al., 2012;Rull et al., 2012;Asadpor et al., 2013;Pluchino et al., 2014;Dean et al., 2019). However, the cause and pathogeneses of approximately 40-50% of RSA cases remain unknown (Stephenson, 1996;Weimar et al., 2012); these are termed unexplained RSA (URSA). URSA poses a major obstacle for many couples who desire to conceive a child. Recent studies have shown a correlation between noncoding RNAs (ncRNAs) and the pathogenesis of URSA (Wang et al., 2012;Dong et al., 2014). Circular RNAs (circRNAs) a class of ncRNAs, are characterized by a covalent bond linking the 3' and 5' ends generated by back splicing. With the advent of high-throughput RNA sequencing technology, circRNAs that were considered as nonfunctional (Capel et al., 1993), were found to be involved in many diseases, such as atherosclerosis, neurological diseases, diabetes, and cancer, all of which involve processes similar to those of embryo implantation (Chen L. et al., 2015;Qu et al., 2015). In recent years, some studies have investigated the role of circRNAs in URSA. Li et al. compared the decidua tissues of the early recurrent miscarriage group with those of the control group using circRNA microarray analysis and identified differentially expressed (DE) circRNAs (Li C. et al., 2020). Another study used an RSA mouse model to demonstrate a molecular mechanism involving migration and invasion of circ-ZUFSP regulating trophoblast cells (Li Z. et al., 2020). This provides a new indicator for early diagnosis and potential treatment of RSA. However, a comprehensive understanding of the processes involved in URSA is lacking, and the potential role of circRNAs, miRNAs, and mRNAs in URSA remains unresolved.
The purpose of the current study was to gain an overview of the expression profiles of circRNAs, miRNAs, and mRNAs in URSA via RNA sequencing. In addition, because the role of circRNAs as miRNA sponges that regulate placental behavior has been confirmed (Ou et al., 2019), we constructed a circRNA-miRNA-mRNA interaction network to investigate the potential functions of these DE circRNAs. Our findings provide new evidence that may help better understand the molecular mechanisms underlying the role of circRNAs, miRNAs, and mRNAs in the pathogenesis of URSA.

Sample Collection
All participants were recruited from the Outpatient Department of Guangdong Provincial Fertility Hospital (Guangzhou, China), during the period between June 2016 and December 2019.
Patients included in the current study provided written consent before surgery, and the Ethical Committee of the Guangdong Provincial Fertility Hospital approved the study protocol. The decidua tissues were obtained from 9 women (age: 27 ± 6.38 years and gestational age: 8.36 ± 0.36 weeks) via dilation and curettage following spontaneous abortion. Normal decidua tissues were obtained from 9 women (age: 26.67 ± 3.64 years and gestational age: 7.88 ± 0.8 weeks) during early pregnancy (6-12 weeks) via dilatation and curettage after the termination of pregnancy. The clinical data are presented in Table 1. Inclusion criteria for URSA were three or more consecutive spontaneous abortions and no history of successful pregnancies. Exclusion criteria for URSA were infections, endocrine or metabolic disorders, anatomic abnormalities, autoimmune diseases, paternal or maternal chromosomal abnormalities, and abnormal karyotyping of the miscarriage product. Inclusion criteria for the normal controls were 1) elective termination of normal early pregnancy, 2) at least one live birth, and 3) no history of miscarriage, preeclampsia, ectopic pregnancy, preterm delivery, and systemic diseases. After dilation and curettage, decidua tissues were immediately separated from the products of conception, washed thoroughly with sterile normal saline, and stored in liquid nitrogen for future use.

RNA Library Preparation and Sequencing
A total amount of 5 μg RNA per sample was used as input material for the RNA sample preparations. Firstly, ribosomal RNA was removed by Epicentre Ribo-zero ™ rRNA Removal Kit (Epicentre, United States), and rRNA free residue was cleaned up by ethanol precipitation. The sequencing libraries were generated by NEBNext ® Ultra ™ Directional RNA Library Prep Kit for Illumina ® (NEB, United States) following manufacturer's recommendations. In addition, A total amount of 3 μg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, United States) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. The clustering of the index coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3 cBot HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50bp single-end reads were generated. The raw sequencing data presented in the study has been deposited in the Sequence Read Archive (SRA) database (SRA: PRJNA819210, PRJNA819201).

Protein-Protein Interaction Network Construction, Hub Gene Selection, and Module Analysis
The protein level interactions of DE mRNAs were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/cgi/input.pl). Interactions with a combined score ≥0.4 were considered statistically significant. Cytoscape software (version 3.7) was applied to visualize the PPI network. The proteins in the central node (hub proteins) that were highly connected to a series of genes in the network were considered most closely associated with pathogenesis. The top 10 hub genes, ranked by degree of connectivity, were identified using the cytoHubba APP plug-in of Cytoscape. Significant modules in the PPI networks were determined using the Molecular Complex Detection (MCODE) APP (v1.5.1) plug-in of Cytoscape.

Real-Time Quantitative PCR (RT-qPCR) Validation
Total RNAs were extracted using a Qiagen RNeasy mini kit (Qiagen, Hilden, Germany) and reverse transcribed using SuperScript III First-Strand Synthesis system (Life Technologies, Carlsbad, CA, USA). Real-time quantitative PCR was carried out using TB Green Premix EX Taq II (Takara, Beijing, China) and analyzed using the AB Step One Plus System (Applied Biosystems, AB, United States). Convergent primers were designed for linear mRNA production while divergent primers targeting the predicted back-spliced region were designed for the circular form. The specific quantitative primers designed using oligo7 software are listed in Supplementary Table S1. Primers of β-actin (for circRNA and mRNA) and U6 (for miRNA) were designed as an endogenous control. Each experiment was performed in triplicate.

GO and KEGG Analyses and GSEA
GO, which describes genes from any organism (http://www. geneontology.org), covering the domains of Biological Process (BP), Cellular Component (CC), and Molecular Function (MF), was used. KEGG pathway analysis (http://www.genome.jp/kegg) was used to analyze the functions of mapped genes. Fisher's exact test, the chi-square test, and the false discovery rate (FDR) were used for level of significance detection, where the p-value denoted the significance of the GO term and pathways correlated with the condition. A small FDR indicated a small error in judging the p-value. Gene Set Enrichment Analysis (GSEA) was performed using GSEA-P software to evaluate the RNA-Seq data at the gene level, with gene sets from the Molecular Signatures Database (MSigDB).

Construction of a circRNA-miRNA-mRNA Network
The expression levels of circRNAs, miRNAs, and mRNAs that showed significant differences between URSA and NEP were analyzed. We used miRanda (http://www.microrna.org/ microrna/) to predict miRNA binding seed sequence sites, where overlapping of the same miRNA binding sites on both circRNAs and mRNAs indicated circRNA-miRNA-mRNA interactions. Finally, we constructed the network using Cytoscape for Windows (National Institute of General Medical Sciences) software (Supplementary Figure S1).

Western Blot Analysis
Briefly, whole protein samples were extracted, and their concentrations were measured using a Pierce BCA Protein Assay Kit (Thermo Fisher Scientific), following which 20 μg of protein was loaded and separated on a 10% sodium dodecyl sulfate (SDS) polyacrylamide gel. The membranes were incubated overnight at 4°C with specific primary antibodies MFGE8 (Abcam, dilution in 1:500) and GAPDH (Proteintech, dilution in 1:1000), followed by incubation with goat anti-rabbit IgG conjugated to horseradish peroxidase (Santa Cruz, dilution in 1:5000) for 1 h at room temperature. The membranes were then visualized via a film following incubation with ECL solution (Cell Signaling Technology).

Statistical Analysis
Statistical analyses were performed using SPSS 20.0 software. All data were expressed as mean ± SD. Statistical significance was set at p < 0.05. The Student's t test was used to compare RT-qPCR results.

Identification of DE circRNAs, miRNAs, and mRNAs
The expression of circRNAs, miRNAs, and mRNAs in the URSA and NEP control groups were analyzed using RNA sequencing. After processing the raw data, we identified 5126 circRNAs, 1540 miRNAs, and 16065 mRNAs that were aberrantly expressed in URSA. Subsequently, we used p < 0.05 as the selection criteria and determined that 139 circRNAs were DE in the decidua of RSA patients, of which 65 DE circRNAs were upregulated and 74 DE circRNAs were downregulated ( Figure 1A; Table 2). Further, differences in miRNA profiles of URSA and control groups in response to changes in the expression profiles of circRNAs were investigated. Of these, 88 were DE miRNAs, including 51 upregulated DE miRNAs and 37 downregulated DE miRNAs ( Figure 1B; Table 3). Besides these, 355 and 195 DE mRNAs were upregulated and downregulated, respectively ( Figure 1C;  Figure 1G).

GO and KEGG Pathway Analysis of DE mRNA Combined With GSEA
GO, and KEGG analyses were conducted to predict the potential functions of the DE mRNAs. DE mRNAs in the BP group were mainly enriched in immune response, immune system process, and regulation of response to stimuli ( Figure 2A). Furthermore, evaluation of RNA-seq data at the genome-wide level via GSEA confirmed and supplemented the GO results. Using BP GO annotations, GSEA showed that the activation of immune response, positive regulation of interleukin 2 production and positive regulation of CD4 positive alpha-beta T cell differentiation were upregulated  in the decidua of URSA ( Figure 2C). To determine the signaling pathways associated with DE-mRNAs, KEGG pathway analysis was performed. The DE mRNAs were primarily enriched in antigen processing and presentation, Toll-like receptor signaling pathway, and Natural killer cell mediated cytotoxicity ( Figure 2B). GSEA results annotated via the KEGG pathway further showed that gene sets regulating antigen processing and presentation, the MAPK signaling pathway, and allograft rejection were upregulated in the decidua of URSA ( Figure 2D).

PPI Analysis of DE mRNAs
To further analyze the functions of DE mRNAs at the protein level and also to explore the core mRNAs involved in the cellular processes of URSA, we used the STRING database to screen for functional genes and constructed an intuitive network for the annotation of PPIs ( Figure 3A). The core mRNAs (known as hub genes) that corresponded to the DEGs   Figure 3B). Module analysis revealed two closely connected regions with the highest scores in the PPI network ( Figures 3C,D). CCL4, DDX58, and CXCL10 in the PPI network were upregulated in decidua of patients with URSA ( Figure 3E).

Distinct circRNA and miRNA Expression Profiles in Decidua Tissues of Patients With URSA
GO analysis showed that the target genes of DE circRNAs in the biological processes mainly regulated response to DNA damage stimuli, T cell activation via T cell receptor contact with antigen bound to MHC molecules on antigen presenting cells, and the Toll-like receptor 7 signaling pathway ( Figure 4A). The DE miRNAs were significantly enriched in cellular developmental processes, cellular component assembly, and developmental processes ( Figure 4C).
KEGG pathway enrichment analysis revealed that the DE circRNAs were mainly enriched in the Rap1 signaling pathway, the TGF-beta signaling pathway, and the tight junction and T cell receptor signaling pathway ( Figure 4B).
The pathways enriched in miRNA included the MAPK signaling pathway, the Ras signaling pathway, and gap junction and focal adhesion ( Figure 4D).

Construction of circRNA-miRNA-mRNA ceRNA Network
To identify the ceRNA network in URSA, we extracted miRNAs, mRNAs, and circRNAs that were abnormally expressed in the triple network and constructed a network based on the ceRNA hypothesis using Cytoscape ( Figures 5A,B). The associated mRNA correlated pathways were complement and coagulation cascades, protein processing in the endoplasmic reticulum, and transcriptional misregulation in cancer ( Figure 5C). We selectively analyzed pairs in which the circRNAs, miRNAs, and target genes were possibly associated with URSA. A subnetwork consisting of circRNAs (hsa_circ_0004558, novel_circ_0006147, hsa_circ_0008546, etc.) included those acting as ceRNAs of hsa-miR-204-5p targeting MFGE8 ( Figure 5B). Subsequently, hsa-miR-204-5p and MFGE8 were selected for RT-qPCR validation, which revealed that the expression trends and amplitudes were consistent with those of RNA-seq ( Figure 5D). To build on these observations, we further investigated MFGE8 protein expression in the decidua of URSA, and the result was in concordance with those of RT-qPCR ( Figure 5E). Additional results are listed in Supplementary Table S2.

DISCUSSION
The current study revealed significant differences between the RNA expression profiles in the decidua tissues of URSA and NEP groups. The pathways of DE mRNAs as well as the potential  target genes of DE miRNAs, and DE circRNAs were mainly involved in antigen processing and presentation, allograft rejection, the T cell receptor signaling pathway, and junction related pathways. A ceRNA network was constructed, which revealed that the pathways associated with the target mRNAs were complement and coagulation cascades and protein processing in the endoplasmic reticulum. Additionally, strict limitative requirements were applied to select possible circRNA-associated-ceRNA networks that were likely involved in the occurrence and development of URSA. Moreover, the potential mechanisms via which circRNAs binding to hsa-miR-204-5p enhanced MFGE8 expression, thereby contributing to URSA were verified.
Our results showed that mRNA regulated URSA via multiple signaling pathways, including antigen processing and presentation, the Toll-like receptor signaling pathway, and Natural killer cell mediated cytotoxicity. Classically, activated macrophages in decidua are mainly involved in antigen processing and presentation (Li et al., 2021) which has been shown to be associated with URSA . Toll like receptors are a family of pattern recognition receptors, which play an essential role in maintaining innate host immunity and Th1/ Th2 balance (Koga and Mor, 2010;Chen et al., 2016), thereby permitting and accepting implanted and growing semi allogenic embryos. Uterine NK (uNK) cells, which form an important part of the innate endometrial immune system, make up approximately 70% of all uterine lymphocytes during pregnancy and play a critical role in establishing and maintaining pregnancy. Numerous studies have investigated the role of uNK cells in adverse pregnancy outcomes, such as recurrent miscarriage (RM) or recurrent implantation failure (RIF), where many fertility clinics have provided immunotherapy to patient groups based on abnormal levels or activity of NK cells (Tenger et al., 2005;Perricone et al., 2006;Moffett and Shreeve, 2015;Mekinian et al., 2016;Plaçais et al., 2020). Using the PPI network, we identified ten hub genes, CCL4, DDX58, CXCL10, CXCL9, MX1, CD44, RPS2, SOCS3, RPS3A, and CXCL11. According to reports, these genes, which are mainly involved in inflammatory functions (Camacho-Arroyo et al., 2021;Jiang et al., 2021), and innate immunity (Nizyaeva et al., 2019), are essential for maintaining normal pregnancy.
Previous studies have indicated that ncRNAs are involved in several stages of implantation and contribute to the dysfunction of trophoblasts in pregnancy-related disorders (He et al., 2013;Chen H. et al., 2015). A study conducted by Dong et al., revealed that miRNAs that were abnormally and significantly expressed in villi and decidua, were involved in several biological processes, including the regulation of transcription and protein amino acid phosphorylation, while KEGG pathway analysis indicated that the predicted target genes were related to adherence junctions, apoptosis, and the T cell receptor signaling pathway (Dong et al., 2014). Another study revealed that the target genes of these miRNAs were closely associated with the mitogen-activated protein kinase (MAPK) pathway, as well as with the B cell receptor, and T cell receptor signaling pathways in the URSA development process (Tang et al., 2016). Similar to these results, the results of the current study indicated that 139 circRNAs were significantly differentially expressed in human decidua, compared to NEP, and that these were primarily involved in Rap1 signaling pathways, the TGF-beta signaling pathway, and the tight junction and T cell receptor signaling pathways. Rap1 signaling pathway is implicated in several basic cellular functions, including cell-cell interactions, cell-matrix adhesion, proliferation, and regulation of cellular polarity during fetal development. Magdalena et al., highlighted the key role of Rap1 signaling in maintaining epithelial and endothelial cell junction integrity that helps maintain fetus viability (Chrzanowska-Wodnicka et al., 2015). TGF-beta, which is important in ensuring maternal support for embryo development, modulates the extent of decidua invasion (Jones et al., 2006), which finding was substantiated by our circRNA pathway analysis. Both circRNAs and miRNAs were enriched in junction related pathways, including tight junction, gap junction and focal junction, which are essential for normal pregnancy (Fukuda and Sugihara, 2008;, and in T cell adhesion to extracellular matrix proteins, which have been described as determinants of successful or failing pregnancy (Chaiworapongsa et al., 2002). Tang et al., suggested that downregulation of miR-181d (Tang et al., 2016) which contains relevant binding sites for upregulated circRNAs, may contribute to the inhibition of adhesion and cell proliferation in RSA NK cells (Qian et al., 2017). But, hsa-miR-141-3p which was top 10 of differentially expressed microRNAs was low expression in the present study, while high expression in NK cells isolated from the decidua of URSA (Li and Li, 2016). A previous study in mouse showed that embryo implantation sites were significantly decreased (Liu et al., 2013), following uteri horn injection of inhibitors or mimic of miR-141. The inconsistency between the results of the present study and previous studies may be due to the existence of a balance between the differentiation and apoptosis functions performed by hsa-miR-141-3p. Considered together, our findings further demonstrated the vital role of ncRNAs in URSA development.
The functions of circRNAs include sponging miRNAs, competing with endogenous RNAs, regulating gene transcription, and interacting with RNA-binding proteins (RBP) (Hansen et al., 2013). In particular, our results confirmed that the upregulated circRNA (has_circ_0004558, novel_circ_0006147, etc.) acted as ceRNAs of hsa-miR-204-5p, targeting MFGE8 in patients with URSA. MFGE8 is a glycoprotein which has been widely studied by many engaged in physiological and pathological research, especially regarding its role in the immune system (Matsuda et al., 2011;Oba et al., 2011). Previous studies have suggested that it participates in the implantation process, because it was found upregulated during the window of implantation (Mirkin et al., 2005;Franchi et al., 2011;Schmitz et al., 2014). On the other hand, MFGE8 is also known to be involved in the inflammatory process (Komura et al., 2009), and is regulated by TNF-α in the human endometrium . It has been reported that overexpression of MFGE8 may impair physiological processes (Yamamoto et al., 2014;Zhao et al., 2015). In keeping with this notion, Schmitz et al. (Schmitz et al., 2017) reported that MFGE8 expression in the glandular epithelium of patients with endometriosis and infertility was Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 858641 8 significantly increased, when compared to that of healthy fertile patients. The present study also demonstrated that the expression levels of MFGE8 mRNA and protein were increased in women with URSA. In summary, these results indicated that circRNAs played a crucial role in the occurrence of URSA, and that the sponging function of miRNAs may be the main pathway involved. Our study was affected by several limitations which should be acknowledged. Firstly, the small cohort studied was a major limitation that prevented us from drawing clear conclusions from our results. Secondly, the present study was unable to substantiate direct interaction between ceRNAs.
Collectively, we found that the expression levels of circRNAs, miRNAs, and mRNAs in the URSA differed significantly from those in the NEP group. In addition, our results enabled us to propose probable pathways via which circRNAs function in URSA. Furthermore, circRNA-associated ceRNAs were constructed and MFGE8 was identified as a critical component of the circRNA-miRNA network. These novel networks may lead to the discovery of potential biomarkers or therapeutic targets in URSA.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://dataview.ncbi.nlm. nih.gov/object/PRJNA819201 and https://dataview.ncbi.nlm.nih. gov/object/PRJNA819210.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethical Committee of the Guangdong Provincial Fertility Hospital. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual (s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
WQ and LH designed the study and coordinated and participated in every part of the study. XL and JW participated in the design of the study and helped in the study and drafting of the manuscript. GS participated in the design and statistical analyses. XZ and HN participated in the design of the study and performed the PCR analyses. All of the authors read and approved the final manuscript.