High-throughput small RNA and transcriptome sequencing reveal reproduction-related microRNAs and mRNA in ovary of Geese in Counter-Season Production

Background Broodiness is a phenomenon that occurs in most avian species and significantly reduces productivity. Several genes are known to play an important role in regulating the progress of reproduction, but the molecular regulation mechanism of broodiness remains unclear. In the current study, via high throughput sequencing, we identified and explored the differentially expressed miRNAs and mRNAs involved in ovarian atrophy. We identified a total of 901 mRNAs and 50 miRNAs that were differentially expressed in egg-laying and atrophic ovaries. Among them, numerous DEGs transcripts and target genes for miRNAs were significantly enriched in reproductive processes, cell proliferation, and apoptosis pathways. In addition, a miRNA- gene-pathway network was constructed by considering target relationships and correlation of the expression levels between ovary development-related genes and miRNAs. Conclusions We discovered mRNA and miRNAs transcripts that are candidate regulators of ovary development in broody geese. Our findings expanded our understanding of the functional of miRNAs in ovarian atrophy and demonstrated that RNA-Seq is a powerful tool for examining the molecular mechanism in regulating broodiness. of counter-season egg-laying and brooding geese by high-throughput technologies.


Abstract Background
Broodiness is a phenomenon that occurs in most avian species and significantly reduces productivity.
Several genes are known to play an important role in regulating the progress of reproduction, but the molecular regulation mechanism of broodiness remains unclear. In the current study, via high throughput sequencing, we identified and explored the differentially expressed miRNAs and mRNAs involved in ovarian atrophy.

Results
We identified a total of 901 mRNAs and 50 miRNAs that were differentially expressed in egg-laying and atrophic ovaries. Among them, numerous DEGs transcripts and target genes for miRNAs were significantly enriched in reproductive processes, cell proliferation, and apoptosis pathways. In addition, a miRNA-gene-pathway network was constructed by considering target relationships and correlation of the expression levels between ovary development-related genes and miRNAs.

Conclusions
We discovered mRNA and miRNAs transcripts that are candidate regulators of ovary development in broody geese. Our findings expanded our understanding of the functional of miRNAs in ovarian atrophy and demonstrated that RNA-Seq is a powerful tool for examining the molecular mechanism in regulating broodiness.

Background
Goose is a seasonally bred poultry, which laying eggs in a specific season and month, this feature restricts the development of goose industry, resulting in uneven production throughout the year [1].
The laying performance of goose mainly depends on the growth and development of follicles. The follicles establish a strict hierarchical system through mechanisms such as development, selection and atresia. The more follicles that are selected for grade development, the longer laying period and the higher laying performance of egg-laying poultry [2][3]. In poultry, there is abundant evidence indicating that follicular development and grading is a complex biological process regulated by many factors, which depends on the regulation of hypothalamus-pituitary-ovary axis [4]. HGC axis regulates ovarian follicle growth and ovulation for egg production by coordinating several reproductive-related endocrine hormones and growth factors [5][6].
In the process of broodiness, decreased gonadotrophin-releasing hormone (GnRH) and increased vasoactive intestinal polypeptide (VIP) release from the hypothalamus induced production of prolactin (PRL) [7][8]. GnRH synthesized and secreted by the hypothalamus and stimulate the synthesis and secretion of gonadotropins by binding with its receptor GnRHR [9]. During the reproductive cycle of birds, PRL cause atresia of lower-grade follicles and broody, VIP directly acts on the anterior pituitary to stimulate PRL secretion [10]. In addition to GnRH, VIP and PRL, several genes such as luteinizing hormone (LH) [11], dopamine (DA) [12] and neuropeptide Y (NPY) [13], have been chosen as candidates to analyze their genetic effect on egg production and broodiness traits. However, the molecular mechanisms and associated signaling pathways of these candidate genes remain poorly understood.
MicroRNAs (miRNAs) are endogenous non-coding small molecule single-stranded nucleic acids with a length of 20-24 bases that function in post-transcriptional regulation of gene expression by recognizing and combining the 3'UTR of the target genes [14][15]. MiRNAs are involved in many cellular processes including development, immunity, cell cycle control, metabolism and disease [16].
Sirotkin et al have demonstrated that miRNAs can regulate proliferation and apoptosis of granulosa cells by controlling steroidogenesis in the human ovary [17]. And there have been some research focused on the regulation of miRNAs in animals' reproduction, such as mouse [18], sheep [19] and fish [20].Thus, we speculate that miRNAs may play important roles in the process of broodiness of poultry.
Xingguo grey goose is a famous native breed distributing in Jiangxi province of China and has a strong seasonal reproduction characteristic. The objective of this experiment was thus to identify mRNA and miRNAs differentially expressed in the ovary of counter-season egg-laying and brooding geese by high-throughput technologies.

Results
Overview of RNA sequencing After performing transcriptome sequencing quality control, we obtained a total of 110,176,986 and 110,458,247 raw reads and 106,311,727 and 107,125,914 clean reads in egg-laying and atrophic ovaries, respectively (Table 1). Additionally, we identified 33,089 protein-coding transcripts including 25,975 known and 7,114 novel protein-coding transcripts. A boxplot of FPKM was shown in Figure   S1A, the expression levels of transcripts in R and NR groups were similar. Furthermore, 1,019,077 and 1,002,561 SNPs, and 23,058 and22, 003 alternative splicing events were found in egg-laying and atrophic ovaries, respectively. In addition, SE (skipped exon) accounts for the majority of these alternative splicing events types and RI (retained intron) is the least ( Figure S2).

Analysis of differentially expressed mRNAs
Using gene expression profiling and comparing the egg-laying ovaries group with the atrophic ovaries group, a total of 901 DEGs (|log2 FC| ≥ 1.3, with p < 0.05) were screened, of which 492 were upregulated and 409 were down-regulated ( Figure 1A). Genes showing differential expression were subjected to hierarchical clustering analysis (Figure 2A). A total of 2,212 DEGs were annotated and included in three main GO categories i.e., biological process, cellular component and molecular function. ( Figure 3A). Among these, the most important molecular functions involved androgen receptor activity, calcium-dependent phospholipid binding, DNA N-glycosylase activity, protein disulfide oxidoreductase activity and acid-ammonia (or amide) ligase activity. The biological processes consisted of reproductive structure development, reproductive system development, gonad development, androgen receptor signaling pathway and aromatic compound catabolic process.
However, cellular process was not significant and frequent enrichment term. In the KEGG pathway analysis, a total of 541 DEGs were assigned to 121 KEGG pathways. Further analysis showed that their main pathways of enrichment included Wnt signaling pathway, Toll−like receptor signaling pathway, TGF−beta signaling pathway, Protein processing in endoplasmic reticulum, Phagosome, Lysosome, Endocytosis, Apoptosis and Adipocytokine signaling pathway ( Figure 3B).
Overview of small RNA sequencing A total of 12,813,559 and 13,457,346 raw reads and 12,650,516 and 13,250,406 high-quality (clean) reads were obtained from small RNA libraries of egg-laying and atrophic ovaries, respectively ( Table  2). After alignment with small RNAs in GenBank, Rfam and the reference genome, we identified 64.32% mature miRNAs, and the remaining small RNAs (35.68%) included rRNA, snRNA, snoRNA, tRNA, simple repeats, exon sense/antisense and intron sense/antisense ( Figure S3). We identified a total of 321 miRNAs including 272 known and 49 novel miRNAs in six small RNA libraries, and the length of most miRNAs were 18~25 nucleotides ( Figure S4). Total 281 and 299 miRNAs were detected from R and NR, respectively. Additionally, a boxplot of TPM was shown in Figure S1B, the expression levels of miRNAs in R and NR groups were similar.

Analysis of differentially expressed miRNAs
There were 50 significantly differentially expressed miRNAs (DEMs) including 31 that were upregulated and 19 that were down-regulated in egg-laying and atrophic ovaries (Fig. 1B). Moreover, the miRNAs with differential expression were clustered by hierarchical clustering analysis ( Figure 2B).
A total of 7,204 target protein-coding genes for DEMs were identified. Among them, 2396 genes had transcripts identified as DEGs and also significantly negatively correlated with miRNA expression, and were assigned as intersection genes, which were more likely to be predicted miRNA target genes [21].
Additionally, we performed GO and pathway enrichment analysis of intersection genes to identify biological functions of miRNAs. The results were shown in Figure 4A, the most important molecular functions involved organic substance transport, single-organism localization, protein targeting to nucleus and nuclear transport. The biological processes consisted of small molecule binding, phosphorelay sensor kinase act, GTPase activity, cysteine-type peptidase activity, catalytic activity and hydrolase activity. In cellular component, the important GO terms were viral capsid, Membrane, actin cytoskeleton, virion, viral envelope and virion part. The KEGG results revealed that the enriched pathways involved Endocytosis, Lysosome, Protein processing in endoplasmic reticulum, Cysteine and methionine metabolism, RNA degradation and Toll−like receptor signaling pathway ( Figure 4B).

Construction of miRNA-gene-pathway relationship network
We analyzed several important intersection genes and their corresponding DEMs. Among them, 18 intersection genes involved in reproductive endocrine system, hormonal regulation and proliferation and apoptosis corresponded with 29 DEMs were shown in Table 3. Phospholipase C beta 1 (PLCB1), TNF receptor superfamily member 1A (TNFRSF1A), and protein phosphatase 3 catalytic subunit beta (PPP3CB) genes play roles in multiple pathways. To further understand and visualize the interactions and investigate the function of corresponding DEMs, a miRNA-gene-pathway network was constructed ( Figure 5) using the data from Table 3. Through the interaction analysis, we identified potential functions of a few miRNAs in the ovary. The results showed that gga-miR-22-3p and gga-miR-365-2-5p were closely associated with reproductive processes; gga-miR-125b-3p and gga-miR-30d may play important roles in cell proliferation and apoptosis; gga-miR-32-3p and gga-miR-92-3p were mainly involved in hormone regulation. In addition, gga-miR-15b-5p, gga-miR-15c-5p, gga-miR-29a-3p and gga-miR-29c-3p were frequently involved in multiple biological processes.

RT-qPCR Validation
Nine mRNAs and nine miRNAs were selected for validation of the RNA-seq results using real time quantitative PCR (RT-qPCR). Our RT-qPCR results showed that the trend of differential expression of mRNAs and miRNAs in egg-laying and atrophic ovaries were consistent with the differential expression patterns observed in RNA-seq data ( Figure 6). Even though there were some differences in the numerical, which may be ascribed to the fact that the sensitivity of high-throughput sequencing and that of qRT-PCR detection method for specific mRNAs or miRNAs are relative and inconsistent. In general, the results indicates that our sequencing data were reliable.

Discussion
Egg production is an important economic trait of the goose. However, Xingguo grey goose has a strong seasonal reproductive characteristic, which seriously affects its egg production. A large number of studies on the mechanism of reproductive traits in geese are mainly focused on hormones and genes, but few researches have focused on the miRNAs involved in broodiness. In this study, a new insight into geese (Anser anser) ovary transcriptome profiles of egg-laying and broody were provided through deep sequencing.
The hypothalamic pituitary gonadal (HPG) axis regulates ovarian development through reproductive endocrine hormones including PRL, GnRH, LH, and FSH [1]. In this study, we demonstrated that the expressions of hormone-related genes were altered in the broody ovary. Follicle stimulating hormone receptor (FSHR) [22], growth hormone inducible transmembrane protein (GHITM) [23] and NPY [13] genes, which encode important receptor or signaling factors in intracellular signaling of hormones regulating folliculogenesis, were decreased in broody geese. Phospholipase C beta 1 (PLCB1), as an important enzyme for oocyte activation, promote the release of prolactin by mediating GnRH [24].
And DNA polymerase lambda (POLL) is a crucial enzyme in ovarian steroidogenesis [25]. They were also significantly lower expressed in broody geese. Moreover, lymphocyte antigen 96 (LY96) could lead to lower aromatase expression and reduced estradiol secretion in granulosa cells. Our results showed that it was significantly higher expressed in the broody geese. These results indicate that hormones played an important role in regulating broodiness in geese.
The transformation of primary follicles into secondary follicles is a key step in the development of early stages of folliculogenesis [26]. In our research, in comparison to the normal ovary with many primary follicles and secondary follicles in egg-laying geese, the atrophic ovaries of broody geese had numerous primary follicles but few secondary follicles, suggesting that they lacked the primary to secondary follicle transition. Consistent with this, in the atrophic ovaries there was higher expression of TNF receptor superfamily member 1A (TNFRSF1A) and secreted frizzled related protein 4 (SFRP4) genes, which promote oocyte apoptosis [27], inhibit angiogenesis in reproductive system [28], and then inhibit follicular development and ovulation.
Cell growth and death involve many responsible cellular processes, including cell proliferation, differentiation, transformation, survival and apoptosis. Tuberous sclerosis 2(TSC2) is phosphorylated by akt, which eliminates the inhibition of Rheb and activates mTOR to regulate cell growth, cell cycle and other physiological functions [29]. Zinc finger FYVE-type containing 9 (ZFYVE9) is a FYVE domain protein first identified as a binding partner for SMAD2/3, which mediated signal is essential for the proliferation of hepatic cells during the expansion of liver bud [30]. In the current study, these genes were significantly lower expressed in the atrophic ovaries. In addition, retinoic acid receptor beta (RARB), associates with the viability, apoptosis, differentiation and the iodine uptake function of cancer cell lines by regulating MAPK signaling pathway [31]. RB transcriptional corepressor like 1 (RBL1) mainly plays an important role in the context of G1-S transition and as an apoptotic trigger [32]. These apoptosis-related genes were significantly overexpressed in the atrophic ovaries. Our results suggested that the broody geese have relatively weaker cell growth and development but strong cell apoptosis in the ovary.
Compared to the normal ovaries in egg-laying geese, the absence of secondary follicles in broody geese were closely related to the reduced expression of genes significantly enriched in many signal transduction pathways including p53 signaling pathway, mTOR signaling pathway, calcium signaling pathway and TGF-beta signaling pathway, resulting in relatively weaker cell growth and development.
Newton suggested that p53 functions as a transcription factor and is responsible for the transactivation and repression of key genes involved in cell growth, apoptosis and the cell cycle [33].
Li et al. pointed out that mTOR is regarded as a central controller of cell growth in response to growth factors, cellular energy and nutrient levels [34].Some researchers figured out that one function of calcium signal is to activate the immediate early genes responsible for inducing resting cells (G0) to re-enter the cell cycle and may also promote the initiation of DNA synthesis at the G1/S transition [35]. According to multiple studies, the TGF-beta pathway promotes granulosa cell proliferation and follicular growth in the antral follicle stage [36][37]. In addition, we also have abundant genes that are highly expressed in atrophic ovaries, mainly involved in Wnt signaling pathway related to apoptosis.
These results indicated that broodiness in geese affected proliferation of granulosa cells and oocyte maturation mediated by the above pathways. found 116 differentially expressed microRNAs in chicken ovaries at broody and reproductive stages [40]. In this work, we detected 50 DEMs, many of which are related to reproduction. Our results are consistent with many studies, that is, miRNAs highly expressed in chicken ovaries is also abundant in goose ovaries, which indicates that these miRNAs are the basis of ovarian function and play an important role in its function. For example, miR-15 was shown to promote apoptosis through suppression of bcl-2 and other pro-proliferating and antiapoptotic substance and probably prevents proliferation of germ cell [41]. MiR-101 can inhibit the transformation of epithelial cells into mesenchymal cells by targeting ZEB1 and ZEB2 genes in ovarian cancer [42]. MiR-200 family was implicated in ovarian cancer initiation and progression via stage-specific regulation [43]. And these miRNAs were mainly expressed in atrophic ovaries in our study.
By targeting the RNA, miRNAs regulate and stabilize its translation, thus affecting the protein level of the RNA translation [44]. In this study, we selected 18 important intersection genes and 29 their corresponding DEMs to construct the miRNA-gene-pathway relationship network. This result improves the molecular regulatory network of ovarian growth and development in geese, and provides an effective reference for the study of miRNAs related to ovarian atrophy in geese. However, due to the target sequence of microRNAs may not be conservative, how these DE miRNAs and mRNAs interact with each other to regulate essential biological functions in ovary during counter-season production, remain unrequited and await further elucidation.

Conclusion
In the present study, mRNA and miRNA transcript profiles of the goose ovary were identified by highthroughput sequencing. Many of DEGs functioned predictably on the pathways of hormone response, follicular development, and apoptosis. Also a regulatory network of the molecular mechanisms of ovarian atrophy was constructed. Overall, this study expanded our understanding of the functional of miRNAs in ovarian atrophy and demonstrated that RNA-Seq is a powerful tool for examining the molecular mechanism in regulating broodiness.

Animals
For the purpose of this study, 100 newborn healthy laying geese aged 50 weeks (in late period of egg production) were randomly divided into two groups, one group used forced melting method and short illumination to regulate the reproduction cycle (counter-season egg-laying group), and another group was fed normally(broody group). At 70 weeks of age, three geese were selected from each group to collect ovaries. All animal care and experimental procedures were approved by the Institutional

Identification of miRNAs and Target Gene Prediction
Raw data (raw reads) of fastq format were firstly processed through custom perl and python scripts.
In this step, clean datas(clean reads) were obtained by removing reads containing ploy-N, with 5' adapter contaminants, without 3' adapter or the insert tag, containing ploy A or T or G or C and low quality reads from raw data. At the same time, Q20, Q30, and GC-content of the raw datas were calculated. Chose a certain range of length from clean reads to do all the downstream analyses.
Differential expression analysis of two groups was performed using the DESeq R package (1.16.1, EMBL Heidelberg, Germany). The P-values was adjusted using the Benjamini&Hochberg method.
Corrected P-value of 0.05 was set as the threshold for significantly differential expression by default.
Then, using BLAST software, we compared the predicted target gene sequences to the NR, Swiss-Prot, GO, COG, KEGG, KOG, Pfam database, and subsequently, the target gene annotation information was obtained. Finally, GO and KEGG analyses of predicted target genes were conducted.

Transcriptome sequencing and Analysis
A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations.
Sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Clean reads were generated by removing reads with adapter and lowquality sequences and mapped to the reference chicken genome and genes (Gallus gallus, Galgal4; available at https://www.ncbi.nlm.nih.gov/assembly/GCF_000002315.3) using TopHat 1.3.2(https://ccb.jhu.edu/software/tophat). Differential expression analysis of two groups was performed using the DESeq2 R package (1.16.1, EMBL Heidelberg, Germany). The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate.
Genes with an adjusted P-value <0.05 found by DESeq2 were assigned as differentially expressed. GO In order to verify the accuracy of high throughput sequencing data, we randomly selected and validated the differentially expressed miRNAs (n = 9) and mRNAs (n = 9) using qRT-PCR. Specific primers were designed by Primer 5.0 software dependent on GeneBank sequences (Table S1). The

Competing interests
All authors declared that no conflict of interest exists. manuscript. JJL and PR worked at the Illumina sequencing data generation and data preprocessing. JH and YQT conceived and designed the study, participated in coordination and contributed to draft manuscript writing and revisions. ZQL and XXZ carried out the RT-qPCR. YPL and SML granted, concept, designed the experiment and revised, given final approval version of the manuscript to be published. All authors read and approved the final manuscript.

Figure 5
The miRNA-gene-pathway network between eighteen intersection genes, and their corresponding pathways and differentially expressed miRNAs. RT-qPCR validation of differentially expressed mRNAs and miRNAs.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.