Zhuangbiao Zhang1†
Jishun Tang1,2†
Ran Di1Qiuyue Liu1Xiangyu Wang1
Shangquan Gan3Xiaosheng Zhang4Jinlong Zhang4
Mingxing Chu1*Wenping Hu1*- 1Key Laboratory of Animal Genetics and Breeding and Reproduction of Ministry of Agriculture and Rural Affairs, Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing, China
- 2Institute of Animal Husbandry and Veterinary Medicine, Anhui Academy of Agricultural Sciences, Hefei, China
- 3State Key Laboratory of Sheep Genetic Improvement and Healthy Production, Xinjiang Academy of Agricultural and Reclamation Sciences, Shihezi, China
- 4Tianjin Institute of Animal Sciences, Tianjin, China
Early studies have provided a wealth of information on the functions of microRNAs (miRNAs). However, less is known regarding their functions in the hypothalamus involved in sheep reproduction. To explore the potential roles of hypothalamic messenger RNAs (mRNAs) and miRNAs in sheep without FecB mutation, in total, 172 and 235 differentially expressed genes (DEGs) and 42 and 79 differentially expressed miRNAs (DE miRNAs) were identified in polytocous sheep in the follicular phase versus monotocous sheep in the follicular phase (PF vs. MF) and polytocous sheep in the luteal phase versus monotocous sheep in the luteal phase (PL vs. ML), respectively, using RNA sequencing. We also identified several key mRNAs (e.g., POMC, GNRH1, PRL, GH, TRH, and TTR) and mRNA–miRNAs pairs (e.g., TRH co-regulated by oar-miR-379-5p, oar-miR-30b, oar-miR-152, oar-miR-495-3p, oar-miR-143, oar-miR-106b, oar-miR-218a, oar-miR-148a, and PRL regulated by oar-miR-432) through functional enrichment analysis, and the identified mRNAs and miRNAs may function, conceivably, by influencing gonadotropin-releasing hormone (GnRH) activities and nerve cell survival associated with reproductive hormone release via direct and indirect ways. This study represents an integral analysis between mRNAs and miRNAs in sheep hypothalamus and provides a valuable resource for elucidating sheep prolificacy.
Introduction
Reproduction, one of the major factors significantly affecting the sheep industry, is a complicated but important physiological process. The success of reproduction is mainly dependent on the release of hormones, including gonadotropin-releasing hormone (GnRH) released from the hypothalamus, follicle-stimulating hormone (FSH) and luteinizing hormone (LH), which are both secreted from the pituitary (Cao et al., 2018a). Following the release of hormones, a series of events associated with reproduction, such as ovulation and fertilization, could occur.
It is well known that reproductive traits, such as litter size, are controlled by minor polygene. Researchers have found several major fecundity genes which considerably influence sheep prolificacy, such as bone morphogenetic protein receptor IB (BMPRIB), bone morphogenetic protein 15 (BMP15) (Chu et al., 2007), and growth differentiation factor 9 (GDF9) (Chu et al., 2011). FecB is a mutation in BMPRIB occurring in base 746 from A to G. This base change further results in changes in protein function due to a key amino acid transition from glutamine to arginine (Fogarty, 2009). Sheep with one copy of the FecB mutation can experience significant increase in litter size, by 0.67, while this increase is about 1.5 when there are two mutated copies (Liu et al., 2014). Moreover, this mutation was also detected in diverse sheep species, such as Booroola Merino sheep (Mulsant et al., 2001) (Australia), Garole sheep (Polley et al., 2010) (India), Hu sheep (Davis et al., 2006) (China), and Small Tail Han sheep (STH sheep; China) (Davis et al., 2006). STH sheep, an indigenous species in China, has attracted much attention for its excellent traits (Liu et al., 2016; Chao et al., 2017), especially the higher prolificacy (Davis et al., 2006). Furthermore, STH sheep can be divided into three genotypes based on the effects of FecB mutation, better known as FecB BB (with two-copy FecB mutations), FecB B+ (with one-copy FecB mutation), and FecB++ (with no FecB mutation). Usually, compared to sheep with the other two genotypes, STH sheep with FecB++ show a monotocous phenomenon. However, the fact is that there are STH sheep with FecB++ and which show a polytocous phenomenon (Davis et al., 2006), and how this mechanism was established remains largely unclear.
With advances in sequencing, the application of RNA sequencing (RNA-seq) in animals, including sheep (Jiang et al., 2014; Zhang et al., 2019a; Zhang et al., 2019b), mice (Beck et al., 2018), and cattle (Correia et al., 2018), enables integral analysis of the expression profiling of mRNA and miRNAs. Therefore, RNA-seq has been widely used to understand some complex traits. Regarding the generation of miRNA, precursor miRNA is transcribed mainly by RNA polymerase II, then processed into mature miRNA (Gebert and Macrae, 2019). Significantly, miRNAs play pivotal roles in life processes, such as muscle growth (Cao et al., 2018c), fleece and hair development (Liu et al., 2018), and neural development (Schratt et al., 2006). Additionally, reproduction is an extremely complex process, and the use of RNA-seq may contribute to enhancing our understanding of sheep fecundity. By comparing the mRNA and miRNA expression patterns in European mouflon and sheep, a research (Yang et al., 2018) found several key mRNAs, such as INHBA, SPP1, and ZP2, and miRNAs, such as miR-374a and miR-9-5p, which may be responsible for the success of female sheep reproduction. Pokharel et al. (2018) detected and characterized some key miRNAs and mRNAs in sheep ovary which may be responsible for sheep prolificacy. Thereby, the identification and functional analysis of mRNAs and miRNAs and characterization of their mutual interaction through sequencing technology may provide new insights into the prolific mechanism in STH sheep with the FecB++ genotype, which has so far been difficult to elucidate using standard approaches.
Therefore, in the present study, we applied transcriptomics analysis in PF vs. MF and PL vs. ML to identify DEGs and DE miRNAs and analyze their potential functions, expecting to elucidate the potential prolific mechanism in sheep with the FecB++ genotype and act as a reference for other female mammals.
Material and Methods
Preparation of Animals
First, the TaqMan probe (Liu et al., 2017) was applied to genotype the sheep population (n = 890). Then, 12 sheep with no significant differences in sheep age, weight, height, body length, chest circumference, and tube circumference were selected from 142 STH sheep with the FecB++ genotype and grouped into the polytocous group (n = 6, litter size ≥2) and monotocous group (n = 6, litter size = 1) according to their litter size records. Additionally, all the sheep were bred under the same conditions, with free access to water and feed, in a sheep farm of the Tianjin Institute of Animal Sciences.
All selected sheep were processed by estrus synchronization with Controlled Internal Drug Releasing Device (CIDR; progesterone 300 mg; Zoetis Australia Pty. Ltd., NSW, Australia) for 12 days. The six sheep, comprising three polytocous sheep and three monotocous sheep, were slaughtered within 45–48 h after CIDR removal (follicular phase), the remaining six sheep were slaughtered on day 9 after CIDR removal (luteal phase). Finally, the selected sheep were divided into four groups, including polytocous sheep in the follicular phase (PF), polytocous sheep in the luteal phase (PL), monotocous sheep in the follicular phase (MF), and monotocous sheep in the luteal phase (ML), on the basis of their littering record and estrous cycle.
Preparation of Tissues, RNA Extraction, and Sequencing
Hypothalamic tissues were collected from 12 killed sheep and immediately stored at −80°C until being used. Then, total RNA was isolated using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) under the manufacturer’s instructions, and the quality and integrity of isolated RNA were assessed by an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and electrophoresis. The high-quality RNA of 3 μg of each sample was used to build the mRNA library using a NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, USA), which has been described in our previous work (Zhang et al., 2019b). All the sequencing works were conducted in Annoroad Gene Technology Co., Ltd. (Beijing, China).
The fragments with lengths of 18–30 nt, which were obtained from total RNA through the gel separation technique, were used as templates to synthesize the first strand of complementary DNA (cDNA). The second strand of cDNA was also synthesized in the presence of deoxynucleoside triphosphates (dNTPs), ribonuclease H, and DNA polymerase I. Then the obtained double-stranded cDNA was processed with end-repair, the addition of base A and sequencing adaptors, and uracil-N-glycosylase (UNG) enzyme digestion. Finally, polymerase chain reaction was conducted to build the miRNA library.
In addition, a paired-end sequencing approach for mRNAs and miRNAs was conducted using an Illumina HiSeq 2500.
Quality Control, Mapping and Assembly
Raw reads were filtered using in-house software of fqtools_plus-v2.0.0 according to strict criteria, including removing reads with adaptor contaminants, low-quality reads, and reads with N bases accounting for more than 5%. Then, HiSAT2 (Kim et al., 2015) was used to map the cleaned reads to the reference genome (Oarv3.1), and both the sheep reference genome and genome annotation file were downloaded from ENSEMBL (http://www.ensembl.org/index.html). Subsequently, StringTie 1.3.2d (Pertea et al., 2015) was used to assemble transcripts of mRNAs.
Several criteria were also implemented to generate clean miRNA reads, including removing reads without a 3′ adapter, reads without insert fragment, reads with lengths beyond the normal range, raw reads containing too much A/T, and some low-quality reads using in-house scripts. Furthermore, the cleaned data of miRNA were matched against the sheep reference genome (Oarv3.1) by Bowtie v1.1.2 (Langmead et al., 2009).
Differential Expression and Functional Enrichment Analysis of mRNAs
To validate the expression level of mRNAs, the fragments per kilobase per million mapped reads (FPKM) values (Trapnell et al., 2010) were calculated to represent the gene expression level, and DESeq 2-1.4.5 (Wang et al., 2010) was also used to detect the DEGs between two comparisons based on FPKM values. Additionally, a gene with fold change >1.5 and p < 0.05 was considered as a DEG in PF vs. MF and PL vs. ML. In addition, we also performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. We first downloaded the Uniprot database, where each sequence contains the GO annotation and KEGG annotation species (sheep) of the sequence as well as gene and protein names. All genes of sheep to be analyzed were compared with the Uniprot database by blast (NCBI-blast 2.2.28) to find the best alignment result for each sequence, and corresponding to GO and KEGG annotation results. Then, we also downloaded the corresponding relationship between the entry name and number provided on the websites of GO and KEGG, as well as the classification hierarchy file, and summarize the GO and KEGG classification of the genes we obtained. Lastly, a particular GO term or KEGG pathway with a hypergeometric p value < 0.05 was thought to indicate significant enrichment.
Differential Expression Analysis and Prediction of Target Genes of miRNAs
The miRDeep v2.0.0.8 (Friedländer et al., 2012) was applied to identify the known and novel miRNAs by mapping clean reads and hairpins to mature miRNAs recorded in the miRbase database (Griffiths-Jones, 2006). In addition, transcripts per million (TPM) were calculated to represent miRNA expression levels on the basis of the reads number. DESeq2-1.4.5 (Wang et al., 2010) was also applied to identify DE miRNAs in PF vs. MF and PL vs. ML, and the threshold of fold change >1.5, p < 0.05 was considered to indicate differential expression. Furthermore, miRanda v3.3a (Enright et al., 2004) was used to predict the target genes of miRNAs.
Integral miRNA–mRNA Networks Analysis
To precisely identify key DE miRNAs and DEGs associated with reproduction, a network containing DE miRNAs and DE mRNAs, on the basis of miRNA functions (Gebert and Macrae, 2019), was built using Cytoscape_v3.5.0 (Shannon et al., 2003), and only mRNAs exhibiting negative relationship with miRNAs were included in miRNA–mRNA interaction networks.
Data Validation
In order to validate the accuracy of sequencing data, four DEGs, including CRH, FOXG1, TTR, and POMC, and four DE miRNAs, including oar-miR-433-3p, oar-miR-495-3p, oar-miRNA-16b, and oar-miR-143, were selected for data validation. First, the primers of DEGs and DE miRNAs were synthesized by Beijing Tianyi Huiyuan Biotechnology Co., Ltd. (Beijing, China) (Supplementary Table 1) for subsequent reverse transcription, which was performed using PrimeScript™ RT reagent kit (TaKaRa) for mRNAs and miRcute Plus miRNAs First-Strand cDNA Kit (TIANGEN, Beijing, China) for miRNA. Furthermore, quantitative PCR (qPCR) was conducted with the SYBR Green qPCR Mix (TaKaRa, Dalian, China) for mRNAs and miRcute Plus miRNA qPCR Kit (TIANGEN, Beijing, China) for miRNAs using a RocheLight Cycler®480 II system (Roche Applied Science, Mannheim, Germany). In addition, β-actin (for mRNA) and U6 small nuclear RNA (snRNA; for miRNA) were utilized as reference gene/miRNA to calculate the relative expression level with the method of 2-ΔΔct (Livak and Schmittgen, 2001). The qPCR for mRNAs was conducted in the following procedure: initial denaturation at 95°C for 5 minutes, followed by 40 cycles of denaturation at 95°C for 5 s, then annealing at 60°C for 30 s. While the qPCR for miRNA was conducted in the following procedure: initial denaturation at 95°C for 15 minutes, followed by 40 cycles of denaturation at 94°C for 20 s, then annealing at 60°C for 34 s. All the qPCR results were presented as the mean ± SD.
Results
mRNA and miRNA Profiling
To fully characterize the globally hypothalamic mRNA and miRNA expression differences between sheep with the same genotype but different litter sizes, RNA-seq was used to detect their expression profile in the hypothalamus. In total, RNA-seq for mRNA generated approximately 1,519 million raw reads and 1,460 million clean reads (Supplementary Table 2) after data filtering. Overall, 21,221 mRNAs were identified (Supplementary Table 3) after mapping to sheep genome, and our results also suggested that many mRNAs were located in the intergenic region (nearly 45%), followed by the intron (about 35%) and exon (more than 20%) regions (Figure 1A and Supplementary Table 4).
Figure 1 Mapping region and chromosome distribution of identified mRNAs. (A) Mapping region of identified genes at the reference genome in polytocous sheep in the follicular phase (PF) (a), polytocous sheep in the luteal phase (PL) (b), monotocous sheep in the follicular phase (MF) (c), and monotocous sheep in the luteal phase (ML) (d). (B) Chromosome distribution of identified genes from the hypothalami in the PF, MF, PL, and ML.
Regarding the expression level of mRNAs, our results showed that the FPKM of those genes obtained from RNA-seq at <50 constituted nearly 90%, and the high-expression genes, i.e., those with FPKM >500, constituted about 0.5% (Supplementary Table 3), which suggested that the data obtained from the hypothalamus via RNA-seq were relatively reasonable. Furthermore, the chromosome distribution of mRNAs indicated that chromosome 3 contains 9.79% of the genes identified from the hypothalamus, followed by chromosome 1 (9.55%) and chromosome 2 (7.22%) (Figure 1B and Supplementary Table 5). Additionally, the number of DEGs identified from PF vs. MF (Figure 2A and Supplementary Table 6) and PL vs. ML (Figure 2B and Supplementary Table 6) were 172 and 235, respectively. Among these DEGs, 79 and 90 were upregulated, while 93 and 145 were downregulated in PF vs. MF and PL vs. ML, respectively. In addition, the expression density of DEGs displayed obviously different expression patterns between PF and MF, and between PL and ML (Figures 2C, D).
Figure 2 Differentially expressed genes (DEGs) analysis. (A) Volcano plot of identified genes in PF vs. MF, where red and green represent up- or downregulation, respectively, same below. (B) Volcano plot of identified genes in PL vs. ML. (C) Heat maps showing the expression intensity of 794 DEGs in the follicular phase, including PF and MF. (D) Heat maps showing the expression intensity of 1,044 DEGs in the luteal phase, including PL and ML.
Regarding miRNAs, RNA-seq generated approximately 315 million raw reads and 267 million clean reads (Supplementary Table 7) with lengths ranging from 18 to 30 nt (Figure 3A) after removing low-quality reads. Overall, 623 miRNAs were detected (Supplementary Table 8). In addition, the chromosome distribution of identified miRNAs was also determined. As Figure 3B shows, the chromosome distribution of miRNAs from 1 to X varies (Supplementary Table 9), and most of the identified miRNAs were located at chromosome 3 (nearly 40%), followed by chromosome 9 (nearly 15%) and chromosome 18 (nearly 9%). Interestingly, chromosome 3 also contains the most mRNAs (Figure 3B). Also, a diversity of non-coding RNAs (ncRNAs), including transfer RNAs (tRNAs), snRNAs, miRNAs, etc., were also identified (Figure 3C and Supplementary Table 10), and the known miRNAs account only for a small part of all the identified ncRNAs. In addition, the target genes of miRNAs in PF vs. MF and PL vs. ML were predicted to be 1,611 and 2,120, respectively (Supplementary Table 11).
Figure 3 Characterization of microRNA (miRNA) profiling and the percentage of detected miRNAs from ncRNAs. (A) Length distribution of clean reads from identified miRNA fragments. (B) The chromosome distribution of identified miRNAs from hypothalami. (C) Categories of identified non-coding RNAs (ncRNAs) via sequencing in PF (a), PL (b), MF (c), and ML (d).
Additionally, the DE miRNAs identified from PF vs. MF and PL vs. ML were 42 and 79, respectively. Of these DE miRNAs, 20 and 23 were upregulated, while 22 and 56 were downregulated, respectively (Figure 4A and Supplementary Table 12). In addition, the expression density of DEGs displayed obviously different expression patterns between PF and MF, and between PL and ML (Figures 4B, C).
Figure 4 Differentially expressed (DE) microRNA (miRNA) analysis. (A) DE miRNAs in PF vs. MF and PL vs. ML. Heat maps showing the expression intensity of 42 and 79 DE miRNAs in the follicular phase including PF and MF (B) and the luteal phase including PL and ML (C), the names of miRNAs were also labeled.
GO and KEGG Enrichment Analysis of DEGs
To better understand the potential functions of the DEGs, GO term and KEGG pathway analyses were performed. In GO analysis, the most enriched term in PF vs. MF was the MHC protein complex (GO:0042611). Other GO terms related to the MHC protein were also enriched, such as MHC class II protein complex binding (GO:0023026) and MHC protein complex binding (GO:0023023), indicating the crucial role of the MHC protein in the hypothalamic functions (Figure 5A and Supplementary Table 13). Regarding PL vs. ML, the top 2 enriched terms were the immune system process (GO:0002376) and immune response (GO:0006955). In addition, some GO terms associated with chemokine receptors, including CXCR3 chemokine receptor binding(GO:0048248) and chemokine receptor binding (GO:0042379), were also highly enriched, suggesting the important roles of the immune system and chemokine receptors in the hypothalamus at the luteal phase (Figure 5A and Supplementary Table 13).
Figure 5 Functional enrichment analysis of DEGs. (A) Top enriched GO terms at the biological process, molecular function, and cellular component level in PF vs. MF and PL vs. ML, in addition, the gray represents no enrichment, same below. (B) Top enriched KEGG pathways in PF vs. MF and PL vs. ML.
KEGG analysis in PF vs. MF (Figure 5B and Supplementary Table 14) showed that the most enriched pathway was type I diabetes mellitus (map04940). In addition, other metabolic pathways, such as alpha-linolenic acid metabolism (map00592) and arachidonic acid metabolism(map00590), were also enriched. Regarding PL vs. ML, the top enriched pathways were cytokine–cytokine receptor interaction (map04060). A pathway named the Jak-STAT signaling pathway (map04630), which has been found to participate in the reproductive process (Ko et al., 2018), was also enriched.
Analysis of Integrated miRNA–mRNA Co-Expression Network
To fully understand the potential reproductive roles of miRNAs, we built interactome networks using DE miRNAs and their targets (DEGs). In total, 42 DE miRNAs (novel miRNAs) in PF vs. MF were predicted to target 1,611 genes (Supplementary Table 15). The number of overlapped genes, which means the target genes were also DEGs, was 8 (Figure 6A and Supplementary Table 16). An mRNA–miRNA co-expression network was then constructed, where 5 DEGs were targeted by 3 novel miRNAs (Figure 6B). Regarding PL vs. ML, 38 known and 41 novel DE miRNAs were predicted to target 1,747 and 1,659 genes (Supplementary Table 15), and the numbers of overlapped genes were 179 and 9, respectively (Figures 6C, D and Supplementary Table 16). The main upregulated miRNA–mRNA co-expression network suggested that 55 DEGs were targeted by 11 DE miRNAs containing the top 10 upregulated known miRNAs and one novel miRNA (Figure 6E). The main downregulated miRNA–mRNA co-expression network suggested that 33 DEGs were targeted by 11 DE miRNAs containing the top 10 downregulated known miRNAs and one novel miRNA (Figure 6F).
Figure 6 Overview of mRNA–miRNA networks. (A) Overlapped genes in PF vs. MF between DEGs and miRNA-targeted genes. (B) Hypothalamus network in PF vs. MF of four miRNAs containing three upregulated miRNAs and one downregulated miRNA and five target genes, where large and small ellipses represent miRNAs and DEGs, respectively, in addition, red and green represent up- or downregulation, respectively; same below. (C) Overlapped genes in PL vs. ML between DEGs and predicted target genes negatively modulated by known miRNAs. (D) Overlapped genes in PL vs. ML between DEGs and novel miRNA-targeted genes. (E) Main hypothalamic upregulated network in PL vs. ML containing the top 10 upregulated miRNAs and one novel miRNA and 55 target genes. (F) Main hypothalamic downregulated network in PL vs. ML containing the top 10 downregulated miRNAs and one novel miRNA and 33 target genes.
Data Validation
>In order to assess the accuracy of sequencing, qPCR was applied to verify the RNA-seq data. The results indicated that both mRNAs and miRNAs in sheep hypothalamus displayed expression patterns similar to the sequencing results (Figure 7), demonstrating the reliability of the data generated from RNA-seq.
Figure 7 Data validation of mRNAs (A) and miRNAs (B) by qPCR in PF, PL, MF and ML, meanwhile ** represents p < 0.01, while * represents p < 0.05. FOXG1: Forkhead box L1, CRH: corticotropin-releasing hormone, TTR: transthyretin, POMC proopiomelanocortin.
Discussion
In this study, we initially identified 172 and 235 DEGs, and 42 and 79 DE miRNAs in two comparisons (PF vs. MF and PL vs. ML) through RNA-seq. Of these DE miRNAs, miRNA family members including the let-7 and oar-miRNA-200 family exhibited differential expression levels. Furthermore, one study detecting 48 DE miRNAs from sheep ovary, including the let-7 and oar-miRNA-200 family members, suggested that those identified miRNAs were differentially expressed in seasonal and non-seasonal sheep breeds (Zhai et al., 2018). Therefore, some miRNAs, such as let-7 and oar-miRNA-200 family members may not be only species-specific but also phase- or fecundity-specific in sheep. In addition, some miRNAs, including miRNA-138 and miRNA-212, were detected in rat hypothalamus (Amar et al., 2012), which differed significantly from miRNAs identified in sheep hypothalamus (both miRNA-138 and miRNA-212 in our results failed to be detected). Besides, several miRNAs, such as miRNA-200 family members, were conserved in the hypothalamus of mice (Choi et al., 2008; Crépin et al., 2014), rat (Sangiaoalvarellos et al., 2014), and zebrafish (Garaffo et al., 2015), as well as sheep (our results). In summary, we confirmed that several miRNAs are conserved in many animals, but there were also miRNAs that showed a species-specific distribution in the hypothalamus, which means those differences may be responsible for the differences between sheep and rats, and even other non-mammals.
Functional Analysis of DEGs in PF vs. MF
In the functional enrichment analysis of DEGs in PF vs. MF, several key genes, including prolactin (PRL), proopiomelanocortin (POMC), and gonadotropin releasing hormone 1 (GNRH1), were found to participate in the reproductive process. Some researchers have proven that PRL and E2 could respond rapidly to stimulation in the arcuate nucleus (ARC) of rat hypothalamic slices (Nishihara and Kimura, 1989). Araujo-Lopes et al. (2014) revealed that PRL could regulate the activities of GnRH through modulating kisspeptin neurons in the ARC of female rats and inhibit LH secretion, causing a series of alterations in the estrous cycle. Our results indicated that the expression of PRL in PF was more than three times that of PRL in MF. Therefore, coupled with the inhibitory role of PRL on LH, we speculate that PRL may affect LH or FSH activities by influencing the pulsatile GnRH wave in the hypothalamus.
POMC neurons, as a key upstream factor affecting hypothalamic hormone release, were found to be sensitive to metabolic hormones such as leptin (Wilson and Enriori, 2015) and enhance kisspeptin neuron activities in rodents, resulting in increased GnRH secretion (Muroi and Ishii, 2016). Leptin can act in the hypothalamus directly, eliciting the release of GnRH (Guzmán et al., 2019), and promoting the expression of POMC (Perello et al., 2007). Although the stimulatory effects of POMC on kisspeptin have been known for a long time, how this signaling is established remains poorly understood (Saedi et al., 2018). Significantly, our results indicated that the expression of POMC in PF was relatively lower than in MF, while GNRH1, which has been reported to play a key role in determining sheep litter size (An et al., 2013), displayed a reverse expression pattern between PF and MF. Therefore, we hypothesized that a negative regulatory relationship between POMC and GNRH1 may exist in sheep hypothalamus.
Functional Analysis of DEGs in PL vs. ML
In functional enrichment analysis of DEGs in PL vs. ML, some pathways including the Jak-STAT signaling pathway (PRL, GH, CRLF2, ENSOARG00000007618, ENSOARG00000016231, and IL2RB) were highly enriched. The current study argued that the Jak-STAT signaling pathway in mice was involved in GnRH activities (Ko et al., 2018). PRL, as mentioned above, plays an important role in GnRH activities (Araujo-Lopes et al., 2014). The expression of PRL was detected not only in the follicular phase but also in the luteal phase, and interestingly, there was a reverse expression pattern of PRL between PF vs. MF and PL vs. ML, suggesting its crucial roles in reproduction. The effects of leptin on GnRH release have been revealed (Guzmán et al., 2019), and the infusion of leptin into the arcuate nucleus in rats could cause PRL release (Watanobe, 2010), which suggested that PRL can be a downstream factor activated by leptin to function in GnRH activities. In addition, the overexpression of growth hormone (GH) could disrupt the state of reproduction, mainly through mediating leptin activities (Chen et al., 2018). Additionally, estrogen could play an inhibitory role on GH in vivo (Leung et al., 2003). Collectively, considering the effects of PRL and GH on leptin, we speculated that GH, leptin, and PRL may coordinate to inhibit GnRH release.
The Regulatory Network of miRNA–mRNA After Transcription in PF vs. MF
To better understand the functions of miRNAs, a negative interactome containing 5 mRNAs and 4 miRNAs in PF vs. MF was built. Cyclin-dependent kinase 3 (CDK3), targeted by Novel_237, was reported that the downregulation of activities of CDK3-related kinase could promote cell apoptosis in the rat (Braun et al., 1998). Immediate early response 3 (IER3), targeted by Novel_327, was also involved in enhancing (Zhou et al., 2017) or mediating (Jin et al., 2015) cell apoptosis. Polycystic kidney and hepatic disease gene 1 (PKHD1), targeted by Novel_401, has been discovered to induce cell apoptosis, after being downregulated through the PI3K and NF-κB pathways (Sun et al., 2011). Furthermore, our sequencing data indicated that CDK3 and IER3 were downregulated while PKHD1 was upregulated in PF vs. MF. All in all, we hypothesized that more nerve cell apoptosis occurred in MF than PF, which may further influence hormone activities associated with reproduction and may lead to the final observed litter size differences.
The Regulatory Network of miRNA–mRNA After Transcription in PL vs. ML
The regulatory network of miRNA–mRNA after transcription in PL vs. ML was divided into two main negative networks: the main upregulated and the main downregulated network. In the main upregulated network, thyrotropin-releasing hormone (TRH), co-regulated by oar-miR-379-5p, oar-miR-30b, oar-miR-152, oar-miR-495-3p, oar-miR-143, oar-miR-106b, oar-miR-218a, and oar-miR-148a, has been reported to function in GnRH release (see below). Triclosan in mice was found to reduce the production of TRH and thyroid-stimulating hormone (TSH), and this decreased effect could further cause hyperprolactinemia. Hyperprolactinemia was suggested to cause a suppressive effect on kisspeptin expression, resulting in deficits in reproductive and endocrine function (Cao et al., 2018b). In addition, TRH can not only stimulate PRL release but also inhibit LH release, and this inhibitory effects may occur through prohibiting the release of GnRH (Araujo-Lopes et al., 2014). Collectively, TRH in the hypothalamus may be responsible, at least in part, for the suppression of GnRH activities.
In the main downregulated network of miRNAs, transthyretin (TTR) was reversely regulated by oar-miR-432. The expression level of TTR’s in rats could be enhanced by progesterone via progesterone receptors both in vitro and in vivo (Quintela et al., 2011), and a similar upregulated effect of TTR caused by progesterone in mouse uterus was also observed (Diao et al., 2010). Furthermore, TTR could drive the nuclear translocation of insulin-like growth factor 1 receptor (IGF-1R) (Vieira et al., 2015), which could lead to functional changes in insulin-like growth factor 1 (IGF1). Interestingly, the stimulatory effect of IGF1 on GnRH release has been discovered (Hiney et al., 2009). Therefore, we speculated that the negative feedback effects of progesterone on GnRH release may be mediated by TTR, which reduces the binding probability between IGF1 and its receptor, further resulting in a suppression of GnRH activities.
All results indicated that several key DEGs and DE miRNAs in the hypothalamus directly or indirectly participate in hormone activities associated with reproduction, and further studies involving gene/miRNA knockout or overexpression could help us to understand their real functions in female reproductive traits.
Conclusion
As far as we know, this study provides the first integral mRNA–miRNA interactome in sheep without FecB mutation from the perspective of the hypothalamus. We identified several DEGs (e.g., POMC, GNRH1, PRL, TRH, and TTR) and mRNA–miRNA pairs (e.g., TRH coagulated by oar-miR-379-5p, oar-miR-30b, oar-miR-152, oar-miR-495-3p, oar-miR-143, oar-miR-106b, oar-miR-218a and oar-miR-148a and PRL regulated by oar-miR-432) from the RNA-seq data obtained from sheep hypothalamus, which may function through influencing the activities of GnRH. Our results provide novel insights into the prolificacy mechanism of sheep, which may facilitate the discovery of novel major genes and a deeper understanding of female sheep reproduction.
Data Availability Statement
All the data obtained from RNA-seq has been deposited in the Sequence Read Archive database under the bioproject numbers PRJNA529384 and PRJNA532808.
Ethics Statement
The animal study was reviewed and approved by the Science Research Department (in charge of animal welfare issues) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China). Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
WH and MC designed the research. ZZ wrote the paper. JT, ZZ, WH, SG, XZ, and JZ collected the data. ZZ performed the study. ZZ and JT analyzed data. MC and WH revised the final manuscript. All authors reviewed the manuscript and approved the final version.
Funding
This work was supported by the National Natural Science Foundation of China (31501941, 31772580 and 31472078), the Genetically Modified Organisms Breeding Major Program of China (2016ZX08009-003-006 and 2016ZX08010-005-003), the Earmarked Fund for China Agriculture Research System (CARS-38), the Central Public-Interest Scientific Institution Basal Research Fund (2018-YWF-YB-1, Y2017JC24, 2017ywf-zd-13), the Agricultural Science and Technology Innovation Program of China (ASTIP-IAS13, CAAS-XTCX2016010-01-03, CAAS-XTCX2016010-03-03, CAAS-XTCX2016011-02-02), the China Agricultural Scientific Research Outstanding Talents and Their Innovative Teams Program, the China High-level Talents Special Support Plan Scientific and Technological Innovation Leading Talents Program (W02020274), and the Tianjin Agricultural Science and Technology Achievements Transformation and Popularization Program (201704020), Joint Funds of the National Natural Science Foundation of China and the Government of Xinjiang Uygur Autonomous Region of China (U1130302). The APC was funded by the National Natural Science Foundation of China (31772580).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Annoroad Gene Technology (Beijing) Co., Ltd. for RNA-sequencing.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.01296/full#supplementary-material
References
Amar, L., Benoit, C., Beaumont, G., Vacher, C. M., Crepin, D., Taouis, M. (2012). MicroRNA expression profiling of hypothalamic arcuate and paraventricular nuclei from single rats using Illumina sequencing technology. J. Neurosci. Methods 209, 134–143. doi: 10.1016/j.jneumeth.2012.05.033
An, X. P., Hou, J. X., Zhao, H. B., Li, G., Bai, L., Peng, J. Y., et al. (2013). Polymorphism identification in goat GNRH1 and GDF9 genes and their association analysis with litter size. Anim. Genet. 44, 234–238. doi: 10.1111/j.1365-2052.2012.02394.x
Araujo-Lopes, R., Crampton, J. R., Aquino, N. S., Miranda, R. M., Kokay, I. C., Reis, A. M., et al. (2014). Prolactin regulates kisspeptin neurons in the arcuate nucleus to suppress LH secretion in female rats. Endocrinology 155, 1010–1020. doi: 10.1210/en.2013-1889
Beck, A. S., Wood, T. G., Widen, S. G., Thompson, J. K., Barrett, A. D. T. (2018). Analysis by deep sequencing of discontinued neurotropic yellow fever vaccine strains. Sci. Rep. 8, 13408. doi: 10.1038/s41598-018-31085-2
Braun, K., Hölzl, G., Pusch, O., Hengstschläger, M. (1998). Deregulated expression of CDK2- or CDK3-associated kinase activities enhances c-Myc-induced apoptosis. DNA Cell Biol. 17, 789–798. doi: 10.1089/dna.1998.17.789
Cao, C., Ding, Y., Kong, X., Feng, G., Xiang, W., Chen, L., et al. (2018a). Reproductive role of miRNA in the hypothalamic-pituitary axis. Mol. Cell. Neurosci. 88, 130–137. doi: 10.1016/j.mcn.2018.01.008
Cao, X. Y., Hua, X., Xiong, J. W., Zhu, W. T., Zhang, J., Chen, L. (2018b). Impact of triclosan on female reproduction through reducing thyroid hormones to suppress hypothalamic kisspeptin neurons in mice. Front. Mol. Neurosci. 11, 6. doi: 10.3389/fnmol.2018.00006
Cao, Y., You, S., Yao, Y., Liu, Z. J., Hazi, W., Li, C. Y., et al. (2018c). Expression profiles of circular RNAs in sheep skeletal muscle. Asian-australas. J. Anim. Sci. 31, 1550–1557. doi: 10.5713/ajas.170563
Chao, T., Wang, G., Ji, Z., Liu, Z., Hou, L., Wang, J. M., et al. (2017). Transcriptome analysis of three sheep intestinal regions reveals key pathways and hub regulatory genes of large intestinal lipid metabolism. Sci. Rep. 7, 5345. doi: 10.1038/s41598-017-05551-2
Chen, J., Cao, M., Zhang, A., Shi, M., Tao, B., Li, Y., et al. (2018). Growth hormone overexpression disrupts reproductive status through actions on leptin. Front. Endocrinol. 9, 131. doi: 10.3389/fendo.2018.00131
Choi, P. S., Zakhary, L., Choi, W. Y., Caron, S., Alvarez-Saavedra, E., Miska, E. A., et al. (2008). Members of the miRNA-200 family regulate olfactory neurogenesis. Neuron 57, 41–55. doi: 10.1016/j.neuron.2007.11.018
Chu, M. X., Liu, Z. H., Jiao, C. L., He, Y. Q., Fang, L., Ye, S. C., et al. (2007). Mutations in BMPR-IB and BMP-15 genes are associated with litter size in Small Tailed Han sheep (Ovis aries). J. Anim. Sci. 85, 598–603. doi: 10.2527/jas.2006-324
Chu, M. X., Yang, J., Feng, T., Cao, G. L., Fang, L., Di, R., et al. (2011). GDF9 as a candidate gene for prolificacy of Small Tail Han sheep. Mol. Biol. Rep. 38, 5199–5204. doi: 10.1007/s11033-010-0670-5
Correia, C. N., McLoughlin, K. E., Nalpas, N. C., Magee, D. A., Browne, J. A., Rue-Albrecht, K., et al. (2018). RNA Sequencing (RNA-Seq) reveals extremely low levels of reticulocyte-derived globin gene transcripts in peripheral blood from horses (Equus caballus)) and cattle (Bos taurus). Front. Genet. 9, 278. doi: 10.3389/fgene.2018.00278
Crépin, D., Benomar, Y., Riffault, L., Amine, H., Gertler, A., Taouis, M. (2014). The over-expression of miR-200a in the hypothalamus of ob/ob mice is linked to leptin and insulin signaling impairment. Mol. Cell. Neurosci. 384, 1–11. doi: 10.1016/j.mce.2013.12.016
Davis, G. H., Balakrishnan, L., Ross, I. K., Wilson, T., Galloway, S. M., Lumsden, B. M., et al. (2006). Investigation of the Booroola (FecB) and Inverdale (FecX(I)) mutations in 21 prolific breeds and strains of sheep sampled in 13 countries. Anim. Reprod. Sci. 92, 87–96. doi: 10.1016/j.anireprosci.2005.06.001
Diao, H., Xiao, S., Cui, J., Chun, J., Xu, Y., Ye, X. (2010). Progesterone receptor-mediated up-regulation of transthyretin in preimplantation mouse uterus. Fertil. Steril. 93, 2750–2753. doi: 10.1016/j.fertnstert.2010.01.009
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., Marks, D. S. (2004). MicroRNA targets in Drosophila. Genome Biol. 5, R1. doi: 10.1186/gb-2003-5-1-r1
Forgaty, N. M. (2009). A review of the effects of the Booroola gene (FecB) on sheep production. Small Ruminant Res. 85, 75–84. doi: 10.1016/j.smallrumres.2009.08.003
Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Garaffo, G., Conte, D., Provero, P., Tomaiuolo, D., Luo, Z., Pinciroli, P., et al. (2015). The Dlx5 and FOXG1 transcription factors, linked via miRNA-9 and -200, are required for the development of the olfactory and GnRH system. Mol. Cell. Neurosci. 68, 103–119. doi: 10.1016/j.mcn.2015.04.007
Gebert, L. F. R., Macrae, I. J. (2019). Regulation of microRNA function in animals. Nat. Rev. Mol. Cell Biol. 20, 21–37. doi: 10.1038/s41580-018-0045-7
Griffiths-Jones, S. (2006). MiRBase : the MicroRNA Sequence Database. Methods Mol. Biol. 342, 129–138. doi: 10.1385/1-59745-123-1:129
Guzmán, A., Hernández-Coronado, C. G., Rosales-Torres, A. M., Hernández-Medrano, J. H. (2019). Leptin regulates neuropeptides associated with food intake and GnRH secretion. Ann. Endocrinol. 80, 38–46. doi: 10.1016/j.ando.2018.07.012
Hiney, J. K., Srivastava, V. K., Pine, M. D., Les Dees, W. (2009). Insulin-like growth factor-I activates KISS-1 gene expression in the brain of the prepubertal female rat. Endocrinology 150, 376–384. doi: 10.1210/en.2008-0954
Jiang, Y., Xie, M., Chen, W., Talbot, R., Maddox, J. F., Faraut, T., et al. (2014). The sheep genome illuminates biology of the rumen and lipid metabolism. Science 344, 1168–1173. doi: 10.1126/science.1252806
Jin, H., Suh, D. S., Kim, T. H., Yeom, J. H., Lee, K., Bae, J. (2015). IER3 is a crucial mediator of TAp73β-induced apoptosis in cervical cancer and confers etoposide sensitivity. Sci. Rep. 5, 8367. doi: 10.1038/srep08367
Kim, D., Langmead, B., Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth3317
Ko, E. K., Chorich, L. P., Sullivan, M. E., Cameron, R. S., Layman, L. C. (2018). JAK/STAT signaling pathway gene expression is reduced following Nelf knockdown in GnRH neurons. Mol. Cell. Endocrinol. 470, 151–159. doi: 10.1016/j.mce.2017.10.009
Langmead, B., Trapnell, C., Pop, M., Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. doi: 10.1186/gb-2009-10-3-r25
Leung, K. C., Doyle, N., Ballesteros, M., Sjogren, K., Watts, C. K., Low, T. H., et al. (2003). Estrogen inhibits GH signaling by suppressing GH-induced JAK2 phosphorylation, an effect mediated by SOCS-2. Proc. Natl. Acad. Sci. U. S. A. 100, 1016–1021. doi: 10.1073/pnas.0337600100
Liu, Q., Pan, Z., Wang, X., Hu, W., Di, R., Yao, Y., et al. (2014). Progress on major genes for high fecundity in ewes. Front. Agric. Sci. Eng 1, 282–290. doi: 10.15302/J-FASE-2014042
Liu, Z. Z., Ji, Z. B., Wang, G. Z., Chao, T. L., Hou, L., Wang., J. M. (2016). Genome-wide analysis reveals signatures of selection for important traits in domestic sheep from different ecoregions. BMC Genomics 17, 863. doi: 10.1186/s12864-016-3212-2
Liu, Q., Hu, W., He, X., Pan, Z., Guo, X., Feng, T., et al. (2017). Establishment of high-throughput molecular detection methods for ovine high fecundity major gene FecB and their application. Acta Veterinaria Et Zootechnica Sin. 48, 39–51. doi: 10.11843/j.issn.0366-6694.2017.01.005
Liu, Y., Zhang, J., Xu, Q., Kang, X., Wang, K., Wu, K., et al. (2018). Integrated miRNA-mRNA analysis reveals regulatory pathways underlying the curly fleece trait in Chinese tan sheep. BMC Genomics 19, 360. doi: 10.1186/s12864-018-4736-4
Livak, K. J., Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔct method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Mulsant, P., Lecerf, F., Fabre, S., Schibler, L., Monget, P., Lanneluc, I., et al. (2001). Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola Mérino ewes. Proc. Natl. Acad. Sci. U. S. A. 98, 5104–5109. doi: 10.1073/pnas.091577598
Muroi, Y., Ishii, T. (2016). A novel neuropeptide Y neuronal pathway linking energy state and reproductive behavior. Neuropeptides 59, 1–8. doi: 10.1016/j.npep.2016.09.002
Nishihara, M., Kimura, F. (1989). Postsynaptic effects of prolactin and estrogen on arcuate neurons in rat hypothalamic slices. Neuroendocrinology 49, 215–218. doi: 10.1159/000125117
Perello, M., Stuart, R. C., Nillni, E. A. (2007). Differential effects of fasting and leptin on proopiomelanocortin peptides in the arcuate nucleus and in the nucleus of the solitary tract. Am. J. Physiol. Endocrinol. Metab. 292, 300–313. doi: 10.1152/ajpendo.004662006
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt3122
Pokharel, K., Peippo, J., Honkatukia, M., Seppälä, A., Rautiainen, J., Ghanem, N., et al. (2018). Integrated ovarian mRNA and miRNA transcriptome profiling characterizes the genetic basis of prolificacy traits in sheep (Ovis aries). BMC Genomics 19, 104. doi: 10.1186/s12864-017-4400-4
Polley, S., De, S., Brahma, B., Mukherjee, A., Vinesh, P. V., Batabyal, S., et al. (2010). Polymorphism of BMPR1B, BMP15 and GDF9 fecundity genes in prolific Garole sheep. Trop. Anim. Health Prod. 42, 985–993. doi: 10.1007/s11250-009-9518-1
Quintela, T., Gonçalves, I., Martinho, A., Alves, C. H., Saraiva, M. J., Rocha, P., et al. (2011). Progesterone enhances transthyretin expression in the rat choroid plexus in vitro and in vivo via progesterone receptor. J. Mol. Neurosci. 44, 152–158. doi: 10.1007/s12031-010-9398-x
Saedi, S., Khoradmehr, A., Mohammad Reza., J. S., Tamadon, A. (2018). The role of neuropeptides and neurotransmitters on kisspeptin/kiss1r-signaling in female reproduction. J. Chem. Neuroanat. 92, 71–82. doi: 10.1016/j.jchemneu.2018.07.001
Sangiaoalvarellos, S., Penabello, L., Manfredilozano, M., Tenasempere, M., Cordido, F. (2014). Perturbation of hypothalamic microRNA expression patterns in male rats after metabolic distress: impact of obesity and conditions of negative energy balance. Endocrinology 155, 1838–1850. doi: 10.1210/en.2013-1770
Schratt, G. M., Tuebing, F., Nigh, E. A., Kane, C. G., Sabatini, M. E., Kiebler, M., et al. (2006). A brain-specific microRNA regulates dendritic spine development. Nature 439, 283–289. doi: 10.1038/nature04367
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Sun, L., Wang, S., Hu, C., Zhang, X. (2011). Down-regulation of PKHD1 induces cell apoptosis through PI3K and NF-κB pathways. Exp. Cell. Res. 317, 932–940. doi: 10.1016/j.yexcr.2011.01.025
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt1621
Vieira., M., Gomes, J. R., Saraiva, M. J. (2015). Transthyretin induces insulin-like growth factor i nuclear translocation regulating its levels in the hippocampus. Mol. Neurobiol. 51, 1468–1479. doi: 10.1007/s12035-014-8824-4
Wang, L., Feng, Z., Wang, X., Wang, X., Zhang, X. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. doi: 10.1093/bioinformatics/btp612
Watanobe, H. (2010). Leptin directly acts within the hypothalamus to stimulate gonadotropin-releasing hormone secretion in vivo in rats. J. Physiol. 545, 255–268. doi: 10.1113/jphysiol.2002.023895
Wilson, J. L., Enriori, P. J. (2015). A talk between fat tissue, gut, pancreas and brain to control body weight. Mol. Cell. Endocrinol. 418, 108–119. doi: 10.1016/j.mce.2015.08.022
Yang., J., Li, X., Cao, Y. H., Pokharel, K., Hu, X. J., Chen, Z. H., et al. (2018). Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries) provides novel insights into the genetic mechanisms for female reproductive success. Heredity 122, 172–186. doi: 10.1038/s41437-018-0090-1
Zhai, M., Xie, Y., Liang, H., Lei, X., Zhao, Z. (2018). Comparative profiling of differentially expressed microRNAs in estrous ovaries of Kazakh sheep in different seasons. Gene 664, 181–191. doi: 10.1016/j.gene.2018.04.025
Zhang, Z. B., Tang, J. S., Di, R., Liu, Q. Y., Wang, X. Y., Gan, S. Q., et al. (2019a). Comparative transcriptomics reveal key sheep (Ovis aries) hypothalamus lncRNAs that affect reproduction. Animals 9, 15. doi: 10.3390/ani9040152
Zhang, Z. B., Tang, J. S., He, X. Y., Zhu, M. X., Gan, S. Q., Guo, X. F., et al. (2019b). Comparative transcriptomics identify key hypothalamic circular RNAs that participate in Sheep (Ovis aries) reproduction. Animals 9, 557. doi: 10.3390/ani9080557
Keywords: hypothalamus, mRNAs, miRNAs, GnRH, reproduction, sheep
Citation: Zhang Z, Tang J, Di R, Liu Q, Wang X, Gan S, Zhang X, Zhang J, Chu M and Hu W (2020) Integrated Hypothalamic Transcriptome Profiling Reveals the Reproductive Roles of mRNAs and miRNAs in Sheep. Front. Genet. 10:1296. doi: 10.3389/fgene.2019.01296
Received: 15 May 2019; Accepted: 25 November 2019;
Published: 15 January 2020.
Edited by:
Robert J. Schaefer, University of Minnesota Twin Cities, United StatesReviewed by:
Ikhide G. Imumorin, Georgia Institute of Technology, United StatesZhibin Ji, Shandong Agricultural University, China
Copyright © 2020 Zhang, Tang, Di, Liu, Wang, Gan, Zhang, Zhang, Chu and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mingxing Chu, mxchu@263.net; Wenping Hu, huwenping@caas.cn
†These authors have contributed equally to this work.