Insight Into Pituitary lncRNA and mRNA at Two Estrous Stages in Small Tail Han Sheep With Different FecB Genotypes

The pituitary is a remarkably dynamic organ with roles in hormone (FSH and LH) synthesis and secretion. In animals with the FecB (fecundity Booroola) mutation, the pituitary experiences hormone fluctuations during the follicular–luteal transition, which is implicated in the expression and regulation of many genes and regulators. Long non-coding RNAs (lncRNAs) are a novel type of regulatory factors for the reproductive process. Nevertheless, the expression patterns of lncRNAs and their roles in FecB-mediated follicular development and ovulation remain obscure. Thus, we profiled the pituitary transcriptome during the follicular (F, 45 h after evacuation vaginal sponges) and luteal (L, 216 h after evacuation vaginal sponges) phases in FecB-mutant homozygous (BB) and wild-type (WW) Small Tail Han sheep. We identified 78 differentially expressed genes (DEGs) and 41 differentially expressed lncRNAs (DELs) between BB_F and BB_L, 32 DEGs and 26 DELs between BB_F and WW_F, 16 DEGs and 29 DELs between BB_L and WW_L, and 50 DEGs and 18 DELs between WW_F and WW_L. The results of real-time quantitative PCR (RT-qPCR) correlated well with the transcriptome data. In both the follicular and luteal phases, DEGs (GRID2, glutamate ionotropic receptor delta type subunit 2; ST14, ST14 transmembrane serine protease matriptase) were enriched in hormone synthesis, secretion, and action. MSTRG.47470 and MSTRG.101530 were the trans-regulated elements of ID1 (inhibitor of DNA binding 3, HLH protein) and the DEG ID3 (inhibitor of DNA binding 3, HLH protein), and EEF2 (eukaryotic translation elongation factor 2), respectively; these factors might be involved in melatonin and peptide hormone secretion. In the FecB-mediated follicular phase, MSTRG.125392 targeted seizure-related 6 homolog like (SEZ6L), and MSTRG.125394 and MSTRG.83276 targeted the DEG KCNQ3 (potassium voltage-gated channel subfamily Q member 3) in cis, while MSTRG.55861 targeted FKBP4 (FKBP prolyl isomerase 4) in trans. In the FecB-mediated luteal phase, LOC105613905, MSTRG.81536, and MSTRG.150434 modulated TGFB1, SMAD3, OXT, respectively, in trans. We postulated that the FecB mutation in pituitary tissue elevated the expression of certain genes associated with pituitary development and hormone secretion. Furthermore, this study provides new insights into how the pituitary regulates follicular development and ovulation, illustrated by the effect of the FecB mutation.

Introduction

The FecB (fecundity Booroola) mutation was first detected in high-fecundity Australian Booroola Merino sheep (1). Mulsant et al. (2) Souza et al. (3), and Wilson et al. (4) revealed that the effect of FecB mutation was due to the substitution of a base (A to G at 746 bp, A746G) in the BMPR1B (bone morphogenetic protein receptor 1B) coding region on ovine chromosome 6, resulting in the corresponding substitution of one amino acid (glutamine to arginine, Q249R). Through an impairment of the activity of the BMPR1B receptor kinase, FecB mutation weakened the inhibitory effect of the BMPR1B gene on steroid synthesis in granulosa cells, leading to the maturation and ovulation of ovulatory follicles at significantly small diameters and fewer granulosa cells in homozygous FecB sheep than those in non-FecB sheep (2, 3, 5–7). In our previous study, litter size was significantly correlated with the FecB genotype in high-prolificacy Chinese Small Tail Han sheep (8). The gonadotropic hormones follicle-stimulating hormone (FSH) and luteinizing hormone (LH) are the dominant hormones in the coordinate regulation of follicle development and antral follicle ovulation (9). A previous study predicted that FecB might impact the pituitary gland through the regulation of hormone release during one estrus cycle in ewes (10). At a particular physiological time during the follicular phase, the plasma concentration of FSH in FecB homozygous ewes was higher than that in non-FecB ewes (11–13). The follicular–luteal transition is a complex and elaborate process that involves pituitary hormones and several cell-specific gene products that are under strict developmental regulation. The identification of the key regulators of these developmental processes could be beneficial for understanding the mechanisms of how FecB regulates the follicular–luteal transition at the molecular level and for screening candidate protein-coding genes with impacts on reproductive ability in ewes.

Previous transcriptomics studies have focused solely on protein-coding genes (14). Nevertheless, transcriptomic analyses and chromatin signatures showed a plethora of non-coding transcripts produced by eukaryotic cells (15, 16). Of these, long non-coding RNAs (lncRNAs) are a heterogeneous group of non-protein-coding transcripts of greater than 200 nucleotides in length with unknown coding potential. They include circular RNAs, long intergenic RNAs, pseudogenes, and antisense RNAs (17). Accumulating evidence has implicated lncRNAs in a wide array of cellular processes, such as transcriptional regulation, cellular pluripotency and reprogramming, and differentiation (reviewed in 18, 19). In terms of transcriptional regulation, lncRNAs act in cis or trans through the recruitment of RNAs and proteins. Cis-acting lncRNAs activate, repress, or modulate gene expression in a manner dependent on the location of their own transcription sites (20). Unfortunately, it has been difficult to accurately identify lncRNA sequences from transcribed sequences. Trans-acting lncRNAs affect nuclear structure and organization and interact with and modulate the activity of RNAs and proteins (21). The cytoplasmic localization of lncRNAs and the biological activities they participate in suggest that some lncRNAs can regulate the secretion of reproductive hormones.

Pituitary gonadotropins form the center of the endocrine axis involved in the reproductive process. At puberty, hypothalamic gonadotropin-releasing hormone (GnRH) stimulates the pituitary tissue to synthesize FSH and LH. Both the FSH and LH receptors are G protein-coupled receptors and have cascade amplification effects on the estrous or menstrual cycle. Gonadotropins are produced at demanded levels and specific times in response to multiple hormonal signals; thus, regulatory elements (lncRNAs) are expected to play key roles in both the stage specificity and hormonal responsiveness of their expression. Several lncRNAs have been proposed to affect the signal transduction of FSH and LH. In rats, the pituitary intergenic lncRNA m433s1 promoted the expression level of FSHB and the secretion of FSH by functioning as a competing endogenous RNA of miR-433 and reducing its inhibition of FSHB messenger RNA (mRNA) (22). In humans, silencing the pituitary lncRNA UCA1 significantly suppressed the secretion of prolactin secretion and the growth of pituitary cancer cells (23). In ram pituitary cells, short-interfering RNA knockdown of lncRNA-TCONS_00066406 reduced the mRNA levels of FSHB and LHB (24). FecB mutation has been reported to cause putative transcriptional changes in the pituitary. Pituitary MSTRG.259847.2 was upregulated in FecB homozygous ewes compared with FecB heterozygous ewes. A subsequently constructed lncRNA–mRNA interaction network highlighted Smad2, and experiments at the cell level showed that MSTRG.259847.2 can inhibit Smad2 and LHB (25). However, FecB mutation reduced the expression levels of some genes in the pituitary, including BMP2, TGFB2, TGFB3, ID1, and ID3 (26). Studies identifying lncRNAs associated with the follicular–luteal transition have been performed mainly on ovarian tissues, and investigation of lncRNAs at the pituitary level in ewes with different FecB status and during different stages of estrus has been rarely reported.

Small Tail Han sheep are an indigenous Chinese high-fecundity sheep breed with year-round estrus and a high ovulation rate (27). The estrus cycle consists of the follicular and luteal phases, accompanied by changes in the concentrations of hormones such as FSH and LH. Our previous study observed that FecB-mutant ewes regulate the expression levels of the BMPR1B, BMP15, and GDF9 genes, which are relevant to hormone and signal transduction (28)—that is, FecB mutation might affect the follicular–luteal transition to further influence sheep fertility. In this study, we collected pituitary tissues to analyze how the pituitary regulates follicular development, maturation, and ovulation in the context of FecB mutation. Our data provided a useful resource for further studying the functional roles of these genes and lncRNAs in the follicular development and ovulation of BB genotype ewes.

Materials and Methods

Ethical Consideration

All experimental ewes were housed and fed by the Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agriculture Sciences (IAS-CAAS), and ethical approval was received from the Animal Ethics Committee of the IAS-CAAS (no. IAS 2019-49).

Experimental Design and Pituitary Collection

Based on TaqMan genotyping (29), six FecB-mutant homozygous (BB, litter size = 2.2 ± 0.6) and six FecB wild-type (WW, litter size = 1.0 ± 0.0) Small Tail Han ewes were selected from a nucleus herd (Shandong, China). The selected ewes were approximately 3 years old and were similar in weight (Table 1). All non-pregnant animals had free access to food and water under natural temperature and lighting conditions. We divided genotype BB and WW ewes into a follicular phase group and a luteal phase group. Estrus synchronization was performed in September, within the usual breeding period (from August to November) described for this breed in Tianjin (117.2 E, 39.13 N). Firstly, ewes were treated with vaginal sponges (300 mg progesterone; InterAg Co., Ltd., Hamilton, New Zealand), followed by injection of vitamin AD to protect the endometrium. The vaginal sponges were evacuated after 12 days, and the removal time was set as 0 h. In accordance with previous reproductive trait analyses, six ewes (three BB and three WW ewes) were euthanized at 45 h (follicular phase) and six ewes (three BB and three WW ewes) at 216 h (luteal phase) (30, 31). Pituitary samples were dissected immediately after euthanasia, frozen in liquid nitrogen, and stored at −80°C for further analysis.

TABLE 1

Table 1 Litter size and body condition characteristics of genotype BB and WW ewes.

RNA Isolation, Library Construction, and Sequencing

Total RNA was extracted from the pituitary using the TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The purity, concentration, and integrity of RNA were checked with a Nano Photometer^® spectrophotometer (IMPLEN, Westlake Village, CA, USA), Qubit^® RNA Assay kits (Thermo Fisher Scientific, Waltham, MA, USA), and RNA Nano 6000 Assay, respectively. The RNA integrity number (RIN) value of all the samples was greater than 7.

Three micrograms of total RNA was used as the starting amount for the construction of an lncRNA library. Ribo-Zero™ Gold Kits (Epicenter, Madison, WI, USA) were applied to remove ribosomal RNA (rRNA) from the total RNA. Following the instructions for the NEB Next^® Ultra Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA), 12 RNA sequencing libraries were constructed for paired-end sequencing. Pooled libraries were sequenced on the Illumina HiSeq X platform (Illumina, San Diego, USA) using a chain-specific library construction strategy and the number and structure of the mRNA, known lncRNA, and novel lncRNA transcripts. All the experimental sequencing data were generated by Annoroad Gene Technology Co., Ltd. (Beijing, China).

Data Quality Control and Assembly

To ensure the quality of further analytical data, raw reads were filtered using in-house Perl scripts (Annoroad Gene Technology Co., Ltd, Beijing, China) to obtain high-quality reads. The filtering step involved deleting reads with low quality (Phred quality score <5%), adapter contamination, matches to rRNA sequence, or a rate of ambiguous bases higher than 5% (32). The Phred quality score refers to the rate of sequencing errors for a given base; for example, Q30 indicates that the base sequencing error rate was 0.1%. We aligned paired-end clean reads to the oar4.0 sheep reference genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_000298735.2) using HISAT2 (v.2.0.5) (33) with the parameters “–rna-strandness RF” and “–dta -t -p 4”. Only the uniquely mapped reads were assembled, and the expression levels were predicted using String Tie (v.1.3.2d) (33) with the parameter “-G ref.gtf -rf -1”. Thereafter, we calculated the number and ratio of the uniquely mapped reads within the three gene functional elements: exons, introns, and intergenic elements. Homogeneity analysis was subsequently performed to guarantee that the sequencing results did not impact further transcriptome analysis. Thus, mRNA and known lncRNA transcripts were identified.

Identification of Potential lncRNA Candidates

To decrease the false-positive rate, the screening of candidate lncRNAs is divided into two steps: one is based on the position of the coding reads [intronic lncRNA, long intergenic non-coding RNA (lincRNA), and antisense lncRNA] and the other is to count the encoding potential (34). The following were applied: 1) Transcripts with lengths longer than 200 bp, exon numbers greater than 2, and read coverage greater than 5 were retained. 2) The GffCompare program (v.0.10.1) was used to annotate the assembled transcripts and discard known mRNAs and other non-coding RNAs (rRNAs, tRNAs, snoRNAs, etc.) (35). 3) Transcripts with class_code “i” (potential intronic lncRNA), “x” (potential antisense lncRNA), “u” (potential intergenic lncRNA), “j” (potential novel isoform with more than one splice junction of a reference transcript), or “o” (genetic exonic overlapping with a reference transcript) were retained (36). 4) The Coding–Non-Coding Index (CNCI), Coding Potential Calculator (CPC), the Pfam database, and the Coding Potential Assessment Tool (CPAT) were utilized to predict the coding potential of preliminary candidate lncRNAs. The results predicted by the four aforementioned tools were intersected to obtain potential lncRNA candidates. The CNCI was determined by profiling adjoining nucleotide triplets to distinguish coding/non-coding transcripts independent of known annotations (37), with the parameter “–score 0 –length 199 –exon_num 2”. CPC can access the protein-coding potential of each transcript according to six biologically meaningful sequence features (38). Both CNCI and CPC are dependent on a score <0, thus defining lncRNA transcripts as non-coding RNAs. In Pfam, the profile hidden Markov models (HMMs) were searched for protein domains to obtain transcripts with known protein domains, based on the UniProt Knowledgebase (UniProtKB) (39). LncRNAs could be regarded as non-coding RNAs when the E-value ≥0.001. CPAT (v.1.2.1) works on a logistic regression model with four sequence features: Fickett TESTCODE statistic, open reading frame (ORF) size, ORF coverage, and hexamer usage bias (40).

Differential Expression Analysis of lncRNAs and mRNAs and Clustering

To verify the rationality of the samples and the reliability of the experiment, the correlation of the expression levels between samples was examined using Pearson’s correlation methods. Gene expression levels are generally measured on the basis of the amount of mRNA transcribed from a gene. Based on the reference GTF file and the HISAT BAM files, the HTSeq Python package (v.0.6.1) was used to calculate read counts, with the parameters “-I gene_id -f bam -s” and “reverse -a 10 -q”. Because three biological replicates were performed for each experimental group, the DEseq2 package (v.1.28.1) was used for the analysis of differential expression to calculate the log₂(fold change) and p-value based on the negative binomial distribution model (41). Meanwhile, log2(fold change) > 1 and p_adj < 0.05 were regarded as the essential criteria to identity the mRNAs and lncRNAs that were significantly differentially expressed in BB_F vs. BB_L (different estrus stages under genotype BB ewes, named gene set 1), BB_F vs. WW_F (different FecB genotypes ewes in the follicular phase, named gene set 2), BB_L vs. WW_L (different FecB genotypes ewes in the luteal phase, named gene set 3), and WW_F vs. WW_L (different estrus stages under genotype WW ewes, named gene set 4) individually. The intersection of BB_F vs. BB_L and WW_F vs. WW_L was named gene set 5; genes in this set were considered highly conserved during the follicular–luteal transition and to be potentially affected at the transcriptional levels by the FecB mutation. The intersection of BB_F vs. BB_L and BB_F vs. WW_F was named gene set 6; these genes were implicated in the process of follicular development and ovulation modulation by FecB. StringTie (v1.3.2d) was adopted to count the fragments per kilobase of transcripts per million mapped reads (FPKM). Subsequently, systematic clustering analysis based on the log₂(FPKM) value of each mRNA and lncRNA was performed with pheatmap (v.1.0.2) to analyze the similarities and relationships between the different libraries (42). The analyses included Pearson’s correlation and Euclidean distance methods.

Target Gene Prediction for Differentially Expressed lncRNAs and Network Construction

The relationships between lncRNAs and potential protein-coding genes were predicted on the basis of their distances and expression correlations, including the cis- and trans-acting relationships. Protein-coding genes adjacent to the lncRNAs (100 kb upstream and downstream) were identified as potential cis-target genes (43). Potential trans-target genes were predicted using Spearman’s correlation coefficients between the expressions of mRNAs and lncRNAs (r² ≥ 0.95) (43). According to the cis- or trans-acting relationship between the differentially expressed lncRNAs (DELs) and target genes, we visualized the lncRNA–mRNA networks with Cytoscape (v.3.8.0), which was run with layout = “attribute circle layout”.

Functional Annotation of Candidate Genes

To reveal the potential roles of DELs, we performed Gene Ontology (GO) category and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on their target genes using the clusterProfiler package (v.3.16.0) with the false discovery rate (FDR) pAdjustMethod. GO analysis (http://geneontology.org) provides structured and quantifiable knowledge concerning the functions of genes and gene products, and the molecular function ontology represents the activity of the gene products, especially focusing on transcription regulator activities (44). KEGG systematically analyzes gene functions, linking genomic information with functional information (45). Furthermore, KEGG analysis can predict the protein interaction networks involved in various cellular processes. A GO term or a KEGG pathway with corrected p < 0.05 was considered significantly enriched. We mainly focused on the pathways related to reproduction and pituitary function.

RT-qPCR Validation

Six mRNAs and six lncRNAs were randomly selected for validation via real-time quantitative PCR (RT-qPCR). β-actin was used as a control. The primers were designed with Primer Premier 5 and synthesized by Qine Ke Biotech (Beijing, China) (Supplementary Table S6). For RT-qPCR analysis, 1 μg of total RNA was reverse transcribed to complementary DNA (cDNA) following the instructions in the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Beijing, China). The genes and lncRNAs were diluted fourfold before use for RT-qPCR according to the instructions in TB Green^® Premix Ex Taq™ II (Takara, Beijing, China) and analyzed with QuantStudio^® 3 (ABI, Foster City, CA, USA). All genes and lncRNAs were amplified in triplicate. The data of each gene and lncRNA were analyzed with the 2^−ΔΔCt method (46). The correlation between the RT-qPCR data and the FKPM value from the RNA sequencing (RNA-seq) data was calculated using Pearson’s correlation method in R (v.4.0.2).

Results

Overview of the Sequencing Data in Small Tail Han Sheep Pituitary Tissue

The raw RNA-seq data from the Small Tail Han sheep were analyzed for quality control before aligning reads to the reference genome. More than 97.50 Gb of clean read sequences were generated from each library, with a corresponding rate of clean Q30 bases higher than 91.24%, indicating that the filtered sequence was appropriate for mapping. All the high-quality clean reads were mapped to the oar4.0 sheep reference genome using HISAT2. The percentage of the total mapped and uniquely mapped reads in these samples reached 92.66% and 88.27%, respectively. The rate of multi-mapping reads from the 12 libraries was less than 4.62% (Supplementary Table S1).

Identification of lncRNAs and mRNAs in Small Tail Han Sheep Pituitary Tissue

According to their locations, the unique reads from 12 samples were categorized into three types: intergenic, intronic, and exon (Figure 1A). Approximately 15.17% of the lncRNAs were distributed in intergenic regions, 42.85% were introns overlapping lncRNAs, and 41.98% were exons overlapping lncRNAs. After screening preliminary candidate lncRNAs, the intersecting results of the CNCI, CPC, Pfam, and CPAT tools were used to obtain potential lncRNA candidates, shown in Figure 1B. Furthermore, we identified the essential features of the lncRNAs and contrasted them with those of mRNAs. The genome distribution of most mRNAs and lncRNAs was randomly assigned to the 26 autosomes and the X chromosome; 480 mRNAs (2.3%) and 431 lncRNAs (3.8%) not aligned with any chromosomal location, and 10 mRNAs and 1 lncRNA were found to align to a mitochondrial location (Figure 1C). The length distribution of the lncRNAs was mainly appropriately 2,900–3,000 bp, similar to that of the mRNA transcripts in pituitary tissue (Figure 1D). Most lncRNAs possessed two exons, significantly fewer than the average number of mRNA transcripts (Figure 1E). The expression abundance of lncRNAs was lower than that of mRNAs in the 12 samples (Figure 1F).

FIGURE 1

Figure 1 Identification of long non-coding RNAs (lncRNAs) and messenger RNAs (mRNAs) in Small Tail Han sheep pituitary. (A) Classification of the uniquely mapped read locations, including exon, intron, and intergenic regions. (B) Venn diagram analysis showing the number of common and unique novel lncRNAs identified by four methods: Coding–Non-Coding Index (CNCI), Coding Potential Calculator (CPC), the Pfam database, and the Coding Potential Assessment Tool (CPAT). (C) Distribution of lncRNAs and mRNAs on chromosomes. (D) Lengths of lncRNAs and mRNAs. (E) Exon numbers of lncRNAs and mRNAs. (F) FPKM value of lncRNAs and mRNAs.

Profiling of DE lncRNA and DEGs in Small Tail Han Sheep Pituitary

The correlation coefficient of the gene expression levels between samples was higher than 0.80, indicating that the sample selection was consistent and reliable. Under the criteria of log2(fold change) > 1 and p_adj < 0.05, 78 differentially expressed genes (DEGs; 37 upregulated and 41 downregulated) and 41 DELs (26 upregulated and 15 downregulated) were screened between BB_F and BB_L (Figure 2A and Supplementary Table S2), 32 DEGs (20 upregulated and 12 downregulated) and 26 DELs (15 upregulated and 11 downregulated) were screened between BB_F and WW_F (Figure 2B and Supplementary Table S2), 16 DEGs (9 upregulated and 7 downregulated) and 29 DELs (9 upregulated and 20 downregulated) were screened between BB_L and WW_L (Figure 2C and Supplementary Table S2), and 50 DEGs (19 upregulated and 31 downregulated) and 18 DELs (9 upregulated and 9 downregulated) were screened between WW_F and WW_L (Figure 2D and Supplementary Table S2). The mRNA and lncRNA expression patterns were similar in both the follicular and luteal phases. Nevertheless, a different expression profile was observed at the follicular–luteal transition, which may indicate that changes in transcript expression accompanied FecB mutation (Figures 2E, F). Furthermore, we compared the quantity and distribution of the DEGs and DELs in each of the four pairwise comparisons. The DEGs and DELs from the intersection of BB_F vs. BB_L and WW_F vs. WW_L were involved in follicular–luteal transition and conserved across the different genotypes, including GRID2, LOC105610376, MSTRG.163804, and MSTRG.81824. The intersection of BB_F vs. BB_L and BB_F vs. WW_F contained vital genes and lncRNAs, including LRRC7, INSM2, SEZ6L, BMX, KCNQ3, GRM8, LDB3, CDH12, LOC101109936, MSTRG.125394, MSTRG.125388, MSTRG.125392, MSTRG.83276, and MSTRG.81824, which are involved in the process by which FecB mutations affect follicular development. The intersection of BB_F vs. BB_L and BB_L vs. WW_L included LOC105610376, MSTRG.26398, MSTRG.127921, and MSTRG.134082, which are involved in ovulation in FecB mutant ewes (Figure 2G).

FIGURE 2

Figure 2 Analysis of differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs). (A–D) Volcano plots showing the upregulated and downregulated genes in BB_F vs. BB_L (A), in BB_F vs. WW_F (B), in BB_L vs. WW_L (C), and in WW_F vs. WW_L (D). (E) Hierarchical clusters of DEGs, where all the fragments per kilobase of transcripts per million mapped reads (FPKM) values of genes were normalized by log₁₀(FPKM). (F) Hierarchical clusters of DELs, where all the FPKM values of known and novel lncRNAs were normalized by log₁₀(FPKM). (E) Venn diagram showing the analysis of unique and shared DEGs and DELs between the BB and WW groups and between the follicular and luteal groups. (G) Venn diagram showing the analysis of unique and shared DEGs and DELs between the BB and WW groups and between the follicular and luteal groups. F, follicular phase; L, luteal phase, BB, FecB-mutant homozygous genotype; WW, FecB wild-type genotype.

GO and Pathway Enrichment Analysis of DEGs

The described methods were employed to predict the functions of the DEGs; that is, the enriched GO functional terms and KEGG pathways were considered to be the predicted functional terms and pathways for the protein-coding gene and lncRNA (47). In the comparison of the follicular and luteal phases, GO analysis of the DEGs showed that prolactin secretion, regulation of mitotic cell cycle, and amino acid metabolism were among the top 50 enriched terms (Supplementary Table S3). In addition, 125 pathways were annotated to hormone synthesis, secretion, and action in follicular–luteal transition, including steroid hormone, growth hormone, thyroid hormone, estrogen, and insulin (Figure 3A and Supplementary Table S3).

FIGURE 3

Figure 3 Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to pituitary function and reproductive process of the differentially expressed genes (DEGs) and the target genes of differentially expressed lncRNAs (DELs). (A, C) KEGG enrichment pathways for DEGs (A) and the target genes of DELs (C) in the follicular phase and the luteal phase. (B, D) KEGG enrichment pathways for DEGs (B) and the target genes of DELs (D) in the BB (FecB-mutant homozygous) and WW (FecB wild type) genotypes.

DEGs between the different FecB genotypes were enriched in pituitary hormone signal transduction processes (in the top 50 enriched terms), such as neurotransmitter metabolic process, choline metabolic process, and regulation of synapse process (Supplementary Table S3). KEGG pathway analysis identified 65 enriched pathways, including neuroactive ligand–receptor interaction, oxidative phosphorylation, dopaminergic synapse, melanogenesis, gap function, Wnt, mTOR, MAPK, PI3K–Akt, and Hippo signaling pathways (Figure 3B and Supplementary Table S3).

Screening of Potential Functional LncRNAs Involved in Small Tail Han Sheep Reproduction

To further reveal the potential biological function of pituitary lncRNAs in Small Tail Han sheep, an interaction network of lncRNAs and their corresponding target genes was constructed. In the comparison BB_F vs. BB_L, 12 known lncRNAs corresponded to 86 target genes, and 28 novel lncRNAs corresponded to 189 target genes (Figure 4A and Supplementary Table S5). In the comparison WW_F vs. WW_L, two known lncRNAs corresponded to eight target genes, and 15 novel lncRNAs corresponded to 75 target genes (Figure 4B and Supplementary Table S5). In the comparison of the follicular and luteal phases, the top 50 enriched terms of the target genes included tight junction, steroid hormone binding, and transcription elongation factor complex (Supplementary Table S4). In addition, 144 pathways were related to circadian entrainment, oocyte meiosis and maturation, and hormone signaling pathways, including thyroid hormone, oxytocin, and growth hormone (Figure 3C and Supplementary Table S4). These follicular–luteal transition DELs might play similar roles in pituitary hormone synthesis and secretion. In Table 2, four DEGs and six genes targeted by DELs are shown in detail. The targets of lncRNAs MSTRG.125392, MSTRG.125394, and MSTRG.83276 were cis-regulated, while others were trans-regulated. The expression patterns of most lncRNAs were similar to those of their corresponding target genes during the follicular–luteal transition.

FIGURE 4

Figure 4 Interactions of the differentially expressed lncRNAs (DELs) with target genes form a network. (A) In the comparison BB_F vs. BB_L, 40 DELs acting in cis or in trans with 275 messenger RNAs (mRNAs) formed the interactive network. Green and red colors denote upregulated and downregulated long non-coding RNAs (lncRNAs), respectively. Quadrilaterals and boxes denote lncRNAs and target genes, respectively. Dotted and straight lines represent cis-acting and trans-acting, respectively. (B) In the comparison WW_F vs. WW_L, 17 DELs acting in cis or in trans with 83 mRNAs formed the interactive network. (C) In the comparison BB_F vs. WW_F, 23 DELs acting in cis or in trans with 117 mRNAs formed the interactive network. (D) In the comparison BB_L vs. WW_L, 27 DELs acting in cis or in trans with 151 mRNAs formed the interactive network. F, follicular phase; L, luteal phase; BB, FecB-mutant homozygous genotype; WW, FecB wild-type genotype.

TABLE 2

Table 2 Details of the differentially expressed genes (DEGs) and genes targeted by differentially expressed lncRNAs (DELs).

In the comparison BB_F vs. WW_F, two known lncRNAs corresponded to four target genes, and 21 novel lncRNAs corresponded to 113 target genes (Figure 4C and Supplementary Table S5). In the comparison BB_L vs. WW_L, six known lncRNAs corresponded to 19 target genes, and 21 novel lncRNAs corresponded to 132 target genes (Figure 4D and Supplementary Table S5). In the comparison of the FecB genotypes, adenosine receptor binding, postsynaptic cytosol, kinase activity, and neurohypophyseal hormone activity were enriched in the top 50 terms (Supplementary Table S4). In addition, 208 pathways were identified, which included tight junction, TGF-beta, Hippo, Wnt, circadian rhythm, melanogenesis, MAPK, and PI3K–Akt signaling pathways (Figure 3D and Supplementary Table S4). Enrichment analyses of the DELs and DEGs showed similarities, indicating that DELs might be involved in the process by which FecB mutation had an effect on pituitary function. Two DEGs and five genes targeted by DELs are shown in detail (Table 2).

Validation of RNA Sequencing Using RT-qPCR

To validate the RNA-seq data, eight genes and eight lncRNAs were examined with RT-qPCR, and the relative gene expression was calculated on the basis of the 2^−ΔΔCt values (Figure 5). The results indicated that the data on the mRNA and lncRNA expression levels were reliable. Meanwhile, the results of Pearson’s correlation analysis of all genes showed a strong positive correlation between the RT-qPCR and RNA-seq data (r² > 0.93, p < 0.05).

FIGURE 5

Figure 5 Validation of the RNA sequencing (RNA-seq) data by real-time quantitative PCR (RT-qPCR). The RT-qPCR data are presented as relative gene expression. RNA-seq data are presented as fragments per kilobase of transcripts per million mapped reads (FPKM). (A–D) Selected genes and long non-coding RNAs (lncRNAs) from BF and BL (A, B) and from WF and WL (C, D) were validated by RT-qPCR and RNA-seq, respectively. (E–H) Selected genes and lncRNAs from BF and WF (E, F) and from BL and WL (G, H) were validated by RT-qPCR and RNA-seq, respectively. BF, FecB-mutant homozygous genotype ewes in the follicular phase; BL, FecB-mutant homozygous genotype ewes in the luteal phase; WF, FecB wild-type genotype ewes in the follicular phase; WL, FecB wild-type genotype ewes in the luteal phase.

Discussion

Pituitary gonadotropin is essential for the activation of reproduction through the synthesis and secretion of gonadotropic hormones (FSH and LH) at the required levels. Given the unique function of the pituitary, i.e., the production of FSH and LH at specific times and appropriate levels, and its response to diverse hormone signals, it is expected that some regulatory elements play dominant roles in this cell-specific expression and hormonal responsiveness. In our study, 8,562 novel lncRNAs identified from the sheep pituitary shared many characteristics with those from other mammals, such as goats, mice, and pigs. In these mammals, most lncRNAs possessed two exons (followed by three exons), and lncRNAs were shorter and lower than mRNAs and were expressed at lower levels, reinforcing the legitimacy of the transcripts identified here (48, 49). Fewer lncRNAs were identified here than in the hypothalamus and pituitary gland of sheep, which had 34,293 and 19,672 lncRNAs, respectively (25, 50), likely reflecting tissue-specific expression patterns and characteristics and the different sampling times employed. The most lncRNAs in the present study were detected on ovine chromosomes 1–3, which was perhaps unsurprising as these chromosomes are the largest in the ovine genome. This implies a relationship between these chromosomes and pituitary function in the context of the FecB mutant. Similar chromosome locations of lncRNAs were found in previous studies on Hu sheep (chromosomes 1, 2, and 3) (25). In other words, the total distribution of lncRNAs was consistent between Hu sheep and Small Tail Han sheep, indicating that the functions of these lncRNAs might be conserved. Meanwhile, 20,248 mRNAs were identified, including OXTR, which is associated with the neurohypophysis, and hormone genes (FSHB, LHB, PRL, and INHA) and receptors (FSHR, ESR1, and ESR2) associated with the adenohypophysis (51, 52).

To understand the regulatory mechanism of the follicular–luteal transition, a previous study focused mainly on the microRNA (miRNA)–mRNA interactions in ovarian tissue and revealed that some miRNAs, such as miR-200a, miR-200b, and miR-200c, appeared to be critical regulators of ovarian follicular and luteal development (53, 54). However, pituitary hormones (FSH and LH) are required for the transition from the follicle to the corpus luteum, which involves receiving feedback signals from gonadal steroids and activin/inhibin signaling pathways (55, 56). In our study, DEGs from the follicular–luteal transition of wild-type ewes were mainly enriched in the G2/M transition of mitotic cell cycle, response to BMP, and Wnt signaling pathway of the GO terms and tight junction, TGF-beta signaling pathway, circadian entrainment, gap junction, and hormone-related signaling pathways of the KEGG pathways. IGFBP3, PLCE1, and GRM1 were annotated to hormone-related signaling pathways, which might belong to the adenohypophysis. Moreover, GRM1 was annotated as a gap junction gene; gap junctions mediate intercellular communication or cytosolic connections to regulate the amount of gonadotropin produced (57). With the FecB mutation, the number of DEGs increased and were annotated to amino acid metabolic process, prolactin secretion, and melanin biosynthetic process of the GO terms and the neuroactive ligand–receptor interaction, amino acid metabolism, hormone biosynthesis, and synapse-related pathways of the KEGG pathways. RLN3, GRM8, GRID2, and VIP were enriched in neuroactive ligand–receptor interaction, and GRM8, KCNQ3, and MAOB were enriched in synapse-related pathways, suggesting that these genes might be part of the neurohypophysis. Among these DEGs, GRID2 and ST14 were downregulated in the follicular phase compared with the luteal phase. A previous study revealed that GRID2 is the receptor for glutamate that controls anterior pituitary hormone secretion and reproduction (58, 59); that is, GRID2 may promote LH secretion to be involved in ovulation like glutamate receptor AMPA 1 (60). ST14 reduced cell growth by downregulating the cell cycle-related proteins, and ST14 was detected in the pituitary, indicating that it played an important role in pituitary development (61). Meanwhile, we identified some RAB proteins, among which RAB1A and RAB26 had higher relative expressions. RAB proteins are master regulators of vesicle transport, which is an essential process for secretion (62). Reductions in RAB1B, RAB3B, RAB6, and RAB11 may disrupt pituitary secretion of FSH and LH (63). The current results suggested that RAB1A and RAB26 might be new candidates involved in the secretory pathway in the pituitary gland. Moreover, some regulatory elements are expected to play dominant roles in follicular–luteal transition. The lncRNA MSTRG.81824 was upregulated in the comparison of WW_F with WW_L, but downregulated in the comparison of BB_F with BB_L, indicating that FecB could decrease the expression level of MSTRG.81824 in the follicular phase. MSTRG.81824 may have a trans effect on CDK2 (cyclin dependent kinase 2) and affect follicle maturation (53, 64, 65). Other lncRNAs (MSTRG.47470 and MSTRG.101530) worked in trans when they affected the target genes ID1 and ID3 (a DEG), and EEF2, respectively. ID1 and ID3, and EEF2 in the pituitary induced melatonin and peptide hormone secretion, indicating that MSTRG.47470 and MSTRG.101530 are associated with hormone secretion (66, 67).

Mulsant et al. reported that Small Tail Han sheep were characterized by a high ovulation rate and litter size due to the action of the major gene FecB (2). It was suggested that, in ewes with the BB genotype, the BMPR1B mutation would advance granulosa cell differentiation and ovulatory follicle maturation, resulting in differential gene expression in the pituitary gland (68). In our study, the BMPR1B gene did not reach a significant difference in the different genotypes. However, estrogen functions in ovarian tissue and is involved in the process by which it downregulates the production and secretion of pituitary gonadotropin via negative feedback to the pituitary (69). Comparison of BB_F vs. WW_F showed that the FecB mutation significantly increased the expression of the WNT11 (Wnt family member 11) gene in the follicular phase. The expression levels of the three DELs (MSTRG.125392, MSTRG.125394, and MSTRG.83276) and the two target genes (SEZ6L and KCNQ3) in the pituitary of BB genotype ewes were significantly higher than those of WW genotype ewes, and the expression trends of lncRNAs and their corresponding cis-target genes were consistent. Moreover, cis-regulatory lncRNAs showed enhancer-like activity and promoted the expression of neighboring genes. The DEG SEZ6L was cis-regulated by MSTRG.125392 and MSTRG.125394 and has been reported to influence dendritic morphology and synapse numbers (70). Previous conservation genomic analysis identified SEZ6L as a candidate gene associated with reproduction and system development in Black Slavonian pigs through the modulation of key oncogenic pathways such as Wnt (71, 72). Meanwhile, the DEG KCNQ3, a cis-regulatory element of MSTRG.83276, was important in the regulation of NPY/AgRP neuron excitability, which promoted hormone communication to maintain homeostasis. These results indicated that these three DELs might regulate the transcription of the target genes by promoting the enhancer and the promoter of the transcriptional machinery molecule, similar to the manner in which classical lncRNAs (HOTTIP and HOTAIRM1) promoted the expression of the adjacent target gene HoxA by cis regulation (73). FKBP4, an element that is trans-regulated by MSTRG.55861, is relevant to steroid hormone receptor binding and transportation and to reproductive and neurological diseases (74). Our study showed that FKBP4 was somewhat differentially expressed and might induce steroid hormone synthesis in pituitary tissue. As an important paralog of this gene, FKBP5 (FKBP prolyl isomerase 5) might produce more FSH in the hypothalamus tissue (75); that is, FKBP4 and FKBP5 are closely tied to hormone synthesis and secretion in the hypothalamic–pituitary–ovarian (HPO) axis, and FecB affects their expression trend. In the luteal phase in different FecB genotypes, the target genes TGFB1 (transforming growth factor beta 1), SMAD3 (SMAD family member 3), and OXT (oxytocin/neurophysin I prepropeptide) were trans-elements of the lncRNAs LOC105613905, MSTRG.81536, and MSTRG.150434, respectively. It was suggested that, with FecB mutation, lncRNAs in the follicular phase might be involved in pituitary development and steroid hormone communication, while lncRNAs in the luteal phase probably participate in the TGFβ/SMAD signaling pathway and in oxytocin secretion.

Overall, our results suggested that some genes were differentially expressed in the pituitary tissue between the follicular and luteal phases via the regulation of reproductive hormone synthesis and secretion. During the follicular phase, the expression of some genes in the pituitary tissue was elevated in the FecB mutant in comparison with the wild-type genotype. Whether the change in the gene expression levels affected the corresponding reproductive pathways and processes remains to be further elucidated.

Conclusions

This pituitary transcriptome analysis showed the expression profiles of the pituitary tissue with the FecB mutation in the follicular–luteal transition. The comparison of the follicular phase and the luteal phase suggested that lncRNAs might be implicated in hormone secretion. In ewes with the FecB mutation, lncRNAs in the follicular phase (MSTRG.125392, MSTRG.125394, MSTRG.83276, and MSTRG.55861) might be associated with pituitary development and steroid hormone communication, while those in the luteal phase (LOC105613905, MSTRG.81536, and MSTRG.150434) probably participated in the TGFβ/SMAD signaling pathway and in oxytocin secretion. Moreover, the data obtained in this study should be useful for improving breed conservation and for better exploiting the genetic resources of sheep.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, accession ID: PRJNA782215.

Ethics Statement

All experimental ewes were supported by the Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agriculture Sciences (IAS-CAAS). Ethical approval was in compliance with the Animal Ethics Committee of the IAS-CAAS (no. IAS 2019-49).

Author Contributions

XG, XW, and MC conceived and designed the experiments. XG and XW performed the experiments. SC analyzed the data and wrote the manuscript, with input from the other authors. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by The National Natural Science Foundation of China (32172704 and 31902150), the Earmarked Fund for China Agriculture Research System of MOF and MARA (CARS-38), the Central Public-Interest Scientific Institution Basal Research Fund (Y2017JC24 and 2018ywf-yb-2), the Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13), the Natural Science Foundation of Tianjin (20JCQNJC00630), the Natural Science Foundation of Jilin Province (20210101376JC), and China Postdoctoral Science Foundation (2021M703202).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2021.789564/full#supplementary-material

Supplementary Table S1 | Alignment of statistical results of reads.

Supplementary Table S2 | The results of differentially expressed mRNAs and lncRNAs.

Supplementary Table S3 | DEGs GO and KEGG enrichment analysis.

Supplementary Table S4 | Target genes of DELs GO and KEGG enrichment analysis.

Supplementary Table S5 | The lncRNAs-mRNAs network

Supplementary Table S6 | Primer sequence and product size of mRNAs and lncRNAs for RT-qPCR.

Supplementary Figure S1 | The Pearson correlation between samples.

References

1. Davis G, Montgomery G, Allison A, Kelly R, Bray A. Segregation of a Major Gene Influencing Fecundity in Progeny of Booroola Sheep. N Z J Agric Res (1982) 25:525–9. doi: 10.1080/00288233.1982.10425216

CrossRef Full Text | Google Scholar

2. Mulsant P, Lecerf F, Fabre S, Schibler L, Monget P, Lanneluc I, et al. Mutation in Bone Morphogenetic Protein Receptor-IB Is Associated With Increased Ovulation Rate in Booroola Merino Ewes. Proc Natl Acad Sci (2001) 98:5104–9. doi: 10.1073/pnas.091577598

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Souza C, MacDougall C, Campbell B, McNeilly A, Baird D. The Booroola (FecB) Phenotype Is Associated With a Mutation in the Bone Morphogenetic Receptor Type 1 B (BMPR1B) Gene. J Endocrinol (2001) 169:R1. doi: 10.1677/joe.0.169r001

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Wilson T, Wu X-Y, Juengel JL, Ross IK, Lumsden JM, Lord EA, et al. Highly Prolific Booroola Sheep Have a Mutation in the Intracellular Kinase Domain of Bone Morphogenetic Protein IB Receptor (ALK-6) That Is Expressed in Both Oocytes and Granulosa Cells. Biol Reprod (2001) 64:1225–35. doi: 10.1095/biolreprod64.4.1225

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Driancourt M, Cahill L, Bindon B. Ovarian Follicular Populations and Preovulatory Enlargement in Booroola and Control Merino Ewes. Reproduction (1985) 73:93–107. doi: 10.1530/jrf.0.0730093

CrossRef Full Text | Google Scholar

6. Kumar S, Rajput PK, Bahire SV, Jyotsana B, Kumar V, Kumar D. Differential Expression of BMP/SMAD Signaling and Ovarian-Associated Genes in the Granulosa Cells of FecB Introgressed GMM Sheep. Syst Biol Reprod Med (2020) 66:185–201. doi: 10.1080/19396368.2019.1695977

PubMed Abstract | CrossRef Full Text | Google Scholar

7. McNatty K, Henderson K, Lun S, Heath D, Ball K, Hudson N, et al. Ovarian Activity in Booroola× Romney Ewes Which Have a Major Gene Influencing Their Ovulation Rate. Reproduction (1985) 73:109–20. doi: 10.1530/jrf.0.0730109

CrossRef Full Text | Google Scholar

8. Chu M, Liu Z, Jiao C, He Y, Fang L, Ye S, et al. Mutations in BMPR-IB and BMP-15 Genes Are Associated With Litter Size in Small Tailed Han Sheep (Ovis Aries). J Anim Sci (2007) 85:598–603. doi: 10.2527/jas.2006-324

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Richards JS, Pangas SA. The Ovary: Basic Biology and Clinical Implications. J Clin Invest (2010) 120:963–72. doi: 10.1172/JCI41350

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Campbell B, Marsters P, Baird D, Walkdenbrown S, Werf J, Nimbkar C, et al. The Mechanism of Action of the FecB (Booroola) Mutation. Use FecB Gene Sheep-breeding Program (2009) 46:46–56.

Google Scholar

11. Isaacs K, McNatty K, Condell L, Shaw L, Heath D, Hudson N, et al. Plasma FSH, LH and Immunoreactive Inhibin Concentrations in FecBB/FecBB and FecB+/FecB+ Booroola Ewes and Rams From Birth to 12 Months of Age. Reproduction (1995) 103:89–97. doi: 10.1530/jrf.0.1030089

CrossRef Full Text | Google Scholar

12. McNatty K, Hudson N, Lun S, Heath D, Shaw L, Condell L, et al. Gonadotrophin-Releasing Hormone and the Control of Ovulation Rate by the FecB Gene in Booroola Ewes. Reproduction (1993) 98:97–105. doi: 10.1530/jrf.0.0980097

CrossRef Full Text | Google Scholar

13. Phillips D, Hudson N, McNatty K. Effects of Ovariectomy and Genotype on Bioactive FSH in Plasma and Pituitary of Booroola Ewes. Reproduction (1993) 98:559–65. doi: 10.1530/jrf.0.0980559

CrossRef Full Text | Google Scholar

14. Laiho A, Kotaja N, Gyenesei A, Sironen A. Transcriptome Profiling of the Murine Testis During the First Wave of Spermatogenesis. PloS One (2013) 8:e61558. doi: 10.1371/journal.pone.0061558

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, et al. Chromatin Signature Reveals Over a Thousand Highly Conserved Large Non-Coding RNAs in Mammals. Nature (2009) 458:223–7. doi: 10.1038/nature07672

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Mercer TR, Gerhardt DJ, Dinger ME, Crawford J, Trapnell C, Jeddeloh JA, et al. Targeted RNA Sequencing Reveals the Deep Complexity of the Human Transcriptome. Nat Biotechnol (2012) 30:99–104. doi: 10.1038/nbt.2024

CrossRef Full Text | Google Scholar

17. Chan JJ, Tay Y. Noncoding RNA: RNA Regulatory Networks in Cancer. Int J Mol Sci (2018) 19:1310. doi: 10.3390/ijms19051310

CrossRef Full Text | Google Scholar

18. Samata M, Akhtar A. Dosage Compensation of the X Chromosome: A Complex Epigenetic Assignment Involving Chromatin Regulators and Long Noncoding RNAs. Annu Rev Biochem (2018) 87:323–50. doi: 10.1146/annurev-biochem-062917-011816

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Sherstyuk VV, Medvedev SP, Zakian SM. Noncoding RNAs in the Regulation of Pluripotency and Reprogramming. Stem Cell Rev Rep (2018) 14:58–70. doi: 10.1007/s12015-017-9782-9

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Gil N, Ulitsky I. Regulation of Gene Expression by Cis-Acting Long Non-Coding RNAs. Nat Rev Genet (2020) 21:102–17. doi: 10.1038/s41576-019-0184-5

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kopp F, Mendell JT. Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell (2018) 172:393–407. doi: 10.1016/j.cell.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Han D-X, Sun X-L, Wang C-J, Yu Z-W, Zheng Y, Huang Y-J, et al. Differentially Expressed lncRNA-M433s1 Regulates FSH Secretion by Functioning as a miRNA Sponge in Male Rat Anterior Pituitary Cells. Biol Reprod (2019) 101:416–25. doi: 10.1093/biolre/ioz100

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Liu G, Wang L, Li Y. Inhibition of lncRNA-UCA1 Suppresses Pituitary Cancer Cell Growth and Prolactin (PRL) Secretion via Attenuating Glycolysis Pathway. Vitro Cell Dev Biol Anim (2020) 56:642–9. doi: 10.1007/s11626-020-00494-x

CrossRef Full Text | Google Scholar

24. Yang H, Ma J, Wang Z, Yao X, Zhao J, Zhao X, et al. Genome-Wide Analysis and Function Prediction of Long Noncoding RNAs in Sheep Pituitary Gland Associated With Sexual Maturation. Genes (2020) 11:320. doi: 10.3390/genes11030320

CrossRef Full Text | Google Scholar

25. Zheng J, Wang Z, Yang H, Yao X, Yang P, Ren C, et al. Pituitary Transcriptomic Study Reveals the Differential Regulation of lncRNAs and mRNAs Related to Prolificacy in Different FecB Genotyping Sheep. Genes (2019) 10:157. doi: 10.3390/genes10020157

CrossRef Full Text | Google Scholar

26. Tian D, Liu S, Tian F, Ding N, Li X, Zhao K. Comparative Transcriptome of Reproductive Axis in Chinese Indigenous Sheep With Different FecB Genotypes and Prolificacies. Anim Reprod Sci (2020) 223:106624. doi: 10.1016/j.anireprosci.2020.106624

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Di R, He J, Song S, Tian D, Liu Q, Liang X, et al. Characterization and Comparative Profiling of Ovarian microRNAs During Ovine Anestrus and the Breeding Season. BMC Genomics (2014) 15:1–15. doi: 10.1186/1471-2164-15-899

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tang J, Hu W, Di R, Liu Q, Wang X, Zhang X, et al. Expression Analysis of the Prolific Candidate Genes, BMPR1B, BMP15, and GDF9 in Small Tail Han Ewes With Three Fecundity (FecB Gene) Genotypes. Animals (2018) 8:166. doi: 10.3390/ani8100166

CrossRef Full Text | Google Scholar

29. Liu Q, Hu W, He X, Pan Z, Guo X, Feng T, et al. Establishment of High-Throughput Molecular Detection Methods for Ovine High Fecundity Major Gene FecB and Their Application. Acta Vet Et Zootechnica Sin (2017) 48:39–51.

Google Scholar

30. Guo X, Wang X, Di R, Liu Q, Hu W, He X, et al. Metabolic Effects of FecB Gene on Follicular Fluid and Ovarian Vein Serum in Sheep (Ovis Aries). Int J Mol Sci (2018) 19:539. doi: 10.3390/ijms19020539

CrossRef Full Text | Google Scholar

31. Wang X, Guo X, He X, Liu Q, Di R, Hu W, et al. Effect of FecB Mutation on Estrus, Ovulation and Endocrine Characteristics in Small Tail Han Sheep. Front Vet Sci (2021) 1305. doi: 10.3389/fvets.2021.709737

CrossRef Full Text | Google Scholar

32. Zhu W, Tian L, Yue X, Liu J, Fu Y, Yan Y. LncRNA Expression Profiling of Ischemic Stroke During the Transition From the Acute to Subacute Stage. Front Neurol (2019) 10:36. doi: 10.3389/fneur.2019.00036

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-Level Expression Analysis of RNA-Seq Experiments With HISAT, StringTie and Ballgown. Nat Protoc (2016) 11:1650. doi: 10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ran M, Chen B, Li Z, Wu M, Liu X, He C, et al. Systematic Identification of Long Noncoding RNAs in Immature and Mature Porcine Testes. Biol Reprod 94 (2016) 77:1–9. doi: 10.1095/biolreprod.115.136911

CrossRef Full Text | Google Scholar

35. Pertea G, Pertea M. GFF Utilities: GffRead and GffCompare. FResearch (2020) 9:304. doi: 10.12688/f1000research.23297.2

CrossRef Full Text | Google Scholar

36. Ouyang Q, Hu S, Wang G, Hu J, Zhang J, Li L, et al. Comparative Transcriptome Analysis Suggests Key Roles for 5-Hydroxytryptamlne Receptors in Control of Goose Egg Production. Genes (2020) 11:455. doi: 10.3390/genes11040455

CrossRef Full Text | Google Scholar

37. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing Sequence Intrinsic Composition to Classify Protein-Coding and Long Non-Coding Transcripts. Nucleic Acids Res (2013) 41:e166–6. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, et al. CPC: Assess the Protein-Coding Potential of Transcripts Using Sequence Features and Support Vector Machine. Nucleic Acids Res (2007) 35:W345–9. doi: 10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, et al. Pfam: The Protein Families Database. Nucleic Acids Res (2014) 42:D222–30. doi: 10.1093/nar/gkt1223

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang L, Park HJ, Dasari S, Wang S, Kocher J-P, Li W. CPAT: Coding-Potential Assessment Tool Using an Alignment-Free Logistic Regression Model. Nucleic Acids Res (2013) 41:e74–4. doi: 10.1093/nar/gkt006

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Love MI, Huber W, Anders S. Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Kolde R, Kolde MR. Package ‘Pheatmap’. R Package (2015) 1:790.

Google Scholar

43. Xiang S, Li Z, Weng X. The Role of lncRNA RP11-154D6 in Steroid-Induced Osteonecrosis of the Femoral Head Through BMSC Regulation. J Cell Biochem (2019) 120:18435–45. doi: 10.1002/jcb.29161

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Consortium GO. The Gene Ontology Resource: 20 Years and Still GOing Strong. Nucleic Acids Res (2019) 47:D330–8. doi: 10.1093/nar/gky1055

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Arocho A, Chen B, Ladanyi M, Pan Q. Validation of the 2-ΔΔct Calculation as an Alternate Method of Data Analysis for Quantitative PCR of BCR-ABL P210 Transcripts. Diagn Mol Pathol (2006) 15:56–61. doi: 10.1097/00019606-200603000-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Guttman M, Rinn JL. Modular Regulatory Principles of Large Non-Coding RNAs. Nature (2012) 482:339–46. doi: 10.1038/nature10887

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gao X, Ye J, Yang C, Luo L, Liu Y, Ding J, et al. RNA-Seq Analysis of lncRNA-Controlled Developmental Gene Expression During Puberty in Goat & Rat. BMC Genet (2018) 19:1–9. doi: 10.1186/s12863-018-0608-9

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Gong Y, Zhang Y, Li B, Xiao Y, Zeng Q, Xu K, et al. Insight Into Liver lncRNA and mRNA Profiling at Four Developmental Stages in Ningxiang Pig. Biology (2021) 10:310. doi: 10.3390/biology10040310

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Zhang Z, Tang J, Di R, Liu Q, Wang X, Gan S, et al. Comparative Transcriptomics Reveal Key Sheep (Ovis Aries) Hypothalamus lncRNAs That Affect Reproduction. Animals (2019) 9:152. doi: 10.3390/ani9040152

CrossRef Full Text | Google Scholar

51. Fletcher PA, Smiljanic K, Maso Prévide R, Iben JR, Li T, Rokic MB, et al. Cell Type-and Sex-Dependent Transcriptome Profiles of Rat Anterior Pituitary Cells. Front Endocrinol (2019) 10:623. doi: 10.3389/fendo.2019.00623

CrossRef Full Text | Google Scholar

52. St-Amand J, Yoshioka M, Tanaka K, Nishida Y. Transcriptome-Wide Identification of Preferentially Expressed Genes in the Hypothalamus and Pituitary Gland. Front Endocrinol (2012) 2:111. doi: 10.3389/fendo.2011.00111

CrossRef Full Text | Google Scholar

53. Donadeu FX, Schauer S, Sontakke S. Involvement of miRNAs in Ovarian Follicular and Luteal Development. J Endocrinol (2012) 215:323. doi: 10.1530/JOE-12-0252

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Wang H, Li X, Zhou R, Xi J, Wei Y, Li L, et al. Genome-Wide Transcriptome Profiling in Ovaries of Small-Tail Han Sheep During the Follicular and Luteal Phases of the Oestrous Cycle. Anim Reprod Sci (2018) 197:212–21. doi: 10.1016/j.anireprosci.2018.08.031

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Refael T, Melamed P. Enhancing Gonadotrope Gene Expression Through Regulatory lncRNAs. Endocrinology (2021) 162(8):1–16. doi: 10.1210/endocr/bqab116

CrossRef Full Text | Google Scholar

56. Shalev D, Melamed P. The Role of the Hypothalamus and Pituitary Epigenomes in Central Activation of the Reproductive Axis at Puberty. Mol Cell Endocrinol (2020) 518:111031. doi: 10.1016/j.mce.2020.111031

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Edwards BS, Clay CM, Ellsworth BS, Navratil AM. Functional Role of Gonadotrope Plasticity and Network Organization. Front Endocrinol (2017) 8:223. doi: 10.3389/fendo.2017.00223

CrossRef Full Text | Google Scholar

58. Elegheert J, Kakegawa W, Clay JE, Shanks NF, Behiels E, Matsuda K, et al. Structural Basis for Integration of GluD Receptors Within Synaptic Organizer Complexes. Science (2016) 353:295–9. doi: 10.1126/science.aae0104

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Todman M, Han S-K, Herbison A. Profiling Neurotransmitter Receptor Expression in Mouse Gonadotropin-Releasing Hormone Neurons Using Green Fluorescent Protein-Promoter Transgenics and Microarrays. Neuroscience (2005) 132:703–12. doi: 10.1016/j.neuroscience.2005.01.035

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Sugimoto M, Sasaki S, Watanabe T, Nishimura S, Ideta A, Yamazaki M, et al. Ionotropic Glutamate Receptor AMPA 1 Is Associated With Ovulation Rate. PloS One (2010) 5:e13817. doi: 10.1371/journal.pone.0013817

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Wang Y, Rathinam R, Walch A, Alahari SK. ST14 (Suppression of Tumorigenicity 14) Gene Is a Target for miR-27b, and the Inhibitory Effect of ST14 on Cell Growth Is Independent of miR-27b Regulation. J Biol Chem (2009) 284:23094–106. doi: 10.1074/jbc.M109.012617

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Takai Y, Sasaki T, Matozaki T. Small GTP-Binding Proteins. Physiol Rev (2001) 81:153–208. doi: 10.1152/physrev.2001.81.1.153

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Ren J-C, Zhu Q, LaPaglia N, Emanuele NV, Emanuele MA. Ethanol-Induced Alterations in Rab Proteins: Possible Implications for Pituitary Dysfunction. Alcohol (2005) 35:103–12. doi: 10.1016/j.alcohol.2005.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Shabbir S, Boruah P, Xie L, Kulyar M-e, Nawaz M, Yousuf S, et al. Genome-Wide Transcriptome Profiling Uncovers Differential miRNAs and lncRNAs in Ovaries of Hu Sheep at Different Developmental Stages. Sci Rep (2021) 11:1–12. doi: 10.1038/s41598-021-85245-y

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Sugiura K, Naito K, Tojo H. Cdk2 Activity Is Essential for the First to Second Meiosis Transition in Porcine Oocytes. J Reprod Dev Biol (2005) 51:143–9. doi: 10.1262/jrd.51.143

CrossRef Full Text | Google Scholar

66. Argüelles S, Muñoz MF, Cano M, Machado A, Ayala A. In Vitro and In Vivo Protection by Melatonin Against the Decline of Elongation Factor-2 Caused by Lipid Peroxidation: Preservation of Protein Synthesis. J Pineal Res (2012) 53:1–10. doi: 10.1111/j.1600-079X.2011.00961.x

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Konishi H, Ogawa T, Nakagomi S, Inoue K, Tohyama M, Kiyama H. Id1, Id2 and Id3 Are Induced in Rat Melanotrophs of the Pituitary Gland by Dopamine Suppression Under Continuous Stress. Neuroscience (2010) 169:1527–34. doi: 10.1016/j.neuroscience.2010.06.030

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Regan S, McFarlane JR, O’Shea T, Andronicos N, Arfuso F, Dharmarajan A, et al. Flow Cytometric Analysis of FSHR, BMRR1B, LHR and Apoptosis in Granulosa Cells and Ovulation Rate in Merino Sheep. Reproduction (2015) 150:151–63. doi: 10.1530/REP-14-0581

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Goodman R, Gibson M, Skinner D, Lehman M. Neuroendocrine Control of Pulsatile GnRH Secretion During the Ovarian Cycle: Evidence From the Ewe. Reprod Suppl (2002) 59:41–56.

PubMed Abstract | Google Scholar

70. Qiu WQ, Luo S, Ma SA, Saminathan P, Li H, Gunnersen JM, et al. The Sez6 Family Inhibits Complement at the Level of the C3 Convertase. bioRxiv (2020). doi: 10.1101/2020.09.11.292623

CrossRef Full Text | Google Scholar

71. Lukić B, Ferenčaković M, Šalamon D, Čačić M, Orehovački V, Iacolina L, et al. Conservation Genomic Analysis of the Croatian Indigenous Black Slavonian and Turopolje Pig Breeds. Front Genet (2020) 11:261. doi: 10.3389/fgene.2020.00261

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Martel-Martel A, Gutiérrez-Cerrajero C, Vallejo-Fuente S, Casado-Garcia A, Postigo S, Sanchez Tapia E, et al. Study of SEZ6L Molecular Alterations in Solid Tumors. Eur J Hum Genet (2019) 27(Suppl 2):1568–9. doi: 10.1038/s41431-019-0494-2

CrossRef Full Text | Google Scholar

73. Shi T, Guo D, Xu H, Su G, Chen J, Zhao Z, et al. HOTAIRM1, an Enhancer lncRNA, Promotes Glioma Proliferation by Regulating Long-Range Chromatin Interactions Within HOXA Cluster Genes. Mol Biol Rep (2020) 47(4):1–11. doi: 10.1007/s11033-020-05371-0

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Hausch F. FKBPs and Their Role in Neuronal Signaling. Biochim Biophys Acta -General Subj (2015) 1850:2035–40. doi: 10.1016/j.bbagen.2015.01.012

CrossRef Full Text | Google Scholar

75. Chen S, Guo X, He X, Di R, Zhang X, Zhang J, et al. Transcriptome Analysis Reveals Differentially Expressed Genes and Long Non-Coding RNAs Associated With Fecundity in Sheep Hypothalamus With Different FecB Genotypes. Front Cell Dev Biol (2021) 9:633747. doi: 10.3389/fcell.2021.633747

CrossRef Full Text | Google Scholar