Key Gene and Functional Pathways Identified in Unexplained Recurrent Spontaneous Abortion Using Targeted RNA Sequencing and Clinical Analysis

Identifying the mechanisms underlying unexplained recurrent spontaneous abortion (URSA) can help develop effective treatments. This study provides novel insights into the biological characteristics and related pathways of differentially expressed genes (DEGs) in URSA. Nineteen patients with URSA and three healthy fertile women with regular menstruation (control group) were recruited. RNA was extracted from the two groups to determine the differential expression of immunoregulatory gene sequences. Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analyses were used to identify the biological functions and pathways of the identified DEGs. A protein-protein interaction (PPI) network was constructed using the STRING database. Furthermore, qRT-PCR and ELISA were performed to validate the differential expression of the hub genes. We also explored the regulatory mechanism of Th1/Th2 imbalance. A total of 99 DEGs were identified, comprising 94 upregulated and five downregulated genes. Through GO analysis, nine immune cell function-related clusters were selected, and genes with significant differential expression were primarily enriched in eight immune regulatory functions related to the KEGG signalling pathway. Subsequently, five hub genes (TLR2, CXCL8, IFNG, IL2RA, and ITGAX) were identified using Cytoscape software; qRT-PCR confirmed the differential expression among the hub genes, whereas ELISA revealed a significant difference in extracellular IFN-γ and IL-8 levels. The levels of Th1 (IFN-γ) and the Th1/Th2 ratio were higher in the peripheral blood of URSA patients than in control group patients. These findings suggest that the occurrence of URSA may be associated with the abnormal expression of some specific immunoregulatory genes involved in T-cell activation and differentiation. Among the identified DEGs, IFNG may play a key role in regulating maternal immune response. Although further validation is required, our data provide an important theoretical basis for elucidating the pathogenesis of recurrent spontaneous abortion.


INTRODUCTION
Recurrent spontaneous abortion (RSA), one of the most common complications of pregnancy, refers to two or more consecutive spontaneous abortions with the same partner (1). It occurs in up to 5% of women of reproductive age (2). RSA has a complex aetiology which includes genetic factors, endocrine disorders, anatomical and structural abnormalities, thyroid dysfunction, infectious diseases, coagulation mechanism disorders, and immune factors (3)(4)(5)(6)(7). Embryonic chromosomal abnormalities are a major cause of spontaneous abortion during the first trimester; multiple studies have shown that at least 50% of first trimester miscarriages are associated with embryonic chromosomal abnormalities (8)(9)(10). However, there remain 40% to 50% of patients with unexplainable miscarriage, referred to as unexplained recurrent spontaneous abortion (URSA) (11). Previous studies have suggested that RSA may be associated with uncharacteristic chromosomal microstructure, gene expression, and noncoding RNA expression variation (12)(13)(14).
Mammalian fertilisation as well as embryo implantation and developmental processes are regulated by hundreds of genes (15)(16)(17). Several genes involved in the immune response (IFNG, IL10, KIR2DS2, KIR2DS3, KIR2DS4, MBL, and TNF), coagulation (F2, F5, PAI-1, and PROZ), metabolism (GSTT1 and MTHFR), and angiogenesis (NOS3 and VEGFA) have been shown to be pathogenetically associated with RSA (18). Most of these genes are associated with overactive immune and inflammatory responses, hypercoagulability, and disturbed metabolic regulation that may contribute to the pathogenesis of RSA. Previous studies have shown that multiple factors involved in immune regulation play a key role in maintaining balance in immune maternal foetal tolerance (19). Indeed, a balance in the immune system influences successful pregnancies, and cytokines secreted by immune cells play important roles in that balance at different stages of implantation. Therefore, in this study, we investigate the potential factors associated with miscarriage in patients with URSA, such as abnormal expression of certain genes that are associated with immune dysfunction. We applied a targeted RNA sequencing chip to screen differentially expressed genes (DEGs) in the peripheral blood of URSA patients, connected to the clinical characteristics, trace the action pathways of the DEGs, and elucidate the underlying immunological mechanisms of URSA at the gene expression level.

URSA Group
In total, 19 patients with URSA attending the Guangdong Family Planning Hospital from March 2020 to September 2020 were enrolled in the study. The inclusion and exclusion criteria for the URSA group were as follows: (1) two or more consecutive pregnancy losses before 12 weeks of gestation; (2) the examination of tissues from aborted embryos suggested a normal karyotype; (3) both biological parents had normal chromosomal karyotypes; (4) no genital tract malformations; (5) no reproductive tract infections; (6) menstrual cycle, basic sex hormones, thyroid function, and fasting blood sugar were normal; (7) autoantibodies (antinuclear antibodies, anti-thyroid autoantibodies, anti-phospholipid antibodies) were negative; (8) routine semen examination of the father was normal, excluding those with anatomical, microbial, viral, hormonal, or genetic disease in both couples.

Control Group
Throughout the same period, three women who visited the hospital for routine physical examination were selected as the control group. The inclusion and exclusion criteria for these individuals were: (1) given birth at least once; (2) no history of adverse pregnancy events such as spontaneous abortion, stillbirth, or premature delivery; (3) no pregnancy complications such as gestational diabetes or preeclampsia; (4) normal menstrual cycle; (5) no personal or family history of autoimmune system diseases or metabolic diseases.
In this study, peripheral blood samples (10 ml) were collected in the morning during nonpregnancy and nonmenstrual periods after at least 8 h of fasting from all subjects, and 5 ml blood was added to anticoagulant tubes containing 2% ethylenediaminetetraacetic acid (EDTA), which were subsequently stored at −80°C. Another 5 ml blood was centrifuged at 1,760g for 10 to 15 min, and the separated serum was stored in tubes at −80°C.

Preparation and Hybridisation of Immunomodulatory Gene Chips
Total RNA was extracted with an QIAamp RNA Blood Mini Kit (Qiagen, Cat. No. 52304) following the manufacturer's protocol and its quality was first examined using an Agilent Bioanalyzer 2100 (Agilent, USA). The purified RNA was used as a template to synthesise cDNA using oligo (dT) as a primer. The experimental procedure was carried out in accordance with the instructions of the reverse transcription kit (Takara, Cat. No. RR036A). After the determination of nucleic acid concentration by using Qubit 3.0 fluorometer (ABI, USA), the targeted gene region transcriptome RNA expression sequencing hybridisation chip was used for RNA capture and target enrichment to the formation of sequence libraries, sequencing was performed on Ion PGM platform (ABI, USA). Further, 212 target genes (detailed in supplementary table) that were mainly related to biological processes-including proliferation, apoptosis, migration, and differentiation of T and B cells and natural killer cells (NK cells), as well as apoptosis and migration of macrophages-were custom designed on a hybridisation chip by Agilent. The entire process was performed in the Guangzhou Ardent Clinical Laboratory.

DEG Screening
Raw data were trimmed of adapter and primer sequences in reads hereafter low-quality bases using fastp software (version 0.19.7) (20) to obtain clean reads. Subsequently, aligned the clean reads to the reference human genome hg38 using HISAT software (version 2.1.0) (21). Genes expression was calculated using SringTie software (version 1.3.5) (21) based on the alignment. The expression levels of different samples were then normalised using DESeq2 (version 1.22.2) (22) and subjected to differential analysis; DEGs were screened using the "limma" package (23), and a volcano map of DEGs was drawn using the "ggplot2" package (24) to show the differential expression of DEGs; genes with a p-value ≤ 0.05, and an absolute value of log2 fold-change greater than 1 were screened as DEGs.

GO and KEGG Enrichment Analysis
To explain the biological process and molecular mechanism of DEGs, the biological functions and pathway enrichment were analysed using Gene Ontology (GO) and KEGG enrichment analyses. The enrichment and a series of functional annotation were performed on an online bioinformatics resource website with the name of the Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifcrf.gov; version 6.8). The functions of genes and proteins related to three main categories, biological process (BP), cytological component (CC), and molecular function (MF), are defined as described among a standard semantic vocabulary as a GO database, whereas biological interpretation of genome sequences and other high throughput data were contained in the KEGG integrated database (25). Cluster Profiler (version 3.10.1) (26) was used to perform GO and KEGG enrichment analyses and unsupervised hierarchical clustering of DEGs. Heat maps were used to show the expression patterns of DEGs between the case and control groups. A p-value < 0.05 was considered statistically significant in the enrichment in characteristic biological functions in GO annotation and among specific potential pathways in KEGG annotation.

PPI Network Analysis
Protein-Protein interaction (PPI) regulates several cellular processes including replication, transcription, translation, splicing, secretion, cell cycle, signal transduction, and intermediate metabolism. In this study, the PPI network of DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org; version 11.0b) (27) with an interaction score ≥ 0.4.

Hub Gene Selection and Analysis
The top five DEGs were defined as hub genes based on the ranking of nodes' score calculated with the maximum correlation criteria with the topological analysis method of maximal clique centrality algorithm via the CytoHubba plugin of Cytoscape (28,29).

Real-Time Quantitative PCR
Total RNA was extracted using a QIAamp RNA Blood Mini Kit (Qiagen, Cat. No. 52304) following the manufacturer's protocol. The purified RNA was used as a template to synthesise cDNA by using oligo (dT) as a primer. The experimental procedure was carried out in accordance with the instructions of the reverse transcription kit (Takara, Cat. No. RR036A). To validate the expression differences of hub genes between the URSA case and normal control groups, according to the results of the previous transcriptome expression sequencing of targeted genes, qRT-PCR was conducted using TB Green Premix Ex Taq ™ (Takara, Cat. No. RR820A) to detect the hub genes in a StepOnePlus ™ Real-Time PCR System with Tower (ABI, USA). Three assays were carried out for each sample. Data were analysed with the 2 −DDCt method using ACTB as the internal control. The primers used in qRT-PCR are listed in Table 1.

Enzyme Linked Immunosorbent Assay
All peripheral blood samples were collected with suited EDTA anticoagulant tube, stewing at 37°C for 30 min and centrifuged for 15 min at the 1,760g to remove supernatant.
ELISA kits were used to detect the concentrations cytometric factors (IL-8, IFN-g, IL2Ra and TLR2) via double antibody sandwich method. Purified capture antibodies of the cytometric factor were coated on a microtiter plate and formed a solid-phase antibody. The revealed cytometric factors in the plasma were captured using HRP-labelled (horseradish peroxidase-labelled) solid-phase antibody to form the antibody-antigen-enzyme labelling-antibody complex and coloured by adding substrate of TMB (3,3′,5,5′-Tetramethylbenzidine) after thorough washing. TMB was transformed into the colour blue under the catalysis of HRP enzyme, and yellow under the effect of acid.
The destination intensity of the colour is positively correlated with the concentration of the cytometric factors in the plasma.
We detected the absorbance (OD value) using enzyme-labelled instrument of Infinite M Plex microplate reader (Tecan, Männedorf, Switzerland) at wavelength of 450 nm to calculate the concentration through standard curve. ELISA kits (MLBio, China) used in this study provided the R value of the correlation coefficient between the linear regression of the sample and the expected concentration above 0.95, whereas the intra-assay and inter-assay coefficient of the variation is less than 10% and 15% respectively.

Retrospective Cytokine Data Analysis
Th1 (IFN-g), Th2 (IL-4), and Th1/Th2 ratio data were retrospectively analysed from 227 URSA patients and 18 healthy fertile women attending the Guangdong Family Planning Hospital from January 2014 to December 2018, using the same inclusion and exclusion criteria listed in 2.1.
The levels of Th1 (IFN-g),Th2 (IL-4), and Th1/Th2 cytokines were detected in fresh peripheral blood within 8 h using BD FACSCalibur flow cytometry (BD Biosciences, San Jose, CA, USA). Briefly, 150 ml of fresh peripheral blood sample within 8 h was placed in a flow tube and 150 ml of RPMI1640 (without calf serum FBS), 10 ml of 1 mg/ml PMA, 10 ml of 50 mg/ml ionomycin, and 10 ml of 0.5 mg/ml BFA were added to the working solution. The plates were incubated at 37°C in 5% CO 2 for 3 to 6 h. After addition of 10.0 ml CD3percp and 2 ml CD8apc, the mixture was incubated for 15 min in the dark at room temperature, followed by membrane disruption and the addition of IFN-g/IL-4 antibody. After incubation in the dark at room temperature for 20 min, and the supernatant was discarded after centrifugation and used for detection.

Statistical Analysis
SPSS 26 software was used for data analysis. The data of general characteristics and laboratory measurement results were tested for normality, and data with a normal distribution were tested for homogeneity of variance. Comparison of data between two groups was conducted using Student's t-test, and those that did not conform to normal distribution were tested with the rank test. Statistical significance was set at p < 0.05. GraphPad Prism6 software was used to generate histograms depicting the Th1 (IFN-g), Th2(IL-4), and Th1/Th2 ratios.

General Characteristics and Laboratory Measurements
A comparison of the basic clinical characteristics of the two groups is shown in Table 2. The ages of the subjects ranged from 27 to 35 years, with no significant differences in the mean age, weight, and various hormonal contrasts between the two groups (p > 0.05).

Identification of DEGs
The differential expression of immune regulated genes was selected by calculating the p value and fold-change value of the chip probes between the case and control groups. A total of 99 DEGs, including 94 upregulated (log2 fold-change > 1, p < 0.05), and five downregulated (log2 fold-change < -1, p-value < 0.05) DEGs were filtered, and the results are illustrated in Figure 1 as volcano plots.

GO and KEGG Enrichment Analyses of DEGs
GO analysis indicated that the 99 DEGs were mainly enriched in processes involved in the T cell activation and proliferation and regulation of leukocyte differentiation ( Table 3). The results of KEGG pathway enrichment analysis demonstrated that the DEGs were significantly enriched in RSA and immune regulated-associated signalling pathways of NK cell-mediated cytotoxicity and cytokine-cytokine receptor interaction ( Table 4).

PPI Network and Hub Gene Analysis
PPI network was analysed using the STRING database; 81 nodes and 205 edges with an average node degree of 5.06 were established ( Figure 2). Five hub genes (TLR2, CXCL8, IFNG, IL2RA and ITGAX) out of the 99 included were identified using CytoHubba according to the nodes' score ranking using the maximum correlation criteria and the topological analysis method of maximal clique centrality algorithm via the CytoHubba plugin of Cytoscape. The PPI network was analysed using the STRING database; the abridged version with the hub gene indicated is shown in Figure 3.

qRT-PCR Validation
Based on the results of the previous targeted transcriptome expression sequencing, RT-PCR was used to validate the high/ low expression of hub genes and housekeeping genes differentially expressed in the peripheral blood samples of cases versus controls. Differences in TLR2, CXCL8, IFNG, IL2RA, and ITGAX genes between the URSA and control groups were statistically significant (p < 0.05), which were consistent with the RNA sequencing results (Figure 4).

ELISA Validation
ELISA validation experiment was performed on revealing cytokines and receptors including IFN-g, IL-2Ra, IL-8, and TLR2 among the results of hub genes analysis to validate the revealing protein

Retrospective Analysis of Cytokines
Th1 (IFN-g), Th2(IL-4), and Th1/Th2 cell ratio data from 227 URSA patients who attended the Guangdong Family Planning Hospital from January 2014 to December 2018 and 18 healthy female controls showed that the Th1 (IFN-g) and Th1/Th2 ratios were significantly higher in the case group (p < 0.05), whereas the Th2 (IL-4) ratios were not significantly different ( Table 5 and Figure 6).

DISCUSSION
Among all RSA causes, the well-defined causes include genetic factors, endocrine factors, thrombosis, immune factors, and genital malformations, whereas at least 50% of RSA remain unexplained and categorised as URSA. It has been suggested that aberrant gene expression is a major cause of URSA and an important indicator of pregnancy disorders (30,31). The regulation of genes in the embryo includes the self-regulation of genes critical for embryonic development and the biparental genetic influence on embryonic development (32). During early pregnancy, the expression of maternal genes plays a crucial role in embryonic development. The involvement of genetic clusters in the occurrence of miscarriage warrants further investigations to improve our understanding of this complicated condition. It has been shown that immune factor dysfunction plays an important role in the pathogenesis of URSA (33). Utilising gene variation and expression to determine the differences in the expression of immune regulatory genes in the pathogenesis of URSA may help to accurately identify the potential risk and associated factors of URSA to ameliorate pregnancy outcomes.
In this study, 99 DEGs between the URSA and control groups were significantly enriched in T-cell activation and proliferation, and regulation of leukocyte differentiation. Meanwhile, most genes in the relevant clusters were significantly upregulated in differentially expression states, suggesting that significant upregulation of genes in relevant immunoregulatory clusters may induce immune dysfunction in URSA patients. According to PPI analysis, TLR2, CXCL8, IFNG, IL2RA, and ITGAX were among the top five immunoregulatory hub genes of DEGs.
KEGG pathway analysis revealed that the main DEGs in URSA patients, compared with normal pregnant women, were enriched in several signalling pathways associated with immune regulatory functions, including natural killer cell-mediated cytotoxicity, cytokine receptor interaction, and other signalling pathways. NK cells are present in the endometrium and synergise with T lymphocytes to generate immune tolerance at the maternal foetal interface, which is important for a successful pregnancy (34,35). Abnormal gene expression of cytokines and cytokine receptor signalling pathway suggest immune cell dysregulation in women with recurrent pregnancy loss.
The protein encoded by TLR2 is a member of the TLR family, which plays an important role in pathogen recognition and innate immune activation and is widely expressed during the maternal foetal interface, including immune cells, trophoblasts, and decidual cells. Multiple studies have shown that TLRs play a crucial role in the pathogenesis of autoimmune diseases, since TLR stimulation has been reported to induce foetal resorption during early pregnancy loss (36,37). TLR2, an initiation factor of the TLR signalling pathway, can increase the release of IL-8 through signalling factors such as MyD88 and is critical in recognising microbial infections and mediating innate immune responses (38).
IL-8, encoded by the CXCL8 gene, is an important mediator in the acute immune response, stimulation of chemoattraction, and promotion of angiogenesis, that may play an important role in URSA by regulating the maternal immune response (39). IL8 recruits neutrophils to the surrounding tissues thereby promoting the spread of acute inflammation and influences the levels of IL6 and IL1 expressed in the foetal membranes of women with chorionicity. Furthermore, IL-8 may control the chemotaxis and migration of immune cells at the site of inflammation and plays a significant role in embryo implantation and the establishment of pregnancy through a network cooperating with other cytokines (40). Antagonising CXCL8 and its receptor could reduce the inflammatory response and improve the symptoms of inflammation-related diseases (41). However, the role of IL-8 and other cytokines in pregnancy and their prognostic value for pregnancy outcome remains to be fully elucidated.
The IFNG gene is located on 12q15 and encodes IFN-g, which is a Th1 proinflammatory cytokine with multifaceted regulatory effects on immune responses in the body (42). Nakagawa et al.  showed that Th1 and Th2 cells play important roles in the immune response, especially in immune rejection and tolerance, and the imbalance of Th1/Th2 ratio is not conducive to maintaining normal pregnancy (43). Tangri et al. (44) demonstrated that Th1 type cytokines, such as IFN-g and TNF-a, are highly expressed in the placentas of pregnant mice with propensity to miscarriage, but not in normal pregnant mice. The immune response is biased toward a Th2 type during normal pregnancy, which protects the embryo from rejection by the maternal immune system, favouring the implantation of fertilised eggs and the development of the foetus (45). If a Th1/ Th2 imbalance occurs during pregnancy, a Th1 type bias with increased secretion of Th1 type cytokines and enhanced cellular immunity may damage the placental trophoblasts and foetus, thereby inducing URSA (46). Animal studies have shown that the key genes in the regulation of maternal immune responses and the maintenance of normal pregnancy is important, since IFNG can affect the immune response by regulating the expression of several immune-related genes, including CIITA, KYNU, IDO1, WARS, and MHC (47). Previous studies have reported the imbalance of Th1/Th2 in URSA females, besides, IFN-g, an indicator of Th1 bias, encoding genes have been previously screened in animal transcriptome studies. Conversely, in this study, we considered combining the imbalance of Th1/Th2 in URSA and the key genes in the regulation of maternal immune responses identified in previous animal studies to illustrate the significance of the potential indicators in the peripheral blood among the Th1/ Th2 differentiation and other signalling pathways.
In this study, TLR2, CXCL8, IFNG, IL2RA, and ITGAX genes were upregulated. To further verify the results of the RNA expression, we confirmed the high/low expression of the above genes using qRT-PCT which was consistent with the RNA sequencing results. We also detect these gene-regulated cytokines using ELISA. The levels of IFN-g and IL-8 were significantly higher in the URSA group than in control.
The levels of IL-2R and IL2Ra in the URSA group were higher than in control, but not significantly. Using flow cytometry, we showed that IFN-g was significantly increased in the peripheral blood of URSA patients. Our results also showed that the high expression of IFNG can promote the secretion of IFN-g both inside and outside of the cell. IFN-g, which is an important Th1 cytokine and an initiator of Th17 cell differentiation and the Th1/Th2 cell differentiation signalling pathway, regulates T cell activation through IFNG/STAT1/T-bet signalling as shown in Figure 7 (48). Meanwhile, T-bet factors within the nuclear envelope act through RUNX1-RORgt, which can inhibit the RUNX1 factor, a positive regulator of Th17 cell differentiation (49,50), by blocking Th17 cell differentiation and resulting in a Th17/Treg imbalance. Thus, the IFNG gene may cause maternal foetal Th1/Th2/Th17/Treg imbalance by regulating the release of IFN-g, thereby increasing the risk of miscarriage. To further confirm the contribution of the IFNG gene to the release of IFN-g and modulation of the Th1/Th2 balance, we investigated the expression of Th1 (IFN-g), Th2 (IL-4), and Th1/Th2 ratio in the peripheral blood of 227 URSA patients and 18 normal controls. The results show that the Th1 (IFN-g) and Th1/Th2 ratios were significantly higher in the URSA group, which suggests that the main cause of the Th1/Th2 imbalance is an increase in Th1 (IFN-g), possibly associated with the regulation of the IFNG gene. Our results suggest that IFNG may play an important role in regulating the maternal immune function.
Qiu et al. have reported that the imbalance of GATA-3/T-bet transcription factors may interfere with type 1 helper T cell (Th1)/type 2 helper T cell (Th2) differentiation at the foetalmaternal interface and possibly cause URSA via situ hybridisation (51). Furthermore, Wu et al. observed a Th1/Th2 imbalance and an increase in sTim-3 and Galectin-9 expression in the patients with URSA, which may be involved in the regulation of immunity during pregnancy (52). Recently, Chen et al. screened out ATP6V1G3, a key gene in RSA patients, via RNA-seq and qRT-PCR validation (25). Therefore, we have performed target region RNA-seq to analyse the expression pattern of the specific immunoregulatory related genes in the peripheral blood of URSA patients, and verified the results of RNA-seq by qRT-PCR. Based on the output of analysis, we summarised the clinical features of immunological disorders in URSA patients for comparison. The data from analysis of bioinformatics and clinical immunology showed that the CXCL8 and IFNG gene we screened out have been confirmed to play an important role in immune imbalance, supporting the reliability of our research. Furthermore, we used ELISA to evaluate the cytokines and their associated receptors from the results of the hub genes analysis to confirm the serum protein concentration between the control and case groups, which provided further evidence of indicators and sensors present in the peripheral blood among the Th1/Th2 differentiation and other signalling pathways. Our results identified a significant indicator, i.e., IFN-g and IL-8, and additional indicators that are not significant, i.e., IL-2Ra and TLR2 as a reference for future research.

CONCLUSIONS
In this study, several immune-related genes involved in T cell activation and proliferation and cytokine-cytokine receptor interaction were screened from URSA patients. However, gene expression was controlled by multiple factors, such as RNA splicing, regulation of transcription factors, ethylation of the promoter region, and histone modification. Therefore, determining the underlying mechanism of DEGs may contribute to the development of therapeutic strategies against URSA. Another possible limitation of our study is the relatively small sample size, and hence the conclusions cannot be generalised.
In summary, we explored the pathogenesis of URSA at the genetic level and screened the DEGs enriched in immune-related gene clusters while mapping specific immune-associated signalling pathways. The key genes and functional pathways identified in this study provide new insights into the molecular mechanisms involved in URSA pathogenesis and provide potential diagnostic and therapeutic targets. In vivo and in vitro validation of our findings and elucidation of the specific mechanisms of this cluster in URSA are warranted.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the Guangdong Family Planning Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
All authors contributed to the article and approved the submitted version.