Identification and Regulatory Network Analysis of Genes Related to Reproductive Performance in the Hypothalamus and Pituitary of Angus Cattle

In this study, we explored the gene expression patterns of the pituitary gland and hypothalamus of Angus cows at different growth and developmental stages by deep sequencing and we identified genes that affect bovine reproductive performance to provide new ideas for improving bovine fertility in production practice. We selected three 6-month-old (weaning period), three 18-month-old (first mating period), and three 30-month-old (early postpartum) Angus cattle. The physiological status of the cows in each group was the same, and their body conformations were similar. After quality control of the sequencing, the transcriptome analyses of 18 samples yielded 129.18 GB of clean data. We detected 13,280 and 13,318 expressed genes in the pituitary gland and hypothalamus, respectively, and screened 35 and 50 differentially expressed genes (DEGs) for each, respectively. The differentially expressed genes in both tissues were mainly engaged in metabolism, lipid synthesis, and immune-related pathways in the 18-month-old cows as compared with the 6-month-old cows. The 30-month-old cows presented more regulated reproductive behavior, and pituitary CAMK4 was the main factor regulating the reproductive behavior during this period via the pathways for calcium signaling, longevity, oxytocin, and aldosterone synthesis and secretion. A variant calling analysis also was performed. The SNP inversions and conversions in each sample were counted according to the different base substitution methods. In all samples, most base substitutions were represented by substitutions between bases A and G, and the probability of base conversion exceeded 70%, far exceeding the transversion. Heterozygous SNP sites exceeded 37.68%.


Introduction
As the demand for beef increases, beef cattle are distributed worldwide. Angus cattle are a superior breed of beef cattle owing to their excellent meat performance and they are considered to be a typical specialized breed of beef cattle worldwide. In addition to their delicious meat, Angus cattle exhibit good lactation ability. Because of the significant advantages of Angus cattle, research on this breed continues. Angus cattle growth and development has been well studied [1]. Angus cattle are normally weaned at 6 months old and gradually reach sexual maturity at 12-14 months old. Primary mating can occur at 18 months of age, and the primiparity usually occurs at 30 months of age. Factors affecting the pregnancy rates and lactation of dairy cows and the relationship between the two have been discussed, and various methods for improving reproductive performance have been proposed [2,3]. However, for Angus cattle, the genes which cause the differences in reproductive performance at three different growth and development stages are unknown.
The pituitary-hypothalamus-gonadal (HPG) axis stems from interaction between nerve and hormone signals from three main sources. Through the release of gonadotropinreleasing hormone (GnRH) and gonadotropin-inhibitory hormone (GnIH), the hypothalamus acts on the pituitary gland to regulate gonadotropins and secrete the pituitary hormone, luteinizing hormone (LH). prolactin (PRL) and follicle-stimulating hormone (FSH) are involved in oocyte maturation, corpus luteum formation, and ovulation. These hormones act on the gonads and regulate steroid hormone production. Interactions among these three signals ensures that the HPG axis regulates reproductive activity. Ovulation, fertilization, and lactation are the keys to successful reproduction, and sexual maturity of the gonad organs helps ensure these behaviors. However, in dairy cattle, substantial genetic antagonism may exist between high yield and fertility in dairy cattle, in other words, high milk production is often accompanied by low fertility [4] Although prostaglandin and GnRH are combined early to induce ovulation to promote pregnancy and alleviate this dilemma, the mechanism behind this remains unclear. Comparing the different sexual development stages of dairy cows may provide clues to understanding this phenomenon. Early research on dairy cows showed that body quality scores during calving and early lactation are related to milk quantity. Few studies on beef cattle have been published in this area, and insufficient evidence exists to show a connection; however, studies have found that milk from beef cattle influenced these mechanisms differently from that of other cows [5].
With the development of high-throughput sequencing technology, transcriptomics research has become increasingly widely used to systematically study genes involved in biological processes and developmental stages, as well as responses to biological changes [6]. The establishment and continuous improvement of annotation information of the bovine genome [7] has progressed the use of transcriptomics to study cattle-related issues. Nextgeneration RNA sequencing technology applied to cattle will affect their short-term dietary control, immune responses, neuroprotection, energy homeostasis, and eventually, their estrus cycles. In addition to exploring the external factors that affect bovine estrus, sequencing technology is being used to determine the differences in HPG axis levels before and after estrus. Conjoint analysis of five reproductive tissues (i.e., the hypothalamus, pituitary, ovary, uterus, and endometrium) and some growth-and metabolism-related tissues with genome-wide association studies, transcriptomes, and bovine transcription factors have enabled the study of co-expression networks around puberty [8].
Beef cattle breeding is an important part of the animal husbandry industry, which has a significant impact on the development of the agriculture and animal husbandry economies and the improvement of national diet structure. Beef is favored by consumers because of its high lean meat rate and low fat. Angus cattle are one of the three high quality beef cattle breeds, with good meat production performance, high meat rate of quality characteristics, and is a relatively common breed of beef cattle breeding. We collected the hypothalamus and pituitary tissues from three stages of sexual development and explored the regulatory mechanisms guiding sexual maturation in angus cattle gonadal development via high-throughput RNA sequencing technology to further resolve existing problems in production practice.

Animal Tissue Collection
Nine Angus cattle at the ages of 6, 18, and 30 months from the Changsha Institute of Animal Science were selected and divided into 3 groups according to age, with 3 cows per group. Six-month-old cows are freshly weaned calves, 18-month-old cows are sexually mature cows that meet the requirements for priming, and 30-month-old cows have completed producing their first progeny. These cattle were maintained in similar housing and feeding conditions, the selected cattle were healthy and disease free, and their body weights and physiological states were similar in each group. External observation and rectal ovarian and uterine evaluation before slaughter confirmed that all cattle were not in heat. After evaluation, the cattle were stunned with an electric shock, and then euthanized by carotid bloodletting. After the cattle was completely dead, the upper hypothalamus and pituitary glands of the cattle were collected, frozen in liquid nitrogen until RNA isolation, and numbered according to age and tissue type. Supplementary Table S1 shows the numbering and grouping information. All animal protocols were carried out according to the guidance of the Animal Care Committee at Jilin University. All ethical assessments were carried out in accordance with the requirements of China's Laboratory Animal Management Regulations.

RNA Extraction and cDNA Library Preparation
The collected tissue samples were processed by grinding in liquid nitrogen, and the total RNA was extracted using the TRIzol protocol (Invitrogen, Carlsbad, CA, USA). The RNA concentration was measured using an ultramicro spectrophotometer (Tuohe, Shanghai, China) and the RNA quality was evaluated using the RNA integrity number (RIN) value from an Experion automated electrophoresis system (BioRad, Hercules, CA, USA). To make sure there was enough mRNA to build the cDNA library, the total RNA of each sample needed to be ≥4.0 µg, and the RIN values needed to be ≥8.0. The PolyATract mRNA separation system (Promega, Madison, WI, USA) was used to further isolate the poly-a-containing mRNA from each sample of total RNA. The mRNA was fragmented, then the reverse transcriptase and random primers were used to convert RNA fragments into first-strand cDNA. DNA polymerase I and RNase H were used to synthesize the second-strand cDNA, The cDNA sticky ends were repaired into flat ends using T4 DNA polymerase and Klenow DNA polymerase, and the cDNA 3 ends were added base A, then, the adapters were ligated. Next, the product was purified using an AMPure XP system, and PCR amplification yielded the final cDNA library. An Aligent Technologies 2100 bioanalyzer (Agilent, Palo Alto, CA, USA) was used to assess the quality of the cDNA library.

Quality Control
Raw data of fastq format were firstly processed through in-house perl scripts. In this step, clean reads were obtained by removing reads containing adapter, ploy-N, and low-quality reads from the raw data. At the same time, Q20, Q30, GC content, and sequence duplication level of the clean data were calculated. The quality score was integer mapping of the probability of base identification errors, and the Phred base quality value formula was commonly used: $$Q = −10 × log {10}P$$, P was the probability of base identification error, and the corresponding relationship between base mass value and the probability of identification error was: Q20 = 1/100, Q30 = 1/1000. All the downstream analyses were based on clean data with high quality (Q30 > 90%, Supplementary Table S2). Then, these clean reads were mapped to the reference genome sequence (Bos taurus reference genome, version UMD 3.1.1); only reads with a perfect match or one mismatch were further analyzed and annotated. The Hisat2 tools software was used to map with reference genome [9].

Variable Splicing Event Prediction
As described above, cleaned reads from each sample were mapped to the reference genome with HISAT2, StringTie was used for transcript assembly, the variable splicing types and corresponding expression levels of each sample were obtained using ASprofile software (Supplementary Figure S1) [19]. Variable shear was mainly based on the alternative 5 first exon (TSS) and alternative 3 last exon (TTS), and a few alternative exons ends or skipped exons were present. Such changes may cause premature termination of the protein coding, thereby acting as a molecular switch (Supplementary Table S3).

Differential Expression Analysis
For the samples with biological replicates, a differential expression analysis of two conditions/groups was performed by DESeq2 [20]. DESeq2 provided statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p-value <0.05 found by DESeq2 were assigned as differentially expressed. The gene ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented using the GOseq R packages based on Wallenius noncentral hypergeometric distribution [17], which could adjust for gene length bias in the DEGs. Gene numbers were calculated for every term, significantly enriched GO terms in the DEGs as compared with the genome background were defined using a hypergeometric test. The calculating formula of p-value is: where N is the number of all genes with GO annotation, n is the number of DEGs in N, M is the number of all genes that are annotated to the certain GO terms, and m is the number of DEGs in M. N stands for total background gene (TB gene number), n stands for total significant gene (TS gene number), M stands for background gene (B gene number), and m stands for significant gene (S gene number). GO terms meeting this condition with p < 0.05 were defined as significantly enriched GO terms in the DEGs. This analysis was able to recognize the main biological functions that the DEGs exercise.

KEGG Pathway Enrichment Analysis
The KOBAS software was used to test the statistical enrichment of differentially expression genes in KEGG pathways [16,21].

PPI (Protein-Protein Interaction)
The sequences of the DEGs were blast (blastx) to the Bos taurus reference genome (version UMD 3.1.1). The PPI information of the DEGs was obtained from the STRING database (http://stringdb.org/ (accessed on 14 March 2021)). Then, the PPI of these DEGs were visualized in Cytoscape [21,22].

Real-Time RT-PCR and Statistical Analysis
The RNA extraction method as described above, and a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Tokyo, Japan) was used for reverse transcription. RT-qPCR was implemented according to the SYBR Remix Ex TaqII kit. The reaction was performed on a real-time fluorescence quantitative PCR instrument (Mx3005P, Agilent, Santa Clara, CA, USA). The RT-qPCR conditions were as follows: 95 • C for 5 min, (95 • C for 15 s, 58 • C for 20 s, and 72 • C for 20 s) × 40 cycles, 72 • C for 5 min, finally get the melting curve, 95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s [23]. Total 20 genes were validated, and the genes are listed in Table 1. The results were analyzed by the 2 −∆∆Ct method. The primer sequences of these genes are shown in the Supplementary Table S4; PPP1R11 was the internal reference control gene. All primers were synthesized by the Comate Bioscience Co. Ltd. (Changchun, China). The SPSS version 22.0 (SPSS Inc., Chicago, IL, USA) was used for statistical analysis. The unpaired t-test was used to evaluate the significance of the differences. Statistical significance was indicated at * p < 0.05, ** p < 0.01, and *** p < 0.01. All data were expressed as mean ± SEM.

Overview of Reads
In the transcriptome analysis of 18 samples, an average of 6.30 GB clean reads were obtained, the GC content was >47.61%, and the ratio of the base percentage of Q30 in each sample was >94.19%. The transcriptome data and reference genome sequence comparison showed that the alignment efficiency of the reads and reference genomes of each sample ranged from 95.06% to 96.53% (Supplementary Table S2), indicating that the selected reference genome assembly met the information analysis needs. The mapped reads in the different regions of the specified reference genome (i.e., exons, introns, and intergenic regions) were counted, and reads from mature mRNAs were compared with exon regions.
In this study, 65% of the reads aligned to the exon region, 18.4% of the reads were aligned to the intron region, and because of the imperfect genome annotation, 16.6% of the reads aligned to the intergenic region (Supplementary Figure S2).

SNP/Indel (Insertion and Deletion) Analysis
The GATK software was used to identify base mismatches, insertions, and deletions between the sequenced samples and the reference genome, and then to analyze whether these SNPs affected gene expression levels or protein product types. The SnpEff software was used to annotate and predict the variation effects. Supplementary Tables S2 and S3 list the SNP/Indel site information.
The SNP inversions and conversions in each sample were counted according to the different base substitution methods. In all samples, most base substitutions were represented by substitutions between bases A and G, and the probability of base conversion exceeded 70%, far exceeding the transversion. Heterozygous SNP sites exceeded 37.68% (Supplementary Table S5).
Most SNP densities were concentrated at 0-2 base substitutions per 1000 bp sequence. As the SNP density increases, the number of genes decreases (Supplementary Figure S3a). The SNP/Indel annotation analysis revealed synonymous mutations in each sample, which may be related to altered protein functions (Supplementary Figure S3b,c).

Screening of Gene Expression in the Pituitary Gland and Hypothalamus
Expression level FPKM values ranged from six orders of magnitude 10 −2 to 10 4 , and the logarithm was mainly concentrated below 1 (Supplementary Figure S4a). To eliminate the effect of heterogeneity in co-expression, we used box plots to more intuitively show the FPKM of each sample (Supplementary Figure S4b). Repeated correlation evaluations were performed on all samples to exclude the biological variability on the results. Pearson's correlation coefficient (r) was used as the evaluation index of biological repetitive correlations. An r2 closer to 1 means a stronger correlation, as shown in Supplementary Figure S2c; a gradient from purple to blue represents a correlation coefficient from zero to one. The follow-up analysis was performed without abnormal samples (Supplementary Figure S4c). As listed in Supplementary Table S6, DEseq2 was performed on the hypothalamus and pituitary samples from different growth stages to obtain the differentially expressed gene sets between the nine groups of biological conditions.

New Gene Annotation Information in the Pituitary Gland and Hypothalamus
Based on the selected reference genome sequence, the StringTie software was used to splice and quantify mapped reads, which were compared with the original genome annotation information to find the original unannotated transcription regions, discover new transcripts and new genes of the species, and filter out peptide-coding chains that were too short (<50 amino acid residues) or contained only a single exon sequence. A total of 5999 new genes were discovered. A total of 4072 new genes were annotated; 4053 new genes received annotations in NR, 2933 new genes received annotations in GO, 2585 new genes received annotations in eggnog, 2012 new genes received annotations in KEGG (Supplementary Table S7).

Functional Identification of Differentially Expressed Genes in the Pituitary Gland and Hypothalamus
We performed GO and a KEGG enrichment analysis in the pituitary and hypothalamus tissues. In the three groups of pituitary tissues, six differentially expressed genes significantly enriched biological processes, cellular components, and molecular function. The differentially expressed genes were mainly concentrated in cellular processes (GO:0009987), biological regulation (GO:0065007), membrane portion (GO:0044425), macromolecular complexes (GO:0032991), binding (GO:0005488), and transporter activity (GO:0005215) ( Figure 3A). Next, we investigated the KEGG pathways in which all differentially expressed genes were significantly enriched in the pituitary tissues. Thirty enriched terms were involved in environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems ( Figure 3B). Figure 3C lists the top 20 pathways, including phototransduction(map04744), NOD-like receptor signaling (map04621), inflammatory bowel disease (map05321), and progesterone-mediated oocyte maturation (map04914).  Figure 3D). Next, we investigated the KEGG pathways in which all differentially expressed genes were significantly enriched in the hypothalamus. Thirty-seven enriched terms were involved in environmental information processing, human diseases, metabolism, and organismal systems ( Figure 3E). Figure 3F lists the top 20 pathways, including inositol phosphate metabolism, TGF-β signaling, biosynthesis of unsaturated fatty acids, ovarian steroidogenesis, and GnRH signaling pathways.
In G3 vs. G1, DEGs were mainly enriched in gene information editing, metabolism, environmental information processing, and immune pathways ( Figure 4A,B). In the hypothalamus, G2 vs. G4 involves differential gene functions mainly involved in metabolism, neuroimmunity, environmental information processing, biological systems (i.e., circadian rhythm, PPAR signaling, phosphatidylinositol signaling system, ovarian steroid production, aldosterone synthesis and secretion, GABA synapses, and prolactin signaling) ( Figure 4C,D). In G5 vs. G1, differential gene functions were focused on biological system functions (i.e., phototransduction, aldosterone synthesis and secretion, longevity regulation, oxytocin signaling, and taste transduction), environmental information processing (i.e., the calcium signaling and cAMP signaling pathways), and immune disease pathways ( Figure 4E,F). In G6 vs. G2, the hypothalamic differential gene functions were concentrated on biological system functions (GnRH signaling pathway and taste transduction), immune diseases, and environmental information processing (the calcium signaling and phospholipase D signaling pathways) ( Figure 4G,H).

The Major Genes of the Pituitary Gland and Hypothalamus
According to the KEGG pathway enrichment results, five genes related to reproduction were screened out and listed in Table 2. In pituitary tissues, as compared with G1, CAMK4 was significantly downregulated in G5 which enriched the oxytocin signaling pathway. As compared with G3, HSP90AB1 was significantly upregulated in G5 which enriched progesterone-mediated oocyte maturation and the estrogen signaling pathway. In hypothalamus tissues, as compared with G2, TH, HSPB7, and LOC511936 were significantly upregulated in G4, TH was involved in prolactin pathway, HSPB7 and LOC511936 were engaged in steroid hormones synthesis pathway. As compared with G2, PTK2B was significantly downregulated in G6 which enriched the GnRH signaling pathway. Table 2. Differentially expressed genes related to reproduction (G1, G3, and G5 = 6-month-old, 18-month-old, and 30-month-old Angus cattle pituitary; G2, G4, and G6 = 6-month-old, 18-month-old, and 30-month-old Angus cattle hypothalamus). Other differentially expressed genes are listed in Tables 3 and 4. In pituitary tissue, as compared with G1, TBX2 and SHANK2 were significantly downregulated in G3, PPP1R14A and HAPLN2 were significantly upregulated in G3. As compared with G1, TEX15 and CAMK4 were significantly downregulated in G5, ANGPTL7 and ORM1 were significantly upregulated in G5. As compared with G3, HSP90AB1 was significantly upregulated in G5 (Table 3).

Group
In the hypothalamus, as listed in Table 4, as compared with G2, Rorb, SYNJ2, and TBR1 were downregulated in G4, and PCP4, KCNK3, and KCNK1 were significantly upregulated in G4. As compared with G2, only AGO3 was downregulated in G6, and CHRNA2, ORM1, and RLF were significantly upregulated in G6. As compared wit G4, CENPA and ZNF282 were downregulated in G6.

The Results of RT-qPCR
As shown in Figure 5, real-time quantitative PCR was used to verify the previous transcriptome sequencing results, and a total of 20 genes were detected. In pituitary tissues, as compared with G1, TBX21 and SHANK2 were significantly downregulated in G3, PPP1R14A and HAPLN2 were significantly upregulated in G3, TEX5 and CAMK4 were significantly downregulated in G5, and ANGPTL7 and ORM1 were significantly upregulated in G5. As compared with G3, HSP90AB1 was significantly upregulated in G5. In the hypothalamus, as compared with G2, RORB and TBR1 were significantly downregulated in G4, HSPB7 and NOG were significantly upregulated in G4, PTK2B and AGO3 were significantly downregulated in G6, and ORM1 and RLF were significantly upregulated in G4. As compared with G4, CENPA, LOC506989, and ZNF282 were significantly downregulated in G6. We compared RT-qPCR with the NGS results. As shown in Supplementary Figure S5, these genes showed the same change trend in RT-qPCR and NGS.

Protein Interaction Network Analysis
Fifty-six and 61 nodes and 1540 and 974 edges were included in the protein-protein interaction network of the pituitary gland and hypothalamus based on significant DEGs (minimum required interaction score >0.7). Analyzing the degree of centrality of the protein-protein interaction network showed that all genes were hub genes (degree = 55), and 40 hub genes (degree = 43) were obtained. The PRL family occupied a central position in the pituitary tissue, and the gene function of this family was mainly related to lysosomes. In the hypothalamus, the MRP family was in the core regulatory position, participating in the activation of the TGF-β signaling pathway and lysosome formation (Supplementary Figure S6).

Discussion
Previous studies have shown that at least 22,000 genes are present in the cattle genome, and 14,345 core group orthologs are shared among seven mammal species, but genes associated with lactation and immunity vary widely among species [24]. In our study, without calculating newly annotated genes, we detected 13,280 and 13,318 genes in the pituitary and hypothalamus tissues, therefore, genes expressed in the pituitary and hypothalamus accounted for 57.9% (13,280/22,000) and 58.1% (14,345/22,000) of the total number of genes in cattle.
The comparison of 6-month-old and 18-month-old Angus cattle showed that pituitary and hypothalamic tissue metabolism was active. This may be because calf growth and development in this period inevitably require a large energy supply. In the hypothalamus, TH enriched in the prolactin pathway. In this pathway, changes of TH directly affect tyrosine synthesis. TH is highly co-expressed with KiSS-1 neurons in the anterior ventricleperiventricular zone (AVPV/PeN). The AVPV/PeN TH neurons have dopaminergic activity. Dopamine signal transduction due to AVPV/PeN kisspeptin neurons may directly regulate GnRH secretion, thereby, regulating the neuroendocrine reproductive axis. The evidence suggests that some AVPV TH neurons project to PVN and activate prolactin neurons, thereby, increasing prolactin secretion [25].
In G6 vs. G2, PTK2b was significantly downregulated by Ca 2+ in the calcium signaling pathway. In zebrafish models, PTK2b plays an important role in fertilization of zebrafish mother cells, and PTK2 is activated by fertilization-induced calcium [26]. PTK2 is also important for fertilization in mammals [27]. PTK2B is a downstream factor in this pathway, and its changes cause changes in other signal transductions, such as in the longevity regulation, long-term growth, and GnRH signaling pathways. PTK2b belongs to the protein tyrosine kinase family and plays an important role in various cellular processes.
In the pituitary, the role of CAMK4 in the calcium signaling pathway determines that changes in this pathway will directly affect fertilization, proliferation, learning, and memory. In the present study, CAMK4 was downregulated in G5 vs. G1 which enriched the oxytocin pathway. CAMK4 is an upstream regulator of nitric oxide synthase (NOS). Some studies have revealed the effect of CAMK4 on regulating NOS in neurological diseases [28], but the mechanism of its reproductive regulation remains unclear. CAMK4 downregulation directly caused PGC-1α downregulation, resulting in oxidative stress damage. Previous studies have shown that the oxidative state at the beginning of female reproduction affects the start and quality of the reproductive behavior [29]. Oxidative stress levels in previous reproductive behaviors may affect the quality of the next reproductive behavior. Activating CAMK4 via the calcium signaling pathway can downregulate PGC-1α and regulate oxidative stress in the body.
In G5 vs. G3, HSP90AB1 was upregulated which enriched the estrogen signaling and progesterone-mediated maturation pathway. In pigs, heat stress and lipopolysaccharides specifically regulate the protein abundance of ovarian HSP in the follicular and luteal phases. In cattle production, low estrus rates caused by heat stress may be related to HSP90AB1 [30]. Three differentially expressed genes, i.e., CENPA, LOC506989, and ZNF282, were found in G4 vs. G6. Although these three genes were not enriched in the KEGG pathway, early research has shown that CENPA is highly stable and is not supplemented after birth.
In this study, compared with the 6-month-old calves, the 18-month-old adults showed more enhanced sensory, metabolic, and physiological functions. This period is the main period of bodily development, but the growth rate of these functions gradually stabilizes with age until the body is fully developed. Between 18 months and 30 months of age, Angus cows showed better ovulation, which is a characteristic of sexual maturity. In 6-month-old pituitary tissue, only LOC101904796 was persistently highly expressed in 18-month-old and 30-month-old cows as compared with the 6-month-old cows, and this high expression did not differ between 18 and 30 months of age. LOC101904796 is a heterogeneous ribonucleoprotein responsible for the RNA recognition motif. No related reports have been published on this gene. According to the COG classification, we speculate that this gene is related to transcription and ribosome structure, which requires further research.
In production, factors such as photoperiod, dietary structure, and living environment affect carcass metabolism. Neurons in the brain that control metabolism are connected to reproductive neurons. The metabolic state affects the reproductive state. In our study, different differentially expressed genes were activated in multiple pathways such as metabolism, feeding, neurobehavior, and taste. ORM1 was co-expressed in both tissues, and the difference appeared in 30-month-old cattle as compared with 6-month-old cattle. The endoplasmic reticulum membrane protein family plays a central role in lipid homeostasis and protein quality control. In our study, ORM1 was significantly upregulated in 30-month-old cattle, likely related to the accelerated decomposition of ketone bodies in cows at this time. ORM1 upregulation may regulate lipid deposition and meat properties in cows.
The photoperiod must be controlled to improve bovine reproductive performance. This study showed activation of the circadian rhythm pathway in the hypothalamus in 18-month-old cattle. RORB is highly expressed in various animal models in the areas of the sensory organs, spinal cord, and brain that process circadian rhythm information [31]. In the pituitary glands of the 30-month-old cows, GNAT2 upregulation activated the photoperiod pathway. Although different genes regulate the photoperiod, the photoperiod affects the reproductive state.
Testosterone may play a role in ovulation in cows. Testosterone deficiency can effectively inhibit the surge of progesterone and LH, thereby, preventing ovulation [32]. In our study, TEX15 was significantly downregulated in the 30-month-old pituitary tissues. TEX11, TEX12, TEX14, and TEX15, are germ cell-specific genes expressed in the testes. TEX15 plays an important role in the correct assembly of synaptic complexes and meiosis. Interestingly, however, TEX15 expression was found in both the testes and ovaries [33]. In mouse models, loss of TEX15 gene function leads to early arrest of meiosis in male mice, which prevents meiotic recombination; however, this does not occur in female mice [34]. Downregulation of the TEX genes has been found in patients with azoospermia [35], confirming that this transcriptome is essential in the early stages of spermatogenesis. Although most functions of TEX15 are unknown, evidence suggests that TEX15, as a new susceptibility gene, is associated with breast cancer [36].

Conclusions
In general, we obtained mRNA expression profiles of the pituitary and hypothalamus of Angus cattle at different growth and developmental periods, providing transcriptome data for the study of the hypothalamic pituitary of Angus cattle. Through transcriptome data analysis, we found that in both hypothalamus and pituitary, the number of different expressed genes between 18 months and 6 months was greater than that between 30 months and 18 months. This suggests that gene expression in the hypothalamus and pituitary gland tend to stabilize in Angus cattle after adulthood. Consistent with our conjecture, both in the hypothalamus and pituitary, the differentially expressed genes between 18 months vs. 6 months were enriched in growth, development, and sexual maturity, while the differentially espressed genes between 30 months and 18 months were enriched in milk promotion. In the future, we need to conduct further functional verification of these differentially expressed genes to further understand the regulation mechanism of the hypothalamic pituitary gonad axis on Angus cattle reproduction.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13060965/s1, Figure S1: Statistics for the variable splicing events; Figure S2: Distribution of reads in different regions of the genome; Figure S3: SNP density map and annotation classification diagram; Figure S4: Sample gene expression and overall distribution; Figure S5: The results of RT-qPCR were consistent with the NGS; Figure S6: Pituitary (left) and hypothalamus (right) DEGs network; Table S1: Grouping and numbering information; Table S2: Sequencing data and comparison results for each sample; Table S3: Statistics of the variable shear event structure; Table S4: Primers for real-time quantitative PCR; Table S5: SNP site statistics;  Table S6: Numbers of differential genes; Table S7: New gene function annotation results; Table S8: Reproductive-related GO classification table involving new transcripts between groups; Table S9: Comparison of pituitary hormone genes and related receptor expression levels; Table S10 Institutional Review Board Statement: All animal protocols were carried out according to the guidance of the Animal Care Committee at Jilin University (SY201903016).

Informed Consent Statement: Not applicable.
Data Availability Statement: All data involved in this article are original and available from the corresponding authors on reasonable request.

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