Exploring the Roles of Fecundity-related Long Non-coding RNAs and mRNAs in the Adrenal Gland of Small-tailed Han Sheep

DOI: https://doi.org/10.21203/rs.2.15698/v2

Abstract

Background: Long non-coding RNAs (lncRNAs) can play important roles in uterine and ovarian functions. However, little researches have been done on the role of lncRNAs in the adrenal gland of sheep. Herein, RNA sequencing was used to compare and analyze gene expressions in adrenal tissue between follicular phase and luteal phase in FecBBB (MM) and FecB++ (WW) sheep, respectively, and differentially expressed lncRNAs and genes associated with reproduction were identified.

Results: In MM sheep, 38 lncRNAs and 545 mRNAs were differentially expressed in the adrenal gland between the luteal and follicular phases; In WW sheep, 513 differentially expressed lncRNAs and 2481 mRNAs were identified. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses indicated that differentially expressed lncRNAs and their target genes are mainly involved in the circadian rhythm, the mitogen activated protein kinase, thyroid, ovarian steroidogenesis and transforming growth factor beta signaling pathways. Differentially expressed lncRNAs can regulate reproduction by modulating genes involved in these signaling pathways and biological processes. Specifically, XLOC_254761, XLOC_357966, 105614839 and XLOC_212877 targeting CREB1, PER3, SMAD1 and TGFBR2, respectively, appear to play key regulatory roles.

Conclusion: These results broaden our understanding of lncRNAs in adrenal gland of sheep and provide new insights into the molecular mechanisms underlying sheep reproduction.

Background

Small-tailed Han sheep are famous for their high fertility, precocious puberty, good fur quality and tall body shape in China[1]. The behavior of estrus and mate in Small-tailed Han sheep appear year-round. The lambing rate of primiparous ewes is about 200%, and in produced ewes is higher than 250%[1]. FecBB mutations have huge economic benefits in production, which can significantly increase the number of ovulation and lambs in sheep. According to the previous research report of this team, all three genotypes of FecB (FecBBB, FecBB+ and FecB++) are distributed in Small-tailed Han sheep, and there is a significant correlation between the three genotypes of FecB and the lamb size of ewes[2]. Therefore, Small-tailed Han sheep can be used as an ideal animal model to study the molecular mechanism of FecB gene regulation of reproductive traits.

In recent years, there have arisen many methods for identifying differentially expressed candidate long noncoding RNAs (lncRNAs) and genes using transcriptome sequencing. lncRNAs are RNAs of >200 nucleotides in length. Studies have shown that lncRNAs play important roles in many life activities such as dose-compensation effects, epigenetic regulation, cell cycle regulation, cell differentiation and regulation, and have become a hot research topic in genetics[3-6]. For example, Yang et al. (2018) identified differentially expressed lncRNAs and mRNAs in the testes of prepubertal and mature rams that were enriched in spermatogenetic and male gonadal developmental signaling pathways[5]. Zheng (2019) analyzed the pituitaries of Hu sheep with high and low fertility and found 57 differentially expressed lncRNAs and 298 differentially expressed mRNAs[6]. Miao et al. (2016, 2017) analyzed the ovaries of Small-tailed Han and Dorset sheep strains and found that differentially expressed lncRNAs were significantly enriched in the oxytocin signaling pathway. Methylation of lncRNAs might contribute to improving the reproduction of Small-tailed Han sheep[3,4]. Feng (2018) identified 76 differentially expressed mRNAs and 5 differentially expressed lncRNAs by analyzing the ovaries of Hu sheep with high and low reproduction rate[7].These studies showed that lncRNAs in the pituitary and ovaries of sheep have regulatory functions in reproduction. It is known that the sheep adrenal gland also has an impact on reproduction[8-12], but studies on the functions of lncRNAs in this organ are limited.

In this study, the differentially expressed genes in the adrenal gland between Small-tailed Han sheep with FecBBB(MM) and Small-tailed Han sheep with FecB++ (WW) (hitherto simply MM and WW sheep) were analyzed using RNA sequencing (RNA-Seq). The molecular mechanisms of differentially expressed lncRNAs and genes in the adrenal gland related to reproduction were subjected to a preliminary exploration. Our results provide an effective theoretical basis for studying the molecular mechanisms by which lncRNAs regulate sheep reproduction.

Results

Transcript assembly and quality control

The RNA-seq data of 12 samples were subjected to quality control, and the results are shown in Table 1. The clean reads for each sample ranged from 80 to 100 million and the Q30 values ranged from 92.50% to 95.30%. About 89.47%–91.63% of the clean reads were mapped to the sheep reference genome, and about 80% were mapped uniquely.

LncRNA identification and characterization

A total of 17201 candidate lncRNAs was identified, including 1174 anti-sense lncRNAs, 10636 intronic lncRNAs and 5391 large intergenic (lin) cRNAs (Fig. 1a). As shown in Fig. 1b, most lncRNAs had two exons: significantly fewer than the exons of mRNAs. The expression levels of mRNAs and lncRNAs were further analyzed according to FPKM values, and the boxplot (Fig. 1c) shows that the expression levels of mRNAs in the adrenal tissues were higher than lncRNAs. The distributions of lncRNA and mRNA lengths were consistent..

Expression levels of genes and differentially expressed analysis

The Fig. 1C box plot shows that lncRNA transcript expression levels were all lower than those of mRNAs in adrenal of both MM and WW Small-tailed Han sheep. Based on P values <0.05, differentially expressed gene analysis of adrenal tissues showed that there were 15 lncRNAs upregulated and 23 downregulated in MM sheep; 279 mRNAs were upregulated and 266 downregulated in MM sheep. In WW sheep, 354 lncRNAs were upregulated and 159 downregulated; 1334 mRNAs were upregulated and 1147 downregulated (Supplemental Table S1 & S2).

GO annotation and KEGG enrichment analysis of differentially expressed mRNAs

The GO enrichment analysis of all differentially expressed mRNAs showed that they were categorized into BP, CC, and MF. The differentially expressed mRNAs in MM and WW sheep associated with reproduction were mainly annotated in biological processes, and these were involved in reproduction processes, mating behavior, ovarian follicle cell development, the mitogen-activated protein kinase (MAPK) cascade and steroid metabolic processes in MM sheep (Fig. 2a, Supplemental Table S3). In WW sheep, these were involved in the Wnt receptor signaling pathway, meiosis, post-embryonic development, asexual reproduction and responses to steroid hormone stimuli (Fig. 2b, Supplemental Table S4).

The KEGG enrichment analysis showed that the majority of differentially expressed mRNAs were enriched in the same pathways in both the MM and WW sheep. The differentially expressed mRNAs were enriched in progesterone-mediated oocyte maturation, the MAPK signaling pathway, circadian rhythm, oocyte meiosis and insulin signaling pathway (Fig. 2c & d, Supplemental Table S5).

GO and KEGG enrichment analysis of target genes of lncRNAs

LncRNAs target genes were analyzed according to GO annotation and KEGG enrichment. As shown by the GO annotation results, lncRNA target genes associated with reproduction were also mainly annotated in biological processes in the MM and WW sheep. As shown in Fig. 3, both cis-target genes and trans-target genes associated with reproduction were enriched in biological processes in the MM type and WW sheep. In the MM sheep, cis-targets were mainly enriched for female gamete generation, reproduction and oogenesis, and trans-targets were enriched for the Wnt receptor signaling pathway, reproduction and circadian rhythm (Fig. 3a). Cis-targets and trans-targets were all mainly enriched for reproduction, gamete generation and ovarian follicle cell development (Fig. 3b).

The KEGG enrichment analysis associated with reproduction showed that the majority of cis-target genes and trans-target genes were enriched in the same pathways in MM and WW sheep (Fig. 3). In MM sheep, the lncRNA targets were assigned to 20 reproduction pathways, such as Wnt signaling pathway, transforming growth factor beta (TGFb) signaling pathway, ovarian steroidogenesis, MAPK signaling pathway, circadian rhythm and other (Fig. 3c). The lncRNA targets were enriched in progesterone-mediated oocyte maturation, MAPK signaling pathway, circadian rhythm, oocyte meiosis and insulin signaling pathway (Fig. 3d; Supplementary Table S6).

LncRNA and mRNA co-expression network analysis

In the MM sheep, a lncRNA/mRNA co-expression network was constructed using 46 differentially expressed lncRNAs and 17 target genes involved in reproductive-related pathways. As shown in Figure 4, some lncRNA targets are located in the center of the network, for example BAD, PLCB3, MYL6, MEL6B, CACNA2S, CSNK1A1, PRKG1 and HIF1A. In the WW sheep, a lncRNA-mRNA co-expression network was constructed for reproductive-related pathways using 29 differentially expressed lncRNAs and 12 target genes. As shown in Figure 5, some lncRNA targets are located in the center of the network, for example MAP3K11, ANAPC11, FKBP5, FLNB and PDPK1. The network model shows that reproductive-related lncRNA targets are co-expressed with lncRNAs, indicating that lncRNA and mRNA are mutually regulated during reproduction.

RNA-Seq data validation by qPCR

To further validate the sequencing data, 8 differentially expressed lncRNAs and 7 differentially expressed lncRNA target genes were selected to detect expression levels by RT–qPCR (Supplementary Table S7). As shown in Fig. 6, BAD, XLOC_397965, XLOC_2882304, XLOC_212066 and XLOC_108057 were upregulated in the adrenal glands of the MM sheep and PER3, SMAD1, TGFBR2 were downregulated in MM sheep. NTRK2, SMAD1, CREB1, XLOC_274809, XLOC_374827, XLOC_381122 and XLOC_393803 were upregulated in the WW sheep. These results were consistent with the transcriptome data results.

Discussion

It is known that lncRNAs play an important regulatory role in sheep reproduction. Functional lncRNAs have been identified in the brain, heart, kidney, liver, lung, ovary, skin, white adipose tissue, and pituitary in sheep[13]. In addition, they have also been identified in human, mouse and pig uterine tissue[14,15]. The adrenal gland, hypothalamus and pituitary were concluded in the hypothalamic–pituitary–adrenal axis (HPA), and the HPA axis interacts with the hypothalamic–pituitary–gonadal (HPG) axis at the brain and pituitary levels to maintain a balance between sheep reproduction and survival. However, little research has been done on lncRNAs and their targets in sheep adrenal glands.

GO and KEGG enrichment analysis indicated that the differentially expressed lncRNA targets were mainly involved in the Wnt signaling pathway, TGF-b signaling pathway, ovarian steroidogenesis, MAPK signaling pathway and circadian rhythm. Analysis of differential lncRNA–mRNA co-expression patterns and functions of target genes revealed that lncRNA might affect the fecundity of sheep by modulating genes associated with the above signaling pathways and biological processes. Among the MM sheep, these pathways were enriched for four differentially expressed lncRNAs (XLOC_212877, XLOC_357966, 105614839, XLOC_154131) and four lncRNA targets (TGFBR2, PER3, SMAD1, LHCGR). In the WW sheep, these pathways were enriched in three differentially expressed lncRNAs (XLOC_254761, 105603287, 105615642) and 3 lncRNA targets (CREB1, BRAF, NTRK2).

The HPA axis can be regulated by stress hormone signaling. The paraventricular nucleus in the hypothalamus is activated to release corticotropin releasing hormone (CRH), and CRH stimulates the release of adrenocorticotropic hormone (ACTH). ACTH stimulates the adrenal gland to release glucocorticoids[16,17]. The HPA can interact with the HPG. For example, estradiol secreted by the ovaries can enhance the activity of the HPA axis. At the level of hypothalamus, CRH inhibits the secretion of gonadotropin releasing hormone (GnRH). In the pituitary, luteinizing hormone levels were significantly reduced in mice chronically exposed to corticosterone[18]. At the level of the adrenal gland, women with congenital adrenal hyperplasia can suffer from impaired fertility or impaired steroid secretion caused by increased androgen level[19]. High concentrations of glucocorticoids have an inhibitory effect on the activities of GnRH-secreting neurons, pituitary gonadotropins and gonads[12]. In addition, glucocorticoids inhibit thyroid stimulating hormone secretion and reduce the conversion of T3 (inactive thyroxine) to T4 (effective triiodothyronine) during stress. Endogenous thyrotropin-releasing hormone and TSH secretion levels can be inhibited by CRH-induced somatostatin[9]. In sheep, acute stress inhibits the release of LH from the pituitary by inhibiting the synthesis of GnRH and GnRH receptors and promotes functional enhancement of gonadotropin-inhibiting hormone-secreting neurons[20-22]. Clearly, the adrenal gland plays an important role in the HPG.

Cyclic AMP response element binding protein 1 (encoded by CREB1) is a member of the CREB family, which is present in different ovarian compartments, including follicular granulosa cells[23-25]. CREB protein concentration increases during sexual maturation and ovarian follicular development[26,27]. CREB1 can bind to the estrogen receptor alpha to play a key role in chicken sexual maturation[28]. The gene encoding serine/threonine kinase (BRAF) is a proto-oncogene that has been identified previously as a candidate target gene in human endometriosis. CREB1 is a potential transcription factor of BRAF. Compared with normal endometrial tissue, the mRNA expression levels of BRAF and CREB1 are significantly upregulated in endometrial tissues of patients with endometriosis; thus, CREB1 can increase the expression of BRAF and regulate cell proliferation by binding directly to BRAF in such women[29]. One study determined the role of CREB1 in mouse granulocytes (mGCs) by knocking down the expression of CREB1. It was found that CREB1 could be regulated by steroid synthesis, cell proliferation, cell cycle, apoptosis and other follicular factors as key regulators of mGCs[30]. Moreover, XLOC_254761 can transregulate CREB1, and 10563287 trans-regulated BRAF. CREB1 and BRAF are enriched in progesterone-mediated oocyte maturation. These results indicated that those genes enriched in progesterone-mediated oocyte maturation may promote follicular maturation and regulate reproduction in Small-tailed Han sheep.

Using six estrous Hu sheep ovaries and six non-estrous China merino sheep ovaries to identify genes associated with off-season reproduction, the gene encoding neurotrophic receptor tyrosine kinase 2 (NTRK2) was found to be differentially upregulated in the ovaries[31]. The NTRK2 gene encodes the NTRK2 receptor for neurotrophin-4/5 and brain-derived neurotrophic factor leading to oocyte death in late adolescence, loss of follicular tissue and infertility in early adulthood. The preovulatory gonadotropin surge promotes oocyte survival at the onset of reproductive cyclicity by inducing oocyte expression of NTRK2 in Hu sheep[32]. 105615642 can transregulate NTRK2, which is enriched in the MAPK signaling pathway. The results suggest that NTRK2 might regulate follicular survival.

The expression level of the gene encoding luteinizing hormone receptor (LHCGR) is low in mural granulosa cells and cumulus cells of antral follicles, and LHCGR is induced in granulosa cells by follicle stimulating hormone (FSH). The synthesis of LHCGR is mediated by retinoic acid (RA), the demethylation of its promoter region is a key mechanism regulating cell type-specific differentiation during follicular development[33].

Mothers against decapentaplegic homolog 1 (SMAD1) is a downstream signal transduction element of bone morphogenetic proteins[34,35], and SMAD1 responds to bone morphogenetic protein 15 (BMP15), which is important for proliferation of granulosa cells in sheep[36]. The transforming growth factor-beta (TGFb) subfamily is encoded by TGFB1, TGFB2 and TGFB3. In mammals, there are seven TGF-1 receptor subtypes (encoded by TGFBR1) and five type 2 receptor subtypes (TGFBR2) associated with signal transduction ligand of the TGFβ superfamily[37,38]. The TGFb superfamily ligand initiates an intracellular signaling pathway upon binding to a cell surface receptor complex. In addition, TGFBR2 in the TGF-b signaling pathway were found to be associated with litter size in pigs[39]. There was a significant correlation between litter size in mutant indigenous Chinese Hu sheep and mutations in the three genes encoding the TGFβ superfamily (FecB, GDF9, and TGFBR2); moreover, the FecB, GDF9 and TGFBR2 polymorphisms are considered to be potentially important genetic markers in marker-assisted selection (MAS) strategies to increase litter size in these sheep)[40]. Moreover, here we show that 105614839 and XLOC_212877 cis-regulate SMAD1 and TGFBR2, respectively. These genes are enriched in the TGF-b signaling pathway. The results indicate that SMAD1 and TGFBR2 might regulate reproduction in terms of uterine decidual function and litter size.

Conclusion

In summary, the adrenal gland plays a key role in sheep female reproductive processes. For example, the adrenal gland affects the reproductive of sheep through hypothalamic-pituitary-adrenal axis. These functions of the adrenal gland are achieved by regulating different signaling pathways and related genes. In this study, we screened the differentially expressed lncRNAs (XLOC_212877, XLOC_357966, XLOC_154131, XLOC_254761, 105614839, 105603287 and 105615642) and differentially expressed lncRNA target genes (TGFBR2, PER3, SMAD1, LHCGR, CREB1, BRAF and NTRK2) associated with sheep prolificacy, and constructed networks of interactions between lncRNAs and mRNAs. The lncRNA and mRNAs associated with sheep breeding were enriched by KEGG analysis. The results of the study will help to elucidate the regulatory mechanisms of lncRNAs in sheep reproduction.

Methods

Ethics statement

All the animals were authorized by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS; Beijing, China). In addition, ethical approval of animal survival was given by the animal ethics committee of IAS-CAAS (No. IAAS 2019-449, 18 September 2019).

Animals and sample collection

Animals were from the core group of Small-tailed Han sheep in the Luxi area of Shandong Province, P. R. China. All were free to eat and drink. We chose healthy nonpregnant sheep aged 2–4 years. Jugular vein blood was collected and the FecB mutation was identified by TaqMan assays. Six WW and six MM sheep were used for the experiments. All experiments complied with the rules established by the the animal ethics committee of IAS-CAAS (No. IAAS 2019-449, 18 September 2019).

All sheep were treated with vaginal sponges (InterAg Co., Ltd., New Zealand) (progesterone 300mg) for 12 days to synchronize estrus. Three WW and three MM ewes were euthanized on the 50th hours after removing the vaginal sponges, and the adrenal glands were collected (follicular phase; WF and MF, respectively), and of the remaining six sheep (three in each group) were euthanized on the 7th day after sponge removal (luteal phase; WL and ML, respectively) and the adrenal glands were collected. All samples were stored immediately at -80 °C for total RNA extraction. All experiments complied with the rules established by the the animal ethics committee of IAS-CAAS (No. IAAS 2019-449, 18 September 2019).

RNA isolation, library preparation and sequencing

Total RNA was extracted from the adrenal glands of all 12 ewes using TRIzol (Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer’s instruction. RNA purity was checked using a Nano Photometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA concentrations were measured using Qubit® RNA Assay kits in a Qubit® 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay kits of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

The lncRNA library we generated was chain-specific. The method of reverse transcription synthesis of the first strand of cDNA is the same as the New England Biolabs (NEB) general library construction (NEB, Ipswich, MA, USA). The difference is that when the second strand is synthesized, the dTTP in the dNTPs is replaced by dUTP, followed by cDNA end repair, addition of an A tail, ligation sequencing and length screening. Then we used the USER enzyme (New England Biolabs, Inc., Ipswich, MA, USA) to degrade the second strand of U-containing cDNA and performed polymerase chain reaction (PCR) amplification to obtain a library. After the lncRNA library was established, the library preparations were sequenced on an Illumina Hiseq platform (Illumina, San Diego, CA, USA). Raw data of the performed RNA-seq have been recorded in the SRA public database (Accession number: SRP222893).

Reference genome mapping and transcriptome assembly

Raw reads were obtained by removing reads with an adapter, reads containing > 10% poly-N sequences and low-quality reads. At the same time, Q20, Q30 and GC contents of the clean data were calculated. An index of the reference genome was built using bowtie2 v. 2.2.8 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)[41] and paired-end clean reads were aligned to the reference genome Oar_v3.1 using HISAT2 v. 2.0.4 (https://ccb.jhu.edu/software/hisat2/index.shtml)[42].

The mapped reads of each sample were assembled using StringTie v. 1.3.1 (https://ccb.jhu.edu/software/stringtie/)[42] in a reference-based approach. StringTie uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus.

Identification of potential lncRNA candidates

Potential lncRNA candidates were identified as follows: selecting transcripts with exon number ³2 and transcript length > 200 nt. The transcripts obtained by combining each sample using Cuffdiff v. 2.1.1 (https://software.broadinstitute.org/cancer/software/genepattern/modules/docs/Cuffcompare/8)[43], and transcripts with uncertain chain direction were removed to obtain complete transcriptome information.

The intersection of transcripts with no coding potential in these software analysis results was used as the lncRNA data set. CNCI v. 2.0 (https://github.com/www-bioinfo-org/CNCI)[44], CPC v. 2.81 (http://cpc.cbi.pku.edu.cn/)[45] and PFAM v. 1.3 (https://pfam.xfam.org)[46] were used to evaluate the protein coding potential of the transcript and the results of these software as the end result.

Analysis of differentially expressed genes

Gene expression estimation using the FPKM value (expected number of Fragments Per Kilobase of transcript sequence per million base pairs sequenced) eliminates the effect of sequencing depth and gene length on fragment counts. Different types of transcripts (lncRNA and mRNA) were analyzed for differences using Cuffdiff v. 2.1.1 (https://software.broadinstitute.org/cancer/software/genepattern/modules/docs/Cuffdiff/7)[43]. LncRNA and mRNA with P values < 0.05 were considered as differentially expression between the two groups of data.

GO annotation and KEGG pathway enrichment analysis of differentially expressed genes

The GO is an international standard system including molecular functions, biological processes, and cellular components for classifying gene function. Pathway enrichment analysis can identify major metabolic pathways and signaling pathways enriched by differentially expressed genes. Enrichment analysis was performed on each Pathway in KEGG using a hypergeometric test. The calculated P value and 0.05 being defined as the significant threshold, the genes were screened and enriched for the pathways. Next, the significance of the pathway enrichment analysis was corrected by FDR, and the corrected P-value (Q-value) was obtained. Differentially expressed genes were further studied using the GO and KEGG databases to study the functions of the genes and identify the pathways in which they participate. If a P value was≤0.05, enrichment was considered significant.

GO annotation and KEGG pathway enrichment analysis of differentially expressed lncRNA targets

Differentially expressed lncRNAs regulate the target genes by cis-regulating nearby genes and trans-regulating distal protein-coding genes. Here, protein-coding genes with a distance of less than 100 Kb were assumed to be the cis-target genes, and Pearson correlation coefficients with the lncRNAs of >0.95 were assumed to represent the trans-target genes[47]. Statistical enrichment of differentially expressed lncRNA target genes in the GO annotation and KEGG pathways were evaluated, and P value£ 0.05 was considered significant enrichment.

Construction of lncRNA/mRNA networks

To infer the functions of differentially expressed lncRNAs and their target genes in sheep prolificacy, we constructed a network based on mRNA and lncRNA in Cytoscape v. 3.1.1 (https://cytoscape.org).

Reverse transcription (RT)–qPCR verification

RT–qPCR was used to verify the expression levels of differentially expressed lncRNAs and their targets. About 0.1 mg of RNA was used per sample and this was reverse transcribed into cDNA using RT reagents (Thermo Fisher Scientific, Waltham, MA, USA). All experiments were performed in triplicate, and β-actin was used as an internal reference to normalize target gene expression. The qPCR was performed on a LightCycler 480II (Roche, Basel, Sweden) using SYBR Premix Ex Taq II. The procedure involved 40 cycles of pre-denaturation at 95 °C for 5 s; denaturation at 95 °C for 5 s, 60 °C for 30 s. After the reaction was completed, melting curve analysis was performed. The relative expression level of the target genes was calculated by the 2DDCt method[48]. The lncRNAs and target gene primers are shown in Supplemental Table S1.

Statistical analyses

All data are presented as the mean ± SD. Student’s t tests were performed to compare means, and P < 0.05 was considered statistically significant.

Declarations

Acknowledgements

We thank James Cummins, PhD, from EDANZ Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Authors’ contributions

These studies were designed by QX and MXC, who performed the experimental analyses and prepared the figures and tables. QX and QLL analyzed the data and drafted the manuscript. MXC contributed to revisions of the manuscript. XSZ, JLZ, XFG and MXC assisted in interpreting the results and revised the final version of the manuscript. All authors read and approved the final manuscript for publication.

Funding

This research was funded by the following bodies: National Natural Science Foundation of China (31472078 and 31772580); Earmarked Fund for China Agriculture Research System (CARS-38), Central Public-interest Scientific Institution Basal Research Fund (2018-YWF-YB-1, Y2017JC24, 2017ywf-zd-13); Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13); China Agricultural Scientific Research Outstanding Talents and Their Innovative Teams Program, China High-level Talents Special Support Plan Scientific and Technological Innovation Leading Talents Program (W02020274); Joint Funds of The National Natural Science Foundation of China and The Government of Xinjiang Uygur Autonomous Region of China (U1130302); Tianjin Agricultural Science and Technology Achievements Transformation and Popularization Program (201704020); The Youth Innovative Research and Experimental Project of Tianjin Academy of Agricultural Sciences (201915). The APC was funded by the National Natural Science Foundation of China (31472078). The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

The following information was supplied regarding data availability: Data is available at the Sequence Read Archive (SRA222893):

Ethics approval and consent to participate

All the Small Tail Han sheep were euthanized and were authorized by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS; Beijing, China). In addition, ethical approval of animal survival was given by the animal ethics committee of IAS-CAAS (No. IAAS 2019-449, 18 September 2019).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Key Laboratory of Animal Genetics and Breeding and Reproduction of the Ministry of Agriculture and Rural Affairs, Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, P. R. China

2College of Life Sciences, Langfang Normal University, Langfang, 065000, P.R. China

3State Key Laboratory of Sheep Genetic Improvement and Healthy Production, Xinjiang Academy of Agricultural and Reclamation Sciences, Shihezi, 832000, P. R. China

4Tianjin Institute of Animal Sciences, Tianjin 300381, P. R. China

Abbreviations

LncRNAs: long noncoding; MAPK: the mitogen-activated protein kinase cascade; HPA: the hypothalamic–pituitary–adrenal axis; HPG: the hypothalamic–pituitary–gonadal axis; CRH: corticotropin releasing hormone; ACTH: the release of adrenocorticotropic hormone; GnRH: gonadotropin releasing hormone; CREB1: Cyclic AMP response element binding protein 1; BRAF: serine/threonine kinase; mGCs: mouse granulocytes; NTRK2: neurotrophic receptor tyrosine kinase 2; LHCGR: luteinizing hormone receptor; FSH: follicle stimulating hormone; SMAD1: Mothers against decapentaplegic homolog 1; BMP15: bone morphogenetic protein 15; TGFBR1: TGF-1 receptor subtypes; TGFBR2: type 2 receptor subtypes; MAS: marker-assisted selection.

References

[1] Zhao Y. Sheep Production Science, 3rd Edition. China Agriculture Press, 2011.

[2] Chu M, Jia L, Zhang Y, Jin M, Chen H, Fang L, Di R, Cao G, Feng T, Tang Q, Ma Y, Li K. Polymorphisms of coding region of BMPR-IB gene and their relationship with litter size in sheep. Mol Biol Rep, 2011, 38(6):4071-4076.

[3] Miao X, Luo Q, Zhao H, Qin X. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries. Sci Rep, 2016; 6(1):39398.

[4] Miao X, Luo Q, Zhao H, Qin X. An integrated analysis of miRNAs and methylated genes encoding mRNAs and lncRNAs in sheep breeds with different fecundity. Front Physiol, 2017; 8:1049.

[5] Yang H, Wang F, Li F, Ren C, Pang J, Wan Y, Wang Z, Feng X, Zhang Y. Comprehensive analysis of long non-coding RNA and mRNA expression patterns in sheep testicular maturation. Biol Reprod, 2018;99(3):650-661.

[6] Zheng J, Wang Z, Yang H, Yao X, Yang P, Ren C, Wang F, Zhang Y. Pituitary transcriptomic study reveals the differential regulation of lncRNAs and mRNAs related to prolificacy in different FecB genotyping sheep. Genes (Basel), 2019: doi: 10.3390/genes10020157.

[7] Feng X, Li FZ, Wang F, Zhang GM, Pang J, Ren CF, Zhang TT, Yang H, Wang ZY, Zhang YL. Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. Bioscience Reports, 2018; doi.org/10.1042/BSR20171350.

[8] Jimena P, Castilla JA, Peran F, Ramirez JP, Jr FV, Molina R, Vergara F, Herruzo A. Adrenal hormones in human follicular fluid. Acta Endocrinologica, 1992;127(5):403.

[9] Magiakou MA, Mastorakos G, Webster E, Chrousos GP. The hypothalamic-pituitary-adrenal axis and the female reproductive system. Ann N Y Acad Sci, 2010; 816(1):42-56.

[10] Gao Y, Chen F, Kong QQ, Ning SF, Yuan HJ, Lian HY, Luo MJ, Tan JH. Stresses on female mice impair oocyte developmental potential: effects of stress severity and duration on oocytes at the growing follicle stage. Reprod Sci, 2016;23(9):1148-1157.

[11] Yuan HJ, Han X, Nan H, Guo LW, Shuai G, Juan L, Min G, Tian JH. Glucocorticoids impair oocyte developmental potential by triggering apoptosis of ovarian cells via activating the Fas system. Sci Reprod, 2016; 6:24036.

[12] Whirledge S, Cidlowski JA. Glucocorticoids and reproduction: traffic control on the road to reproduction. Trends Endocrinol Metab, 2017; 28(6):399-415.

[13] Bakhtiarizadeh MR, Hosseinpour B, Arefnzhad B, Shamabadi N, Salami SA. In silico prediction of long intergenic non-coding RNAs in sheep. Genome, 2016;59(4):263-275.

[14] Zhou C, Zhang T, Liu F, Zhou J, Ni X, Huo R, Shi R. The differential expression of mRNAs and long noncoding RNAs between ectopic and eutopic endometria provides new insights into adenomyosis. Mol Biosyst, 2015;12(2):362-370.

[15] Wang Y, Xue S, Liu X, Hu T, Qiu X, Zhang J, Lei M. Analyses of long non-coding RNA and mRNA profiling using RNA sequencing during the pre-implantation phases in pig endometrium[J]. Sci Reprod, 2016;6:20238.

[16] Herbison AE, Robert P, Jean-Remi P, Mora JM, Hurst PR. Gonadotropin-releasing hormone neuron requirements for puberty, ovulation, and fertility. Endocrinology, 2008;149(2):597-604.

[17] Lehman MN, Hileman SM, Goodman RL. Neuroanatomy of the kisspeptin signaling system in mammals: comparative and developmental aspects. Advances in Experimental Med Biol, 2015;784:27-62.

[18] Rivier C, Vale W. Corticotropin-releasing factor (CRF) acts centrally to inhibit growth hormone secretion in the rat. Endocrinology, 1984;114(6):2409-2411.

[19]Gomes LG, Bachega TASS, Mendonca BB. Classic congenital adrenal hyperplasia and its impact on reproduction. Fertil and Steril, 2019;111(1):7-12.

[20] Dobson H, Fergani C, Routly JE, Smith RF. Effects of stress on reproduction in ewes. Animal Reprod Sci, 2012;130(3-4):135-140.

[21] Ciechanowska M, Lapot M, Antkowiak B, Mateusiak K, Paruszewska E, Malewski T, Paluch M, Przekop F. Effect of short-term and prolonged stress on the biosynthesis of gonadotropin-releasing hormone (GnRH) and GnRH receptor (GnRHR) in the hypothalamus and GnRHR in the pituitary of ewes during various physiological states. Anim Reprod Sci, 2016;174:65-72.

[22]Clarke IJ, Bartolini D, Conductier G, Henry BA. Stress increases gonadotropin inhibitory hormone cell activity and input to gnRH cells in ewes. Endocrinology, 2016;157(11):en20161513.

[23] Gubbay O, Rae MT, McNeilly AS, Donadeu FX, Zeleznik AJ, Hillier SG. cAMP response element-binding (CREB) signalling and ovarian surface epithelial cell survival. J Endocrinol, 2006;191(1):275-285.

[24] Pei JH, Fujimoto Y, Yamauchi N, Hattori MA. Real-time monitoring of cAMP response element binding protein signaling in porcine granulosa cells modulated by ovarian factors. Mol Cell Biochem, 2006;290(1-2):177-184.

[25] Suman R, Androulla E, Zara J, Laura P, Maon HD. Metformin inhibits follicle-stimulating hormone (FSH) action in human granulosa cells: relevance to polycystic ovary syndrome. J Clinical Endocrinol Metab, 2013;98(9):E1491-E1500.

[26] Natalie YO, Noa S, Naomi MB, Sarah E, Moriah K, Manna PR, Stocco DM, Joseph O. Transcription of steroidogenic acute regulatory protein in the rodent ovary and placenta: alternative modes of cyclic adenosine 3', 5'-monophosphate dependent and independent regulation. Endocrinology, 2009;150(2):977-989.

[27] Sirotkin AV, Gupta K, Kapoor R, Dwivedi A. MicroRNA-145 targets Smad1 in endometrial stromal cells and regulates decidualization in rat. J Med, 2019;97(4):509-522.

[28] Guo M, Li Y, Chen Y, Guo X, Yuan Z, Jiang Y. Genome-wide mapping of estrogen receptor α binding sites by ChIP-seq to identify genes related to sexual maturity in hens. Gene, 2018, 642:32-42.

[29] Lv X, Wang D, Ma Y, Long Z. Analysis of the oncogene BRAF mutation and the correlation of the expression of wild-type BRAF and CREB1 in endometriosis. Int J Mol Med, 2018;41(3):1349-1356.

[30] Zhang P, Wang J, Lang H, Wang W, Liu X, Liu H, Tan C, Li X, Zhao Y, Wu X. Knockdown of CREB1 promotes apoptosis and decreases estradiol synthesis in mouse granulosa cells. Biomed Pharmacother, 2018;105:1141-1146.

[31] Lei C, Liu K, Zhao Z, Blai HT, Peng Z, Li D, Ma RZ. Identification of sheep ovary genes potentially associated with off-season reproduction. J Genet Genomics, 2012;39(4):181-190.

[32] Dorfman MD, Cecilia GR, Zefora A, Bredford K, Alejandro L, Dissen GA, Juan MC, David GG, Francisco G, Baoji X. Loss of Ntrk2/Kiss1r signaling in oocytes causes premature ovarian failure. Endocrinology, 2014;155(8):3098-3111.

[33] Kawai T, Richards JS, Shimada M. The cell type–specific expression of Lhcgr in mouse ovarian cells: evidence for a DNA-demethylation–dependent mechanism. Endocrinology, 2018;159(5):2062-2074.

[34] Caterina C, Tripurani SK, Large MJ, Edson MA, Creighton CJ, Hawkins SM, Ertug K, Vasa K, Lydon JP, Pangas SA. Activin-like kinase 2 functions in peri-implantation uterine signaling in mice and humans. PLoS Genet, 2013;9(11):1245-1245.

[35] Kim MR, Park DW, Lee JH, Choi DS, Hwang KJ, Ryu HS, Min CK. Progesterone-dependent release of transforming growth factor-beta1 from epithelial cells enhances the endometrial decidualization by turning on the Smad signalling in stromal cells. Mol Hum Reprod, 2005; 11(11):801-808.

[36] Moore RK, Otsuka F, Shimasaki S. Molecular basis of bone morphogeneric protein-15 signaling in gramulosa cells. J Biol Chem, 2003;278(1):304-310.

[37] Claire G, Kemp CF, Knight PG. Bone morphogenetic protein (BMP) ligands and receptors in bovine ovarian follicle cells: actions of BMP-4, -6 and -7 on granulosa cells and differential modulation of Smad-1 phosphorylation by follistatin. Reproduction, 2004;127(2):239-254.

[38] Knight PG, Glister C. TGF-beta superfamily members and ovarian follicle development. Reproduction, 2006;132(2):191-206.

[39]Li XJ, Ye JW, Han X, Qiao RM, Li XL, Lv G, Wang KJ. Whole-genome sequencing identifies potential candidate genes for reproductive traits in pigs. Genomics, 2019; doi.org/10.1016/j.ygeno.2019.01.014.

[40] Wang Q,Wang N, Cai R, Zhao F, Xiong Y, Li X, Wang A, Lin P, Jin Y. Genome-wide analysis and functional prediction of long non-coding RNAs in mouse uterus during the implantation window. Oncotarget, 2016;8(48):84360-84372.

[41] Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods, 2012;9(4):357-359.

[42] Pertea M, Kim D, Pertea GM, Leek JT, Salzerg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc, 2016;11(9):1650-1667.

[43] Trapnell C, Wiliams BA, Pertea G, Mortazavi A, Kwan G, ven Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol, 2010;28(5):511-515.

[44] Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res, 2013;41(17):e166-e166.

[45]Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, Gao G. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res, 45(Web Server issue): 2017; W12–W16.

[46] Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, Griffiths-Jones S, Howe KL, Marshall M, Sonnhammer EL. The Pfam protein families database. Nucleic Acids Res, 2002, 30(1):276-280.

[47] Alessandro F, Irene B. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet, 2014;15(1):7-21.

[48] Schmittgen TD, Libak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc, 2008;3(6):1101-1108.

Tables