circRNA landscape of non-pregnant endometrium during the estrus cycle in dairy goats

Liu Xiao-rui,ZhANg Lei,Cui Jiu-zengYANg Li-chunhAN Jin-chengChe Si-chengCAo Bin-yun,Li guang,SoNg Yu-xuan

1 College of Animal Science and Technology,Northwest A&F University,Yangling 712100,P.R.China

2 Shaanxi Provincial Key Laboratory of Animal Genetics,Breeding and Reproduction,College of Animal Science and Technology,Northwest A&F University,Yangling 712100,P.R.China

Abstract Endometrial development is a complicated process involving numerous regulatory factors.Circular RNAs (circRNAs)have been known as a member of the naturally occurring non-coding RNA family,and are reportedly crucial for a variety of physiological processes.This study investigated the circRNA landscape of non-pregnant endometrium of dairy goats during estrus.Non-pregnant endometrial samples of goats at estrus day 5 (Ed5) and estrus day 15 (Ed15) were used to methodically analyze the circRNA landscape using strand-specific Ribo-Zero RNA-Seq.A total of 2 331 differentially expressed (P<0.05) circRNAs (DEciRs) between Ed5 and Ed15 were discovered in the goat endometrium.It was found that Nipped-B-like (NIPBL) and calcium responsive transcription factor (CARF) may participate in the development of the endometrium by decreasing (P<0.05) the levels of their circRNA-transcript forms.Furthermore,Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of DEciR host genes (hgDEciRs) revealed that tight junctions and GTPases may be involved in endometrial development during the estrus cycle.A total of 2 331 DEciRs were discovered in the endometrium at Ed5 and Ed15.Based on GO and KEGG enrichment analyses,it could be inferred that tight junctions and GTPases are likely to play an important role in the development of goat endometrium during the estrus cycle.This circRNA study greatly enhances our knowledge of global trends in the development of non-pregnant endometrium during the estrus cycle in goats;these results help us to better understand the molecular regulation of endometrial development in dairy goats.

Keywords:circRNAs,Ribo-Zero RNA-Seq,non-pregnant endometrium,dairy goats

1.introduction

During breeding seasons,the mammalian uterus grows and is substantially remodeled for implantation of the embryo and to accommodate the development of the conceptus(Igwebuike 2009).The endometrium also undergoes defined morphological changes to attain a receptive state(Gonget al.2014) owing to the integration and synchronized interaction of numerous genes and factors (Zhang Set al.2013).Transcriptome analysis at different physiological periods is pivotal for fully understanding the relevant gene regulatory networks,including those involved in endometrial development in mammals (Zhanget al.2015,2017a;Wanget al.2016).Moreover,such studies on the human (Huet al.2014;Lira-Albarránet al.2019),pig (Ma Xet al.2019),cattle(Mazzoniet al.2020),and goat endometrium (Zhanget al.2015,2017a) have been reported.

Increasing evidence indicates that many genes expressed in the endometrium are regulated through epigenetic pathways,including those controlled by miRNA(Munroet al.2010;Zhanget al.2017c,e).For example,in dairy goats,miR-26a induces apoptosis in ESCs(endometrial stromal cells) and enhances the proliferative capacity of EECs (endometrial epithelium cells) through the PTEN/PI3K/AKT axis (Zhanget al.2017b).By targeting the 3´-untranslated region (3´-UTR) of Pleiotrophin (PTN),miR-182 induces apoptosis in EECs (Zhanget al.2017a)and regulates HOXA10 expression through the 3´-UTR(Zhanget al.2017e).

RNA-Seq shows potential for simultaneous transcript discovery and abundance estimation (Mortazaviet al.2008;Trapnellet al.2010);this technology has been widely used to identify novel and conservative circular RNAs (circRNA)in multifarious species,and several potential circRNAs have been identified using this RNA-Seq technology (Danget al.2016;Wanget al.2016;Zhaoet al.2017;Zhanget al.2018).circRNA has the potential effect of competing with miRNA sponges or endogenous RNA (Hansenet al.2013).circRNAs from back-spliced exons have been confirmed as belonging to the non-coding RNA family naturally found in eukaryotes (Jeck and Sharpless 2014;Lasda and Parker 2014).Recently,circRNAs have become a focus of study globally,and several circRNAs have been reported in different cell lines and species (Guoet al.2014;Rybak-Wolfet al.2015;Youet al.2015).

Goats are spontaneously ovulating,polyestrous animals,with gestational days 5 and 15 being the most important time points for embryo implantation (Igwebuike 2009;Songet al.2015;Zhanget al.2015).Moreover,in our previous study,well-formed pinopodes were found on the surface of the endometrium on day 15 with high E2 and P4 concentrations in the serum (Zhanget al.2017d).Expression levels of some marker genes of the receptive endometrium (e.g.,osteopontin (OPN),vascular endothelial growth factor (VEGF),leukemia inhibitory factor (LIF),and prolactin receptor (PRLR)) increased on day 15 compared to those on day 5,indicating that the endometrium of dairy goats attains a receptive state at gestational day 15 (Zhanget al.2017d).circRNA-3175 is known to be differentially expressed in receptive endometrium compared to that in pre-receptive endometrium and displays a sponging effect for miR-182 in EECs during embryo implantation (Zhanget al.2018).In addition,circRNA-8073 promotes EECs proliferationviathe PI3K/AKT/mTOR pathway by sponging miR-449a,suggesting that circRNA-8073 plays a regulatory role in the endometrial receptive development of dairy goats (Liuet al.2018).Furthermore,miR-34a/c induces EECs apoptosis by regulating circRNA8073/CEP55viathe RAS/RAF/MEK/ERK and PI3K/AKT/mTOR pathways(Liuet al.2020).circRNA-8073 also functions as a ceRNA that sequesters miR-181a,thereby protecting neurotensin transcripts from miR-181a-mediated suppression in EECs(Zhanget al.2019).

However,it is not known whether circRNAs participate in the development of the non-pregnant endometrium in goat estrus.Considering the universal expression of circRNAs and their key role in regulating gene expression,we hypothesized that there would be differentially expressed circRNAs (DEciRs) in the non-pregnant endometrium of dairy goats at estrus days 5 (Ed5) and 15 (Ed15),and these circRNAs (DEciRs) would play important roles in the development of goat endometriumviatheir host genes and/or miRNA during the estrus cycle.To verify these hypotheses,we first investigated the circRNA expression profile in the non-pregnant endometrium at Ed5 and Ed15 using strand-specific Ribo-Zero RNA-Seq,and then analyzed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) of the DEciR host genes(hgDEciRs) and the interaction network of circRNA-miRNA using bioinformatics analysis.

2.Materials and methods

2.1.experimental design and sample collection

A completely randomized design was adopted in this study.A total of 10 healthy and multiparous dairy goats(24-mon-old) (Xinong Saanen) were induced to estrous synchronization (Zhanget al.2017d;Liuet al.2019).On day 1,the goats were injected with prostaglandin F2α(0.2 mg) (Ningbo Pharmaceutical Co.,Ltd.,China) and then treated with intravaginal sponges which contained MPA (medroxyprogesterone acetate;60 mg).The goats were then injected twice with FSH (20 IU,follicle-stimulating hormone) (Ningbo Pharmaceutical Co.) on day 10 at 8:00 h and 19:00 h.On day 11,PMSG (200 IU,pregnant mare’s serum gonadotropin) (Ningbo Pharmaceutical Co.,Ltd.) and prostaglandin F2α (0.1 mg) were intramuscularly injected into the goats.The goats were observed three times a day for signs of estrus (using vasectomized teaser goats).None of the goats mated on the first day of estrus,which was defined as estrus day 0 (Ed0) in this study.Five goats were euthanized by intravenous injection of barbiturate(30 mg kg–1) on day 5 (Ed5) and five goats were euthanized on day 15 (Ed15) of estrus;endometrium samples were obtained from the anterior wall of the uterine cavity (Songet al.2015).All samples were washed briefly with ice-cold phosphate buffered saline (PBS) and immediately frozen in liquid nitrogen.Thus,there were two groups (Ed5 and Ed15) with five replicates in each group and one dairy goat in each replicate in this study.

For each experimental goat,a total of 2.8 kg basal diet(2.45 kg DM,Appendix A) was supplied each day.The basal diet was formulated to meet the nutrient requirements of dairy goats specified in the NRC (2007),and all goats were fed in the semi-open shed of Xinlongmen Milk Goat Breeding Farm,Xi’an City,Shaanxi Province,China.

2.2.Library construction and RNA-Seq

Total RNA was extracted from goat endometrium using Trizol reagent (Invitrogen,CA,USA) according to the manufacturer’s instructions (Zhanget al.2017a).RNA amount and purity of each sample were quantified using NanoDrop ND-1000 (NanoDrop,Wilmington,DE,USA).RNA integrity was assessed using Agilent 2100 with RIN number>7.0 (Songet al.2019).Ed5 and Ed15 libraries were constructed by pooling five homogenized total RNAs from the endometrium on Ed5 and Ed15,respectively.Approximately 5 μg of total RNA was used to deplete ribosomal RNA according to the instructions in the Ribo-Zero™ rRNA Removal Kit (Illumina,San Diego,USA).The remaining RNAs were treated with RNase R (Epicentre Inc,Madison,WI,USA) to remove linear RNAs and enrich circRNAs.After the removal of ribosomal and linear RNAs,the enriched circRNAs were fragmented into small pieces using divalent cations at high temperatures.Cleaved RNA fragments were reverse-transcribed to create cDNA,which were then used to synthesize U-labeled second-stranded DNAs withEscherichia coliDNA polymerase I,RNase H,and dUTP.An A-base was then added to the blunt ends of each strand to prepare them for ligation to indexed adapters.Each adapter contained a T-base overhang for ligating it to an A-tailed DNA fragment.Single-or dual-index adapters were ligated to the fragments and size selection was performed using AMPureXP beads.After the heat-labile UDG enzyme treatment of the U-labeled second-stranded DNAs,the ligated products were amplified using PCR under the following conditions:initial denaturation at 95°C for 3 min,eight 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 final extension at 72°C for 5 min (Songet al.2019).The average insert size for the final cDNA library was 250 bp (±50 bp).Lastly,we performed a paired-end sequencing using an Illumina Hiseq 2500 (LC Bio.,China)following the vendor’s recommended protocol.

2.3.RNA-Seq reads mapping

Before the goat reference genome map was drawn,raw reads (raw data) were processed to remove low-quality reads and all downstream analyses were based on clean data with high quality.The clean reads were mapped to the goat reference genome (http:/www.gencodegenes.org/,GRCh38) using a publicly available program (TopHat 2)(Kimet al.2013),which aligned the reads to genomes using ultra-high-throughput short-read aligner Bowtie2(Langmead and Salzberg 2012).The mapping data were then analyzed to identify splice junctions.Unmapped reads were mapped with TopHat-Fusion on the same chromosome in a noncollinear ordering and were extracted as candidate back-spliced junction reads.Lastly,according to the structure and splice site characters of circRNA,the circRNAs were identified using the back-splice junction reads.

2.4.Goat circRNA quantification

In this study,we used Cufflinks Software to assemble circRNAs,assess their abundance,and detect their differential expression.The unit of measure was Fragment Per Kilobase of exon per Million fragments (FPKM) (Trapnellet al.2010).Only comparisons with“q-value<0.01”and status flagged as“OK”in Cuffdiff output were considered to show differential expression.

2.5.go and Kegg analyses of circRNA-hosting genes

The hypergeometric test was used to map all differentially expressed genes to the appropriate terms in a GO Database(Consortium 2004),which is an international standard gene functional classification system (Ashburneret al.2000).A corrected“P-value<0.05”was used as the threshold to find significantly enriched GO terms in the input list of hgDEciRs compared to their genomic background.

KEGG helps researchers to better understand the biological functions of genes based on large-scale molecular datasets (Guoet al.2014).In this study,a correctedP-value<0.05 was used as the threshold (He and Liu 2013)to find significantly enriched KEGG terms compared to their genomic background in the DEG input list (Kanehisaet al.2008).

2.6.The interaction network analysis of circRNAmiRNA

circRNAs are believed to contribute substantially to studies on competitive endogenous RNA networks,and such studies have reported that they acted as miRNA sponges to regulate gene expression (Hansenet al.2013;Jens 2013).In this study,we used Targetscan 7.0 and miRanda,two publicly available programs,to predict circRNA-miRNA interactions.The 436 goat miRNAs used in the analysis were obtained from the miRBase Public Databases (http://www.mirbase.org/,release 21).

2.7.RNA preparation and RT-qPCR

Trizol (TaKaRa,Dalian,China) was used to extract the total RNA from goat endometrium samples.The extracts were then incubated with or without (mock) RNase R (3 U μg–1) (Epicenter Biotechnologies,Chicago,USA) at 37°C for 20 min.cDNA was synthesized using the Prime Script RT Reagent Kit with gDNA Eraser (TaKaRa,Dalian,China).RTqPCR was performed in a 20 μL reaction using SYBR Green PCR Master Mix (TaKaRa,Dalian,China).Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal reference.The relative expression of circRNAs was calculated using N=2–ΔΔCt.

2.8.Data statistical analyses

SPSS (SPSS statistical software version 17.0;SPSS Inc.,Chicago,USA) was used in this study to analyze the relative expression of circRNAs in the non-pregnant endometrium at Ed5 and Ed15 in dairy goats.Thet-test was used to compare the differences between these two groups,and the differences were considered significant atP<0.05 and very significant atP<0.01.Each replicate served as the experimental unit for the statistical analyses.

3.Results

3.1.Transcript expression in goat endometrium

From each endometrium library,120 million paired-end reads ((250±50) bp,18 Gb) were generated,representing approximately seven-fold coverage of the goat genome(2.66 Gb).There were 13.91 and 14.51 G of valid data and the percentages were 77.29 and 80.58% in the Ed5 and Ed15 libraries,respectively (Table 1).Statistical analysis of the distribution of mapped reads in each library is shown in Appendix B and the TopHat Fusion mapped reads were grouped (into exon group,intron group,and intergenic group) according to gene location information (Fig.1-A).

The numbers of circRNAs and corresponding circRNA host genes in each library are displayed in Table 1 and a total of 10 216 circRNAs were identified (Appendix C).Subsequently,these circRNA candidates were annotated according to the corresponding circRNA host gene location;and we then mapped 3 967 corresponding circRNA host genes in the goat genome (Appendix C).Further investigation revealed 7 253 circRNAs from Ed5 and 1 226 circRNAs from Ed15 (Fig.1-B),suggesting a high stage-specificity of circRNAs in the goat endometrium.Simultaneously,3 692 circRNA host genes were identified at Ed5 and 1 840 genes at Ed15.

3.2.Differential expression of circRNAs (DeciRs) in goat endometrium between ed5 and ed15

Fold change (log2(Ed15/Ed5)) and the corresponding significance threshold of theP-value were estimated according to normalized circRNAs expression levels.Based on the expression levels,significantly different DEciR values between Ed5 and Ed15 were identified usingP-value<0.05 as the threshold.In total,2 331 circRNAs were found to be significantly different (P<0.05) between Ed5 and Ed15 samples,and the full list of DEciRs is provided in Appendix D.Overall,1 077 differentially expressed circRNAs were downregulated;however,1 254 circRNAs were upregulated(P<0.05) on Ed15 compared to Ed5 in the goats (Fig.2-A).Further analysis showed that 262 DEciRs were co-expressed during the two physiological stages of the goat endometrium,whereas 1 022 and 1 047 DEciRs were specifically expressed(P<0.05) on Ed5 and Ed15,respectively (Fig.2-B).To verify whether the back-spliced events were indicative of true circular molecules,we examined the physical properties of 10 DEciRs using RNase R treatment and RT-qPCR,followed by standard Sanger sequencing.Ten DEciRs were found to be resistant to RNase R,and the sequences of the 10 circRNAs were consistent with the sequencing data of rRNA deletions (Fig.3).

Table 1 Number of reads and circular RNAs (circRNAs) in goat endometrium

The 10 DEciRs with the highest FPKM values at Ed15 are shown in Table 2;the FPKM value of circRNA1460(FPKM value=298.80) was the highest among them and it was specifically expressed at Ed15,followed by circRNA279(FPKM value=208.08).The 10 DEciRs with the highest FPKM values at Ed5 are shown in Table 3.The highest among them was circRNA2045 (FPKM value=2 369.17)at Ed5,followed by circRNA8694 (FPKM value=1 650.11),which was specifically expressed at Ed5.

3.3.go and Kegg analyses of circRNA host genes(hgDeciRs)

A total of 1 952 DEciRs were annotated to 1 374 hgDEciRs,revealing that a single gene might have generated more than one circRNA.For example,CRIM1 produced 10 circRNAs and TGFBR3 produced 15 circRNAs in the goat endometrium (Appendix E).Considering the fact that a given circRNA was spliced from a given linear transcript,the function of circRNAs might be associated with the function of the host gene (Youet al.2015).Therefore,circRNA host genes were queried against the GO Database.

We focused on the GO terms that met the designated criteria of aP-value<0.05 and categorized the hgDEciRs into 27 prominent GO terms (Table 4).Thirteen terms were those enriched in biological processes (Fig.4-A),including protein ubiquitination which was referred to as ubiquitindependent protein catabolic processes (GO:0042787),tight junction assembly (GO:0070830),and multicellular organism development (GO:0007275).In the GO category of cellular compartment,nine terms were markedly enriched (Fig.4-B).The most notably enriched terms werecis-Golgi network(GO:0005801) with seven hgDEciRs annotated,followed by the mitochondrial outer membrane (GO:0005741),lysosome (GO:0005764),and coated pit (GO:0005905).Five terms were prominently enriched in molecular functions(Fig.4-C),of which the most notably enriched GO term was small GTPase regulator activity (GO:0005083) with six hgDEciRs annotated,followed by GTPase activator activity(GO:0005096),protein binding (GO:0005515),growth factor binding (GO:0019838),and Rac GTPase binding(GO:0048365).

Fig.1 The detected reads and circular RNAs (circRNAs) in goat endometrium at Ed5 and Ed15.A,genomic features of the detected reads in goat endometrium using Tophat Fusion.B,Venn diagrams of detected transcripts in goat endometrium.Ed5 and Ed15,estrus day 5 and 15,respectively.

Fig.2 Differentially expressed circRNAs (DEciRs) between Ed5 and Ed15.A,scatter plot of the correlation of DEciRs abundance between Ed5 and Ed15.B,Venn diagrams of EDciRs between Ed5 and Ed15 in goat endometrium.Ed5 and Ed15,estrus day 5 and 15,respectively.

Fig.3 The expression levels of 10 differentially expressed circRNAs (DEciRs) as determined by RT-qPCR.GAPDH was used as the internal control gene for normalization in our experiments.Values are means±SD in estrus day 15 (Ed15,n=5) that were relative to estrus day 5 (Ed5,n=5),whose values were all set to 1 in this study because of the calculation step of 2–ΔΔCt.** indicates P<0.01 and * indicates P<0.05.

Several genes expressed their respective biological functions.KEGG analysis of the hgDEciRs was conducted to further understand the biological functions of the DEciRs.In this study,the KEGG pathway annotation showed that 700 hgDEciRs annotated to 220 biological processes(Appendix F).Moreover,four KEGG pathways met the criterion ofP-value<0.05 (Table 5;Fig.5),suggesting that these pathways might play important roles in the development of endometrial receptivity.The pathways that showed the highest levels of significance were those of fructose and mannose metabolism (ko00051,P=0.001)with eight hg-DEciRs enriched,followed by endocytosis(ko04144,P=0.002),autophagy regulation (ko04140,P=0.01),and p53 signaling pathway (ko04115,P=0.041).Our findings indicated that these pathways are active in endometrial development and that a variety of cellular activities are involved in the developing endometrium.

3.4.Target miRNA predictions for identified circRNAs in goat endometrium

Some exonic circRNAs function as miRNA sponges,whichregulate linear protein-coding RNA product expression through“mRNA trap”.In our study,we used Targetscan 7.0 and miRanda to predict circRNA-miRNA interactions.In summary,10 142 circRNAs were predicted as possible targets for 436 goat miRNAs,indicating that a single miRNA could be regulated by more than 23 circRNAs (Appendix G).

Table 2 The 10 differentially expressed circRNAs (DEciRs) with the highest FPKM values at Ed151)

Table 3 The 10 differentially expressed circRNAs (DEciRs) with the highest FPKM values at Ed51)

Table 4 The 27 significant GO Term annotations for the DEciR host genes (hgDEciRs)

4.Discussion

The hypotheses that there would be differentially expressed circRNAs (DEciRs) in the non-pregnant endometrium of dairy goats at estrus days 5 (Ed5) and 15 (Ed15),and that these DEciRs would play important roles in the development of goat endometriumviatheir host genes and/or miRNA have been supported by the results of this study.The present study has demonstrated that there were 2 331 DEciRs between Ed5 and Ed15 in the goat non-pregnant endometrium,while GO and KEGG analyses of the DEciR host genes (hgDEciRs) revealed that tight junctions and GTPases might be involved in endometrial development during the estrus cycle,and some circRNAs were predicted as possible targets for goat miRNAs using bioinformatics analysis.Therefore,these DEciRs might play an important role in the development of the goat endometrium during the estrus cycle.These results enhance the existing knowledge on the expression dynamics of circRNAs in goat endometrium during the estrus cycle and help us better understand the molecular regulation of endometrial development in dairy goats.

Multiple circRNAs are generated by alternative circularization,including exons (Salzmanet al.2013;Zhanget al.2014) and introns (ciRNAs) (Zhang Let al.2013),and a majority (more than 96% in this study) of circRNAs consist of protein-coding exons,in accordance with several reports (Zhenget al.2016).Ours is the first report showing circRNAs in goat non-pregnant endometrium during the estrus cycle.

In recent years,an increasing number of circRNAs have been identified in goats.About 13 320 circRNAs were identified in the skin of Liaoning cashmere goats (Zhenget al.2020),37 818 in the mammary tissue of Laoshan dairy goats (Ma Det al.2019),9 090 in the skeletal muscles of Anhui white goat (Linget al.2020),and 13 950 in preovulatory ovarian follicles of Macheng black goats and Boer goats (Taoet al.2018).There are 92 061 human circRNAs in the circBase (Petaret al.2014) and 65 731 circRNAs have been identified in human brain tissues (Rybak-Wolfet al.2015).In addition,more than 27 296 circRNA candidates have been detected in six normal and seven cancer tissues(Zhenget al.2016).We identified fewer circRNAs in this study than are often found in humans.This is mainly because more studies have been conducted on human circRNAs than other species.The human genome has more detailed information available than the genomes of other species and is better annotated.Moreover,human circRNAs have been identified in various tissues such as lung,brain,heart,colon,stomach,liver,and some tumor tissues (Zhenget al.2016).The identification of circRNAs in the present study was restricted to goat endometrium at two physiological stages.Furthermore,circRNA expression is regulated by the development and changes which occur in different tissues (Salzmanet al.2013;Fanet al.2015;Youet al.2015).In the present study,10 216 circRNAs were identified in the non-pregnant endometrium of dairy goats.

Fig.4 The GO enrichment analyses of DEciR host genes(hgDEciRs) with P-values less than 0.05.A,biological processes.B,cellular components.C,molecular functions.

Our study first aimed to identify circRNAs in the endometrium of dairy goats on Ed5 and Ed15 of their estrus cycle using RNA-Seq.A total of 2 331 circRNAs differed significantly in terms of FPKM levels between Ed5 and Ed15;1 254 circRNAs were upregulated and 1 077 circRNAs were downregulated on Ed15 compared to their levels at Ed5.Further analysis revealed that 262 of the 2 331 DEciRs were expressed in both physiological stages of the goat endometrium,1 023 DEciRs were specifically expressed at Ed5,and 1 047 DEciRs were specifically expressed at Ed15,suggesting that the specifically expressed circRNAs might play an important role in the development of the goat endometrium.

To verify whether the back-spliced events were indicative of true circular molecules,we examined the physical properties of the products using RT-qPCR.Firstly,circRNAs are known to be strongly resistant to exonuclease RNase R compared to linear RNAs.The circRNA candidates here were resistant to RNase R treatment and there was a reduction in the abundance of linear RNAs,thereby verifying the circularized features of the molecules.Secondly,since circRNAs have a distinctive head-to-tail junction site,which is different from typical linear splicing transcripts,we used qPCR with outward-facing primers to verify the circRNA candidates through targeted amplification of head-totail junction,followed by standard Sanger sequencing.Our results showed that the levels of 10 DEciRs were in accordance with the sequencing data.

The DEciR with the highest FPKM value at Ed15 was circRNA1460 (FPKM value 298.80) which was specifically expressed at Ed15,and the DEciR with the highest FPKM value at Ed5 was circRNA2045 (FPKM value=2 369.17),which was downregulated 0.08-fold at Ed15,followed by circRNA8694 (FPKM value=1 650.11),which was specifically expressed at Ed5.The function of any given circRNA is likely related to the function of the host genebecause a circRNA is spliced to generate a given linear transcript (Youet al.2015).

Table 5 The four KEGG functional annotations for the DEciR host genes (hgDEciRs)

Fig.5 Scatter plot of the KEGG pathway enrichment analysis of DEciR host genes (hgDEciRs) with P-values less than 0.05.The x-axis indicates the number of unique sequences assigned to a specific pathway,and the y-axis represents the KEGG pathway.Fisher’s exact test was used to generate pathway enrichment statistics.P<0.05 indicates the significant pathways.

Nipped-B-like (NIPBL),the host gene of circRNA1460,is the most frequently mutated gene (60% of cases,OMIM 122470) found in the human developmental disorder Cornelia de Lange syndrome (CdLS,1 of 10 000–30 000 live births) (Krantzet al.2004;Huismanet al.2013;Moradet al.2014).This syndrome is characterized by upper limb malformations,craniofacial anomalies including hair growth and mental retardation,and other systemic abnormalities(Zuinet al.2017).Hence,NIPBL has evolved a complex response to damaged DNA,thereby revealing its dynamic role in maintaining genomic stability (Botet al.2017).

Calcium is crucial for signal transduction.It initiates the intracellular signaling cascade that transmits signals to the nucleus and activates multiple transcription factors in the nucleus;and these transcription factors combine the regulatory elements of activity-regulated genes (Lyonset al.2012).The role of calcium in endometrium growth and development in humans and animals has been verified(Pietras and Szego 1975;Burnset al.1998;Zhanget al.2011).The host gene of circRNA8694,calcium responsive transcription factor (CARF),may also be used as a basic inhibitor because within a group of genes that CARF binds to their promoters,equal numbers of genes are increased and decreased compared to their control group in CARF knockout neurons (Pfenninget al.2010).Although the mechanisms by which CARF activates and inhibits genes are not entirely clear,the short variants encoded by a CARF gene contain a DNA binding domain but lack the transcription activation domain (Taoet al.2002).As a whole,our results indicated that CARF primarily acts on activity-regulated gene promoters and as a basal activator or repressor,but not as a calcium-inducible activator (Lyonset al.2016).The results of the present study suggested that NIPBL and CARF may participate in the development of endometrium by decreasing the circRNA form;however,this hypothesis requires further investigation.

In this study,1 952 DEciRs were annotated to 1 374 hgDEciRs,suggesting that one gene might produce more than one circRNA,and this was consistent with a previous report (Zhanget al.2014).There were 27 significant GO terms meeting the criterion ofP-value<0.05.Among them,tight junction assembly (GO:0070830) and tight junction (GO:0005923) are interesting because the tight junction structure of the endometrial epithelium is regulated during the menstrual cycle,suggesting that the integrity of the epithelial barrier may be important in preparing the endometrium for implantation (Rogers and Murphy 1992).Tight junctions are the only type of intercellular junctions that have been identified between microvascular endothelial cells in the endometrium (Karakotchian and Fraser 2007).Tight junctions prevent uncontrolled leak flux and sustain mucosal monolayer integrity,which are beneficial for endometrial health (Wiraet al.2010;Kiatprasertet al.2015).Besides,studies have shown that GTPases have significant effects on endometrium development.For example,Rho GTPases in human endometrial stromal cells play opposing roles in regulating embryo invasion and focal adhesion transformation (Grewalet al.2010),and Rho GTPases in mouse uterus endometrium are involved in the process of implantation (Ma and Zhang 2006).In this study,GTPase activator activity (GO:0005096),small GTPase regulator activity (GO:0005083),and Rac GTPase binding activity(GO:0048365) were significant GO terms,suggesting that GTPases have important effects on goat endometrium development;however,this hypothesis requires further study for corroboration.

Enriched pathways might suggest that organisms respond to the conditions in a particular endometrial stage during the estrus cycle.In the present study,the hgDEciRs were markedly enriched in four KEGG pathways,revealing that a few genes were referred to several cell proliferationrelated pathways and regulated endometrium developmentviatheir related circRNAs.

Previous studies have demonstrated that several circRNAs bind competitively with miRNAs to release and increase the levels of miRNA target genes,thereby exerting their biological effects (Ebertet al.2007;Franco-Zorrillaet al.2007;Polisenoet al.2010).In this study,circRNAs(10 142) were predicted as possible targets for the goat miRNAs (436),revealing that one miRNA target was likely to be regulated by two or more circRNAs (Zhenget al.2016).Further analysis showed that a single circRNA was used as a sponge for numerous miRNAs.In addition,the sponge effect of exon circRNAs on miRNAs may be rare or does not require a mass of miRNA binding sites (Jeck and Sharpless 2014).Our results suggested that circRNAs may play significant roles in endometrium development by impairing the effect of miRNAs.However,considering that animal miRNAs were not completely complementary to their targets,it was more challenging to identify targets and determine prediction accuracy.Therefore,the actual regulatory relationships between these circRNAs and miRNAs need further study for verification.

5.Conclusion

In this study,the application of RNA-Seq markedly enhanced the current understanding of transcriptional events in goat endometrium during the estrus cycle.A total of 2 331 DEciRs were discovered in the endometrium at Ed5 and Ed15.Based on GO and KEGG enrichment analyses,tight junctions and GTPases likely play an important role in the development of goat endometrium during the estrus cycle.Thus,the findings of this study helped us to better understand the molecular regulation of endometrial development in goats.

Acknowledgements

This work was supported by the China Postdoctoral Science Foundation (2019M653776,2020M673516),the National Key Research and Development Program of China (2016YFD0500508),the Natural Science Foundation of Shaanxi Province,China (2020JQ-265),and the PhD research startup foundation of Northwest A&F University,China (00400/Z109021811).

Declaration of competing interest

The authors declare that they have no conflict of interest.

ethical approval

All animals in our study were raised in accordance with the No.5 announcement of the former Ministry of Agriculture of China.All animal experimental procedures were approved by the Review Committee for the Use of Animal Subjects of Northwest A&F University,China.

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm