This article is mentioned in:

Introduction

Currently, increasing families are facing infertility problems. It has been estimated that 186 million people are affected by infertility, imposing a heavy burden on both families and society (1). With the development of assisted reproductive technology (ART), especially IVF-embryo transfer (IVF-ET), the pregnancy rate has significantly increased during the past decades (2). However, there are still many factors that lead to implantation failure.

To establish a successful implantation, a developmentally competent embryo and a receptive endometrium are both essential (3). Endometrial receptivity (ER) changes with the menstrual cycle, and an embryo can only be implanted during a short period of time called the ‘window’, when ER is appropriate for implantation (4). ER is a critical factor of gestation and closely related to infertility. Lack of receptivity or an embryo transfer that is not synchronized with a receptive endometrium represents a major reason for implantation failure and the subsequent failure of assisted reproduction (5).

Circular RNAs (circRNAs) are a class of non-coding RNAs generated by back-splicing, forming a covalently closed loop (6). They are widely expressed in mammals, and many circRNAs have been identified (7). Recently, circRNAs were revealed to be modulators in numerous physiological and pathological processes. Specifically, circRNA-9119 was reported to regulate the development of ER in dairy goat by sponging miR-26a (8). However, the role of circRNAs in human infertility remains unclear and requires further investigation.

Acupuncture and moxibustion, which are important in traditional Chinese medicine, have gained increased popularity for treating infertility. A prevalence study that included eight community and academic infertility practices in the United States reported that 22% of their patients had tried acupuncture treatment for infertility (9). Acupuncture combined with embryo transfer can increase pregnancy and live birth rates in women undergoing in vitro fertilization. The pregnancy rate was 42.5% (34 out of 80 patients) in the acupuncture treatment patients, while only 26.3% (21 of 80 patients) pregnancy rate in patients without acupuncture (10). However, the molecular mechanism of acupuncture in treating infertility has not been determined.

The aim of the present study was to verify that acupuncture and moxibustion can be used to treat infertility by altering the expression of circRNAs. CircRNA sequencing and reverse transcription-quantitative (RT-q)PCR were conducted to examine circRNA differences between patients treated with and those treated without acupuncture and moxibustion. Then, the functions and effects of circRNAs were analyzed by bioinformatics methods, including Gene Ontology (GO) analysis and pathway analysis.

Materials and methods

Patients

From January 2016 to June 2017, 48 infertile women of childbirth age treated at the Reproductive Medical Center of Nanjing Jinling Hospital, were recruited into this study. These women suffered from various diseases that lead to infertility and were planning on receiving IVF-ET. To prepare the endometrium for IVF-ET, the women were randomly assigned into two groups: One received a regular endometrial preparation program and the other received acupuncture and moxibustion in addition to the regular endometrial preparation program, with 24 women in each group. All women were monitored for the development of follicles by transvaginal ultrasound. Endometrial samples were extracted on day luteinizing hormone (LH) +7. Acupuncture was applied at the points of guanyuan, zigong, zusanli, sanyinjiao, shenshu, and ciliao. Acupuncture and moxibustion were administered starting the first day of the two natural menstrual cycles before the in vitro embryo transfer cycle up to the implantation day, one time on every other day. Patients received acupuncture for 1 h at a time, including half an hour on the front and half an hour on the back. The clinical characteristics of the women in this study are presented in Table I. The follow-up data revealed that the pregnancy rate in the women treated with acupuncture and moxibustion (58.33%) was higher than that in the women who did not receive this treatment (41.67%). The present study was approved by the Ethics Committee of Nanjing Medical University (No. (2014)204). All women provided written informed consent.

Table I. General characteristics of the patients in this study.

Table I.

General characteristics of the patients in this study.

Terms	Case group (n=24)	Control group (n=24)	P-value
Age, years	32.75±0.40	31.46±0.84	0.17
Infertility years	5.75±0.45	5.125±0.45	0.33
BMI, kg/m2	21.21±0.48	20.93±0.32	0.62
Age of menarche, years	13.83±0.24	13.17±0.32	0.10
Menstrual cycle, days	51.78±16.28	36.22 ±4.39	0.37
Menstrual duration	5.67±0.44	5.56±0.53	0.87

[i] Differences between the two groups were compared using Student's t-test. P<0.05 indicated statistical significance. BMI, body mass index.

RNA extraction, library construction, and circRNA sequencing analysis

Six women in each group were used for circRNA sequencing, and 24 women in each group were used for real time PCR. Total RNA of endometrial samples from the case group and control group were extracted using TRIzol (Invitrogen; Thermo Fisher Scientific, Inc.), according to the manufacturer's instructions. RNA quality and quantity were evaluated with a Nanodrop2000 (Thermo Fisher Scientific, Inc.). Qualified RNA was processed for library construction using a TruSeq Stranded mRNA sample preparation kit and sequenced on a Hiseq™ 2500 platform (both from Illumina, Inc.) with a paired-end 150 bp read run.

Bioinformatics analysis of sequencing data

The FASTQ files generated by RNA sequencing were aligned to the USCS hg38.0 Homo sapiens reference genome with HISAT2 (11). The abundance of gene expression was calculated using fragments per kilobase of exon per million fragments mapped (FPKM). Differentially expressed (DE) mRNAs and circRNAs were screened using the DeSeq package (12) and EBSeq package (13) based on |log2 fold change (FC)|>1 and a false discovery rate (FDR) <0.05. Functional and pathway enrichment analyses for DEmRNAs and for the host genes of DEcircRNAs were predicted using the DAVID and GO Consortium platforms with a threshold of FDR <0.05.

miRNA prediction and construction of the CeRNA network

The miRNA candidates that target the DEGs were predicted using miRanda and RNAhybrid (14,15). The predicted diameters were set to energy <-30 and score >160 for miRanda prediction and energy <-30 for RNA hybrid prediction. Only the overlapped miRNAs predicted by the two databases were retained for further analysis. Cytoscape V3.0 (https://cytoscape.org) was used to construct the CeRNA network.

RT-qPCR for validation of DEmRNAs and DEcircRNAs

Total RNA of endometrial samples from 24 women in each group were extracted using TRIzol. The concentration and purity of RNA were assessed using a microspectrophotometer (Tiangen Biotech Co., Ltd.). First-strand cDNA was synthesized using the RevertAid First Strand cDNA synthesis kit (Thermo Fisher Scientific, Inc.). FastStart Universal SYBR Green Master mix was used to amplify cDNA on the QuantStudio 6 Flex Real-Time PCR System (Thermo Fisher Scientific, Inc.). The primers are presented in Table II. The 2−ΔΔCq method was used to determine the relative gene expression with the internal control of B2M (16).

Table II. The primers used for real-time PCR.

Table II.

The primers used for real-time PCR.

Name	Sequence
circ-RP11-446H18.5-F	5′-CCTGCCTTCCCTAGATATTGTCCT-3′
circ-RP11-446H18.5-R	5′-TCCTAATCTAGCCTCTCCTGGG-3′
circ-MIB1-F	5-AGTATCCAAGTGGCAATAGGCA-3′
circ-MIB1-R	5′-TAGTAATTCGGTAAAAGCGATGTCT-3′
circ-MLH3-F	5′-TGAGGACATTCAGGCTGCTT-3′
circ-MLH3-R	5′-ATGGTAGGTAGAAAGATGGTGAGA-3′
circ-TET2-F	5′-TGTCAAGACTCAATATGATTTCCCA-3′
circ-TET2-R	5′-GGCAAATGGGAGGTGATGGT-3′
circ-SFMBT2-F	5′-AATCCCTTATTGATGCTGCCA-3′
circ-SFMBT2-R	5′-TCGAACCAGTCAAGTCACGTATG-3′
circ-RSF1-F	5′-AAGGAGGTCAACAAGAACAAGGA-3′
circ-RSF1-R	5′-CATCAGAGGAGGGCGAAGG-3′
circ-BABAM1-F	5′-AGATGAGTTGGAAGCCAGCAGA-3′
circ-BABAM1-R	5′-TCCGTCAAGGAGAACTGGGG-3′
circ-BACH1-F	5′-AATCAGAAATTGAGAAGCTGGTTG-3′
circ-BACH1-R	5′-GTTGTCGGGAAGTTCAGTGGA-3′
circ-LPAR1-F	5′-TGTAATTTCACAGCCCCAGGT-3′
circ-LPAR1-R	5′-CAGATCGAAGTCATGCTAGGAGA-3′

[i] F, forward; R, reverse.

Statistical analysis

Real-time PCRs were repeated in triplicate. Data are expressed as the mean ± standard deviation (SD) and compared using Student's t test. P<0.05 was defined as the level of statistical significance.

Results

Several circRNAs are differentially expressed after acupuncture and moxibustion

Using high-throughput sequencing, 2,653 circRNAs were identified, and 86 circRNAs were significantly differentially expressed (|log2 FC|>1, FDR <0.05) between the acupuncture group and the control group (Fig. 1). Moreover, among the DEcircRNAs, 57 circRNAs were upregulated, and 29 were downregulated.

Figure 1. Cluster heatmap. The differential expression profile of circRNAs between the two study groups and the homogeneity within each group are revealed based on a hierarchical clustering analysis.

GO and pathway analysis of differential circRNAs

The roles of circRNAs were investigated by GO analysis within three different domains [biological process (BP), cellular component (CC), and molecular function (MF)]. The top five terms of each domain, ranked by -log10(P-value), are presented in Fig. 2A and Table III. The identified BP terms were chromatin modification, positive regulation of transcription from RNA polymerase II promoter involved in unfolded protein response, oxidative DNA demethylation, regulation of transcription from RNA polymerase II promoter in response to hypoxia, and regulation of smooth muscle cell differentiation. The identified CC terms were nucleoplasm, nucleolus, nucleus, histone acetyltransferase complex, and annulate lamellae. The identified MF terms were methylcytosine dioxygenase activity, chromatin binding, zinc ion binding, histone binding, and protein binding. The significantly enriched KEGG pathways, ranked by -log10(P-value), are shown in Fig. 2B and Table IV. The identified pathways were protein processing in endoplasmic reticulum, degradation of aromatic compounds, shigellosis, mTOR signaling pathway, bacterial invasion of epithelial cells, and prostate cancer.

Figure 2. GO and pathways related to differentially expressed circRNAs. (A) The top 15 terms of BP, CC, and MF. (B) The significantly enriched pathways (-log10P>1.301). GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

Table III. The top five (ranked by -log10P-value) BP terms, CC terms, and MF terms.

Table III.

The top five (ranked by -log10P-value) BP terms, CC terms, and MF terms.

GO ID	-log10(P-value)	Description	GO terms
GO:0016568	6.1914023	Chromatin modification	BP
GO:0006990	3.9159859	Positive regulation of transcription from RNA polymerase II promoter involved in unfolded protein response	BP
GO:0035511	3.6954187	Oxidative DNA demethylation	BP
GO:0061418	3.5660861	Regulation of transcription from RNA polymerase II promoter in response to hypoxia	BP
GO:0051150	3.5206085	Regulation of smooth muscle cell differentiation	BP
GO:0005654	4.7545709	Nucleoplasm	CC
GO:0005730	3.212923	Nucleolus	CC
GO:0005634	2.7622145	Nucleus	CC
GO:0000123	2.3794772	Histone acetyltransferase complex	CC
GO:0005642	2.0758729	Annulate lamellae	CC
GO:0070579	3.9552169	Methylcytosine dioxygenase activity	MF
GO:0003682	3.0894796	Chromatin binding	MF
GO:0008270	3.0872473	Zinc ion binding	MF
GO:0042393	2.8177814	Histone binding	MF
GO:0005515	2.8066047	Protein binding	MF

[i] GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

Table IV. The significantly enriched KEGG pathways [ranked by -log10(P-value)].

Table IV.

The significantly enriched KEGG pathways [ranked by -log10(P-value)].

PathwayID	-log10 (P-value)	Pathway term
PATH:04141	3.3411215	Protein processing in endoplasmic reticulum
PATH:01220	1.9245391	Degradation of aromatic compounds
PATH:05131	1.6033722	Shigellosis
PATH:04150	1.5776843	mTOR signaling pathway
PATH:05100	1.4296001	Bacterial invasion of epithelial cells
PATH:05215	1.3156425	Prostate cancer

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes.

circRNAs associated with infertility are differentially expressed after acupuncture and moxibustion

Seven DEcircRNAs (circ-SFMBT2, circ-MLH3, circ-RSF1, circ-BABAM1, circ-TET2, circ-BACH1, and circ-LPAR1) associated with infertility were selected to verify the sequencing results by RT-qPCR. The relative expression levels of the seven circRNAs are presented in Fig. 3. The variation tendencies of circ-SFMBT2, circ-MLH3, circ-RSF1, circ-BACH1, and circ-LPAR1, but not those of circ-BABAM1 and circ-TET2, were in accord with those observed in high-throughput sequencing. The expression levels of circ-SFMBT2, circ-BACH1, circ-RSF1, and circ-LPAR1 were significantly upregulated and that of circ-MLH3 was downregulated in samples from patients treated with acupuncture and moxibustion compared to levels in the control group.

Figure 3. Relative expression levels of the circRNAs. Expression of seven differentially expressed circRNAs was validated by reverse transcription- quantitative PCR. n=24, *P<0.05 vs. respective C. C, control; T, tumor.

Construction of the circRNA-miR-mRNA interaction network

To evaluate the potential functions of the DE circRNAs, upregulated circRNAs were selected, and a circRNA-miR-mRNA interaction network was constructed (Fig. 4). circBACH1 was associated with 35 downregulated miRs, followed by circSFMBT2 with 7 downregulated miRs and circLPAR1 with 3 downregulated miRs, indicating that circBACH1 may play crucial roles in the network. Additionally, six genes (ATHL1, HTR6, KIAA1755, PLIN4, PRRT2, SNX29) were related to all three circRNAs.

Figure 4. circRNA-miR-mRNA interaction network. circBACH1 was associated with 35 miRs, circSFMBT2 with 7 miRs, and circLPAR1 with 3 miRs. miRs, microRNAs.

Discussion

With the rapid development of high-throughput sequencing technology, the roles of many circRNAs are being clarified. circRNAs have been confirmed to be vital modulators in a large number of physiological and pathological processes. circFBLIM1 was reported to promote hepatocellular cancer progression by sponging miR-346 (7). circNT5E was reported to promote glioblastoma tumorigenesis (17). However, studies on the functions and biological characteristics have mainly involved cancers, and there are few studies on circRNAs regarding the endometrium. Although circRNA-9119 was confirmed to regulate the development of ER in dairy goat (8), studies on the effects of circRNAs on human ER remain very limited.

Using high-throughput RNA sequencing, 2,653 circRNAs were identified. Among the 86 DEcircRNAs (|log2 FC|>1, FDR <0.05), 57 were upregulated and 29 were downregulated between the acupuncture group and the control group. Nine circRNAs were selected for further verification, and three (circ-SFMBT2, circ-BACH1, and circ-LPAR1) were differentially expressed in the acupuncture and moxibustion group.

The potential regulatory role of circRNAs was investigated by conducting KEGG and GO analyses. The significantly enriched terms embryonic digestive tract morphogenesis (GO:0048557) and regulation of smooth muscle cell differentiation (GO:0051150) indicate the possibility of altering the endometrial thickness and influencing ER. Notably, estrogen and progestogen are two important hormones involved in regulating the growth and development of the endometrium (18,19). Estrogen and progestogen are both aromatic compounds, and the degradation of the aromatic compound pathway (PATH:01220) was significantly enriched. Analogously, protein processing in the endoplasmic reticulum pathway (PATH:04141) was also significantly enriched, and the progestogen receptor has been reported to influence the ER (20). Moreover mTOR was reported to play a part in the tumorigenesis and development of endometrial carcinoma (21), and notably, the mTOR signaling pathway (PATH:04150) was revealed to be significantly enriched. Thus, the degradation of the aromatic compound pathway, protein processing in the endoplasmic reticulum pathway, and the mTOR signaling pathway are most likely involved in the regulation of ER.

Six genes (ATHL1, HTR6, KIAA1755, PLIN4, PRRT2, SNX29) were associated with all three differential circRNAs (circ-SFMBT2, circ-BACH1, and circ-LPAR1). In particular, circBACH1 was associated with the most genes, indicating that it may play a dominant role in the regulatory network as a competitively endogenous RNA (ceRNA) and can serve as a biomarker of ER. Furthermore, PLIN4 was revealed to be associated with steroids (22), which included estrogen, indicating that PLN4 may be important for ER. In future studies, PLN4, can be assessed as a potential target for new medicines and treatments. Our results indicated that acupuncture and moxibustion affect ER via regulation of the expression of circRNAs. circBACH1 may influence ER as a ceRNA of miRs, and PLN4 may play a dominate role in ER.

Acknowledgements

Not applicable.

Funding

The present study was supported by the National Natural Science Foundation of China (81473767, 81873371, 81804179, 81403477, 81603674 and 81403481), the Jiangsu Provincial Science and Technology Project of the Traditional Chinese Medicine Bureau (YB2017005), and the Special Foundation for Clinical Medical Research of the Chinese Medical Association-Research and Development Project of Young Physicians of Reproductive Medicine (17020350704).

Availability of data and materials

The datasets used and analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

JS and LC designed the experiments. JS, LC, JC, XJ, YM, QL, LX, YG and YX performed the experiments. JS, JC and XJ analyzed the data; and JS, LC and JC prepared the manuscript. All authors read and approved the manuscript and agree to be accountable for all aspects of the research in ensuring that the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Ethics approval and consent to participate

The present study was approved by the Ethics Committee of Nanjing Medical University (no. (2014)204). All patients provided written informed consent.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References