Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and miRNA of Yili goose ovary at different egg-laying stages

doi:10.21203/rs.3.rs-1527604/v1

Download PDF

Research Article

Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and miRNA of Yili goose ovary at different egg-laying stages

https://doi.org/10.21203/rs.3.rs-1527604/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background: The expression of ovarian tissue at the transcriptional level of goose at different egg laying stages could provide an important basis for screening and identifying key genes regulating ovarian development.

Results: In this study, 4 ovary tissues at each breeding period of prelaying period (KL), laying period (CL), and ceased period (XL), respectively, with significant morphology difference, were used for RNA extraction and mRNAs, lncRNAs, and miRNAs comparison in Yili geese. CeRNA regulatory network was constructed for key genes screening. A total of 337, 1136, and 525 differentially expressed DE mRNAs, 466, 925, and 742 DE lncRNAs and 258, 1131 and 909 DE miRNAs were identified between KL and CL, between XL and CL, and between XL and KL groups, respectively. Functional enrichment analysis showed that the differentially expressed mRNAs and non-coding RNA target genes were mainly involved in the cell process, cytokine-cytokine receptor interaction, phagosome, calcium signaling pathway, steroid biosynthesis and ECM-receptor interaction. Differential genes and non-coding RNAs, PDGFRB, ERBB4, LHCGR, MSTRG.129094.34, MSTRG.3524.1 and gga-miR-145-5p, related to reproduction and ovarian development were highly enriched. Furthermore, lncRNA-miRNA-mRNA regulatory networks related to ovary development were constructed.

Conclusions: Our study found dramatic transcriptomic differences in ovaries of Yili geese at different egg-laying stages, and a differential lncRNA-miRNA-mRNA regulatory network related to cell proliferation, differentiation and apoptosis and involved in stromal follicle development were established and preliminarily validated, which could be regarded as a key regulatory pathway of ovarian development in Yili goose.

Yili goose

Ovary

Laying period

mRNA

Non-coding RNA

Regulation networks

The development of ovary and follicle is an important factor affecting the laying performance of geese [1]. The ovarian volume of female geese is smaller at non-reproductive period than that of the ovary during the laying period. No hierarchical follicles and many pre-hierarchical follicles could be detected of the ovaries at non-reproductive period. The ovary stay degenerated until the next egg-laying period [2]. The biological process of ovarian development and ovulation is transcriptionally regulated by a large number of key genes under dynamic and stage-specific expression [3]. Therefore, studying the expression characteristics of ovarian tissue at the transcriptional level in goose egg laying stages could provide an important basis for screening and identifying key genes regulating goose ovarian development. Yili goose is one of the typically seasonal breeding poultry in China which provided it an ideal animal model to study the gene expression pattern of geese in each egg laying period.

Transcriptome of Zhedong white geese ovary at laying and nesting period showed that differential genes were mainly involved in norepinephrine metabolism, steroid hormone biosynthesis, Wnt signaling pathway, calcium signaling pathway, GnRH signaling pathway and oocyte meiosis, which were closely related to follicular development[4]. Based on high-throughput sequencing technology, the transcriptome sequencing of the gonadal axis tissues in high and low egg production Xinjiang Yili goose resulted in 30 candidate genes related to egg production[5]. However, the regulation mechanism, especially the non-coding RNA involved regulatory network of reproductive traits in goose is still unclear.

In this study, Yili female geese at prelaying, laying and ceased period were selected for exploring the candidate genes and network that regulate the ovarian development by using RNA-seq and small RNA-seq technology. The results could provide new insights into the molecular regulatory mechanism of non-coding RNA-mediated ovary development in goose.

Ovary histomorphometric analysis at different egg laying periods in Yili goose

During the prelaying period, the ovaries of Yili geese showed primary follicles are numerous, while the ovaries during the laying period showed a hierarchy of follicles, and goose at ceased period showed degenerated ovaries (Fig. 1A-B). Moreover, the ovarian weight showed highly significant differences between periods (P < 0.01) (Fig. 1C).

Overview of sequencing and identification of lncRNA and miRNA

A total of 57,811,186 ~ 85,328,377 clean reads (Supplementary Table S2-1) were produced from the 12 ovary cDNA libraries, of which 80.85%~87.21% sequences were uniquely mapped (Supplementary Table S2-2) to produce 38,296 lncRNA and 33,332 mRNA transcripts (Supplementary Table S2-3).

For miRNA libraries, a total of 69.38%~93.34% clean reads were produced from the 12 ovary cDNA libraries (Supplementary Table S2-4), to produce 2700 known-miRNAs (Supplementary Table S2-5) and 2983 novel-miRNAs (Supplementary Table S2-6), and the proportion of annotated miRNAs was 47.51%-74.58% (Supplementary Table S2-7). Most miRNA fragments are 18-25nt in length, and the number of miRNAs with a length of 22nt is the highest in each sample, with an average of 35.44% (Supplementary Table S2-8).

Analysis of the differentially expressed mRNA, lncRNA and miRNA among different egg laying periods

A total of 337 mRNAs were differentially expressed (286 up-regulated and 51 down-regulated) in the CL as compared with KL (KL vs CL) (Fig. 2A and supplementary Table S3-1), 1136 mRNAs were differentially expressed (582 up-regulated and 554 down-regulated) in the CL as compared with XL (XL vs CL) (Fig. 2A and supplementary Table S3-2), and 525 mRNAs were differentially expressed (220 up-regulated and 305 down-regulated) in the KL as compared with XL (XL vs KL) (Fig. 2A and supplementary Table S3-3). Four differential genes were coexpressed in the three comparison groups (Fig. 2B and supplementary Table S3-4).

A total of 466 lncRNAs were differentially expressed (133 up-regulated and 333 down-regulated) in the CL as compared with KL (KL vs CL) (Fig. 3A and supplementary Table S3-5), 925 lncRNAs were differentially expressed (307 up-regulated and 618 down-regulated) in the CL as compared with XL (XL vs CL) (Fig. 3A and supplementary Table S3-6), and 742 lncRNAs were differentially expressed (347 up-regulated and 395 down-regulated) in the KL as compared with XL (XL vs KL) (Fig. 3A and supplementary Table S3-7). Six differential lncRNAs were coexpressed in the three comparison groups (Fig. 3B and supplementary Table S3-8).

A total of 258 miRNAs were differentially expressed (56 up-regulated and 202 down-regulated) in the CL as compared with KL (KL vs CL) (Fig. 4A and supplementary Table S3-9), 1131 miRNAs were differentially expressed (217 up-regulated and 914 down-regulated) in the CL as compared with XL (XL vs CL) (Fig. 4A and supplementary Table S3-10), and 909 miRNAs were differentially expressed (358 up-regulated and 515 down-regulated) in the KL as compared with XL (XL vs KL) (Fig. 4A and supplementary Table S3-11). Forty-four differential miRNAs were coexpressed in the three comparison groups (Fig. 4B and supplementary Table S3-12).

Functional enrichment analysis of differentially expressed genes and target genes

In the KL vs CL, XL vs CL, and XL vs KL groups, 204, 670, and 294 differential mRNAs were annotated by the GO database, respectively. Most of the annotated differential mRNAs were assigned to biological regulation, cellular process, and metabolic process, and also were annotated the reproductive process and reproduction (Fig. 5 and supplementary Table S4).

Differential mRNAs were significantly enriched in 5, 9 and 3 pathways in KL vs CL, XL vs CL and XL vs KL groups, respectively (q-value < 0.05). In the KL vs CL group, the differentially enriched pathways including phagosome, calcium signaling pathway (Fig. 6A and supplementary Table S5-1). In the XL vs CL group, differentially enriched pathways including steroid biosynthesis, neuroactive ligand-receptor interaction (Fig. 6B and supplementary Table S5-2). In the XL vs KL group, differentially enriched pathways including ECM-receptor interaction, neuroactive ligand-receptor interaction, etc (Fig. 6C and supplementary Table S5-3).

It was found that 38,034 lncRNAs regulate 32,352 mRNAs in their adjacent positions by cis mode; 158 lncRNAs have potential targeting regulatory relationship with 12,939 mRNAs by trans mode. Based on the prediction results of lncRNAs target genes, this study conducted GO enrichment analysis on differentially expressed lncRNA target genes. It was found that 2485, 2993 and 2847 differential lncRNAs target genes in KL vs CL, XL vs CL and XL vs KL groups were annotated by GO database, respectively. Most of the annotated differential lncRNAs target genes were assigned to biological regulation, cellular process, metabolic process, binding, catalytic activity, and also were annotated the reproductive process and reproduction (Fig. 7 and supplementary Table S6).

Differential lncRNAs target genes were significantly enriched in 4, 12 and 10 pathways in KL vs CL, XL vs CL and XL vs KL groups, respectively (q-value < 0.05). In the KL vs CL group, the differentially enriched pathways including neuroactive ligand-receptor interaction, ECM-receptor interaction, etc (Fig. 8A and supplementary Table S7-1). In the XL vs CL group, differentially enriched pathways including neuroactive ligand-receptor interaction, calcium signaling pathway, etc (Fig. 8B and supplementary Table S7-2). In the XL vs KL group, differentially enriched pathways including neuroactive ligand-receptor interaction, calcium signaling pathway (Fig. 8C and supplementary Table S7-3).

In the KL vs CL, XL vs CL and XL vs KL groups, 4129, 8178 and 7738 differential miRNAs target genes were annotated by the GO database, respectively. Most of the annotated differential miRNAs target genes were assigned to cellular process, single-organism process, organelle, membrane, binding and catalytic activity, and also were annotated the reproductive process and reproduction (Fig. 9 and supplementary Table S8).

In the KL vs CL group, the differential miRNA target genes were significantly enriched in the ECM-receptor interaction. No KEGG pathway was significantly enriched in the XL vs CL and XL vs KL groups. However, miRNAs target genes were involved in calcium signaling pathway, Notch signaling pathway and GnRH signaling pathway (Fig. 10 and supplementary Table S9).

Construction of lncRNA-miRNA-mRNA networks and functional analysis

The network of KL vs CL composed of 92 nodes and 1018 edges, and the nodes included 38 lncRNAs, 11 miRNAs and 43 mRNAs (Fig. 11A and supplementary Table S10-1). The network of XL vs CL composed of 68 nodes and 496 edges, and the nodes included 32 lncRNAs, 14 miRNAs and 22 mRNAs (Fig. 11B and supplementary Table S10-2). The network of XL vs KL composed of 149 nodes and 3882 edges, and the nodes included 55 lncRNAs, 33 miRNAs and 61 mRNAs (Fig. 11C and supplementary Table S10-3).

MicroRNA of gga-miR-145-5p, eca-miR-145, novel_miR_4043, novel_miR_560, and novel_miR_72 in the ceRNA network were participated in two groups. In XL vs KL group, 22 and 21 lncRNAs were less expressed while miRNAs of gga-miR-145-5p and eca-miR-145 were highly expressed, respectively. In XL vs CL, the prediction suggested that 32 and 25 lncRNAs were binding to gga-miR-145-5p and eca-miR-145, respectively, and showed negative regulation of target genes. Furthermore, functional enrichment analyses of mRNA involved in ceRNA network were deeply analyzed (Fig. 12 and supplementary Table S11). Genes were mostly enriched in the calcium signaling pathway, ECM-receptor interaction, oxytocin signaling pathway and oocyte meiosis pathways.

Validation of RNA-Seq Results

In the lncRNA-miRNA-mRNAs interaction network, the expression levels of 9, 7, and 6 key nodes in the KL vs CL, XL vs CL, and XL vs KL groups, respectively, were verified by RT-qPCR. (Fig. 13). The results showed that the expression trends of key node genes in each group were consistent with the RNA-seq data, indicating that the RNA-seq data were reliable. And the expression trend of miRNA in each group was opposite to that of its target lncRNA and mRNA. The above results showed that the gene expression relationship of these 22 key nodes was in line with the ceRNA hypothesis. These key node genes may be one of the important reasons affecting the ovarian development of Yili goose.

MSTRG.5970.28 Binds to gga-miR-145-5p

MSTRG.5970.28 was upregulated in XL vs CL; Targetscan and miRanda software predicted that MSTRG.5970.28 was targeted to gga-miR-145-5p. Therefore, mutant vectors were constructed to verify the specific binding sites (Fig. 14A). The results showed that compared with the NC group, gga-miR-145-5p significantly decreased the expression of luciferase in MSTRG.5970.28 wt (P < 0.01), and no significant effect on luciferase activity when co-transfected with mutated plasmids (Fig. 14B). It shows that there is a binding effect between the two in this experiment.

The ovary is an important reproductive organ of female animals, its main function is to produce and discharge healthy eggs that could be fertilized, and to secret sex hormones to maintain the sexual characteristics and periodic reproductive activities [6]. The development of animal ovary, especially follicle development and ovulation biology, is distinctly stage-specific and regulated by their own gene expression [3]. Studies have shown that mRNA and non-coding RNA have important regulatory roles in reproductive-related processes such as gonadal development, hormone regulation, sex determination, and embryo implantation [7–14].

In this study, through RNA sequencing of 12 Yili goose ovarian tissues, 82.79%~89.24% of total reads were mapped to the goose genome (Ans Cyg_PRJNA183603_v1.0), the mapping rates was higher than the reports of ovary [15] and intramuscular precursor adipocytes in chicken [16], suggested that the reference genome selected met the requirements of subsequent data analysis.

MicroRNA, which regulates more than thirty percent of protein coding genes, was regarded as a key gene regulatory factor and participates in a variety of biological and metabolic processes such as growth, development, reproduction and disease [17, 18]. The largely number of obtained miRNA of this research might participated in the regulation of ovary function. The proportion of annotated miRNAs ranged from 47.51–74.58%, slightly lower than that of goose pituitary [19], was still met the requirement of data analysis.

In this study, 337 (286 up-regulated and 51 down-regulated genes), 1136 (582 up-regulated and 554 down-regulated genes), and 525 (220 up-regulated and 305 down-regulated genes) differential genes were screened in the KL vs CL, XL vs CL, and XL vs KL groups, respectively. There were four differentially expressed genes appeared in all the three comparison groups, namely ADRA2A, CP, GPNMB and LOC106033756. These 4 overlapping differential mRNAs were significantly up-regulated when geese enters the laying period from prelaying period, and significantly down-regulated when geese enters the ceased period from the laying period. GPNMB is mainly located in the cytoplasm of ovarian cancer cells, and showed no or less expression in normal ovarian tissues [20], which suggested that GPNMB could promote cell proliferation and thus maintain the follicle development.

According to the GO enrichment, the biological process of reproductive process were enriched, which suggested that reproduction was iniated in egg laying stage. The results of KEGG enrichment analysis showed that cytokine-cytokine receptor interaction, cell adhesion molecules (CAMs), ECM-receptor interaction, Phagosome, Neuroactive ligand-receptor interaction, and Calcium signaling pathway were enriched in two or three comparison groups. Differential genes of CXCR4, HGF, PDGFRB, ERBB4 and LHCGR. Elevated CXCR4 have been reported to promote the occurrence of endometriosis [21], suggested its relation to the regulation of ovary function. HGF plays a promoting role in estrogen-mediated stimulation of ovarian cell growth and differentiation [22]. Elevated PDGFRB may promote steroid hormone synthesis in mouse ovaries [23]. High expression of PTGER3 induces ovarian tumorigenesis by upregulating the Ras-MAPK/Erk-ETS1-ELK1/CFTR1 axis [24]. Loss of ERBB4 gene may cause aberrant ovarian function by affecting reproductive function and metabolic function [25]. Elevated expression of LHCGR can promote ovarian ovulation and luteinization [26]. It is suggested that the differentially expressed genes screened in this study may play a key regulatory role in the ovarian development of Yili geese by affecting hormone secretion and ovarian function.

The regulation of follicle and oocyte maturation is a complex multi-factor regulation process that requires specific genes to be accurately expressed at a specific time[3]. Studies have shown that lncRNA has an important regulatory role in reproduction-related processes such as gonadal development, hormone regulation, sexual determination, meiosis, and embryo implantation [27, 28]. In this study, 466 (133 up-regulated and 333 down-regulated lncRNAs), 925 (307 up-regulated and 618 down-regulated lncRNAs), and 742 (347 up-regulated and 395 down-regulated lncRNAs) differential lncRNAs were screened in the KL vs CL, XL vs CL, and XL vs KL groups, respectively. Six lncRNAs were all differentially expressed in the three comparison groups.The target genes of the 6 differentially expressed lncRNAs were significantly enriched in the calcium signaling pathway, apoptosis and ECM-receptor interaction associated with reproduction [29–31], suggested that the six differential lncRNAs affected the ovarian development of Yili goose by regulating the target genes in the above-mentioned signaling pathways.

The GO term with significant enrichment of differential lncRNAs target genes was mainly related to biological regulation, cellular process, metabolic process. The KEGG enrichment analysis showed that apoptosis, calcium signaling pathway, cell adhesion molecules (CAMs), and ECM-receptor interaction pathways were significantly enriched in two or three comparison groups. During reproduction, preventing cell apoptosis and autophagy could obtain high-quality oocytes in the ovary, maintain ovarian function, and protect female fertility [30]. High expression of CTSS in ovary induced excess lipid deposit, oxidative stress and potential ovulation damage [32]. High level of FAS in antral follicle initiated ovarian cell death [33]. RIPK1, FADD, TRADD and TNFRSF1A were all highly expressed in cell apoptosis [34], which would also promote cell apoptosis in the ovary. In this study, 23, 8, 4, 2, 6, 16 differential lncRNAs targeted CTSS, FAS, RIPK1, FADD, TRADD, and TNFRSF1A (Supplementary Table S12-1), respectively, which were significantly enriched in the apoptosis pathway, indicating that the above differential lncRNAs might inhibit ovarian development of Yili goose. It is well known that intracellular Ca²⁺ concentration is a key signaling molecule to control exocytosis and regulate the release of neurotransmitters and endocrine hormones [34]. Wang et al.[35] found that EGFR and GNAS were involved in the reproductive process of pigs. In this study, 4 and 1 differential lncRNAs were found to be significantly enriched in the calcium signaling pathway by targeting EGFR and GNAS, respectively. It is speculated that the above 5 differential lncRNAs participate in the development of the ovary of Yili goose by regulating the release of neurotransmitters and endocrine hormones. In addition, 23, 8, 7 and 24 differential lncRNAs were found in this study to target HGF, PDGFRB, ERBB4 and LHCGR (Supplementary Table S12-2), respectively, which were reported to regulate ovarian function [22, 23, 25, 26].

Studies have found that miRNA has an important function in the development of the mouse ovary, including meiosis of oocytes, reproductive regulation [56]. In this study, 258 (56 up-regulated and 202 down-regulated miRNAs), 1131 (217 up-regulated and 914 down-regulated miRNAs), and 909 (358 up-regulated and 551 down-regulated miRNAs) differential miRNAs were screened in the KL vs CL, XL vs CL, and XL vs KL groups, respectively. There were 44 differential miRNAs overlapped in the three comparison groups, which were involved in ECM-receptor interaction, FoxO signaling pathway, and oxytocin signaling pathway, etc, all related to reproduction regulation [31, 36, 37].

KEGG enrichment analysis showed that miRNAs target genes in the KL vs CL group were significantly enriched in ECM-receptor interaction, in which THBS1 was involved in mediating IGF-I-induced steroid production and proliferation in granulosa cells during follicular development [31, 36, 37]. TNC defective would cause the failure of ovulation [39]. In this study, novel_miR_1844, novel_miR_1333, and novel_miR_336 target THBS1 and TNC, respectively, might regulate the ovulation process of Yili goose.

In this study, oha-miR-375-3p, tgu-miR-375, and gga-miR-375 were significantly different in the KLvs CL group. Compared with the prelaying period, it was increased by 46.8 times during the laying period. Studies have shown that the Wnt signaling pathway was regulated by miR-375 [40], and the Wnt signaling pathway is crucial in the differentiation, proliferation and metastasis of primordial germ cells and the regulation of stromal cells [41]. This indicated that oha-miR-375-3p, tgu-miR-375, and gga-miR-375 might participate in the regulation of the developmental process of the ovaries from the prelaying period to the laying period through the Wnt signaling pathway. In the XLvs CL group, the top 10 foldchanged novel_miR_1116, novel_miR_1473, novel_miR_2688, novel_miR_3346, novel_miR_3441, novel_miR_3838, novel_miR_4002, novel_miR_4129 and novel_miR_4382 target genes participate mainly in Hedgehog signaling pathway, MAPK signaling pathway, and FoxO signaling pathway. Hedgehog signaling pathway was involved in regulating the proliferation and differentiation of germline stem cells [42]. MAPKs mediate signal transmission related to a variety of cell activities in cells, including cell proliferation, differentiation, survival, death and transformation [43, 44]. The main function of FoxO is to regulate cell cycle, apoptosis, atrophy and energy balance [45]. Therefore, these novel miRNAs might play a key role in regulating the proliferation and differentiation of germline stem cells and signal transmission related to cell activities. In the XL vs KL group, gga-miR-34c-5p and ppa-miR-34c were significantly down-regulated with 71.4-fold and 71.2-fold, respectively. Members of the miR-34 family (miR-34a, miR-34b, and miR-34c) have been widely speculated to be key mediators of p53 pathway [46, 47], which regulated a variety of cellular processes, including apoptosis, senescence, cell cycle arrest, differentiation, and DNA repair and replication [48]. It is speculated that gga-miR-34c-5p and ppa-miR-34c might regulate the development, senescence and apoptosis of ovarian cells.

Through the miRNA binding elements, lncRNA could compete with mRNA to combine with miRNA, so as to realize the regulation of lncRNA on mRNA expression [49]. Despite great progress in understanding human disease using the ceRNA, relatively few studies have been carried out in the regulation of ovarian development and reproductive traits.

In this study, 509 ceRNA (lncRNA-miRNA-mRNA) pairs were constructed in the KL vs CL group, including 38 lncRNA, 11 miRNA and 43 mRNA. novel_miR_336 showed the most nodes in the regulatory network, and its expression was down-regulated in the KL vs CL group. The predicted target genes, PDGFRB and LPAR4, were negatively correlated with the expression of novel_miR_336. PDGFRB was involved in controlling the synthesis of steroid hormones in the ovaries of mice [49]. LPAR4 is expressed at high levels in human ovaries [50]. In this network, there were many nodes of MSTRG.129094.34, and its target genes were involved in a number of pathways related to reproduction, such as MAPK signaling pathway and Calcium signaling pathway. At the same time, MSTRG.129094.34 was predicted to target novel_miR_336, and to be negatively regulated by novel_miR_336. It is speculated that MSTRG.129094.34 might up-regulate PDGFRB and LPAR4 by targeting novel_miR_336, thereby participating in the synthesis of steroid hormones in the ovary. Studies have shown that mmu-miR-34b-5p and mmu-miR-107-3p might be involved in the alternative splicing of kitl pre-mRNA in mouse ovarian granulosa cells [51]. WDFY4 was associated with antigen processing, T cell activation, and immune response in humans, mice, and rats [52]. MSTRG.129094.34 targeted gga-miR-34b-5p which further regulated the expression of WDFY4 to promote the immune response in the ovary, and then affecting the ovarian development process.

In this study, 249 ceRNA (lncRNA-miRNA-mRNA) relationship pairs were constructed in the XL vs CL group, including 32 lncRNA, 14 miRNA and 22 mRNA. novel_miR_814 has the most nodes in the regulatory network, and its expression was up-regulated in the XL vs CL group, which targeted to NRG1 and KIT and showed negative correlation with this two genes. The expression of NRG1 was related to LH-induced ovulation and depends on the ERK1/2-C/EBP regulatory pathway. It is supposed that NRG1 could enhance granulosa cell luteinization and regulate the maturation and development of oocytes [53]. KIT was related to the cleavage of vitelline membrane and the formation of primordial follicles and the number of oocytes [54]. LncRNA of MSTRG.128144.1 and MSTRG.158716.1 targeted to NRG1 and KIT, respectively, through trans-action. It is speculated that MSTRG.128144.1 and MSTRG.158716.1 could release the inhibitory factor on the target genes of NRG1 and KIT, respectively, by adsorbing novel_miR_814. The other two miRNAs of gga-miR-145-5p and pccr-miR-133a-3p have multiple nodes in this network. Overexpression of miR-145-5p could inhibit the activation of Notch signaling pathway and apoptosis signaling pathway [55], which regulated cell proliferation, differentiation and apoptosis, and participated in the development of stromal follicles [56]. Both gga-miR-145-5p and MSTRG.5970.28/MSTRG.3524.1 targeted ERBB4. ERBB4 has been demonstrated to be related to ovarian function. Thus, it is speculated that MSTRG.5970.28/MSTRG.3524.1 might be act as the ceRNA sponges of gga-miR-145-5p and affect the expression of ERBB4. MiR-133 has been demonstrated to regulate oocyte meiosis [57] and suppress ovarian cancer cell proliferation. NYAP2 and MSTRG.24644.1 were significantly down-regulated in the XL vs CL group, and both targeted to ccr-miR-133a-3p. It is speculated that MSTRG.24644.1 affects the expression of NYAP2 by acting as the ceRNA sponges of ccr-miR-133a-3p. Moreover, the results of dual-luciferase reporter detection suggested that MSTRG.5970.28 and gga-miR-145-5p have a targeted binding relationship. It is suggested that the lnRNA regulates gene expression by binding miRNA, and then regulates the growth and development of Yili goose ovary.

In this study, 1941 ceRNA (lncRNA-miRNA-mRNA) relationship pairs were constructed in the XL vs KL group, including 55 lncRNA, 33 miRNA and 61 mRNA. novel_miR_2064 has the most nodes in the regulatory network, and was down-regulated in the XL vs KL group. The novel_miR_2064 and MSTRG.7233.4 was predicted to target gene MYL9, which is enriched in the oxytocin signaling pathway. It is speculated that MSTRG.7233.4 affects the expression of MYL9 by acting as a ceRNA of novel_miR_2064. In addition, oan-miR-143-5p and gga-miR-145-5p in this network have been reported to be related to the reproduction process, in which oan-miR-143-5p belongs to the miR-143 family, while miR-143 is related to steroid hormone synthesis, proliferation and apoptosis of granulosa cells [60, 61]. There were 21 lncRNAs targeted to PDGFB and oan-miR-143-5p, 22 lncRNAs targeted to LOC106032597 and gga-miR-145-5p in this network. Therefore, lncRNAs might affect steroid hormone synthesis, granulosa cell proliferation and apoptosis, and then affect the ovarian development of Yili goose.

In this study, differentially expressed mRNAs, lncRNAs and miRNAs were identified in the ovary of Yili goose at different laying stages. It was found that differential mRNAs, lncRNAs and miRNAs were mainly related to cellular processes, biological regulation, ECM-receptor interaction, neuroactive ligand-receptor interaction, and calcium signaling pathway. A differential lncRNA-miRNA-mRNA regulatory network related to cell proliferation, differentiation and apoptosis and involved in stromal follicle development were established and preliminarily validated, which could be regarded as a key regulatory pathway of ovarian development in Yili goose.

Animals, Sample Collection and Histomorphometric Analysis

Yili goose enters the peak egg laying period at the third year. In this study, 4 geese at 3-year-olds at the same batch and similar body weight (3.2 kg ± 0.2 kg) were selected at the period of prelaying (39 months) (KL), laying (41 months) (CL)and ceased period (46 months) (XL), respectively. All the geese were provided by the National Yili goose breeding farm (Xinjiang, China). All selected geese were euthanized by inhaling carbon dioxide and cervical dislocation, which performed by competent personnel who experienced and correctly applied the technique. Then, the ovary removed immediately after slaughter. The morphological characteristics and weight data of the ovary were recorded quickly. Stroma with cortical follicles < 2 mm in diameter were dissected out of the ovaries, rinsed with PBS buffer and fixed in 4% paraformaldehyde. The samples were then embedded in paraffin, sectioned (5 µm) and mounted on slides, and standard hematoxylin and eosin (H&E) staining was performed. Histological characteristics of the ovarian stroma were observed using Nikon Eclipse Ci-L instrument and Nikon E.Z-MET software. For RNA-seq, collected ovaries were stored in liquid nitrogen and then transferred to -80℃.

RNA extraction, RNA-seq library preparation, and sequencing

Total RNA was extracted from each ovarian tissue using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA degradation and contamination were detected by 1.5% agarose gel. The purity and concentration of RNA were determined from OD260/280 readings using the Nanodrop 2000 (Thermo Fisher Scientific Inc, Waltham, MA, USA). RNA integrity was checked using Agilent 2100 Bioanalyzer (Agilent Technologies,Palo Alto, CA, USA).

Approximately 3 µg RNA per sample was used for whole-transcriptome library. Ribosomal RNA was depleted using the Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, USA). After rRNA depletion, the remaining RNA was purified. The NEBNext® Ultra™ The Directional RNA Library Prep Kit for Illumina® (New England Biolabs; NEB, Ipswich, MA. USA) was used according to manufacturer’s guidelines to construct the cDNA libraries. The cDNA fragments were enriched by PCR amplification. The amplification product was purified using AMPure XP system. PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The 12 ovary tissues cDNA libraries were subjected to 2×150 bp paired-end sequencing using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA).

For small RNA sequencing, the same sample was used to construct Illumina small RNA-Seq (RNA sequencing) library by using the NEBNext®MultiplexSmall RNA Library Prep Set for Illumina kit® (NEB) following the manufacturer’s recommendations. In brief, the 5' SR and 3’ SR adaptor for Illumina was ligated to the small RNA, and first strand cDNA was synthesized. The product was then subjected for the second strand synthesis. The product with 3’SR and 5’SR was then PCR amplified. PAGE gel was used to recover the 140-160bp target fragment, rubber cutting recycling as the pieces get small RNA libraries. At last, PCR products were purified (AMPure XP system) and library quality was assessed. The 12 ovary tissues cDNA libraries were subjected to 1×50bp single-end sequencing using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA).

Quality control

Raw data (raw reads) of fastq format were firstly processed through self perl scripts. Clean data(clean reads) were obtained by removing reads containing adapter, ploy-N and low quality reads. Further, the reads of miRNAs were trimmed by removing the sequences smaller than 18 nt or longer than 30 nt. The content of Q20, Q30 and GC and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Analysis of RNA-Seq data

The clean reads were aligned to the geese reference genome (Ans Cyg_PRJNA183603_v1.0) by Hisat2(v2.0.4) [62], the mapped reads of each sample were assembled by StringTie (v1.3.1) [63]. StringTie༈v1.3.1༉[63] was used to assess expression levels of mRNAs and lncRNAs by calculating FPKM [64](Fragments Per Kilobase of transcript per Million fragments mapped). DESeq2 [65] was used for differentially expressed gene (DEG) analysis and Benjamini and Hochberg's [66, 67] method for FDR (False Discovery Rate) correction. Differential expression of mRNAs and lncRNA were defined according to the following criteria: |log2(Fold Change)| ≥1 and adjusted P-value (FDR) < 0.05.

Analysis of microRNA-Seq data

The Clean Reads were aligned with Silva database, GtRNAdb database, Rfam database, and Repbase database, respectively, to filter ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), repeat sequences, and other ncRNA using Bowtie tools (v1.0.0) [68], then unannotated reads containing miRNA were acquired. Unannotated reads were aligned to the geese reference genome (Ans Cyg_PRJNA183603_v1.0) using Bowtie(v1.0.0)[68], to obtain the location information on the reference genome (Mapped Reads). The small RNA that matched to the reference genome was compared to the miRNAs from all animals in the miRBase(v22) database to identify known miRNA. Novel miRNA prediction was performed using miRDeep2 (v2.0.5) [69]. For both known miRNAs and novel miRNAs, the miRNA expression level was calculated and normalized by transcripts per million (TPM) [70]. DESeq2(v1.10.1) [65] was used for differentially expressed miRNA analysis and Benjamini and Hochberg's [66, 67] method for FDR correction. Differential expression of miRNAs were defined according to the following criteria: |log2(Fold Change)|≥0.58 and P-value<0.01.

lncRNAs identification

Transcript class_codes with "i", "x", "u", "o", "e", length ≥ 200bp and Exon number ≥ 2, were selected, FPKM ≥ 0.1, CPC(Coding Potential Calculator) [71], CNCI(coding-noncoding-index) [72], CPAT (coding potential assessment tool) [73] and Pfam-scan [74] were used to distinguish mRNAs from lncRNAs.

GO and KEGG enrichment analyses

Coding genes located within 100 kb upstream and downstream of lncRNA were regarded as the cis-target genes of the lncRNA. Pearson correlation coefficient method was used to analyze the correlation between lncRNA and mRNA and the genes with |r|> 0.95 and P < 0.01 were considered as trans-target gene of the lncRNA. According to the gene sequence information of known miRNAs and novel miRNAs, miRanda[75] and targetscan [76] software were used to predict miRNA target genes. ClusterProfiler(v3.10.1) [77] and KOBAS (v2.0) [78] was used for GO(Gene Ontology) and KEGG(Kyoko Encyclopedia of Genes and Genomes) pathway analyses of differentially expressed mRNA, lncRNA target genes and miRNA target genes. GO enrichment analysis includes biological process (BP), cellular component (CC), and molecular function (MF). GO terms or KEGG pathways with corrected p-value (q-value) < 0.05 were considered to be significantly enriched.

Construction of lncRNA-miRNA-mRNA networks

In this study, lncRNA-miRNA-mRNA network was constructed according to the ceRNA method with the following step, (1) The correlation of differentially expressed miRNA and mRNA was used to predict the negative interaction of miRNA-mRNA. (2) The predicted miRNA binding sites for the differentially expressed lncRNAs were identified using the miRanda [75], to construct a differentially expressed lncRNA-miRNA negative regulatory network. (3) Baes on the predicted miRNA-mRNA and lncRNA-miRNA regulatory network, the lncRNA-miRNA-mRNA network (ceRNA network) was constructed using Cytoscape(v3.8.0) [79] software.

Reverse transcription real-time quantitative PCR (RT-qPCR)

The remaining RNA from RNA-seq library was used for RT-PCR quantification. Based on the lncRNA–miRNA–mRNA correlation networks, specifically, several interaction nodes were validated by RT-qPCR, including 9 mRNA, 6 lncRNA and 7 miRNA(Table S1). For mRNA, lncRNA and miRNA, RNA was reverse transcribed using RevertAid Reverse Transcriptase (Thermo Fisher Scientific) following the manufacturer’s protocol. Quantitative RT-PCR (RT-qPCR) was performed using AceQ Universal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China) and ABI StepOnePlus machine (ABI, Foster City, CA, USA). Primers were designed using Primer Premier 5.0. All primer sequences, including selected genes, miRNAs and internal control genes (GAPDH and U6 snRNA), were listed in Table S1. Relative expression levels of genes and miRNAs were calculated by the2^−ΔΔCt [80]method.

Dual-Luciferase reporter assays

PsiCheck2-MSTRG.5970.28-wt/mut vectors were synthesized by He Fei Yuanen Biotechnology Company (Hefei, Anhui, China). 293T cells were seeded in a 24-well plates (5×105 cells/well). Cells were transfected with psiCheck2-MSTRG.5970.28-wt/mut, gga-miR-145-5p mimic and NC mimic using Lip 3000 (Invitrogen). At 48 h after transfection, cell lysates were prepared and dual luciferase reporter assay kit (RG027, Beyotime Institute of Biotechnology) was used to measure luciferase activities following the manufacturers instructions. The relative luciferase activities were calculated by comparing the Firefly/Renilla luciferase activity ratio.

Statistical analysis

Statistical analysis was performed using the software IBM SPSS Statistics version 22. GraphPad Prism 8 was applied for making graph. It was considered to be statistically significant when p-value < 0.05.

KL: Prelaying period; CL: Laying period; XL: ceased period; LncRNA: Long non-coding RNA; miRNA: MicroRNA; DEG: Differentially expressed gene; PDGFRB: Platelet derived growth factor receptor beta; PDGFRB: Platelet derived growth factor receptor beta; ERBB4: Erb-b2 receptor tyrosine kinase 4; LHCGR: Luteinizing hormone/choriogonadotropin receptor; RNA-seq: RNA sequencing; rRNA: ribosomal RNA; tRNA: Transfer RNA; snoRNA: Small nucleolar RNA; FPKM: Fragments per kilobase of transcript per million fragments mapped; FDR: False discovery rate; GO: Gene Ontology; KEGG: Kyoko Encyclopedia of Genes and Genomes; BP: Biological process; CC: Bellular component; MF: Molecular function; RT-qPCR: Reverse transcription real-time quantitative PCR; PGR: Progesterone receptor; FDFT1: Farnesyl-diphosphate farnesyltransferase 1; ADRA2A: Adrenoceptor alpha 2A; CP: Ceruloplasmin; GPNMB: Glycoprotein nmb; CXCR4: C-X-C motif chemokine receptor 4; HGF: Hepatocyte growth factor; PTGER3: Prostaglandin E receptor 3; CTSS: Cathepsin S; FAS: Fas cell surface death receptor; RIPK1: Receptor interacting serine/threonine kinase 1; FADD: Fas associated via death domain; TRADD: TNFRSF1A associated via death domain; TNFRSF1A: TNF receptor superfamily member 1A; EGFR: Epidermal growth factor receptor; GNAS: GNAS complex locus; THBS1: Thrombospondin 1; TNC: Tenascin C; LPAR4: Lysophosphatidic acid receptor 4; WDFY4: WDFY family member 4; NRG1: Neuregulin 1; KIT: KIT proto-oncogene: receptor tyrosine kinase; NYAP2: Neuronal tyrosine-phosphorylated phosphoinositide-3-kinase adaptor 2; MYL9: Myosin light chain 9.

Ethics approval and consent to participate

Our study was carried out in compliance with the ARRIVE guidelines (AVMA Guidelines for the Euthanasia of Animals: 2020 Edition). All experimental animal procedures complied with the laboratory animal management and welfare regulations approved by the Ethics Committee of Xinjiang Agricultural University, Urumqi, Xinjiang, China. This study was approved by the Ethics Committee of Xinjiang Agricultural University, Urumqi, Xinjiang, China (Permit No. 2019004, date of approval: 22 March, 2019).

Consent for publication

Not applicable.

Availability of data and materials

The raw sequences were deposited into Sequence Read Archive (SRA) database with the BioProject accession number PRJNA825140. (https://www.ncbi.nlm.nih.gov/bioproject/825140).

Competing interests

The authors declare that they have no competing interests.

Funding

This research was funded by the National Natural Science Foundation of China (31960648) and the the Xinjiang Uygur Autonomous Region Graduate Research and Innovation Project (XJ2021G177).

Authors’ contributions

XYZ and HYL designed the project. XYZ, YPW, YC and ZYM managed the experimental geese. XYZ, YPW, JHL and YYY collected the samples. XYZ and YPW performed bioinformatic analysis. XYZ, YPW, YC and ZYM participated in data analysis and experimental validation. HYL supervised this study. The manuscript was prepared by XYZ and substantially revised by HYL. All authors read and approved the final manuscript for publication.All authors read and approved the final manuscript for publication.

Acknowledgments

We are grateful to Prof. Xingyong Chen for suggestions to design the paper and for writing assistance.

1. Wang Y, Chen Q, Liu Z, Guo X, Du Y, Yuan Z, Guo M, Kang L, Sun Y, Jiang Y: Transcriptome Analysis on Single Small Yellow Follicles Reveals That Wnt4 Is Involved in Chicken Follicle Selection. FRONT ENDOCRINOL 2017, 8:317.

2. Ren J, Sun C, Chen L, Hu J, Huang X, Liu X, Lu L: Exploring differentially expressed key genes related to development of follicle by RNA-seq in Peking ducks (Anas Platyrhynchos). PLOS ONE 2019, 14(6):e209061.

3. Bonnet A, Dalbies Tran R, Sirard MA: Opportunities and challenges in applying genomics to the study of oogenesis and folliculogenesis in farm animals. Reproduction (Cambridge, England) 2008, 135(2):119-128.

4. Xu Q, Zhao W, Chen Y, Tong Y, Rong G, Huang Z, Zhang Y, Chang G, Wu X, Chen G: Transcriptome profiling of the goose (Anser cygnoides) ovaries identify laying and broodiness phenotypes. PLOS ONE 2013, 8(2):e55496.

5. Wu Y, Zhao X, Chen L, Wang J, Duan Y, Li H, Lu L: Transcriptomic Analyses of the Hypothalamic-Pituitary-Gonadal Axis Identify Candidate Genes Related to Egg Production in Xinjiang Yili Geese. Animals : an open access journal from MDPI 2020, 10(1):90.

6. Zhao Z, Wang L, Sun X, Zhang J, Zhao Y, Na RS, Zhang J, Leung FCC: Transcriptome analysis of the Capra hircus ovary. PLOS ONE 2015, 10(3):e121586.

7. Zhang Q, Sun H, Jiang Y, Ding L, Wu S, Fang T, Yan G, Hu Y: MicroRNA-181a suppresses mouse granulosa cell proliferation by targeting activin receptor IIA. PLOS ONE 2013, 8(3):e59667.

8. Xu Q, Zhang Y, Chen Y, Tong Y, Rong G, Huang Z, Zhao R, Zhao W, Wu X, Chang G et al: Identification and differential expression of microRNAs in ovaries of laying and Broody geese (Anser cygnoides) by Solexa sequencing. PLOS ONE 2017, 9(2).

9. Rui C, Jun WW, Long ZX, Peng X, Yi W, Lin LH: Expression and preliminary functional profiling of the let-7 family during porcine ovary follicle atresia. MOL CELLS 2015, 38(4).

10. Jianfang Z, Xiaowen J, Doudou Z, Yueqin L, Junkai L, Jiali L, Haoshu L, Sheng C: miR-143 is critical for the formation of primordial follicles in mice. Frontiers in bioscience (Landmark edition) 2013, 18.

11. C BS, A SC, N RK, J DT, H SA, V TML: Manipulation of estrogen synthesis alters MIR202* expression in embryonic chicken gonads. BIOL REPROD 2011, 85(1).

12. V SA, Dmitriy O, Roland G, Marcela L, Milos M: Identification of microRNAs controlling human ovarian cell steroidogenesis via a genome-scale screen. J CELL PHYSIOL 2009, 219(2).

13. Liu G, Liu S, Xing G, Wang F: lncRNA PVT1/MicroRNA-17-5p/PTEN Axis Regulates Secretion of E2 and P4, Proliferation, and Apoptosis of Ovarian Granulosa Cells in PCOS. Molecular Therapy - Nucleic Acids 2020, 20:205-216.

14. Yang G, Lu X, Yuan L: LncRNA: a link between RNA and cancer. Biochimica et biophysica acta 2014, 1839(11):1097-1109.

15. Shen X, Bai X, Xu J, Zhou M, Xu H, Nie Q, Lu X, Zhang X: Transcriptome sequencing reveals genetic mechanisms underlying the transition between the laying and brooding phases and gene expression changes associated with divergent reproductive phenotypes in chickens. MOL BIOL REP 2016, 43(9):977-989.

16. Zhang T, Zhang X, Han K, Zhang G, Wang J, Xie K, Xue Q, Fan X: Analysis of long noncoding RNA and mRNA using RNA sequencing during the differentiation of intramuscular preadipocytes in chicken. PLOS ONE 2017, 12(2):e172389.

17. Bushati N, Cohen SM: microRNA functions. Annu.rev.cell Dev.biol 2007, 23(1):175.

18. ALVAREZGARCIA, Ines, MISKA, Eric A: MicroRNA functions in animal development and human disease. DEVELOPMENT 2005, 132(21):4653-4662.

19. Ouyang Q, Hu S, Li L, Ran M, Zhu J, Zhao Y, Hu B, Hu J, He H, Li L et al: Integrated mRNA and miRNA transcriptome analysis provides novel insights into the molecular mechanisms underlying goose pituitary development during the embryo-to-hatchling transition. POULTRY SCI 2021, 100(9):101380.

20. Ma, R.Q.; Ye X.; Cheng H.Y.; Chen J.; Cui H.; Wei L.H.; Chang X.H.. The expression and cliniacl significance of GPNMB in epithelial ovarian cancer. Chinese Journal of Clinical Obstetrics and Gynecology. 2013, 14, 388-391.

21. Moridi I, Mamillapalli R, Cosar E, Ersoy GS, Taylor HS: Bone Marrow Stem Cell Chemotactic Activity Is Induced by Elevated CXCl12 in Endometriosis. REPROD SCI. 2017, 24(4):526-533.

22. Zhou HY, Pon YL, Wong AST: HGF/MET Signaling in Ovarian Cancer. CURR MOL MED 2008, 8(6):469-480.

23. Schmahl J, Rizzolo K, Soriano P: The PDGF signaling pathway controls multiple steroid-producing lineages. Cold Spring Harbor Laboratory Press 2008, 22(23):3255-3267.

24. Rodriguez-Aguayo C, Bayraktar E, Ivan C, Aslan B, Mai J, He G, Mangala LS, Jiang D, Nagaraja AS, Ozpolat B et al: PTGER3 induces ovary tumorigenesis and confers resistance to cisplatin therapy through up-regulation Ras-MAPK/Erk-ETS1-ELK1/CFTR1 axis. EBIOMEDICINE 2019, 40:290-304.

25. Ville V, Nsrein A, Milena D, Antti K, Ilkka M, Claes O, Matti P, Nafis R, Klaus E, J VS et al: Erbb4 regulates the oocyte microenvironment during folliculogenesis. HUM MOL GENET 2020, 29(17):2813-2830.

26. Menon B, Gulappa T, Menon KMJ: Molecular regulation of LHCGR expression by miR-122 during follicle growth in the rat ovary. MOL CELL ENDOCRINOL 2017, 442:81-89.

27. Xu F, Fengzhe L, Feng W, Guomin Z, Jing P, Caifang R, Tingting Z, Hua Y, Ziyu W, Yanli Z: Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. BIOSCIENCE REP 2018, 38(2):R20171350.

28. Hua Y, Feng W, Fengzhe L, Caifang R, Jing P, Yongjie W, Ziyu W, Xu F, Yanli Z: Comprehensive Analysis of long non-coding RNA and mRNA Expression Patterns in Sheep Testicular Maturation. BIOL REPROD 2018, 99(3):650-661.

29. J MTF: Tuning exocytosis for speed: fast and slow modes. Biochimica et biophysica acta 2003, 1641(2-3):157-165.

30. Yadav PK, Tiwari M, Gupta A, Sharma A, Prasad S, Pandey AN, Chaube SK: Germ cell depletion from mammalian ovary: possible involvement of apoptosis and autophagy. J BIOMED SCI 2018, 25(1):36.

31. Wanqiu Z, Taoyan Y, Yan F, Dong N, Weihu C, Li C, Lizhi L: Seasonal differences in the transcriptome profile of the Zhedong white goose (Anser Cygnoides) pituitary gland. POULTRY SCI 2020, 100(2):1154-1166.

32. Wei Z, Li P, Huang S, Lkhagvagarav P, Zhu M, Liang C, Jia C: Identification of key genes and molecular mechanisms associated with low egg production of broiler breeder hens in ad libitum. BioMed Central 2019, 20(1):408.

33. Fernanda P, Griselda I, Alejandra V, Olga G, Adalí P, Marta T: Gonadotropin-releasing hormone antagonist antide inhibits apoptosis of preovulatory follicle cells in rat ovary. BIOL REPROD 2005, 72(3):659-666.

34. J MM, Emilia S, Leo C, Marion M, W JA, K PK, J SD: Death Processes in Bovine Theca and Granulosa Cells Modelled and Analysed Using a Systems Biology Approach. INT J MOL SCI 2021, 22(9):4888.

35. Wang M, Feng S, Ma G, Miao Y, Zuo B, Ruan J, Zhao S, Wang H, Du X, Liu X: Whole-Genome Methylation Analysis Reveals Epigenetic Variation in Cloned and Donor Pigs. FRONT GENET 2020, 11:23.

36. Manman S, Tingting L, Fuxiang C, Pengfeng W, Ying W, Lan C, Kaizhou X, Jinyu W, Genxi Z: Transcriptomic Analysis of circRNAs and mRNAs Reveals a Complex Regulatory Network That Participate in Follicular Development in Chickens. FRONT GENET 2020, 11:503.

37. Gao Y, Wu F, Ren Y, Zhou Z, Chen N, Huang Y, Lei C, Chen H, Dang R: MiRNAs Expression Profiling of Bovine (Bos taurus) Testes and Effect of bta-miR-146b on Proliferation and Apoptosis in Bovine Male Germline Stem Cells. INT J MOL SCI 2020, 21(11):3846.

38. Grado-Ahuir JA, Aad PY, Ranzenigo G, Caloni F, Cremonesi F, Spicer LJ: Microarray analysis of insulin-like growth factor-I-induced changes in messenger ribonucleic acid expression in cultured porcine granulosa cells: possible role of insulin-like growth factor-I in angiogenesis. J ANIM SCI 2009, 87(6):1921-1933.

39. Takashi O, Kanako O, Shoichiro O: Troponin I controls ovulatory contraction of non-striated actomyosin networks in the C. elegans somatic gonad. J CELL SCI 2010, 123(9):1557-1566.

40. Hu F, Sun B, Xu P, Zhu Y, Meng X, Teng G, Xiao Z: MiR-218 Induces Neuronal Differentiation of ASCs in a Temporally Sequential Manner with Fibroblast Growth Factor by Regulation of the Wnt Signaling Pathway. SCI REP-UK 2017, 7(1):804-810.

41. Dong W, Tan F, Yang W: Wnt signaling in testis development: Unnecessary or essential? GENE 2015, 565(2):155-165.

42. Li S, Wang M, Chen Y, Wang W, Wu J, Yu C, Zheng Y, Pan Z: Role of the Hedgehog Signaling Pathway in Regulating the Behavior of Germline Stem Cells. STEM CELLS INT 2017, 2017:5714608.

43. Satoru T, Takuya Y, Yoshiki T, Eisuke N: ERK MAP kinase in G cell cycle progression and cancer. CANCER SCI 2006, 97(8):697-702.

44. McCubrey JA, Franklin MML, Franklin RA: Reactive oxygen species-induced activation of the MAP kinase signaling pathways. ANTIOXID REDOX SIGN 2006, 8(9-10):1775-1789.

45. Guri T, Melissa D, Gopalakrishnan R: FoxO transcription factors; Regulation by AKT and 14-3-3 proteins. Biochimica et biophysica acta 2011, 1813(11):1938-1945.

46. Chang T, Wentzel EA, Kent OA, Ramachandran K, Mullendore M, Lee KH, Feldmann G, Yamakuchi M, Ferlito M, Lowenstein CJ et al: Transactivation of miR-34a by p53 Broadly Influences Gene Expression and Promotes Apoptosis. MOL CELL 2007, 26(5):745-752.

47. Bommer GT, Gerin I, Feng Y, Kaczorowski AJ, Kuick R, Love RE, Zhai Y, Giordano TJ, Qin ZS, Moore BB et al: p53-Mediated Activation of miRNA34 Candidate Tumor-Suppressor Genes. CURR BIOL 2007, 17(15):1298-1307.

48. Lazo PA: Reverting p53 activation after recovery of cellular stress to resume with cell cycle progression. CELL SIGNAL 2017, 33:49-58.

49. Chiu H, Martínez MR, Bansal M, Subramanian A, Golub TR, Yang X, Sumazin P, Califano A: High-throughput validation of ceRNA regulatory networks. BMC GENOMICS 2017, 18(1):418.

50. Ewa L, Pierrette R, Emmanuelle B, Olivier D, Philippe R, Gilles C: Lysophosphatidic acid signaling during embryo development in sheep: involvement in prostaglandin synthesis. ENDOCRINOLOGY 2009, 150(1):422-434.

51. Wang W, Chen J, Luo L, Li Y, Liu J, Zhang W: Effect of cadmium on kitl pre‐mRNA alternative splicing in murine ovarian granulosa cells and its associated regulation by miRNAs. J APPL TOXICOL 2017, 38(2):227-239.

52. Li H, Rukina D, David Fabrice P A, Li Terytty Y, Oh C, Gao Arwen W, Katsyuba E, Bou Sleiman M, Komljenovic A, Huang Q et al: Identifying gene function and module connections by the integration of multispecies expression compendia. GENOME RES 2019, 29(12):2034-2045.

53. Noritaka N, Ikko K, Heng-Yu F, Youko F, Tomoko K, Yoshinori T, Toshihiro M, S RJ, Masayuki S: LH-induced neuregulin 1 (NRG1) type III transcripts control granulosa cell differentiation and oocyte maturation. Molecular endocrinology (Baltimore, Md.) 2011, 25(1):104-116.

54. Jones RL, Pepling ME: KIT signaling regulates primordial follicle formation in the neonatal mouse ovary. DEV BIOL 2013, 382(1):186-197.

55. Wei B, Liu Y, Guan H: MicroRNA-145-5p attenuates high glucose-induced apoptosis by targeting the Notch signaling pathway in podocytes. EXP THER MED 2020, 19(3):1915-1924.

56. Jing J, Jiang X, Chen J, Yao X, Zhao M, Li P, Pan Y, Ren Y, Liu W, Lyu L: Notch signaling pathway promotes the development of ovine ovarian follicular granulosa cells. ANIM REPROD SCI 2017, 181:69-78.

57. Song Y, Shi L, Liu Z, Qiu G: Global analysis of the ovarian microRNA transcriptome: implication for miR-2 and miR-133 regulation of oocyte meiosis in the Chinese mitten crab, Eriocheir sinensis (Crustacea:Decapoda). BioMed Central 2014, 15(1):547.

58. Luo J, Zhou J, Cheng Q, Zhou C, Ding Z: Role of microRNA-133a in epithelial ovarian cancer pathogenesis and progression. ONCOL LETT 2014, 7(4):1043-1048.

59. Guo J, Xia B, Meng F, Lou G: miR-133a suppresses ovarian cancer cell proliferation by directly targeting insulin-like growth factor 1 receptor. Tumor biology 2014, 35(2):1557-1564.

60. Zhang Z, Chen C, Xu M, Zhang L, Liu J, Gao Y, Jiang H, Yuan B, Zhang J: MiR-31 and miR-143 affect steroid hormone synthesis and inhibit cell apoptosis in bovine granulosa cells through FSHR. THERIOGENOLOGY 2019, 123:45-53.

61. Li Z, XiaoXin Z, Xuejing Z, Yu L, Lei L, Sheng C: MiRNA-143 mediates the proliferative signaling pathway of FSH and regulates estradiol production. The Journal of endocrinology 2017, 234(1):1-14.

62. Kim D, Langmead B, Salzberg SL: HISAT: a fast spliced aligner with low memory requirements. NAT METHODS 2015, 12(4):357-360.

63. Mihaela P, M PG, M AC, Tsung-Cheng C, T MJ, L SS: StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. NAT BIOTECHNOL 2015, 33(3):290-295.

64. Trapnell C, Williams BA, Pertea G, Mortazavi A, 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.

65. Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. GENOME BIOL 2014, 15(12):550.

66. Michael L, Robert T: Transcription regulation and animal diversity. NATURE 2003, 424(6945):147-151.

67. Li C, Guo H, Xu X, Weinberg W, Deng C: Fibroblast growth factor receptor 2 (Fgfr2) plays an important role in eyelid and skin formation and patterning. DEV DYNAM 2001, 222(3):471-483.

68. Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. GENOME BIOL 2009, 10(3):R25.

69. R FM, D MS, Li N, Chen W, Nikolaus R: miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. NUCLEIC ACIDS RES 2012, 40(1):37-52.

70. Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN: RNA-Seq gene expression estimation with read mapping uncertainty. BIOINFORMATICS 2010, 26(4):493-500.

71. Kong L, Zhang Y, Ye Z, Liu X, Zhao S, Wei L, Gao G: CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. NUCLEIC ACIDS RES 2007, 35(suppl 2):w345-w349.

72. 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.

73. Wang L, Park Hyun J, Dasari S, Wang S, Kocher J, Li W: CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. NUCLEIC ACIDS RES 2013, 41(6):e74.

74. D FR, Alex B, Jody C, Penelope C, Y ER, R ES, Andreas H, Kirstie H, Liisa H, Jaina M et al: Pfam: the protein families database. NUCLEIC ACIDS RES 2014, 42(Database issue):D222-D230.

75. Doron B, Manda W, Aaron G, S MD, Chris S: The microRNA.org resource: targets and expression. NUCLEIC ACIDS RES 2008, 36(suppl 1):D149-D153.

76. Lewis BP, Shih I, Jones-Rhoades MW, Bartel DP, Burge CB: Prediction of Mammalian MicroRNA Targets. CELL 2003, 115(7):787-798.

77. Yu G, Wang L, Han Y, He Q: clusterProfiler: an R package for comparing biological themes among gene clusters. Omics : a journal of integrative biology 2012, 16(5):284-287.

78. Chen X, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li C, Wei L: KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. NUCLEIC ACIDS RES 2011, 39(suppl 2):W316-W322.

79. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T: Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Cold Spring Harbor Laboratory Press 2003, 13(11):2498-2504.

80. J LK, D ST: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif.) 2001, 25(4):402-408.

No competing interests reported.

Supplementarymaterials.zip

Download PDF

Editorial decision: Major revision
03 Jun, 2022
Reviews received at journal
02 Jun, 2022
Reviewers agreed at journal
25 May, 2022
Reviewers agreed at journal
21 Apr, 2022
Reviews received at journal
18 Apr, 2022
Reviewers agreed at journal
15 Apr, 2022
Reviewers invited by journal
14 Apr, 2022
Editor assigned by journal
14 Apr, 2022
Editor invited by journal
14 Apr, 2022
Submission checks completed at journal
14 Apr, 2022
First submitted to journal
05 Apr, 2022

You are reading this latest preprint version

Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and miRNA of Yili goose ovary at different egg-laying stages

Status:

Version 1

Abstract

Figures

Background

Results

Ovary histomorphometric analysis at different egg laying periods in Yili goose

Overview of sequencing and identification of lncRNA and miRNA

Analysis of the differentially expressed mRNA, lncRNA and miRNA among different egg laying periods

Functional enrichment analysis of differentially expressed genes and target genes

Construction of lncRNA-miRNA-mRNA networks and functional analysis

Validation of RNA-Seq Results

MSTRG.5970.28 Binds to gga-miR-145-5p

Discussion

Conclusions

Materials And Methods

Animals, Sample Collection and Histomorphometric Analysis

RNA extraction, RNA-seq library preparation, and sequencing

Quality control

Analysis of RNA-Seq data

Analysis of microRNA-Seq data

lncRNAs identification

GO and KEGG enrichment analyses

Construction of lncRNA-miRNA-mRNA networks

Reverse transcription real-time quantitative PCR (RT-qPCR)

Dual-Luciferase reporter assays

Statistical analysis

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1