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.

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

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.

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

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

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

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; Table 4). Results of the cluster analyses of the expression levels of DE circRNAs, microRNAs, and mRNAs are displayed in Figures 1D–F. Six DE RNAs (hsa_circ_0003653 and hsa_circ_0003234, hsa-miR-200c-3p, hsa-miR-141-3p, ITGVA, TP53BP2) were further verified by RT-qPCR, showing a significant difference between URSA and NEP, which was consistent with the RNA-sequencing results (Figure 1G).

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

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 included CCL4, DDX58, CXCL10, CXCL9, MX1, CD44, RPS2, SOCS3, RPS3A, and CXCL11 (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).

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

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.