- Open Access
Analyses of circRNA profiling during the development from pre-receptive to receptive phases in the goat endometrium
Journal of Animal Science and Biotechnologyvolume 10, Article number: 34 (2019)
Recent studies have revealed that noncoding RNAs play important regulatory roles in the formation of endometrial receptivity. Circular RNAs (circRNAs) are a universally expressed noncoding RNA species that have been recently proposed to act as miRNA sponges that directly regulate expression of target genes or parental genes.
We used Illumina Solexa technology to analyze the expression profiles of circRNAs in the endometrium from three goats at gestational day 5 (pre-receptive endometrium, PE) and three goats at gestational day 15 (receptive endometrium, RE). Overall, 21,813 circRNAs were identified, of which 5,925 circRNAs were specific to the RE and 9,078 were specific to the PE, which suggested high stage-specificity. Further analysis found 334 differentially expressed circRNAs in the RE compared with PE (P < 0.05). The analysis of the circRNA-miRNA interaction network further supported the idea that circRNAs act as miRNA sponges to regulate gene expression. Moreover, some circRNAs were regulated by estrogen (E2)/progesterone (P4) in endometrial epithelium cell lines (EECs) and endometrial stromal cell line (ESCs), and each circRNA molecule exhibited unique regulation characteristics with respect to E2 and P4.
These data provide an endometrium circRNA expression atlas corresponding to the biology of the goat receptive endometrium during embryo implantation.
In mammals, endometrial receptivity is essential for successful embryo implantation . The receptive endometrium (RE) is a spatial and temporal phenomenon known as the “window of implantation” [2, 3] that is an absolutely necessary part of the reproductive process. During the “window of implantation”, the endometrium undergoes remarkable changes in construction and function: endometrial stromal cells (ESCs) proliferate, endometrial epithelium cells (EECs) differentiate, and the endometrium acquires adhesive properties that enable embryo adhesion and subsequent invasion [4,5,6]. Many researchers have proposed that the failure of in vitro fertilization embryo transfer using good quality embryos is due to impaired RE [7, 8] or inappropriate epigenetic modifications . Therefore, more studies are urgently needed to explore the molecular regulation mechanisms of RE to improve the success rate of embryo implantation and the fecundity of females, which is one of the most economically important components in animal husbandry [10, 11].
The number of genes and proteins identified from previous studies that are involved in the establishment of RE has exponentially increased [12, 13]. The establishment of RE is a highly dynamic process that is post-transcriptionally regulated, and several epigenetically regulated genes in the endometrium have been identified . Recent studies have demonstrated that noncoding RNAs play important roles in the formation of RE, and we have previously identified 143 differentially expressed miRNAs in the RE compared to the PE . Moreover, another study demonstrated that hsa-miR-30b, hsa-miR-30d, and hsa-miR-494 play important roles in gene reprogramming during the formation of RE in humans .
With potential functions as competing endogenous RNAs or miRNA sponges, circRNAs contain miRNA response elements (MREs) that enable them to act as molecular sponges for miRNAs leading to the derepression of miRNA target genes, thereby influencing post-transcriptional regulation . The circRNAs from back-spliced exons have been identified as a naturally occurring family of noncoding RNAs in eukaryotes [17, 18]. In recent years, circRNAs have been the focus of international research, and a large number of circRNAs have been successfully identified in various cell lines and species [19,20,21,22]. Furthermore, studies have demonstrated that circRNAs are expressed in eutopic and normal endometrium , endometriosis , and endometrial cancer . However, there is currently no information available on the expression and function of circRNAs in the developing normal endometrium, or on the differences in circRNAs expression between the RE and PE.
Considering the universal expression of circRNAs and their key role in regulating gene expression, we hypothesized that circRNAs are expressed in the endometrium of dairy goats as potential regulators of RE at the post-transcriptional level. Focusing on these points, we investigated the circRNA expression profile in the endometrium at gestational day 5 (PE) and gestational day 15 (RE) using Illumina Solexa technology. In addition, we analyzed the differentially expressed circRNAs (DECs) in the RE and PE and conducted GO enrichment and KEGG pathway analyses of the hosting genes of DECs (hg-DECs). Then, we examined the regulatory relationships between circRNAs with miRNAs and corresponding mRNAs.
Ethics statement and sample collection
All experimental goats were provided by XinLongMen milk goat breeding Farm, Xi’an, Shaanxi Province, China. The goats were maintained according to the No. 5 Proclamation of the Ministry of Agriculture, P. R. China, and we confirm that all animal protocols and methods were approved by the Review Committee for the Use of Animal Subjects of Northwest A&F University. A total of 10 healthy, 24-month-old multiparous dairy goats (Xinong Saanen) were treated as be described previously [15, 26]. Endometrium samples from 2 goats at gestational day 5 (pre-receptive endometrium, PE) and 2 goats at gestational day 15 (receptive endometrium, RE) were obtained from the anterior wall of the uterine cavity. What’s more, the other 6 goats were treated to verify the results of RNA-Seq.
Total RNA was isolated and purified using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s procedure. The RNA amount and purity of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by Agilent 2100 with RIN number > 7.0. Approximately 5 μg of total RNA was used to deplete ribosomal RNA according to the manuscript of the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA). After removing ribosomal RNAs, the remaining RNAs were fragmented into small pieces using divalent cations under 94 °C. Then the cleaved RNA fragments were reverse-transcribed to create the cDNA, which were next used to synthesise U-labeled second-stranded DNAs with E. coli DNA polymerase I, RNase H and dUTP. An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR by the following conditions: initial denaturation at 95 °C for 3 min; 8 cycles of denaturation at 98 °C for 15 s, annealing at 60 °C for 15 s, and extension at 72 °C for 30 s; and then final extension at 72 °C for 5 min. The average insert size for the final cDNA library was 300 bp (±50 bp). At last, we performed the paired-end sequencing on an Illumina Hiseq 4000 (LC Bio, China) following the vendor’s recommended protocol. The RNA-Seq data was deposited in GEO (GSE85384).
Sequence and primary analysis
Prior to reference genome mapping, the raw data (raw reads) were processed to remove low quality reads. In this step, clean reads (valid data) were obtained by removing reads containing adapter sequences, the uncertain base number in a read was more than 5%, low quality reads (1, reads containing sequencing adaptors; 2, reads containing sequencing primer; 3, nucleotide with q-value < 20). All of the downstream analyses were based on high quality clean data.
RNA-Seq reads mapping
CircRNAs are thought to be predominantly produced by back-splicing reactions that covalently link the 3′ splice donor (3′ end) of an downstream exon to the 5′ splice acceptor (5′ end) of an upstream exon [27, 28]. Due to the rearranged exon ordering, specific algorithms are required to annotate these back-spliced exon events for circular RNA prediction. Firstly, the clean reads (valid data) were mapped to the goat reference genome (http:/www.gencodegenes.org/, GRCh38) using a publicly available program TopHat 2 (http://ccb.jhu.edu/software/tophat/index.shtml) , what aligns the reads to genomes using the ultra high-throughput short read aligner Bowtie2 , and then analyzes the mapping results to identify splice junctions between exons.
Secondly, the unmapped reads with TopHat but could be mapped with TopHat-Fusion (http://ccb.jhu.edu/software/tophat/fusion_index.shtml) on the same chromosome in a noncolinear ordering (back-spliced ordering) were extracted as candidate back-spliced junction reads. Last, the circRNAs were identified from back-splice junction reads based on structural and splice sites characters of circRNAs: (1) the two ends of splice sites must be GU/AG, (2) mismatch≤2, (3) back-splice junction reads ≥1, (4) if there were more than 2 splice sites, they were on the same chromosome and not too wide apart from each other more than 100 kb.
Identification and quantification of goat circRNAs
Cufflinks software was used to assemble circRNAs, estimate their abundances, and test for differential expression, and the unit of the measurement is Fragment Per Kilobase of exon per Million fragments mapped (FPKM)  in this study. Only the comparisons with “q-value” less than 0.01 and status marked as “OK” in the Cuffdiff output were regarded as showing differential expression.
GO enrichment and KEGG pathway analysis of circRNA-hosting gene
The hypergeometric test was applied to map all differentially expressed genes to terms in the Gene Ontology (GO) database (ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz) , is an international standard gene functional classification system . The corrected P-value< 0.05 was used as the threshold to find significantly enriched GO terms in the input list of hgDECs compared to their genomic background . In addition, GO DAG, what illustrated the relationships of the GO terms , were explained the using cytoscape 2.8.1 version .
KEGG helped researchers to better understanding the biological functions of genes based on large-scale molecular datasets (http://www.genome.jp/kegg/) . In this study, the corrected P-value< 0.05 was used as a threshold  to find significantly enriched KEGG terms in the input list of DEGs compared to their genomic background .
The interaction network analysis of circRNA-miRNA
CircRNAs are believed to contribute substantially to the competing endogenous RNAs network for the studies what reported it acted as miRNA sponges to regulate gene expression [16, 39]. In this study, two publicly available program, Targetscan 7.0 (http://www.targetscan.org/vert_72/ and miRanda (http://www.microrna.org/microrna/home.do), were used to predict the interaction of circRNA-miRNA. The potential regulated relationship between the DECs in the present study and the goat miRNAs were also predicted by Targetscan 7.0 and miRanda.
RNA preparation and RT-qPCR
After sequencing, the other 6 goats (3 goats at day 5, and 3 goats at day 15) were treated to verify the results of RNA-Seq using RT-qPCR. Total RNA from the goat endometrium were extracted using Trizol reagent (TaKaRa, Dalian, China), 1% agarose gel electrophoresis and spectrometry (A260/A280) ratio were detected to judge the RNA quality [40,41,42]. And then the RNA was incubated for 15 min at 37 °C with or without (mock) 3 U/μg of RNase R (Epicenter Biotechnologies, Chicago, USA). To quantify the amount of circRNA, cDNA was synthesized with the Prime Script RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) using random hexamers. In particular, the divergent primers annealing at the distal ends of circRNA and RT-qPCR were used to determine the abundance of circRNA. GAPDH was used as the reference, the relative expression levels of the mRNAs were calculated using the eq. N = 2-ΔΔCt. What’s more, the PCR products from the individual circRNAs were further identified by Sanger sequencing.
Cell culture and treatment
Caprine endometrial epithelial cell (EEC) and stromal cell (ESC) lines were kindly provided by Dr. Y.P. Jin (Northwest A&F University, Yangling, China) . Cells were cultured and maintained in Dulbecco′s Modified Eagle Medium (DMEM)-Hank F12 (Gibco, Invitrogen Corporation, Grand Island, NY, USA) supplemented with 10% charcoal stripped fetal bovine serum (FBS; Gibco, USA) and antibiotics (100 U/mL penicillin and 100 mg/mL streptomycin) at 37 °C in a 5% CO2 atmosphere . The cells were treated with different concentrations of estradiol (E2, Selleck, Catalog No. S1709) and progesterone (P4, Selleck, Catalog No. S1705) respectively (Additional file 1: Table S1), total RNAs were extracted from cells at 24 h. And then, the expression levels of some circRNAs were detected as mentioned above.
Overview of sequencing data
To systematically identify circRNAs and their changes in expression levels between the PE and RE in Xinong Saanen dairy goats, we purified and sequenced RNA using the Illumina paired-end RNA-Seq approach. First, we isolated total RNA from 3 PE tissues and 3 RE tissues, then we ribosomal RNA-depleted and RNase R-treated the total RNA, generating a total of 120 million paired-end reads of 250 (± 50) bp in length. This approach yielded 18 Gb of sequence, representing approximately seven times the size of the genome (2.66 Gb). The distribution and uniformity of the repeated samples were analyzed using principal component analysis (Additional file 2: Figure S1A), and a boxplot was used to identify outliers and evaluate the overall quality (Additional file 2: Figure S1B). The box-plot distribution of the log10 FPKM (Fragment Per Kilobase of exon per Million fragments mapped) values suggested that the median and quartile values among the samples compared in terms of differential expression were almost identical. Thus, these two samples were discarded for poor repeatability, and ~ 90 million clean reads (~14Gb) were obtained after the low-quality data were separated out and discarded. A detailed summary for each sample is provided in Table 1.
The reference genome mapping
The clean reads (valid data, Table 1) of every sample were aligned to the goat reference genome. TopHat-Fusion was used for the nonlinear alignment of the back-splice junction reads (Additional file 3). Moreover, circRNAs can arise from exons or introns, and these distinct species show independent modes of generation. Thus, the reads mapped to the genome were divided into three groups (exonic, intronic, and intergenic) based on their origin (Fig. 1a). Furthermore, approximately 23% circRNAs were more than 1 kb in length (Fig. 1b).
Identification of goat circRNAs
In this study, the circRNA number and the corresponding circRNA-hosting gene number in every sample are shown in Table 2, and a total of 21,813 circRNAs were identified (Additional file 4). More than 95% of the circRNAs (20,720 circRNAs) consisted of protein coding exons, and only 5% of the circRNAs (1,093) identified were circular intronic RNAs (ciRNAs). We annotated these circRNA candidates based on the location of the corresponding circRNA-hosting genes, and 5,457 corresponding circRNA-hosting genes were mapped to the goat genome.
Differential expression of circRNAs in the receptive and pre-receptive endometrium
The expression levels of circRNAs were normalized by FPKM value, and the FPKM distribution of circRNA in each sample is shown in Table 3. More than 80% of the circRNAs were detected at medium expression levels (FPKM), and approximately 2% of the circRNAs were found to be highly expressed in the PE and RE libraries, respectively.
The fold changes (log2(PE/RE)) and corresponding significance threshold of the P-value were estimated according to the normalized circRNAs expression levels. Based on the expression levels, the DECs between the PE and RE were considered significant when P < 0.05. A total of 334 circRNAs significantly differed in terms of FPKM levels between the PE and RE (Additional file 5). Overall, 257 differentially expressed circRNAs were downregulated, whereas 77 circRNAs were upregulated in the RE compared to the PE. Further analysis showed that 30 of the 334 DECs were co-expressed in the two physiological stages of the goat endometrium, whereas 65 DECs were specifically expressed in the RE and 239 were specifically expressed in the PE (Fig. 2, Additional file 6).
The ten most upregulated circRNAs in RE compared to PE with the ten highest FPKM values are shown in Table 4. The FPKM value of circRNA8077 was highest (FPKM value = 57.95) in the RE, followed by circRNA8071 (FPKM value = 53.65). Notably, five of the ten upregulated circRNAs in the RE were derived from cysteine-rich motor neuron 1 (CRIM1), revealing that one gene could produce multiple circRNAs. The downregulated circRNAs with the ten highest FPKM values in the PE are shown in Table 5. The highest was circRNA5540 (FPKM value = 138.64, 15.37-fold decrease) in the PE, followed by circRNA2859 (FPKM value = 86.99, 20.27-fold decrease), and circRNA3537 (FPKM value = 69.31). A heatmap of the Pearson’s correlation values and a dendrogram of the correlation between the 20 DECs are provided in Fig. 3.
GO and KEGG analysis of the circRNA-hosting genes
In this study, 21,813 circRNA were annotated to 5,457 circRNA-hosting genes, suggesting that a single gene might produce more than one circRNA. For example, in the goat endometrium, CRIM1 produced twenty circRNAs, PIK3R3 produced five circRNAs, OGDH produced seven circRNAs, and ALDH18A1 produced three circRNAs (Additional file 4). Considering that a given circRNA is spliced from a given linear transcript, it is likely that the function of the circRNA is related to the function of its hosting gene . Therefore, the circRNA-hosting genes were analyzed by running queries for each gene against the GO database, which provided information related to three ontologies: molecular function, cellular component, and biological process.
Notably, 243 host genes of the 344 DECs (hg-DECs) were annotated to 850 GO terms (Additional file 7). We focused our attention on the GO terms with P < 0.05 and found that the hg-DECs were categorized into 41 significant GO terms (Additional file 8). Out of the 18 terms that were significantly enriched in molecular functions (Fig. 4a), the most significantly enriched GO terms were positive regulation of canonical Wnt receptor signaling pathway (GO:0090263) with 4 hg-DECs annotated, followed by protein polyubiquitination (GO:0000209), and protein auto-ADP-ribosylation (GO:0070213). In the cellular compartment GO category, thirteen terms were significantly enriched (Fig. 4b). The most significantly enriched GO terms were microtubule (GO:0005874) with thirteen hg-DECs annotated, followed by centriole (GO:0005814), nuclear chromosome, telomeric region (GO:0000784), and hemidesmosome (GO:0030056). As for the biological processes, ten GO terms were significantly enriched and were related to various processes (Fig. 4c) such as transcription corepressor activity (GO:0003714), Rac GTPase binding (GO:0048365), protein binding (GO:0005515), and Rac GTPase activator activity (GO:0030675).
Moreover, extracellular matrix structural constituent (GO:0005201) and integrin-mediated signaling pathway (GO:0007229) were interesting because previous studies have proposed that integrins are morphological and biochemical biomarkers of endometrial receptivity that participate in a series of important processes in extracellular matrix (ECM) remodeling in the endometrium . The inclusion of the annotation for phosphatidylinositol-3,4,5-trisphosphate binding (GO:0005547) also interested us because the PI3K/AKT signaling pathway plays important roles in the development of the RE in mice and humans.
Various genes cooperate with each other to carry out biological functions . Accordingly, KEGG analysis of circRNA-hosting genes was conducted to further understand the biological functions of circRNAs. In this study, KEGG pathway annotation showed 2,378 circRNA-hosting genes annotated to 133 biological processes (Fig. 5, Additional file 9). Moreover, the 120 hg-DECs were significantly enriched in 18 KEGG pathways (P < 0.05), suggesting that these pathways may play important roles in the development of endometrial receptivity (Table 6). The results showed the most significant pathways were in cancer (ko05200) with 21 hg-DECs enriched, followed by prostate cancer (ko05215), melanoma (ko05218), pancreatic cancer (ko05212), and renal cell carcinoma (ko05211). These results indicated that cancer-related pathways were active in the development of a RE from the pre-receptive phase.
Target miRNA predictions for differentially expressed circRNAs
There is evidence of functionality based on the finding that some exonic circRNAs act as miRNA sponges that regulate the expression of linear protein-encoding RNA products via “mRNA trap” mechanisms . Two publicly available programs (Targetscan 7.0 and miRanda) were used to predict the circRNA-miRNA interactions in this study. Overall, 21,683 circRNAs were predicted as putative targets for the 436 goat miRNAs using miRBase, suggesting that a single miRNA might be regulated by more than 49 circRNAs (Additional file 10). Moreover, we observed that 308 DECs were predicted as targets of 141 differentially expressed miRNAs between the RE and PE . The competing endogenous RNA (ceRNA) networks were visualized by importing the above interactions into the Cytoscape software (Fig. 6). For example, miR-449a shared MREs with 16 circRNAs (Fig. 6a) and three miRNAs targeted circRNA8077 (Fig. 6b).
Validation of circRNA with RNase R and RT-qPCR
To verify whether the back-spliced events were indicative of true circular molecules, we examined the physical properties of these products using two independent methods. First, circRNAs are known to be strongly resistant to the exonuclease RNase R compared with linear RNAs. The circRNA candidates in this study were also resistant to RNase R treatment and the abundance of linear RNAs decreased, thus confirming the circularized characteristics of these molecules. Second, because the circRNAs have a special head-to-tail junction site that differs from canonical linear splicing transcripts, RT-qPCR was used to validate circRNA candidates via the targeted amplification of the head-to-tail junction region unique to circRNA. To this end, outward-facing primers (Additional file 11) were designed against 11 circRNA candidates with high expression levels, followed by standard Sanger sequencing. Combining the two approaches, the results showed that the eleven circRNAs were consistent with the rRNA-depleted sequencing data, except circRNA8077 (Fig. 7).
Correlation analysis of circRNAs and their hosting genes
To validate the global downregulation of circRNAs in goat endometrial tissues, we selected the upregulated circRNAs with the ten highest FPKM values in the RE and the downregulated circRNAs with the ten highest FPKM values in the PE. To better evaluate the relationship between the circRNAs and their hosting genes, data on mRNA expression in the endometrial tissue from goats at PE and RE phases were used for the correlation analysis. A strong correlation was observed between the circRNAs and their hosting genes in the PE (R = 0.596, P = 0.006, Fig. 8a) and RE (R = 0.764, P < 0.001, Fig. 8b). These results indicated that circRNAs and their host genes may be co-regulated in the goat endometrium.
The effect of E2 and P4 on the expression levels of circRNAs in goat EECs and ESCs
The endometrium is composed mainly of two types of cells, ESCs and EECs, and these cells are regulated by E2 and P4. In EECs, the levels of circRNA3537, circRNA3165, circRNA3314, circRNA8071, circRNA8072, and circRNA8074 were up-regulated by 100 nmol/L E2, while circRNA8077 and circRNA9422 levels were not changed (Fig. 9a). In ESCs, the levels of circRNA3537, circRNA3165, circRNA3314, and circRNA9422 were up-regulated by 1 nmol/L E2, circRNA8072 was down-regulated and circRNA8074 was up-regulated by 10 nmol/L E2, and circRNA8071 and circRNA8077 levels were not changed (Fig. 9b). In EECs, the levels of circRNA3537, circRNA3165, and circRNA3314 were up-regulated by 100 nmol/L P4 (Fig. 9c). The levels of circRNA3537, circRNA3165, circRNA3314 and circRNA8074 were also up-regulated by 1 nmol/L P4 in ESCs. After 100 nmol/L P4 treatment, circRNA8072 levels were up-regulated, circRNA9422 levels were downregulated, and circRNA8071 and circRNA8077 levels were not changed (Fig. 9d).
This study yielded 18 Gb sequence for every sample using RNA-Seq, representing approximately seven times the size of the goat genome. After separating out and discarding low-quality sequences, we obtained ~ 90 million clean reads, and the sequencing depth in this study was sufficient to detect the transcripts expressed at low levels .
A recent study reported 65,731 human circRNAs in brain tissues , and another study further detected at least 27,296 circRNA candidates in six normal tissues and seven cancer tissues . Compared with human studies, fewer circRNAs (21,813 distinct circRNA candidates) were identified in goat. This could be due to the fact that the human genome is better annotated than the goat genome. Moreover, the human circRNAs were derived from various tissues , but the circRNAs in the present study were only obtained from two physiological stages of the goat endometrium. In addition, the expression of circRNAs was developmentally regulated and varied between tissues [22, 47, 48]. Notably, 5,925 circRNAs were specifically identified in the RE, and 9,078 were specifically identified in the PE, but only 6,810 circRNAs were co-expressed at both stages, exhibiting developmental stage-specific expression consistent with previous studies [49, 50]. This indicates that these molecules may play crucial roles in goat endometrium development.
Because a circRNA comes from a given linear transcript, it is likely that the circRNA function is associated with its hosting gene function . CRIM1, the hosting gene of circRNA8077, plays an essential role in activating the expression of vascular endothelial growth factor and ECM , which are important during early pregnancy [52, 53]. CRIM1 was also the hosting gene of circRNA8072, a markedly upregulated circRNA with a high FPKM value. Importantly, we found that half of the upregulated circRNAs with the ten highest FPKM values in the RE were derived from CRIM1. This indicated that multiple circularized exons could be produced from a single gene locus, which is consistent with previous literature [54, 55]. Collectively, these results suggested that CRIM1 might play an important regulatory role in the endometrial receptivity of dairy goats by generating alternative circularization products and extending the complexity of post-transcriptional regulation. However, further study is required as the present study provides no direct evidence that CRIM1 and its circRNAs participated in the formation of the RE.
The downregulated circRNA with the highest FPKM value in the PE was circRNA5540, showing 15.37-fold decreased expression in the RE. Its hosting gene, OGDH, was one of the top ten transcripts with the highest expression values in the bovine endometrium . This gene determines the metabolically active endometrial phenotype, which plays an important role in the formation of RE . In addition, a previous study showed that OGDH mRNA expression increased 2.15-fold  in the RE compared to the PE. Therefore, based on both the previous and present data, we propose that OGDH contributes to the formation of endometrial receptivity by decreasing its noncoding transcript, circRNA5540, thereby increasing its mRNA level. However, this hypothesis should be validated under well-controlled conditions in animal models.
Gene Ontology (GO) provides a simple and quick way to understand the properties of genes and gene products, and it is widely used as a classification system for gene function . Extracellular matrix structural constituent (GO:0005201) and integrin-mediated signaling pathway (GO:0007229) came from the GO analysis of hg-DECs in this study. These pathways interested us because integrins have been implicated as a biomarker of RE and participate in the remodeling of the ECM in the endometrium during early pregnancy . A previous study reported that a number of genes that encoded ECM proteins were enriched in the bovine endometrium during the estrous cycle . Importantly, all of these endometrial events were mediated by mechanisms of proliferation, differentiation, and apoptosis. PTEN (phosphatase and tensin homolog deleted on chromosome ten) antagonized PI3K (phosphatidylinositol 3-kinase) activity and negatively regulated serine/threonine kinase Akt , which modulated the activity of a variety of downstream proteins that were related to cell survival and proliferation . Therefore, the inclusion of the annotation for phosphatidylinositol-3,4,5-trisphosphate binding (GO:0005547) was also interesting to determine the important roles of the PI3K/Akt signaling pathway in the endometrium.
KEGG pathway analysis is widely used to understand the interaction between a set of genes in the context of biological functions [38, 61]. The KEGG pathway analysis of hg-DECs in this study revealed that several genes might be involved in several cancer-related pathways and regulate the development of the endometrium though their cognate circRNAs. There were striking similarities between the behavior of placental cells during the “window of implantation” and that of cancer cells diffusion . An appreciation of the cancer mechanisms may lead to a better understanding of the maternal mechanisms that control embryo implantation. Moreover, a previous study proposed that the cellular mechanisms used by trophoblast cells during implantation are reused by cancer cells to invade and spread within the body . Integrins, matrix metalloproteinase, ECM, and angiogenesis are common features of both implantation and the spread of cancer [63, 64].
Recent studies have confirmed that some circRNAs exert their biological effects by competitively binding miRNA, thus releasing and increasing the expression levels of the miRNA target genes [65,66,67]. This miRNA sponge mechanism has been clearly demonstrated for circRNAs ciRS-7/CDR1as and circRNA Sry [16, 28]. In the present study, 99.41% (21,683) of the identified circRNAs were predicted as putative targets for 436 goat miRNAs, suggesting that a miRNA was regulated by more than one circRNA , and further analysis showed that a circRNA could serve as a sponge for multiple miRNAs. However, further research is needed as effective miRNA sponging by exonic circRNAs may be relatively unusual or may not require a large number of miRNA binding sites . Moreover, we observed that 308 DECs were predicted as targets of 141 differentially expressed miRNAs between the RE and PE , suggesting that these molecules play important roles in the formation of endometrial receptivity. Moreover, 16 DECs were predicted to be regulatory factors of miR-449a that were significantly upregulated in the RE . Surprisingly, circRNA8072, circRNA8073, and circRNA8074 were significantly upregulated in the RE and were also predicted to be targets of miR-449a. Further study showed that circRNA8073 could bind miR-449a via the target site and inhibit miR-449a activity, creating a negative feedback relationship between circRNA8073 and miR-449a in EECs. Furthermore, circRNA8073 could increase the expression levels of centrosomal protein55 (CEP55) by sponging miR-449a in EECs in vitro. Moreover, circ-8073/miR-449a/CEP55 could regulate EECs proliferation and apoptosis via the PI3K/AKT/mTOR pathway. Taken together, circRNA8073 could regulate CEP55 by sponging miR-449a to promote EEC proliferation via the PI3K/AKT/mTOR pathway .
In addition, circRNAs substantially contributed to total RNA  and circular isoforms typically accounted for 5–10% of the total transcripts of their corresponding coding genes. However, certain circRNAs were up to 200 times more abundant than their linear counterparts . In this study, we found no distinct difference in the expression level of circRNAs and their hosting gene mRNA in the PE, but the expression levels of the hosting genes were higher than the levels of circRNAs in the RE. We suspect that this occurs because circRNAs comprise one to several coding exons of linear mRNAs and range between a few hundred and thousands of nucleotides in length . However, circRNAs aligned with antisense, intronic, or intergenic sequences, or untranslated regions have also been described. Nevertheless, circular RNA transcripts were generally expressed at low levels compared with linear RNAs , and our results support this observation.
CircRNAs are RNA molecules with covalently joined 3′ and 5′ ends formed by back-splice events, thus presenting as covalently closed continuous loops . Both mechanisms of “exon skipping” and “direct back-splicing” involve back-splicing by the canonical spliceosome , and a majority of apparent back-splice sequences in RNA-Seq data derive from circRNAs . We validated ten of eleven (91%) circRNA candidates by the targeted amplification of the head-to-tail junction region unique to the circRNAs using reverse transcription-PCR with an outward-facing primer and standard Sanger sequencing. Furthermore, some circRNAs were regulated by E2/P4 in EECs and ESCs, and each circRNA molecule exhibited unique regulation characteristics with respect to E2 and P4 in vitro.
In this study, high-quality circRNA expression profiles were obtained from the endometrial tissues of dairy goats. We identified 21,813 circRNAs, 334 of which were differentially expressed between the PE and RE. GO and KEGG analyses, the interaction network analysis of circRNA-miRNA, and the correlation analysis of circRNAs and their hosting genes could contribute to a better understanding of how circRNAs mediate the regulation of target genes in the development of endometrium receptivity. Additional studies showed that some circRNAs were regulated by E2/P4 in endometrial cells. Taken together, these findings reveal a new level of diversity in the goat transcriptome and introduce new insights into their regulation during the development of endometrial receptivity.
Aldehyde dehydrogenase 18 family member A1
Cysteine-rich motor neuron 1
Directed acyclic graph
Differentially expressed circRNAs
Endometrial epithelium cell lines
Endometrial stromal cell line
Fragment Per Kilobase of exon per Million fragments mapped
Hosting genes of differentially expressed circRNAs
Phosphatase and tensin homolog deleted on chromosome ten
Simmons DG, Kennedy TG. Uterine sensitization-associated gene-1: a novel gene induced within the rat endometrium at the time of uterine receptivity/sensitization for the decidual cell reaction. Biol Reprod. 2002;67:1638–45.
Xia H-F, Jin X-H, Song P-P, Cui Y, Liu C-M, Ma X. Temporal and spatial regulation of miR-320 in the uterus during embryo implantation in the rat. Int J Mol Sci. 2010;11:719–30.
Charlet-Renard C, Goffin F, Foidart M, Geenen V. The implantation window. J Gynecol Obstet Biol Reprod. 2002;31:440–55.
Pan Q, Luo X, Toloubeydokhti T, Chegini N. The expression profile of micro-RNA in endometrium and endometriosis and the influence of ovarian steroids on their expression. Mol Human Reprod. 2007;13:797–806.
Altmäe S, Martinez-Conejero JA, Esteban FJ, Ruiz-Alonso M, Stavreus-Evers A, Horcajadas JA, et al. MicroRNAs miR-30b, miR-30d, and miR-494 regulate human endometrial receptivity. Reprod Sci. 2013;20:308–17.
Revel A, Achache H, Stevens J, Smith Y, Reich R. MicroRNAs are associated with human embryo implantation defects. Human reproduction. 2011;26(10):2830–40.
Macklon NS, Stouffer RL, Giudice LC, Fauser BC. The science behind 25 years of ovarian stimulation for in vitro fertilization. Endocr Rev. 2006;27:170–207.
Sha AG, Liu JL, Jiang XM, Ren JZ, Ma CH, Lei W, et al. Genome-wide identification of micro-ribonucleic acids associated with human endometrial receptivity in natural and stimulated cycles by deep sequencing. Fertil Steril. 2011;96:150–5.
Kao LC, Germeyer A, Tulac S, Lobo S, Yang JP, Taylor RN, et al. Expression profiling of endometrium from women with endometriosis reveals candidate genes for disease-based implantation failure and infertility. Endocrinology. 2003;144:2870–81.
Zhang CY, Chen SL, Li X, Xu DQ, Zhang Y, Yang LG. Genetic and phenotypic parameter estimates for reproduction traits in the Boer dam. Livest Sci. 2009;125:60–5.
Khanum SA, Hussain M, Kausar R. Assessment of reproductive parameters in female dwarf goat (Capra hircus) on the basis of progesterone profiles. Anim Reprod Sci. 2007;102:267–75.
Koler M, Achache H, Tsafrir A, Smith Y, Revel A, Reich R. Disrupted gene pattern in patients with repeated in vitro fertilization (IVF) failure. Hum Reprod. 2009;24(10):2541–8.
Kuokkanen S, Chen B, Ojalvo L, Benard L, Santoro N, Pollard JW. Genomic profiling of microRNAs and messenger RNAs reveals hormonal regulation in microRNA expression in human endometrium. Biol Reprod. 2010;82:791–801.
Munro S, Farquhar C, Mitchell M, Ponnampalam A. Epigenetic regulation of endometrium during the menstrual cycle. Mol Hum Reprod. 2010;16:297–310.
Song Y, An X, Zhang L, Fu M, Peng J, Han P, et al. Identification and profiling of microRNAs in goat endometrium during embryo implantation. PLoS One. 2014;10(4):e0122202.
Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, et al. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495:384–8.
Jeck WR, Sharpless NE. Detecting and characterizing circular RNAs. Nat Biotechnol. 2014;32:453–61.
Lasda E, Parker R. Circular RNAs: diversity of form and function. RNA. 2014;20:1829–42.
P G, P P, N R. circBase: a database for circular RNAs. Rna:new York, Ny. 2014; 20: 1666–70.
Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, et al. Circular RNAs in the mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell. 2015;58:870–85.
Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15:1–14.
You X, Vlatkovic I, Babic A, Will T, Epstein I, Tushev G, et al. Neural circular RNAs are derived from synaptic genes and regulated by development and plasticity. Nat Neurosci. 2015;18:603–10.
Xu XX, Jia SZ, Dai Y, Zhang JJ, Li XY, Shi JH, et al. Identification of circular RNAs as a novel biomarker for ovarian endometriosis. Chinese Med J. 2018;131:559–66. https://doi.org/10.4103/0366-6999.226070.
Shen L, Zhang Y, Zhou W, Peng Z, Hong X, Zhang Y. Circular RNA expression in ovarian endometriosis. Epigenomics. 2018;10:559–72.
Chen BJ, Byrne FL, Takenaka K, Modesitt SC, Olzomer EM, Mills JD, et al. Analysis of the circular RNA transcriptome in endometrial cancer. Oncotarget. 2018;9:5786–96.
Zhang L, An XP, Liu XR, Fu MZ, Han P, Peng JY, et al. Characterization of the transcriptional complexity of the receptive and pre-receptive endometria of dairy goats. Sci Rep. 2015;5:14244.
Starke S, Jost I, Rossbach O, Schneider T, Schreiner S, Hung LH, et al. Exon circularization requires canonical splice signals. Cell Rep. 2015;56:103–11.
Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M, et al. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56:55–66.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:295–311.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.
Consortium GO. The gene ontology (GO) database and informatics resource. Nucleic Acids Res. 2004;32:D258–D61.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Choi H, Jo Y, Lian S, Jo KM, Chu H, Yoon JY, et al. Comparative analysis of chrysanthemum transcriptome in response to three RNA viruses: Cucumber mosaic virus , Tomato spotted wilt virus and Potato virus X. Plant Mol Biol. 2015;88:233–48.
Vinay Sharma MKS, Swami AK, Sarin R. Identification of drought responsive proteins using gene ontology hierarchy. Bioinformation. 2011;8:595–9.
Guo Q, Ma X, Wei S, Qiu D, Wilson IW, Wu P, et al. De novo transcriptome sequencing and digital gene expression analysis predict biosynthetic pathway of rhynchophylline and isorhynchophylline from Uncaria rhynchophylla, a non-model plant with potent anti-alzheimer’s properties. BMC Genomics. 2014;15:676.
He H, Liu X. Characterization of transcriptional complexity during longissimus muscle development in bovines using high-throughput sequencing. PLoS One. 2013;8:e64356.
Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–D4.
Jens M. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495:1439–49.
Zhang Y, Zhang XO, Chen T, Xiang JF, Yin QF, Xing YH, et al. Circular Intronic long noncoding RNAs. Mol Cell. 2013;51:792–806.
Zhang L, Wang Y, Fu M, Li G, An N, Li S, et al. The effects of ovariectomy on meat performance and expression of GH/IGF-I in young goats. Can J Anim Sci. 2014;94:619–26.
Zhang L, Wang Y, Zhou Z, Fu M, Li G, Peng F, et al. Fatty acid composition and mRNA expression levels of lipid-metabolic genes in the muscles of ovariectomised young goats. Anim Prod Sci. 2015;56:1585–92.
Zhang YY, Wang AH, Wu QX, Sheng HX, Jin YP. Establishment and characteristics of immortal goat endometrial epithelial cells and stromal cells with hTERT. J Anim Vet Adv. 2010;9:2738–47.
Qi X, Qu Y, Nan Z, Jin Y, Zhao X, Wang A. Caprine endometrial stromal cell line modulate the effects of steroid hormones on cytokine secretion by endometrial epithelial cells in vitro. Reprod Biol. 2012;12:309–15.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Baren MJV, et al. Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms. Nat Biotechnol. 2010;28:511–5.
Zheng Q, Bao C, Guo W, Li S, Chen J, Chen B, et al. Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs. Nat Commun. 2016;7:11215.
Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO. Cell-type specific features of circular RNA expression. PLoS Genet. 2013;9:119–29.
Fan X, Zhang X, Wu X, Guo H, Hu Y, Tang F, et al. Single-cell RNA-seq transcriptome analysis of linear and circular RNAs in mouse preimplantation embryos. Genome Biol. 2015;16:1–17.
Zhang C, Wu H, Wang Y, Zhu S, Liu J, Fang X, et al. Circular RNA of cattle casein genes are highly expressed in bovine mammary gland. J Dairy Sci. 2016;99:4750–60.
Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, et al. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365:141–8.
Nakashima Y, Takahashi S. Induction of cysteine-rich motor neuron 1 mRNA expression in vascular endothelial cells. Biochem Biophys Res Commun. 2014;451:235–8.
Kaloglu C, Onarlioglu B. Extracellular matrix remodelling in rat endometrium during early pregnancy: the role of fibronectin and laminin. Tissue Cell. 2010;42:301–6.
Castro-Rendón WA, Castro-Alvarez JF, Guzmán-Martinez C, Bueno-Sanchez JC. Blastocyst-endometrium interaction: intertwining a cytokine network. Braz J Med Biol Res. 2006;39:1373–85.
Burd CE, Jeck WR, Liu Y, Sanoff HK, Wang Z, Sharpless NE. Expression of linear and novel circular forms of an INK4/ARF-associated non-coding RNA correlates with atherosclerosis risk. PLoS Genet. 2010;6:e1001233.
Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu J, et al. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19:141–57.
Mesquita FS, Ramos RS, Pugliesi G, Andrade SC, Van Hoeck V, Langbeen A, et al. The receptive endometrial transcriptomic signature indicates an earlier shift from proliferation to metabolism at early diestrus in the cow. Biol Reprod. 2015;93(2):52.
Ling YH, Ren CH, Guo XF, Xu LN, Huang YF, Luo JC, et al. Identification and characterization of microRNAs in the ovaries of multiple and uniparous goats (Capra hircus) during follicular phase. BMC Genomics. 2014;15:339.
Szymanowski K. Apoptosis pattern in human endometrium in women with pelvic endometriosis. European J Obstet Gynecol Reprod Biol. 2007;132:107–10.
Kapucuoglu N, Aktepe F, Kaya H, Bircan S, Karahan N, Çiriş M. Immunohistochemical expression of PTEN in normal, hyperplastic and malignant endometrium and its correlation with hormone receptors, bcl-2, bax, and apoptotic index. Pathol Res Pract. 2007;203:153–62.
Podsypanina K, Ellenson LH, Nemes A, Gu J, Tamura M, Yamada KM, et al. Mutation of Pten/Mmac1 in mice causes neoplasia in multiple organ systems. Proc Natl Acad Sci U S A. 1999;96:1563–8.
Yun KW, Lee JY, Yun SW, Lim IS, Choi ES. Elevated serum level of MicroRNA (miRNA)-200c and miRNA-371-5p in children with Kawasaki disease. Pediatr Cardiol. 2014;35:745–52.
Horikoshi Y, Matsumoto H, Takatsu Y, Ohtaki T, Kitada C, Usuki S, et al. Dramatic elevation of plasma metastin concentrations in human pregnancy: metastin as a novel placenta-derived hormone in humans. J Clin Endocrinol Metab. 2003;88:914–9.
Murray MJ, Lessey BA. Embryo implantation and tumor metastasis: common pathways of invasion and angiogenesis. Semin Reprod Endocrinol. 1999;17:275–90.
Albini A, Benelli R. The chemoinvasion assay: a method to assess tumor and endothelial cell invasion and its modulation. Nature Protocol. 2007;48:504–11.
Ebert MS, Neilson JR, Sharp PA. MicroRNA sponges: competitive inhibitors of small RNAs in mammalian cells. Nat Methods. 2007;4:721–6.
Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39:1033–7.
Poliseno L, Salmena L, Zhang J, Carver B, Haveman WJ, Pandolfi PP. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature. 2010;465:1033–8.
Liu X, Zhang L, Liu Y, Cui J, Che S, An X, et al. Circ-8073 regulates CEP55 by sponging miR-449a to promote caprine endometrial epithelial cells proliferation via the PI3K/AKT/mTOR pathway. Biochim Biophys Acta (BBA) - Mol Cell Res. 2018;1865:1130–47.
Bachmayrheyda A, Reiner AT, Auer K, Sukhbaatar N, Aust S, Bachleitnerhofmann T, et al. Correlation of circular RNA abundance with proliferation--exemplified with colorectal and ovarian cancer, idiopathic lung fibrosis, and normal human tissues. Sci Rep. 2015;5:8057.
The authors would like to acknowledge XinLongMen milk goat Breeding Farm for providing the goats.
This study was supported by PhD research startup foundation of Northwest A&F University (00400/Z109021811), National Key Research and Development Program of China (2016YFD0500508), and Shaanxi Science and Technology Innovation Project Plan (2015KTCQ03–08). The authors declare that they have no competing interests. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and materials
The data generated and transcripts obtained were deposited at NCBI/GenBank as the GEO with accession number GSE85384. The data and material has been also provided in form of Additional files.
Ethics approval and consent to participate
All experimental goats were provided by XinLongMen milk goat Breeding Farm, Lantian County, Xi’an, Shaanxi Province, China. The goats were maintained according to the No. 5 Proclamation of the Ministry of Agriculture, P. R. China, and we confirm that all animal protocols and methods were approved by the Review Committee for the Use of Animal Subjects of Northwest A&F University.
Consent for publication
The authors declare that they have no competing interests. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Table S1. The groups of hormone treatment (nmol/L) (DOCX 15 kb)
Figure S1. The overview of the results of RNA-Seq. (A) Principal component analysis (PCA) of the results of RNA-Seq. The six independent endometrium samples collected from 3 goats at gestational day 5 (PE, GE258-C, GE796-C, GE990-C) and 3 goats at gestational day 15 (RE, GE484-T, GE706-T, GE5928-T), respectively. (B) Boxplot of the results of RNA-Seq for four endometrium samples. Boxplot of the log10 FPKM (Fragments Per Kilobase of exon per million fragments mapped) expression values in four endometrium samples. The Fig. reflects the distribution of FPKM values computed for the circRNAs in each samples from RNA-Seq data, and shows that the median of the expression values across the samples being compared for differential expression are comparable. (TIF 2658 kb)
Table S2. The result of nonlinear alignment using TopHat-Fusion. (XLS) (XLS 24 kb)
Table S3. CircRNA expression in goat endometrium. (XLS) (XLS 9770 kb)
Table S4. Differential expressed circRNA between receptive and pre-receptive endometrium. (XLS) (XLS 1658 kb)
Table S5. CircRNA specifically expressed in the RE or PE. (XLS) (XLS 3082 kb)
Table S6. Enrichment analysis of hgDEC. (XLS) (XLS 420 kb)
Table S7. Sgnificant GO terms of hgDEC. (XLS) (XLS 177 kb)
Table S8. KEGG significant analysis of hgDEC. (XLS) (XLS 58 kb)
Table S9. All predicted miRNAs bind to circRNAs in goat. (XLS) (XLS 9205 kb)
Table S10. The primers used in this study. (XLS) (XLS 26 kb)