Small RNA sequencing revealed aberrant piRNA expression profiles in deciduas of recurrent spontaneous abortion patients

Piwi-interacting RNAs (piRNAs) is a novel class of non-coding RNAs. However, changes in piRNA expression profiles in recurrent spontaneous abortion (RSA) have not yet been investigated. The aim of this study was to identify differentially expressed piRNAs in deciduas of RSA patients. Decidua tissues were collected by curettage from recruited RSA patients and normal early pregnant (NEP) women with their informed consent. Small RNA sequencing was used to evaluate the differences in piRNA expression profiles between RSA and NEP. The present results demonstrated that the counts of total piRNA reads in RSA samples were increased compared with those in NEP samples (0.21% vs. 0.11%). Differential expression analysis identified 29 upregulated piRNAs and 18 downregulated piRNAs in RSA samples. RT-qPCR further confirmed that the expression levels of uniq-109625, uniq-89328, uniq-50651 and uniq-4569 were decreased in 8 RSA tissues, compared with 13 NEP tissues. Otherwise, pi-22628 and uniq-173406 were increased in 8 RSA tissues. Based on GO term and KEGG pathway analysis, we speculate that these piRNAs regulate RSA by targeting extracellular matrix component pathway, cell adhesion pathway and focal adhesion pathway. PiRNAs may be involved in RSA pathogenesis by target genes function on adhesion and extracellular matrix component.


Introduction
Pregnancy is a complex process, comprising fertilization, implantation, organ and tissue differentiation, and fetal growth, which is effectively controlled by a number of both maternal and fetal factors (Qu et al., 2019;Verma et al., 2019a;Verma et al., 2019b). Recurrent spontaneous abortion (RSA) refers to three or more consecutive spontaneous abortions before 20 weeks of pregnancy and affects approximately 1%-2% of couples who attempt to have a child (Coomarasamy et al., 2015;Raj and Lesley, 2006). There are many factors leading to RSA, such as genetic factors, immune and environmental factors, endocrine disorders, infection and trauma, abnormal anatomy of the reproductive organs, cervical dysfunction, age, lifestyle problems and thrombosis tendency (Pei et al., 2019;Pereza et al., 2017;van den Berg et al., 2012). The causes that are already found cannot explain all the cases. The unsatisfactory clinical efficacy of some RSA cases has become a difficult problem to be solved in the field of gynecology and obstetrics. Therefore, more factors of progression of RSA need to be revealed.
Non-coding RNAs (ncRNAs) are generated from noncoding regions that were previously considered to be junk DNA, without the potential of translation into proteins (Goodrich and Kugel, 2009). NcRNAs, such as microRNAs (miRNAs) and long ncRNAs (lncRNAs), have been identified as crucial regulators of gene expression in RSA, endometriosis, pre-eclampsia, infertility and other reproductive system diseases (Dong et al., 2014;Dong et al., 2017;Qin et al., 2016;Sober et al., 2016;Zhang et al., 2019;Zhao et al., 2018). Piwi-interacting RNAs (piRNAs) are small regulatory RNAs expressed in mammalian germ line cells and have been identified as key players in germ line development. These piRNAs molecules, typically of length 26-32 nt, associate with Piwi proteins of the Argonaute family to form the Piwi-interacting RNA complex. Studies had shown that piRNAs generally exist in clusters on chromosomes, ranging in length from 20 kb to 100 kb, and the density of piRNAs ranges from 40 to 4000. PiRNA cluster was transcribed to form a long single-stranded precursor, and then processed to form mature piRNA of about 26-32 nt (Hirakata and Siomi, 2016;Huang et al., 2017). Although piRNAs were initially considered to be only expressed in germ cells, a growing number of studies identified piRNAs are also abnormally expressed in early embryos (Cui et al., 2018;Russell et al., 2016), indicating that piRNAs may be involved in early embryos development. PiRNA pathway plays a role in animal germ cells and controls TE (Transposable Element) activity. The core of piRNA pathway is a ribonucleo-protein complex composed of a small RNA called piRNA and a protein from the Argonaute nuclease PIWI subfamily. The piRNA pathway relies on the specificity provided by piRNAs to identify TES (Transposable Elements) targets, while the effector function is provided by piwi protein. The piwi-piRNA complex silences TEs at both the transcriptional level (by attracting inhibitory chromatin to modify genomic targets) and the post-transcriptional level (by shearing the transcripts in the cytoplasm) (Lee et al., 2011;Ozata et al., 2019). Damage to piRNA pathway leads to overexpression of TEs, significantly impairs genomic structure, and always leads to germ cell death and infertility (Toth et al., 2016).
In the present study, to improve the understanding of the biological function of piRNAs in pregnancy, the differences in piRNA expression profiles between normal control tissues and RSA decidua tissues were compared using second-generation deep sequencing for small RNAs. Subsequently, the results were validated by reverse transcription quantitative polymerase chain reaction (RT-qPCR). The potential clinical utility of pregnancy associated piRNAs was also assessed by analyzing the association of piRNAs with clinic opathological features of patients with RSA. To the best of our knowledge, the present results provided, for the first time, an overall change of piRNA expression profiles in RSA and an outlook into clinical applications of piRNAs as a clinic opathological target.

Patients and samples
In the current study, a total of 8 RSA patients and 13 normal early pregnancy (NEP) patients were recruited at the Family Planning Special Hospital of Guangdong Province between June 2016 and July 2018. Written informed consent was obtained from each RSA patient and control subject before participation in the study. None of the patients had a history of autoimmune conditions, metabolic disorders, hypertension, arterial or venous thrombosis, embryo chromosomal abnormalities, liver or kidney dysfunction, or uterine anomalies. Chromosomal abnormalities were excluded from all couples in the recurrent miscarriage group. 8 RSA patients (age 25.63 ± 4.87 years and gestational age), who had experienced at least three consecutive embryonic losses before the 12 th gestational week, were recruited in the RSA group.
Meanwhile, 13 clinically normal early pregnant women (age 26.46 ± 4.03 years and gestational age), which were terminated for non-medical reasons and undergoing legal abortions around the 6-12 th gestational week were recruited as the control. This study was approved by the Ethics Committee of Family Planning Special Hospital of Guangdong Province [No. 2017(19)]. Chromosome analysis was performed following standard methods for preparing Gbanded metaphases from cultured peripheral blood-derived lymphocytes. The chromosomal abnormalities results showed that there were 6 kinds of chromosome aneuploidy in villi of early pregnancy abortion, including trisomy 13, trisomy 16, trisomy 18, trisomy 21, trisomy 22 and X monomer. Other types of chromosome aneuploidy and chromosome terminal monomer or trisomy were also suggested. Tissues were collected immediately after operations. Tissues were divided into 50 mg in each EP tubes and store in −80°C freezer.
Small RNA library construction, sequencing and data analysis Total RNA used for deep sequencing was extracted from decidua tissues using TRIzol reagent according to the manufacturer's protocol (Invitrogen), respectively. The concentration of the total RNA was measured by NanoDrop (Thermo Scientific, Wilmington, DE, USA), and the RNA integrity was checked on Bioanalyzer2100 (Agilent, Santa Clara, CA). The library was constructed by using Small RNA Sample Pre Kit (Illumina) after the samples were qualified. The library was constructed by using the special structure of 3' and 5' ends of small RNA (complete phosphate group at 5' end and hydroxyl group at 3' end). The two ends of small RNA were directly attached to the adaptor with total RNA as the starting sample, and then reverse transcribed to cDNA. Subsequently, the target DNA fragments were separated by PAGE gel electrophoresis after PCR amplification, and the DNA library was recovered by gel cutting. After the construction of the library, the initial quantification was carried out by using Qubit 2.0 to dilute the library to 1 ng/μL, and then the insert size of the library was detected by Agilent 2100. After the insert size met the expectation, the effective concentration of the library was quantified accurately by RT-qPCR (the effective concentration of the library >2 nM) to ensure the quality of the library. After passing the library inspection, the different libraries were pooling according to the effective concentration and the requirement of the target offline data volume of HiSeq/MiSeq sequencing. PiRNA sequences of six species of humans, mouse, rats, drosophilas, zebra fish and platypuses were included in the piRNA Bank (http:// pirnabank.ibab.ac.in). As the species to be analyzed belongs to one of these six species, we could first conduct known piRNA analysis. The analysis method was to map the known piRNA with the sequenced reads, and then it used unmap reads of AP to predict, otherwise piRNA could be predicted directly from beginning. The method was referent to the literature (Zhang et al., 2011). The new known piRNAs that could not be matched in the bank are named beginning with 'uniq'. If the change was ≥2.0 or ≤0.5, the difference was recognized to be significant; the significance was maintained after correction for multiple testing.
Enrichment analysis of piRNA-generation gene GO term analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of piRNA-generation genes were carried out. Gene Ontology (GO, http://www. geneontology.org/) was the international standard classification system of gene function. GO term analysis method was GO seq. This method was based on Wallenius non-central hyper-geometric distribution, which was different from ordinary HP-geometric distribution in that the probability of extracting an individual from a category was different from that of extracting an individual from another category. This difference in probability was obtained by estimating the preference to gene length, which could more accurately calculate the probability of GO term being enriched by candidate genes.
KEGG enrichment was measured by Rich factor, Q value and the number of genes enriched in this pathway. Rich factor referred to the ratio of the number of differentially expressed genes in the pathway entry to the total number of genes in all annotated genes. Q value was P value after multiple hypothesis test correction. Pathway significant enrichment could determine the main biochemical metabolic pathways and signal transduction pathways in which candidate genes participated (Kanehisa et al., 2008).

RT-qPCR of piRNA
Total RNAs were exacted using QiagenReasy mini kit (Qiagen) and reverse transcribed using a Super Script II kit (Invitrogen) according to the manufacturer's description, respectively. Real-time quantitative PCR was carried out using TB Green Premix EX Taq II (Takara) and analyzed using AB Step One Plus System (Applied Biosystems, AB, USA). The primer pairs used for piRNAs and U6 are shown in Suppl. Table 2. PCR parameters were as follows: 95°C for 20 s, followed by 40 cycles of 95°C for 5 s and 57°C for 30 s. Relative expression level was calculated as 2-(CTU6-CTgene) (Livak and Schmittgen, 2001;Schmittgen and Livak, 2008). The piRNA level was normalized by the housekeeping gene, U6. All experiments were carried out in triplicate. Screen the differential expression of piRNAs following the criteria: P-value ≤ 0.05; fold change ≥ 2.

Statistical analysis
All values are presented as mean ± SD. Statistical comparisons between groups were analyzed by F-test followed by Student's t-test. A P-value < 0.05 was considered to indicate a statistically significant difference. Statistical analysis was performed using SPSS software version 17.0 (SPSS, Inc., Chicago, IL, USA).

Result
Unique and total reads analysis To evaluate the overall differences of small RNAs in RSA and control group (NEP), the common and specific reads were analyzed between the RSA library and NEP library, including the number of unique reads (types of reads) and total reads (Table 1). Small RNA sequencing yielded 1429676 unique reads with 26453565 total reads from these two libraries. Among these unique reads, only 67,866 (15.20%) unique reads were shared by the two libraries, whereas 746290 (47.80%) and 466074 (32.60%) unique reads were specific in the RSA library and NEP library, respectively (Table 2). These specific reads were determined to not be abundant, with only 0.95% of total reads in the RSA library and 0.89% of total reads in the NEP library (Table 2). These data indicate that RSA decidua tissues and NEP decidua tissues have diverse small RNA profiles and that specific small RNAs exhibit low expression levels.
Length distribution of small RNAs In general, the length of small RNAs ranged from 18-35 nt, and the peak of length distribution was beneficial to identify the classes of small RNAs. MiRNAs were concentrated in 21 or 22 nt, whereas piRNAs were concentrated in 26-32 nt . As depicted in Fig. 1, in the RSA and NEP libraries, the lengths of small RNAs were clustered in two ranges: 18-24 and 25-35 nt, and the most abundant cluster was 18-24 nt with a 22 nt peak point, the canonical length of miRNAs. The percentage of reads within the 18-24 nt cluster demonstrated no notable difference between these two libraries (57.77% vs. 58.55%), indicating that the abundance of miRNAs was not significantly changed in RSA. However, an increased 30-33 nt peak, primarily comprised of piRNAs, was determined in the RSA library compared with the NEP library (33.08% vs. 28.8%), indicating that the piRNA pathway is activated in RSA.

piRNA annotation
To analyze the differentially expressed piRNAs between RSA and NEP decidua tissues, piRNAs were first annotated and screened by mapping clean reads to the piRNA database. The results demonstrated that 321 unique reads in the NEP library and 356 unique reads in the RSA library, accounting for 0.06% and 0.04%, respectively, were aligned with piRNA sequences. Although the proportion of unique piRNA reads indicated no notable changes between these two libraries, the counts of total piRNA reads increased from 15265 (0.11%) in the NEP library to 27093 (0.21%) in the RSA library (Table 1), indicating that the overall expression of piRNAs was increased in RSA tissues, compared with NEP tissues, and further implying that piRNAs may be implicated in RSA pathogenesis.
Differentially expressed piRNAs in RSA decidua tissues Comparison of the expression levels of piRNAs in RSA and NEP decidua tissues is beneficial for understanding the roles of piRNAs in the pathogenesis of RSA. As depicted in Fig. 2A total of 47 differentially-expressed piRNAs were identified, of which 29 piRNAs were upregulated in RSA decidua tissues, whereas 18 piRNAs were downregulated in RSA decidua tissues. These differentially expressed piRNAs are detailed in Table 3.

RT-qPCR validation
Top two up-regulated piRNAs (pi-22628 and uniq-173406) and top four down-regulated piRNAs (uniq-109625, uniq-89328, uniq-50651 and uniq-4569) were selected for further validation in 8 RSA and 13 NEP decidua tissues by RT-qPCR. The relative expression levels of these six piRNAs, calculated using the 2 -ΔCq method, are depicted in Suppl. Tables 1 and 2. Consistent with the results from small RNA sequencing, the RT-qPCR results demonstrated that the expression levels of six piRNAs were consistent in RSA decidua tissues, compared with NEP decidua tissues (Fig. 3). Notably, the fold changes from the RT-qPCR results were more than those from small RNA sequencing, which may be due to the increased sensitivity of RT-qPCR.
Enrichment analysis of piRNA-generation genes GO term results showed that five categories had statistical significance. They were protein binding, positive regulation of biological process, anchoring junction, extracellular matrix component, and cell adhesion molecule binding (Fig. 4). Scatter map of candidate target genes was a graphical display of KEGG pathway enrichment analysis results (Fig. 5). We selected 20 pathway entries with the most significant enrichment to display in this figure (Fig. 5). In KEGG pathway enrichment analysis results, 25 genes of our sequencing results were involved in focal adhesion pathway with big rich factor and small Q-value (Q < 0.05). 15 genes were involved in extracellular matrix receptor interaction (ECM-receptor interaction) with the biggest rich factor and small Q-value (Q < 0.05). GO term results had similarity with KEGG pathway enrichment analysis results. Extracellular matrix component and cell adhesion molecule binding in GO term results were similar to focal adhesion and ECM-receptor interaction in KEGG pathway enrichment analysis results. The above results reveal that piRNAs may link to RSA through adhesion and extracellular matrix.

Discussion
It is well known that ncRNAs, as key molecules regulating gene expression, are widely distributed in various tissues (Kaikkonen et al., 2011;Sasaki et al., 2007;Sun et al., 2004). Aberrant expression of ncRNAs is associated with numerous human disorders (de Almeida et al., 2016;Esteller, 2011;Mestdagh et al., 2015). The rapid development of the second-generation deep sequencing technology has provided     an unprecedented platform to comprehensively analyze noncoding transcriptomes as well as reveal a number of novel non-coding transcripts in various tissues, organs and disease models. Furthermore, second-generation sequencing has high sensitivity, which can avoid some minor differentiallyexpressed RNAs being missed. In view of the advantages of second-generation sequencing and our shortage of funds, only 3 pairs of RSA tissues and NEP decidua tissues were used for deep sequencing of small RNAs (≤40 nt). The present data demonstrated that only 15.20% of small RNAs were shared by RSA and NEP decidua tissues, indicating that the expression patterns of small RNAs in RSA and NEP decidua tissues were If the P-value was smaller than 0.05, the difference was recognized to be significant; the significance was maintained after correction for multiple testing. notably different. These observations further implied that small RNAs may be involved in RSA pathogenesis.
To date, the expression patterns of miRNAs and lncRNAs in RSA, endometriosis, pre-eclampsia, infertility and other reproductive system diseases have been extensively investigated (Dong et al., 2014;Dong et al., 2017;Qin et al., 2016;Sober et al., 2016;Zhang et al., 2019;Zhao et al., 2018). However, piRNAs, as a class of newly identified small ncRNAs, and the expression patterns of piRNAs in RSA remain largely unknown, and the field remains in its infancy. RSA is defined as three or more abortions with the same partner in the first 20 weeks of pregnancy . Previous reports indicated that cell adhesion has a critical role in RSA (Agostinis et al., 2012;Quenby et al., 2007;Yurdakan et al., 2008). In addition, different types of recurrent miscarriage are associated with varying patterns of adhesion molecule expression in endometrium (Quenby et al., 2007). In women with RSA, the adhesion of trophoblasts to endothelial cells is inhibited. However, the pathogenetic mechanism which inhibits the adhesion of trophoblasts to endothelial cells remains unclear in women with RSA. In our study, 29 piRNAs were upregulated and 18 piRNAs were down-regulated in RSA, compared with normal controls. We carried out GO and  KEGG pathway analysis for the target genes of piRNAs detected. The most relevant two KEGG pathways were focal adhesion and ECM-receptor interaction (extracellular matrix component, ECM), and GO analysis also indicated that piRNAs regulated cell adhesion and extracellular matrix component. Coincidentally, GO and pathway analysis demonstrated that they also regulated protein binding via target genes. Therefore, we speculate that piRNA could affect cell adhesion and extracellular matrix component, and then regulate the occurrence of RSA.
RNA sequencing of the samples of RSA patients and normal controls provided us some piRNAs that may be differently expression by data analysis. Volcanic map and heat map can be used to infer the overall distribution of differential piRNA clusters, which shows 18 down-regulated cluster and 29 up-regulated clusters with statistical significance. As we used RT-qPCR to confirm the results, a total of six piRNAs were found to be up or down-regulated in the samples of RSA patients compared with those in the samples of NP controls. These results suggested that piRNA may involve in the progression of RSA, but the evidence needs to be further confirmed. PIWI Interacting RNA (piRNA) pathway is a conservative defense system that protects the genomic integrity of animal strains from harmful transposable element (Dennis et al., 2019;Ozata et al., 2019). PiRNA biogenesis and their potential roles as part of an epigenetic network are possibly involved in cancer progression. Moreover, potential strategies based on the use of piRNAs and PIWI proteins as diagnostic and prognostic biomarkers as well as for cancer therapeutics are revealed (Assumpcao et al., 2015). Transposable Elements (TEs) have the ability to replicate and insert new genome locations. In addition, aryl-hydrocarbon receptor deficiency may have different effects on testicular and ovarian development through a process involving piRNA-related proteins, piRNAs and transposon elements (Rico-Leo et al., 2016). Moreover, a growing number of studies have shown that piRNAs and PIWI proteins, which are abnormally expressed in various cancers and other diseases, may serve as novel biomarkers and therapeutic targets for disease diagnostics and treatment .
Collectively, to the best of our knowledge, the present study presented, for the first time, global piRNA expression profiles in RSA and NEP decidua tissues by deep sequencing for small RNAs. Based on the small RNA sequencing data, it was determined that the overall expression levels of piRNAs in RSA decidua tissues were increased or decreased, compared with NEP decidua tissues, implying that piRNAs may be involved in RSA pathogenesis.

Conclusion
In conclusion, the present study is the first to explore the piRNA expression profiles in decidua tissues of RSA. Based on bioinformatics analysis, we speculate that piRNAs may be involved in RSA pathogenesis by target genes function on adhesion and extracellular matrix component. These observations will provide a theoretical basis for piRNAtargeted therapeutic strategies for RSA.
Availability of Data and Materials: The datasets used or analysed during the current study are available from the corresponding author on reasonable request.
Author Contribution: WQ designed the study, drafted the manuscript, and coordinated and participated in every part of the study. JW participated in the design of the study and helped in the study and drafting of the manuscript. LH participated in the design and statistical analyses. XL, HN, YT participated in the design of the study and performed the PCR analyses. LZ, GS and YT helped to write and revise the manuscript and supplied important comments. All of the authors read and approved the final manuscript.
Ethics Approval: This study was performed in line with the principles of the Declaration of Helsinki. Approval was granted by Ethics Committee of Family Planning Research Institute of Guangdong Province (No. 2017-19).