Deep sequencing of the transcriptome in the anterior pituitary of heifers before and after ovulation.

We aimed to determine gene expression patterns in the anterior pituitary (AP) of heifers before and after ovulation via deep sequencing of the transcriptome (RNA-seq) to identify new genes and clarify important pathways. Heifers were slaughtered on the estrus day (pre-ovulation; n=5) or 3 days after ovulation (post-ovulation; n=5) for AP collection. We randomly selected 4 pre-ovulation and 4 post-ovulation APs, and the ribosomal RNA-depleted poly (A)+RNA were prepared to assemble next-generation sequencing libraries. The bovine APs expressed 12,769 annotated genes at pre- or post-ovulation. The sum of the reads per kilobase of exon model per million mapped reads (RPKM) values of all transcriptomes were 599,676 ± 38,913 and 668,209 ± 23,690, and 32.2 ± 2.6% and 44.0 ± 4.4% of these corresponded to the AP hormones in the APs of pre- and post-ovulation heifers, respectively. The bovine AP showed differential expression of 396 genes (P<0.05) in the pre- and post-ovulation APs. The 396 genes included two G-protein-coupled receptor (GPCR) genes (GPR61 and GPR153) and those encoding 13 binding proteins. The AP also expressed 259 receptor and other 364 binding proteins. Moreover, ingenuity pathway analysis for the 396 genes revealed (P=2.4 × 10-3) a canonical pathway linking GPCR to cytoskeleton reorganization, actin polymerization, microtubule growth, and gene expression. Thus, the present study clarified the novel genes found to be differentially expressed before and after ovulation and clarified an important pathway in the AP.

The anterior pituitary (AP) receives signals from both the hypothalamus and various peripheral tissues, and secretes the important hormones that control various important functions in multiple organs. At least 400 bioactive peptides may be synthesized for various purposes, including paracrine, autocrine, and endocrine roles [11] by the heterogeneous secretory cells (corticotrophs, gonadotrophs, lactotrophs, somatotrophs, and thyrotrophs) and non-secretory cells. For these various roles, the AP is expected to express a large number of genes. However, little is known regarding global gene expression in the AP of animals.

Ovulation and fertilization are the critical steps in successful reproduction. Ovulation is controlled by hormones secreted by the AP, namely luteinizing hormone (LH) and follicle-stimulating hormone (FSH). Anterior pituitary hormones also induce important events related to ovulation, such as oocyte maturation, embryo development, and formation and maintenance of the corpus luteum. Peripheral tissues, especially the ovaries, secrete hormones to control gene expression in the AP. On Day 3 after ovulation, when blood concentrations of estradiol and progesterone are low, active pulsatile LH secretion occurs in heifers [22]. Therefore, we can expect that the expression levels of several genes in the AP on Day 3 after ovulation are different from those in the AP before ovulation.

Deep sequencing of the transcriptome (RNA-seq) via next-generation sequencing (NGS) technology is the most recent and high-throughput type of genetic analysis tool that can be used to determine global gene expression. However, only a limited number of previous studies have used this technique for pituitary gene expression. He et al. [18] used this technique on 200 or 300 whole-pituitary samples (and not just the AP), collected from zebrafish before and after sexual maturation. This was because the size of the pituitary gland in zebrafish is too small to analyze the genes in the AP alone. Reads per kilobase of exon model per million mapped reads (RPKM) is a method of quantifying gene expression from RNA-seq data by normalizing for total read length and the number of sequencing reads. Few studies have analyzed global gene expression in APs in post-pubertal animals. To the best of our knowledge, only 2 previous studies utilized oligonucleotide microarrays to compare gene expression in APs in dairy cows in the estrous or luteal phases [25] and in anestrous and cycling postpartum beef cows [38]. However, NGS methods are becoming increasingly common because of certain advantages relative to microarrays, and RNA-seq has the capacity to identify novel transcript variants and is not limited by the potential for cross-hybridization [52].

Therefore, the aims of this study were to determine gene expression patterns in the AP of heifers before and after ovulation using RNA-seq in order to discover new genes involved in this process and to clarify important gene networks in AP before and after ovulation.

MATERIALS AND METHODS

Animals and treatments

All experiments were performed according to the Guiding Principles for the Care and Use of Experimental Animals in the Field of Physiological Sciences (Physiological Society of Japan) and approved by the Committee on Animal Experiments of the School of Veterinary Medicine, Yamaguchi University. Post-pubertal Japanese Black heifers (30-month-old) were housed in a free-stall barn. Their daily diet included Italian ryegrass hay (84.2% dry matter [DM], 2.30 Mcal of metabolisable energy per kilogram of DM [ME·kg–1 DM], 13.3% crude protein [CP]) and concentrate (86.6% DM, 3.82 Mcal ME·kg–1 DM, 21.3% CP), and were given 23.9 Mcal of ME. The heifers received two intramuscular injection of dinoprost (Pronalgon F; Pfizer, Tokyo, Japan) given 11 days apart to control estrous stage. The heifers were slaughtered on the day of estrus (pre-ovulation; n=5) or 3 days after ovulation (post-ovulation; n=5). After the second dinoprost injection, we observed estrous behavior at least thrice a day for 20 min, based on the reported criteria [7] for color change in heat-mount detectors (Kamar Product Inc., Zionville, IN, U.S.A.). Briefly, the pre-ovulation heifers behaved restlessly, roamed, discharged mucous, stood to be mounted, and their vulva were sniffed by other heifers. The pre-ovulation heifers were slaughtered within 3 hr after confirmation of standing heat. In contrast, the post-ovulation heifers did not show estrous behavior on the day of slaughtering. We also macroscopically examined the ovaries and uterus immediately after slaughtering to verify the ovulation stage, based on the reported criteria [32]. The pre-ovulation heifers had a small (5–15 mm) corpus luteum as light yellow to white color and covered with connective tissues, a large (12–20 mm) dominant follicle, pale pink endometrium, and oedema in the stroma and mucus. In contrast, the post-ovulation heifers reached the ovulation point with a blood-colored surface, no dominant follicles, pink to red endometrium, and slight oedema in the stroma. Anterior pituitaries were collected from the head as previously described [40], immediately frozen in liquid nitrogen, and preserved at −80°C until RNA extraction.

Preparation of the transcriptome library for deep sequencing

Total RNA was extracted from the APs of the heifers using an RNeasy Mini kit (Qiagen, Valencia, CA, U.S.A.), and subsequently treated with ribonuclease-free deoxyribonuclease (Qiagen) to eliminate possible genomic DNA contamination. The quality was verified by agarose gel electrophoresis and capillary electrophoresis using a Bioanalyzer (Agilent Technologies, Santa Clara, CA, U.S.A.). The value of RNA integrity number (RIN) in the samples was >8.0, and absence of genome DNA contamination was confirmed by electrophoresis. Poly (A)+RNA was purified from 3 µg total RNA with the Gene Read Pure mRNA kit (Qiagen) and the concentration of poly (A)+RNA was measured using the Qubit RNA Assay kit (Thermo Fisher Scientific, Wilmington, DE, U.S.A.). We randomly selected four pre-ovulation and four post-ovulation heifers. The poly (A)+RNA sample (40.8 ± 4.2 ng) of heifers was fragmented by incubation with RNaseIII for 10 min at 37°C, after which, RNA size was confirmed as being approximately 150 bp using a Bioanalyzer. The purified fragmented RNA was ligated with adaptors using SOLiD Total RNA-Seq kit and reverse-transcribed to cDNA by using ArrayScript Reverse Transcriptase (Thermo Fisher Scientific). After purification with magnetic beads (AMPure, XP, Beckman Coulter Inc., Brea, CA, U.S.A.), cDNA was amplified for 15 cycles using Taq polymerase (AmpliTaq, Thermo Fisher Scientific) with a barcode-labeled 3′-PCR primer (SOLiD RNA Barcoding kit Module 1–16; Thermo Fisher Scientific) to differentiate each cDNA sample. The PCR products were purified by PureLink PCR Nano column (Thermo Fisher Scientific) and their quality was verified using a Bioanalyzer. Each DNA library was determined as approximately 257 bp similar to the molecular mass of the fragmented RNA ligated to the adaptors.

Bead preparation and RNA-sequencing by SOLiD5500

The DNA library was emulsified using SOLiD EZ Bead Emulsifier (Thermo Fisher Scientific) and the products were amplified by emulsion PCR using SOLiD EZ Bead Amplifier (Thermo Fisher Scientific). The bead-combined PCR products were concentrated using SOLiD EZ Bead Enricher (Thermo Fisher Scientific) and bead volumes were measured by NanoDrop. The bead-combined DNA (4.9 × 106 beads) was loaded onto four lanes of a flowchip, with 269,500 beads per lane for the SOLiD5500 system (Thermo Fisher Scientific).

Analysis of read information

Short sequence reads of 75 bp in the XSQ file created by the SOLiD5500 system were assembled and mapped with CLC Genomic Workbench software (ver. 8.5.1, CLC Bio, Aarhus, Denmark) using Bos_taurus UMD3.1 (National Center for Biotechnology Information [NCBI] assembly accession GCA_000003055.3) as a reference sequence. RNA-Seq analysis was carried out using parameters as below. For trimming of reads, the ambiguous limit was two, the minimum number of nucleotides in reads was 25, and the quality limit was 0.05. For reads mapping, the number of mismatch was two, number of insertion was three, number of deletion was three, and number of color error was three. The RPKM value was used to represent the expression level for each gene.

Pathway analysis

Globally expressed genes were compared to elucidate important pathways using Ingenuity Pathways analysis (IPA) build version 377306M (Qiagen) utilizing content version 27216297 (release date 16th March, 2016). The IPA discerns differences between molecular and cellular functions and canonical pathways of heifers based on millions of findings reported in the literature, and the software is updated weekly. The IPA uses a Fisher’s exact test to determine whether the input genes are significantly related to pathways as compared to the complete ingenuity knowledge base. We selected genes (1) with varying expression (P<0.05) between the pre- and post-ovulation APs, and also (2) with the average RPKM equal to or more than 1 in either or both of the pre- and post- ovulation APs. Then, we used the IPA software, employing the fold change of pre-ovulation expression to post-ovulation expression of the selected genes, on 17th June 2016. The parameters were as default. The molecule activity predictor (MAP) function was utilized to predict both upstream and downstream molecules and functions after clarification of canonical pathway.

Real-time PCRs to evaluate differences in the expression of central genes in important pathways

Real-time PCRs were performed to confirm significant differences in the expression of the GPR genes that were differently expressed in the pre- and post-ovulation samples (n=5 for each group) and the genes involved in the important pathway that was elucidated by the IPA. The cDNA was synthesized from 2 µg of the total RNA per AP in 20-µl reaction mixtures containing random hexamer primers using the high-capacity cDNA reverse transcription kit (Thermo Fisher Scientific). Table 1 shows the primers designed for real-time PCR using Primer Express Software V3.0 (Thermo Fisher Scientific) based on the reference sequences. The amount of gene expression was measured in duplicate by real-time PCR analyses with 20 ng cDNA, using CFX96 Real Time PCR System (Bio-Rad, Hercules, CA, U.S.A.) and Power SYBR Green PCR Master Mix (Thermo Fisher Scientific), together with a 6-point relative standard curve, non-template control, and no reverse-transcription control. Standard 10-fold dilutions of purified and amplified DNA fragments were prepared. Temperature conditions for all genes were as follows: 95°C for 10 min for pre-denaturation; five cycles each of 95°C for 15 sec and 66°C for 30 sec; and 40 cycles each of 95°C for 15 sec and 60°C for 60 sec. Melting curve analyses were performed at 95°C for each amplicon and each annealing temperature to ensure the absence of smaller non-specific products such as dimers. To optimize the real-time PCR assay, serial dilutions of a cDNA template were used to generate a standard curve by plotting the log of the starting quantity of the dilution factor against the Cq value obtained during amplification of each dilution. Reactions with a coefficient of determination (R2) >0.98 and efficiency between 95 and 105% were considered optimized. The concentration of PCR products was calculated by comparing Cq values of unknown samples with the standard curve using appropriate software (CFXmanagerV3.1, Bio-Rad). The gene expression levels for each of the eleven genes were normalized to the geometric mean of the expression levels of two house-keeping genes, GAPDH (NCBI reference sequence, NM_001034034) and RANBP10 (NM_001098125). We selected these 2 house-keeping genes from 20 previously reported candidate house-keeping genes [37, 47] because they had the smallest inter-heifer coefficients of variation of RPKM value.

Table 1.

Name and details of the primers used for real-time PCRs

Gene name	Primer	Sequence 5ʹ-3ʹ	Size (bp)
GAPDH	Forward	TGGTGAAGGTCGGAGTGAAC	91
	Reverse	ATGGCGACGATGTCCACTTT
RANBP10	Forward	CCCAGTCCTACCAGCCTACT	133
	Reverse	CCCCCAGAGTTGAATGACCC
GPR61	Forward	CATCAACGTGGAGCGCTACTAT	62
	Reverse	GCGTCATTCGCACCTCATAA
GPR153	Forward	GCTGAGCAACGCCAAGAAG	70
	Reverse	CGACAGGATGAAGGACACCAT
ACTA2	Forward	TGTCCACCTTCCAGCAGATG	58
	Reverse	GATGGCCCGGCTTCGT
CDH10	Forward	ACGGCCCCGGATAACG	56
	Reverse	TGTTCCTTGAGCCTTTCATTGA
ELK1	Forward	AGTACTTCCCCAGCGATGGA	59
	Reverse	AGCAGCTGCAGCAGAAACTG
GNA14	Forward	GGCGGGAGTACCAACTGTCA	62
	Reverse	GGCGATCCGGTCAATGTC
LIMK1	Forward	GTGTTCTCGTTCGGGATCGT	59
	Reverse	TCCGGGTCAGCGTTCACT
LIMK2	Forward	TGCATTCCATGTGTATCATCCA	60
	Reverse	CCAGCTTGATGAGGCAGTTG
RHOC	Forward	CGCCTGCGGCCTCTCT	59
	Reverse	TCGATGGAGAAGCACATGAGAA
SEPT5	Forward	CAGCAGGACCGGGAACTG	57
	Reverse	CGTGTTGCTGCCGATAACAG
WIPF1	Forward	GCTCCCTGCCTGCTTGTG	65
	Reverse	TGGCTGACAAATCGTTAAAGGA

Statistical analysis

The statistical significance of differences in RPKM value of each gene between the pre- and post-ovulation heifers (n=4 for each group) was analyzed by non-paired t-test. We calculated the average of the RPKM values of pre-ovulation APs and divided it by the average of the RPKM values of the post-ovulation APs (Pre/Post ratio) to investigate the differences in the RPKM value of each gene. The statistical significance of the differences between the gene expression levels of the pre- and post-ovulation heifers, as measured by real-time PCR (n=5 for each group), was analyzed by non-paired t-test. The level of significance was set at P<0.05. Data are expressed as mean ± SEM.

RESULTS

Information from NGS

The number of nucleotides detected in the APs from the pre- and post-ovulation heifers was 12,266 ± 737 and 14,163 ± 1,936 Mbases, respectively. After low quality and ambiguous reads were removed from the detected reads, and high quality reads were mapped on the reference sequences of bosTau7 database, 9.7 ± 0.6 and 11.3 ± 1.6 million read numbers were obtained in the pre- and post-ovulation samples, respectively.

Annotated genes expressed in APs in the pre- and post-ovulation periods

The sum of the RPKM values of all transcriptomes (total RPKM values) was 599,676 ± 38,913 and 668,209 ± 23,690 in the APs from the pre- and post-ovulation heifers, respectively. There was no difference between the 2 groups (P>0.1). In total, 12,769 annotated genes (for which the average RPKM was equal to or more than 1 in either or both of the pre- and post- ovulation APs) were expressed in both groups of AP samples put together. There were differences (P<0.05) in the expression level of 396 genes between the pre- and post-ovulation APs. The details of 396 genes are shown in Supplementary Table 1. Among the 396 genes, only 31 were more expressed in the post-ovulation AP than the pre-ovulation AP; the remaining 365 genes showed higher expression in the pre-ovulation AP than the post-ovulation AP. The 31 genes that were expressed more after ovulation included AGXT2 (alanine-glyoxylate aminotransferase 2), ARG2 (arginase 2), CDH23 (cadherin-related 23), DYRK3 (dual specificity tyrosine-Y-phosphorylation regulated kinase 3), PADI1 (peptidyl arginine deiminase, type I), PNP (purine nucleoside phosphorylase), PNRC1 (proline rich nuclear receptor coactivator 1), SPAG8 (sperm associated antigen 8), and TMEM35 (transmembrane protein 35). Among these, 11 genes had only the Ensemble gene identifier (ENSBTAG), and did not have names.

Transcriptomes of pituitary hormone genes

The sum of the RPKM values of the AP hormone genes (PRL, GH1, CGA, LHB, POMC, FSHB and TSHB) accounted for 32.2 ± 2.6 or 44.0 ± 4.4% of the about 630,000 total RPKM values of all transcriptomes, in the pre-ovulation group and the post-ovulation group, respectively. There were no statistically significant differences in the RPKM of the AP hormone genes between the pre- and post-ovulation heifers (Table 2). The RPKM value of PRL (prolactin) was the largest among the all genes, including those encoding AP hormones.

Table 2.

The RPKM values (mean ± SEM) of bovine AP genes in the pre- and post-ovulation periods, arranged as per the Pre/Post ratio (average of the RPKM values of pre-ovulation APs divided by the average of the RPKM values of post-ovulation APs)

Gene name	Description	Accession number	RPKM value		Pre/Post ratio
			Pre-ovulation	Post-ovulation
TSHB	thyroid stimulating hormone beta	NM_174205	4,249 ± 1,564	981 ± 333	4.33
FSHB	follicle stimulating hormone beta	NM_174060	897 ± 519	426 ± 151	2.10
CGA	glycoprotein hormones alpha	NM_173901	25,558 ± 7,763	15,219 ± 2,170	1.68
POMC	proopiomelanocortin	NM_174151	3,055 ± 513	1,889 ± 461	1.62
LHB	luteinizing hormone beta	NM_173930	7,081 ± 2,009	6,783 ± 1,974	1.04
GH1	growth hormone 1	NM_180996	38,182 ± 15,332	56,718 ± 15,351	0.67
PRL	prolactin	NM_173953	114,353 ± 58,951	211,708 ± 31,501	0.54

No significant differences were found between pre- and post-ovulation APs in any of the genes.

Transcriptomes of non-hormonal genes with the highest RPKM values

Supplementary Tables 2 and 3 show the RPKM values of non-hormonal genes with the highest RPKM values. The RPKM value of NNAT (neuronatin) was the largest among the listed genes in both the pre- and post- ovulation APs. Some of the other genes were ribosomal proteins.

Genes of receptors and binding proteins

Bovine APs expressed a total of 259 receptor genes (for which the average RPKM was equal to or more than 1 in all APs in either or both of the pre- and post- ovulation APs). The details of 259 receptor genes are shown in Supplementary Table 4. Table 3 shows all of the 11 receptor genes differently (P<0.05) expressed in the pre- and post-ovulation APs. Two orphan G-protein-coupled receptors (GPCRs), GPR61 and GPR153, are listed among the 11 receptor genes. The real-time PCRs verified the difference in GPR61 and GPR153 expression levels (P<0.01; Fig. 1).

Table 3.

The RPKM values (mean ± SEM) of receptor genes expressed differently (P<0.05) between pre- and post-ovulation periods, arranged as per the Pre/Post ratio

Gene name	Description	Accession number	RPKM value		P-value	Pre/Post ratio
			Pre-ovulation	Post-ovulation
SCARA3	scavenger receptor class A, member 3	XM_002689483	8.06 ± 1.24	3.02 ± 0.50	<0.01	2.67
CRHR1	corticotropin releasing hormone receptor 1	NM_174287	12.85 ± 2.75	4.95 ± 1.02	<0.05	2.60
GPR153	G protein-coupled receptor 153	XM_005217164	8.50 ± 1.16	3.68 ± 0.92	<0.05	2.31
CX3CR1	chemokine (C-X3-C motif) receptor 1	NM_001102558	4.89 ± 0.67	2.56 ± 0.46	<0.05	1.91
P2RY6	pyrimidinergic receptor P2Y, G-protein coupled, 6	NM_001192295	2.20 ± 0.25	1.18 ± 0.04	<0.01	1.86
GPR61	G protein-coupled receptor 61	NM_001038571	2.85 ± 0.38	1.66 ± 0.22	<0.05	1.72
CELSR2	cadherin, EGF LAG seven-pass G-type receptor2	NM_001192931	4.60 ± 0.55	2.80 ± 0.43	<0.05	1.64
RARA	retinoic acid receptor, alpha	NM_001014942	11.13 ± 1.40	7.03 ± 0.76	<0.05	1.58
IL27RA	interleukin 27 receptor, alpha	NM_001098028	37.25 ± 2.95	25.71 ± 1.23	<0.05	1.45
NR1D1	nuclear receptor subfamily 1, group D, member 1	NM_001078100	14.83 ± 0.78	11.55 ± 0.97	<0.05	1.28
OGFR	opioid growth factor receptor	NM_001077019	7.80 ± 0.25	6.23 ± 0.48	<0.05	1.25

The table only shows genes for which the RPKM was equal or more than 1 in all APs in either or both of the pre- and post-ovulation group.

Fig. 1.

The mean ± SEM (n=5 for each stage) of GPR 61 and GPR 153 expression (measured by real time PCR) in bovine AP, before and after ovulation. Gene expression levels were normalized to the geometric means of two housekeeping genes, GAPDH and RANBP10. **P<0.01: significant difference compared to post-ovulation.

There were no significant differences in the RPKM values of steroid hormone receptors and hypothalamic hormone receptors (data not shown).

Supplementary Table 5 shows all the 18 orphan GPCRs expressed in bovine APs and the difference in the expression of each GPCR gene. This table does not contain other GPCRs for which the average RPKM values in both the pre- and post-ovulation APs were less than 1.

The bovine AP expressed receptors for insulin, IGF-1 and IGF-2. The bovine AP expressed also 364 binding proteins including IGF binding protein (IGFBP) types 2, 3, 4, 5, 6 and 7. The details of these and other binding proteins are shown in Supplementary Table 6. Table 4 shows all the 13 binding protein genes differently (P<0.05) expressed in the pre- and post-ovulation APs. This table omitted other binding proteins for which average the RPKM values in both the pre- and post-ovulation APs were less than 1.

Table 4.

The RPKM values (mean ± SEM) of binding protein genes expressed differently (P<0.05) between the pre- and post-ovulation periods, arranged as per the Pre/Post ratio

Gene name	Description	Accession number	RPKM value		P-value	Pre/Post ratio
			Pre-ovulation	Post-ovulation
RLBP1	retinaldehyde binding protein 1	NM_174451	7.39 ± 1.91	2.02 ± 0.73	<0.05	3.66
GNA14	guanine nucleotide binding protein (G protein), alpha 14	NM_174323	2.97 ± 0.61	1.12 ± 0.29	<0.05	2.64
DMTN	dematin actin binding protein	NM_001034431	38.69 ± 5.05	20.58 ± 2.53	<0.05	1.88
RBFOX3	RNA binding protein, fox-1 homolog (C. elegans) 3	NM_001075537	1.59 ± 0.07	0.90 ± 0.01	<0.01	1.78
RAC2	ras-related C3 botulinum toxin substrate 2 (rho family, small GTP binding protein Rac2)	NM_175792	3.04 ± 0.17	1.72 ± 0.16	<0.01	1.77
N4BP3	NEDD4 binding protein 3	NM_001205742	2.27 ± 0.30	1.30 ± 0.18	<0.05	1.75
MAP3K7IP1	TGF-beta activated kinase 1/MAP3K7 binding protein 1	NM_001102057	10.49 ± 1.62	6.27 ± 0.59	<0.05	1.67
CREB3L4	cAMP responsive element binding protein 3-like 4	XM_015462516	1.76 ± 0.15	1.08 ± 0.21	<0.05	1.62
TAPBP	TAP binding protein (tapasin)	NM_001045885	77.93 ± 9.65	50.99 ± 3.10	<0.05	1.53
CIRBP	cold inducible RNA binding protein	NM_001034278	57.72 ± 4.20	39.15 ± 5.35	<0.05	1.47
RBBP9	RB binding protein 9, serine hydrolase	NM_001083424	3.75 ± 0.15	2.62 ± 0.23	<0.01	1.43
MYBBP1A	MYB binding protein (P160) 1a	XM_010815952	11.07 ± 0.76	8.66 ± 0.49	<0.05	1.28
DRG2	developmentally regulated GTP binding protein 2	NM_001014865	11.70 ± 0.55	9.95 ± 0.40	<0.05	1.18

The table only shows genes for which the RPKM was equal or more than 1 in all APs in either or both of the pre- and post- ovulation group.

Results of IPA analysis

The IPA analysis employed the fold change of 396 genes differently (P<0.05) expressed between the pre- and post-ovulation APs. The results clarified a canonical pathway termed as “Signaling by Rho Family GTPases” (Fig. 2). The P-value was 2.4 × 10−3 and the biased Z-score was 2.3. The 236 genes constituted the canonical pathway, and 9 genes [ACTA2 (actin, alpha 2, smooth muscle, aorta), CDH10 (cadherin 10), ELK1 (ELK1, ETS transcription factor), GNA14 (G protein subunit alpha 14), LIMK1 (LIM domain kinase 1), LIMK2 (LIM domain kinase 2), RHOC (ras homolog family member C), SEPT5 (septin5), WIPF1 (WAS/WASL interacting protein family member 1)] belonging to the pathway were identified in the bovine APs used in this study. Table 5 shows the RPKM values of the 9 genes, and all of the 9 genes were higher in the pre-ovulation APs than post-ovulation APs. The results of the real-time PCRs verified these differences (P<0.05; Fig. 3).

Fig. 2.

The subcellular location of genes in the canonical pathway, termed as “Signaling by Rho Family GTPases”, as determined via IPA analysis to compare transcriptomes in the APs of pre- and post-ovulation heifers. Note that the IPA software shows gene names as normal font and not italics. The solid lines and dotted line indicate direct or indirect biological relationship between genes, respectively. The bold grey lines indicate the border between the extracellular space, cytoplasm, and nucleus. The purple indicates the upregulated genes identified in this experiment. This figure includes the results of molecule activity predictor (MAP), and the predicted stimulated genes are shown in orange, and the predicted inhibited genes are shown in blue. In total, 236 genes constituted the canonical pathway, and 9 of these genes (ACTA2, CDH10, ELK1, GNA14, LIMK1, LIMK2, RHOC, SEPT5, and WIPF1) were identified in the bovine APs used in this study. In this figure, F-actin represents ACTA2, Cadherin represents CDH10, Gα represents GNA14, LIMK represents both LIMK1 and LIMK2, Rho represents RHOC, Septin5 represents SEPT5, and WIP represents WIPF1.

Table 5.

The RPKM values (mean ± SEM) of 9 genes consisting the “Signaling by Rho Family GTPases”, arranged in alphabetical order of the gene name

Gene name	Description	Accession number	RPKM value		P-value	Pre/Post Ratio
			Pre-ovulation	Post-ovulation
ACTA2	actin, alpha 2, smooth muscle, aorta	NM_001034502	63.84 ± 17.52	18.14 ± 2.60	<0.05	3.52
CDH10	cadherin 10	NM_001076266	1.03 ± 0.35	0.11 ± 0.03	<0.05	9.21
ELK1	ELK1, member of ETS oncogene family	NM_001191236	7.03 ± 0.53	5.62 ± 0.15	<0.05	1.25
GNA14	guanine nucleotide binding protein (G protein), alpha 14	NM_174323	2.97 ± 0.61	1.12 ± 0.29	<0.05	2.64
LIMK1	LIM domain kinase 1	NM_001206904	10.72 ± 0.91	7.87 ± 0.33	<0.05	1.36
LIMK2	LIM domain kinase 2	NM_001038098	20.99 ± 1.21	16.19 ± 0.98	<0.05	1.30
RHOC	ras homolog family member C	NM_001046138	33.78 ± 3.85	21.38 ± 2.85	<0.05	1.58
SEPT5	septin 5	NM_001076371	32.07 ± 1.20	25.25 ± 1.41	<0.05	1.27
WIPF1	WAS/WASL interacting protein family member 1	NM_001076923	2.88 ± 0.63	1.18 ± 0.20	<0.05	2.43

Fig. 3.

The mean ± SEM (n=5 for each stage) of the expression levels of 9 genes (ACTA2, CDH10, ELK1, GNA14, LIMK1, LIMK2, RHOC, SEPT5 and WIPF1) comprising the “Signaling by Rho Family GTPases” (measured by real time PCR) in bovine AP, before and after ovulation. Gene expression levels were normalized to the geometric means of two housekeeping genes, GAPDH and RANBP10. *P<0.05: significant difference compared to post-ovulation. **P<0.01: significant difference compared to post-ovulation.

The IPA software uses the MAP function to predict both upstream and downstream molecules and their functions after clarifying the canonical pathway. The stimulated functions, as predicted by the MAP function of IPA, were cytoskeleton reorganization, cell trafficking, cytokinesis, cytoskeleton regulation, actin membrane linkage, contraction, membrane ruffling, cell-cell adhesion, microtubule-organizing center orientation, and actin-polymerization. Predicted inhibited functions were actin nucleation, microtubule growth, and actin polymerization. Thus, actin-polymerization was predicted as both stimulated and inhibited function. This canonical pathway also controls gene expression via the MEK1/2, ERK1/2, and Elk1 pathway.

DISCUSSION

The present study showed that bovine APs expressed 12,769 annotated genes. Cattle genome contains a minimum of 22,000 genes, with a core set of 14,345 orthologs shared among seven mammalian species including humans [6]. Therefore, about 58% genes are expressed in bovine AP, suggesting the importance of AP.

Up to 32.2 or 44.0% of total RPKM values corresponded to AP hormones in this study. He et al. [18] reported that a similar percentage (44.1%) of the total RPKM values corresponded to eight major AP hormones in fish (including somatolactin, which is absent in the bovine genome) in the whole pituitary from postpubertal zebrafish. Therefore, AP needs to synthesize several mRNAs for its primary role, namely secreting AP hormones.

The RPKM value of PRL was the highest among the all genes including AP hormone genes in this study. Freeman et al. [15] reported that PRL is expressed in both lactotrophs and mammosomatotrophs, and these cells comprise up to 50% of the cellular population in the AP. Prolactin has a broad range of functions in the body, apart from its defining role in promoting lactation [17]. The present study used non-lactating heifers; therefore, the high RPKM value of PRL suggested multiple important roles of prolactin.

As expected, AP was found to express many receptors and binding proteins in this study. However, the roles of most genes listed in Tables 3 and 4 are not well elucidated. For example, murine gonadotroph-like cells express CRHR1 [42], but the functions of CRH have not been clarified yet. RARA in Table 3 is retinoic acid receptor alpha, and RLBP1 in Table 4 is retinaldehyde binding protein 1. Retinoic acid regulates GnRH release and gene expression in the rat hypothalamic fragments [12]. Vitamin A deficient rats have higher serum FSH and LH concentration than controls [21]. It has been established that vitamin A and retinoic acid are very important for reproduction [13]; however, their roles in the pituitary remain to be clarified. Therefore, further studies are important to clarify the roles of the various genes listed in Tables 3 and 4.

To the best of our knowledge, this is the first report of the two orphan receptors, GPR 61 and GPR153, expressed in the AP of any animal species. Although the ligand(s) and functions of GPR61 are unknown, it is known to associate with the Gs protein [44, 46] and to stimulate ERK signaling in neurons [20]. It was also reported to regulate cyclic AMP in Chinese hamster ovary cells [29]. Both ERK and cAMP pathways are important for GnRH-induced LH secretion in bovine gonadotrophs [33, 34, 41]. These findings implied that GPR61 was expressed on the surface of bovine gonadotrophs. In addition, GPR153 is an orphan receptor that is widely expressed in the brain, including in the hypothalamus arcuate nucleus and pituitary [43]. Thus, further studies are required to clarify whether GPR61 and GPR153 are expressed in gonadotrophs, and whether their expression levels change based on the estrous stage. Furthermore, recent studies have indicated that GnRHR co-localized with insulin and glucocorticoid receptors in lipid rafts on the gonadotroph plasma membrane, where they facilitated downstream signaling [36, 49]. Therefore, further studies are required to clarify whether GPR61 and GPR153 are localized on the lipid raft with the GnRH receptor in gonadotrophs.

The “Signaling by Rho Family GTPases” is a cytoplasmic pathway that mediates for any GPCR bound any ligand. The pathway was not reported in the 2 previous studies that used oligonucleotide microarray analysis for bovine APs [25, 38]. The MAP function of IPA predicted the following stimulated functions: cytoskeleton reorganization, cell trafficking, cytokinesis, cytoskeleton regulation, actin membrane linkage, contraction, membrane ruffling, cell-cell adhesion, microtubule-organizing center orientation, and actin-polymerization. It also predicted the following inhibited functions: actin nucleation, microtubule growth, and actin polymerization. ACTA2 and other actin genes encode the cell cytoskeleton. Both cell cytoskeleton and microtubule are very important parts of gonadotrophs for the secretion of LH and FSH [1, 2, 48]. GnRH regulates the morphology and migration of the immortalized LbetaT2 gonadotroph-like cells via the pathway regulating cytoskeletal reorganization [16]. The protein encoded by RhoC, Ras homolog family member C, may have important roles in cytoskeleton protein reorganization by regulating the polymerized actin to control secretory granule exocytosis in AP cells in response to extracellular signals [14]. The LIM kinases, encoded by LIMK1 and LIMK2, mediate the action of estrogen for regulating the dynamics of actin [53]. WIPF1 is an actin-binding protein that regulates actin polymerization [4]. Although the role of the CDH10 product, cadherin-10, in the pituitary has not yet been studied, it is known to be important for cell-cell adhesion in the construction of the blood-brain barrier [50]. The product of SEPT5, septin 5, is a member of the septin GTPase family and is involved in exocytosis in postmitotic neurons [5]. Therefore, the canonical pathway may control AP hormone secretion.

We next investigated the specific AP cells that contained the “Signaling by Rho Family GTPases” pathway. Because this study used a heterogeneous cell mixture, in which a large number of GPCRs were expressed, it was impossible to further clarify the mechanisms. However, the canonical pathway also includes ELK1. Utilizing ELK1, the activated GnRHR enhances the expression of Egr-1 [30], which binds to the promoter region of LHB and increases LHB expression in gonadotrophs [8, 45]. The GNA14 product, Gα (14), is a G protein that can interact with different classes of receptors to regulate phospholipase C [19], which is important for the release of LH in nonhuman primate pituitary [27]. Therefore, the canonical pathway may exist in the gonadotrophs.

The blood levels of estradiol and progesterone are low on day 3 after ovulation, when active pulsatile LH secretion occurs in heifers [22]. Only 31 genes showed higher expression in the post-ovulation AP than the pre-ovulation AP, whereas about 12 times that number (365) of genes was expressed more in the pre-ovulation AP than the post-ovulation AP. The low concentrations of estradiol and progesterone on day 3 after ovulation may be the direct or indirect cause of this difference. The 31 genes that were expressed more after ovulation included PNRC1, SPAG8, AGXT2, and ARG2. The product of PNRC1, proline rich nuclear receptor coactivator 1, interacts with estrogen receptor α [9], and a single nucleotide polymorphism (SNP) in this gene influences egg weight-related and hatchability traits in ducks [10]. The product of SPAG8, sperm associated antigen 8, enhances the transcription of the cAMP response element modulator (CREM) [51], which may influence the estrogen-mediated negative feedback on GnRH neurons [26]. The product of AGXT2, alanine-glyoxylate aminotransferase 2, was recently found to have several novel functions through genomic and metabolomic studies; it also has a unique role in the intersection of key mitochondrial pathways [39]. The product of ARG2, arginase 2, reduces the activity of nitric oxide synthase (NOS) in genital tissues [24]; NOS is known to control the secretion of LH and FSH in the pituitary gland [31]. Thus, these 4 genes may contribute to the synthesis and secretion of LH and FSH in gonadotrophs. However, the roles of the remaining 27 genes, which are expressed more in post-ovulation AP, are unclear. Furthermore, 11 genes still do not have names. Therefore, further studies are required to clarify the roles of these 31 genes.

The roles of most genes listed in Supplementary Tables 2 and 3 are unclear. For example, the NNAT product, neuronatin, is expressed in prenatal and postnatal rat pituitary [23]. However, the role of neuronatin in the adult pituitary has not been clarified yet. The SCGII product, secretogranin II, may be a paracrine factor secreted from lactotrophs to stimulate gonadotropin release [54]. Moreover, estradiol decreases SCGII expression in rat pituitary cells [3]. Therefore, bovine lactotrophs may express SCGII under the control of estradiol to stimulate gonadotropin secretion. SCG5 encodes neuroendocrine protein 7B2 and is expressed in various types of AP cells, including gonadotrophs [28]. GnRH injection increased blood 7B2 concentration in women, but the role of increased 7B2 levels has not been clarified yet [35]. Therefore, these genes may contribute to AP hormone synthesis and secretion in bovine AP.

In conclusion, the present study clarified the novel genes differentially expressed in the pre- and post-ovulation stages as well as an important pathway in the AP. Silent heat, delayed ovulation, and luteal hypoplasia are common problems in domestic animals; however, little is known regarding the underlying pathophysiological mechanisms. Therefore, further studies are required to clarify the roles of most of these genes that show an increase or decrease in expression after ovulation.