Differential genome-wide gene expression profiling of bovine largest and second-largest follicles: identification of genes associated with growth of dominant follicles.

BACKGROUND: Bovine follicular development is regulated by numerous molecular mechanisms and biological pathways. In this study, we tried to identify differentially expressed genes between largest (F1) and second-largest follicles (F2), and classify them by global gene expression profiling using a combination of microarray and quantitative real-time PCR (QPCR) analysis. The follicular status of F1 and F2 were further evaluated in terms of healthy and atretic conditions by investigating mRNA localization of identified genes.
METHODS: Global gene expression profiles of F1 (10.7 +/- 0.7 mm) and F2 (7.8 +/- 0.2 mm) were analyzed by hierarchical cluster analysis and expression profiles of 16 representative genes were confirmed by QPCR analysis. In addition, localization of six identified transcripts was investigated in healthy and atretic follicles using in situ hybridization. The healthy or atretic condition of examined follicles was classified by progesterone and estradiol concentrations in follicular fluid.
RESULTS: Hierarchical cluster analysis of microarray data classified the follicles into two clusters. Cluster A was composed of only F2 and was characterized by high expression of 31 genes including IGFBP5, whereas cluster B contained only F1 and predominantly expressed 45 genes including CYP19 and FSHR. QPCR analysis confirmed AMH, CYP19, FSHR, GPX3, PlGF, PLA2G1B, SCD and TRB2 were greater in F1 than F2, while CCL2, GADD45A, IGFBP5, PLAUR, SELP, SPP1, TIMP1 and TSP2 were greater in F2 than in F1. In situ hybridization showed that AMH and CYP19 were detected in granulosa cells (GC) of healthy as well as atretic follicles. PlGF was localized in GC and in the theca layer (TL) of healthy follicles. IGFBP5 was detected in both GC and TL of atretic follicles. GADD45A and TSP2 were localized in both GC and TL of atretic follicles, whereas healthy follicles expressed them only in GC.
CONCLUSION: We demonstrated that global gene expression profiling of F1 and F2 clearly reflected a difference in their follicular status. Expression of stage-specific genes in follicles may be closely associated with their growth or atresia. Several genes identified in this study will provide intriguing candidates for the determination of follicular growth.

Background

The final stage of bovine follicular development occurs in a wave-like fashion [1,2]. During a wave, increase of follicle-stimulating hormone (FSH) induces recruitment of a cohort of follicles beyond 4 mm in diameter and usually a single follicle is selected as a dominant follicle (DF) [3,4]. Although the DF continues to grow by transition of gonadotropin dependency from FSH to luteinizing hormone (LH) and secretes large quantities of estradiol (E2), the remaining subordinate follicles (SFs) cease to grow, then undergo atresia [5]. It is well documented that increased expression of LH receptor (LHR) in granulosa cells (GC) and specific changes of intrafollicular factors such as the insulin-like growth factor (IGF) and inhibin-activin-follistatin systems play a critical role in E2 production in the DF [6,7]. Therefore, regulatory mechanisms of follicular development are closely associated with complex interactions between follicular local paracrine/autocrine factors and endocrine hormones.

Increasing evidence using global gene expression analysis such as a DNA microarray, suppression subtractive hybridization and serial analysis of gene expression have identified numerous genes in various aspects of bovine follicular development [8-18]. Some studies compared the gene expression profiles between DF and SF around the time of follicular selection. They showed that DF up-regulates genes regulating E2 synthesis, anti-apoptosis, cell proliferation and gene transcription. Conversely, SF enhanced the expression of genes associated with pro-apoptosis and cell death compared with the DF [8,9,13,14]. Recent studies found that 93 mostly novel genes were differently expressed in the GC of newly selected DF compared with SF and/or growing cohort follicles whereas most of these genes were down-regulated in the GC of preovulatory follicles during final maturation before the LH surge [15,17]. Growth of a DF during 2-5.5 days following follicular wave emergence was associated with a decrease in genes encoding proliferation and pro-apoptotic factors and an increase in genes regulating anti-apoptotic factors [12]. An increase in follicular diameter during follicular growth was accompanied by alteration of gene expression regulating some growth factors and cytokines [16,18]. Ndiaye et al. identified a subset of novel genes down-regulated in preovulatory follicles after human chorionic gonadotropin (hCG) stimulation compared with DF, which may contribute to ovulation and luteinization [11].

These previous studies lead us to suggest that gene expression profiles in individual follicles reflect their developmental status, thus each follicle can be classified by differences in gene expression profiles. On the other hand, details of the genetic processes and biological pathways regulating bovine follicular development still remain to be elucidated. We consider that investigating the global gene expression of follicles after selection can help to understand the molecular mechanisms responsible for the regulation and control of follicular development and atresia. Therefore, in this study, we tried to classify the largest (F1) and second-largest (F2) follicles according to differences in gene expression profiles and to identify differentially expressed genes between the groups using a combination of microarray analysis and quantitative real-time PCR (QPCR) analysis. In addition, spatial expression profiles of several identified genes were investigated using in situ hybridization in healthy and atretic follicles classified based on follicular fluid (FF) concentration of steroids.

Methods

Experiment 1: classification of F1 and F2 and identification of genes by microarray analysis and QPCR analysis

Sample collection and RNA extraction

Paired ovaries were obtained from four pregnant Japanese Black cows in the institute ranch less than 10 min after slaughtering. These cows were pregnant and slaughtered for another study. Both F1 and F2 were dissected from the ovaries. Then, the surrounding stroma and theca externa were removed from the follicular walls. We collected three F1 and three F2 from four cows because two cows had both F1 and F2 collected whereas one cow had only a F1 collected and another cow had only a F2 collected. The follicles were snap-frozen and stored at -80°C until RNA extraction. Total RNA from the follicular wall (i.e., granulosa plus theca interna) was extracted from each follicle using ISOGEN (NipponGene, Tokyo, Japan) according to the manufacturer's instructions. All procedures for animal experiments were carried out in accordance with guidelines approved by the Animal Ethics Committee of the National Institute of Agrobiological Sciences for the use of animals.

Microarray analysis

A custom-made bovine oligonucleotide microarray fabricated by Agilent Technologies (Santa Clara, CA, USA) was used in this study. Sixty-mer nucleotide probes for customized microarray were synthesized on a glass slide. The annotated bovine oligonucleotide array represented 10263 sequences 4466 of which were known bovine genes, 5697 were unknown sequences and possible candidates for novel bovine genes, and 100 internal references.

We performed one-color microarray using five follicles (three F1 and two F2). Fluorescence-labeled (Cy3) cRNA probes were prepared from 150-300 ng of total RNA of each follicle using a Low RNA Input Linear Amplification Kit (Agilent Technologies). Labeled cRNA probes (750 ng each) were hybridized to the customized microarray in hybridization buffer (Gene Expression Hybridization Kit, Agilent Technologies) at 60°C for 17 h. After hybridization, the arrays were washed with 6 × SSC, 0.005% Triton X-102 at room temperature for 10 min, followed by 5-min washes in 0.1 × SSC, 0.005% Triton X-102 at 4°C. Hybridized arrays were blow dried with N2 gas and scanned using an Agilent Microarray Scanner (Agilent Technologies), and Feature Extraction ver. 9.1 (Agilent Technologies) was used for image analysis and data extraction. Gene expression datasets were normalized using the median of the signal intensity for 100 GAPDH genes on a microarray platform as internal control.

After normalization, 3308 genes were left to use for further analysis. The relative abundance of individual genes between follicles was calculated by dividing the normalized value of the genes between each follicle. We used the normalized microarray data of genes that showed an expression level of more than 20-fold between at least two follicles for subsequent hierarchical cluster analysis. The data were transformed log2 values and hierarchical cluster analysis was performed using the TIGR MultiExperiment Viewer 4.0 (MeV 4.0) software program [19]. Two parameters (average linkage and cosine correlation) were selected for constructing the hierarchical tree. Compliance with Minimum Information About a Microarray Experiment (MIAME) [20] was assured by depositing all the data in the Gene Expression Omnibus (GEO) repository [21]. The GEO accession numbers are as follows. Platform: GPL9136; Samples: GSM453634, GSM453635, GSM453636, GSM453637 and GSM453638; Series: GSE18145.

Quantitative real-time RT-PCR analysis

To validate the results of microarray analysis, we confirmed mRNA expression of 16 representative genes using QPCR analysis. All six follicles were used in QPCR analysis. The procedures for QPCR were previously described [22]. Briefly, single-strand cDNA was reverse-transcribed from 50 ng of total RNA using MultiScribe™ reverse transcriptase with a random primer, dNTP mixture, MgCl2 and RNase inhibitor (Applied Biosystems, Foster City, CA, USA). The reverse transcription cycle consisted of 10 min annealing at 25°C, 30 min cDNA synthesis at 48°C and 5 min inactivation at 95°C. The primers were designed using the Primer Express computer software program (Applied Biosystems) based on the bovine sequences. The primer sequences for each gene are given in Table 1. Each QPCR reaction (25 μl) contained 1 μl cDNA template, 0.5 μl forward primer (20 μM), 0.5 μl reverse primer (20 μM), 12.5 μl Power SYBR® Green PCR Master Mix (Applied Biosystems) and 10.5 μl nuclease-free water. The thermal cycling conditions included one cycle at 50°C for 2 min, one cycle at 95°C for 10 min, and 40 cycles at 95°C for 15 s and 60°C for 1 min. Each cDNA template was analyzed for quantitation in duplicate. QPCR and the resulting relative increase in reporter fluorescent dye emission were monitored in real time using an Mx3000P QPCR system (Stratagene, La Jolla, CA, USA). The relative difference in the initial amount of each mRNA species (or cDNA) was determined by comparing the cycle threshold values. To quantify the mRNA concentrations, standard curves for each gene were generated by serial dilution of the plasmid containing its cDNA. The melting curve for detecting the SYBR Green-based objective amplicon were confirmed because SYBR Green also detects double-stranded DNA including primer dimers, contaminating DNA and PCR products from misannealed primers. Contaminating DNA or primer dimers appear as a peak separate from the desired amplicon peak.

Table 1

Details of the primers used for quantitative real-time RT-PCR analysis

Gene name	GeneBank accession number	Primer	Sequences	Position
AMH	NM_173890	Forward	5'-ACACCGGCAAGCTCCTCAT-3'	1647-1665
		Reverse	5'-CACCATGTTTGGGACGTGG-3'	1714-1696
CCL2	NM_174006	Forward	5'-CGCTCAGCCAGATGCAATTA-3'	110-129
		Reverse	5'-GCCTCTGCATGGAGATCTTCTT-3'	186-165
CYP19	NM_174305	Forward	5'-TCCATGGGATTTTCCAGGC-3'	2050-2068
		Reverse	5'-TGGTGGCTTGTCTTTTCCAAC-3'	2123-2103
FSHR	NM_174061	Forward	5'-AATCTACCTGCTGCTCATAGCCTC-3'	1300-1323
		Reverse	5'-TTTGCCAGTCGATGGCATAG-3'	1376-1357
GADD45A	NM_001034247	Forward	5'-CCGCATTCATCACAGTGGAA-3'	592-611
		Reverse	5'-CATCACCGTTCAGGGAGATTAATC-3'	704-681
GPX3	NM_174077	Forward	5'-GCTTCCCCTGCAACCAATT-3'	357-375
		Reverse	5'-TCGAACATACTTGAGGGTGGCT-3'	433-412
IGFBP5	NM_001105327	Forward	5'-ACTGTGACCGCAAAGGGTTCT-3'	682-702
		Reverse	5'-TTCATCCCGTACTTGTCCACG-3'	778-758
PlGF	NM_173950	Forward	5'-TGAATGACTCACTCCCTCCATG-3'	877-898
		Reverse	5'-GGTCTGTCTTCTTTCTCTCACGTTC-3'	957-933
PLAUR	NM_174423	Forward	5'-CGCGGCCCTATGAATCAAT-3'	730-748
		Reverse	5'-CTGATGGTGTAGCTTGGGTTCC-3'	800-779
PLA2G1B	NM_174646	Forward	5'-GGCCTTCATCTGCAACTGTGA-3'	358-378
		Reverse	5'-TGTGCTCCTTGTTGTATGGCA-3'	428-408
SCD	NM_173959	Forward	5'-ATTCCCGACGTGGCTTTTTC-3'	659-678
		Reverse	5'-TTCTTTGACAGCTGGGTGTTTG-3'	729-708
SELP	NM_174183	Forward	5'-GTCAAGCAGGGCCACTGACTAT-3'	1700-1720
		Reverse	5'-TCACTAAGCCTGTTGTACCAGCTG-3'	203-2182
SPP1	NM_174187	Forward	5'-AGCCCTGAGCAAACAGACGAT-3'	304-324
		Reverse	5'-GCGTCGTCGGAGTCATTAGAGT-3'	380-359
TIMP1	NM_174471	Forward	5'-CTATGCTGCTGGTTGTGAGGAAT-3'	508-530
		Reverse	5'-TGAGTGTCGCTCTGCAGTTTG-3'	582-562
TRB2	NM_178317	Forward	5'-GACCTCAAGCTTCGGAAATTCA-3'	525-546
		Reverse	5'-CGTCATCTCCCCGCAGAATAT-3'	621-601
TSP2	NM_176872	Forward	5'-GGAAAACAAGTCATGGCGGA-3'	3845-3864
		Reverse	5'-TTGAGAGAAGACAAACAGACCCAG-3'	3928-3902
GAPDH	U85042	Forward	5'-ACCCAGAAGACTGTGGATGG-3'	444-463
		Reverse	5'-CAACAGACACGTTGGGAGTG-3'	621-602

Experiment 2: localization of characteristic genes identified in experiment 1 in healthy and atretic follicles using in situ hybridization

Sample collection and storage

Ovaries containing follicles more than 8 mm in diameter were obtained from Japanese Black cows at local slaughterhouse. We used only follicles which have a transparent follicular wall and fluid and did not show any aspect of cystic follicles. Eleven follicles were collected and 200 μl of FF was aspirated from each follicle by a syringe fitted with a 27G needle. The FF was snap-frozen and stored at -30°C until hormone determinations. The follicles were dissected from the ovaries and fixed in 10% formalin, embedded in paraffin wax, and stored at 4°C until in situ hybridization.

Steroid hormone determinations

Concentrations of E2 and P4 in the FF samples were determined directly in duplicate using a time-resolved fluoroimmunoasssay (TR-FIA). The TR-FIA for E2 and P4 was performed as previously described by our laboratory [23,24]. The FF samples were diluted to 100-, 2000- and 5000-fold for E2 determination and 25-fold for P4 determination using charcoal-treated plasma (collected from adult Japanese-Black cows). Ranges of the standard curves were 5-200 pg/ml for E2 and 0.33-36 ng/ml for P4. The intra- and interassay coefficients of variation were 8.2 and 11.4% for E2, and 8.5 and 10.5% for P4, respectively.

In situ hybridization

We classified follicles into two groups based on relative levels of FF concentrations of E2 and P4 (E2/P4 ≥ 1: healthy; E2/P4 < 1: atretic). Six representative genes differently expressed between F1 and F2 in experiment 1 were selected for in situ hybridyzation: anti-Mullerian hormone (AMH), cytochrome P450, family XIX (CYP19), growth arrest and DNA-damage-inducible, alpha (GADD45A), IGF binding protein 5 (IGFBP5), placental growth factor (PlGF) and thrombospondin 2 (TSP2). In these genes, CYP19 and IGFBP5 were selected as markers of healthy or atretic follicles since mRNA expression of CYP19 and IGFBP5 were up-regulated in the bovine DF and SF, respectively [25,26].

Digoxigenin (DIG)-labeled antisense and sense cRNA probes were prepared as previously described [27,28]. For hybridization, follicles were sectioned into 7 μm-thick sections. We performed in situ hybridization using an automated Ventana HX System Discovery with a RiboMapKit and a BlueMapKit (Roche Diagnostics, Basel, Switzerland) as previously described by our laboratory [27,28]. Briefly, the sections were hybridized with DIG-labeled probes in RiboHybe (Roche Diagnostics) hybridization solution at 65°C (PlGF) or 61°C (AMH, CYP19, GADD45A, IGFBP5 and TSP2) for 6 hours, then washed for 3 × 6 min in RiboWash (Roche Diagnostics) at 65°C and fixed in RiboFix (Roche Diagnostics) at 37°C, 10 min. The hybridization signals were detected with a rabbit polyclonal anti-digoxin antibody HRP conjugate (Dako Cytomation, Carpinteria, CA, USA) using an AmpMapKit (Roche Diagnostics). The hybridized slides were observed with a Leica DMRE HC microscope (Leica Microsystems, Wetzlar, Germany) and a Nikon Digital Sight DS-Fi1-L2 (Nikon, Tokyo, Japan).

Statistical analysis

In experiment 1, the expression ratio of each gene to GAPDH mRNA was calculated to adjust for variations in the QPCR reaction. The follicular diameter and the QPCR data in experiment 1 and concentrations of E2 and P4 and E2/P4 ratio in FF in experiment 2 were analyzed by Mann-Whitney's U test. Results were presented as the mean ± SEM. Statistical significance was considered to be at P < 0.05.

Results

Mean diameter of F1 and F2 were 10.7 ± 0.7 and 7.8 ± 0.2 mm, respectively (P < 0.05).

Hierarchical cluster analysis of microarray data

The expression level of 76 genes was enhanced between at least two follicles by more than 20-fold. Using the microarray data of these 76 genes, we performed a hierarchical cluster analysis and constructed a cluster heat map (Figure 1). As can be seen from the dendrogram of the sample axis, clustering analysis distinctly separated the five follicles into two clusters (A and B) based on their microarray expression profiles. Cluster A included two follicles that were both F2, whereas cluster B contained the remaining three follicles that were all F1. Cluster analysis also identified two major clusters in the gene axis. One cluster contained 31 genes that were relatively highly expressed in cluster A, while the other contained 45 genes that were relatively highly expressed in cluster B. The details of highly expressed genes in clusters A and B are listed in Table 2 and 3, respectively.

Figure 1

Hierarchical cluster analysis of 76 differentially expressed genes in largest (F1) and second-largest follicles (F2). These genes were enhanced between at least two follicles by more than 20-fold. Red scale indicates relative higher expression level and green scale indicates relative lower expression level. The expression levels were transformed to log2 values. Dendrograms of sample axis (above matrix) and gene axis (to the left of matrix) represent overall similarities in gene expression profiles. Five follicles were classified into two major clusters (A and B). The follicles divided into cluster A were all F2 and the follicles divided into cluster B were all F1. The cluster A was characterized by highly expression of 31 genes, whereas the cluster B was predominately expressed 45 genes.

Table 2

List of differentially expressed genes in cluster A as compared with cluster B.

Accession No.	Gene symbol	Gene name
NM_001101080	ADAMTS1	ADAM metallopeptidase with thrombospondin type 1 motif, 1
NM_206843	AGOUTI	AGOUTI protein
NM_001098982	CCDC80	Coiled-coil domain containing 80
NM_174006	CCL2	Chemokine (C-C motif) ligand 2
NM_001008670	CRABP2	Cellular retinoic acid binding protein 2
NM_174290	CRYAB	Crystallin, alpha B
NM_001034247	GADD45A	Growth arrest and DNA-damage-inducible, alpha
NM_001077112	GSTA3	Glutathione S-transferase, alpha 3
NM_001105327	IGFBP5	Insulin-like growth factor binding protein 5
NM_174357	IL1RN	Interleukin 1 receptor antagonist
NM_001033610	KRT8	Keratin 8
NM_001001138	LOC407171	Fc gamma 2 receptor
NM_174384	LOXL4	Lysyl oxidase-like 4
NM_174132	OLR1	Oxidized low density lipoprotein (lectin-like) receptor 1
NM_176855	OXT	Oxytocin
NM_001101883	PDK4	Pyruvate dehydrogenase kinase, isozyme 4
NM_174141	PENK	Proenkephalin
NM_174423	PLAUR	Plasminogen activator, urokinase receptor
NM_001034681	RTKN	Rhotekin
NM_174183	SELP	Selectin P
NM_174137	SERPINE1	Serpin peptidase inhibitor, clade E (nexin, plasminogen activator inhibitor type 1), member 1
NM_001075764	SFRP4	Secreted frizzled-related protein 4
NM_174601	SLC1A5	Solute carrier family 1 (neutral amino acid transporter), member 5
NM_174187	SPP1	Secreted phosphoprotein 1
NM_003254	TIMP1	TIMP metallopeptidase inhibitor 1
NM_176872	TSP2	Thrombospondin 2
AW430112		Transcription factor B1, mitochondrial
BE721140		Transcribed locus
BP101259		Caldesmon, smooth muscle
XM_587930		Similar to ATP-binding cassette sub-family G member 1 (ABCG1), mRNA.
XM_869699		Similar to tumor necrosis factor receptor superfamily, member 12A

Table 3

List of differentially expressed genes in cluster B as compared with cluster A.

Accession No.	Gene symbol	Gene name
NM_173890	AMH	Anti-Mullerian hormone
NM_174518	CITED1	Cbp/p300-interacting transactivator, with Glu/Asp-rich carboxy-terminal domain, 1
NM_174304	CYP17	Cytochrome P450, subfamily XVII
NM_174305	CYP19	Cytochrome P450, family XIX, aromatase
NM_001034410	CXCL14	Chemokine (C-X-C motif) ligand 14
NM_182786	FOS	V-fos FBJ murine osteosarcoma viral oncogene homolog
NM_174061	FSHR	Follicle stimulating hormone receptor
NM_175801	FST	Follistatin
NM_001083674	GCLC	Glutamate-cysteine ligase, catalytic subunit
NM_174077	GPX3	Glutathione peroxidase 3
NM_174546	GUCA1A	guanylate cyclase activator 1A (retina)
NM_001102365	HSD17B1	Hydroxysteroid (17-beta) dehydrogenase 1
NM_174555	IGFBP2	Insulin-like growth factor binding protein 2, 36 kDa
NM_174094	INHA	Inhibin, alpha
NM_174363	INHBA	Inhibin, beta A (activin A, activin AB alpha polypeptide)
NM_173950	PlGF	Placental growth factor
NM_174646	PLA2G1B	Phospholipase A2, group IB (pancreas)
NM_174753	PTHLH	Parathyroid hormone-like hormone
NM_173957	RGN	Regucalcin (senescence marker protein-30)
NM_173959	SCD	Stearoyl-CoA desaturase (delta-9-desaturase)
NM_174670	SERPINF2	Serpin peptidase inhibitor, clade F (alpha-2 antiplasmin, pigment epithelium derived factor), member 2
NM_001098036	SLC39A14	Solute carrier family 39 (zinc transporter), member 14
NM_001025326	SRGN	Serglycin
NM_001076470	TMEM20	Transmembrane protein 20
NM_001007813	TNFAIP6	Tumor necrosis factor, alpha-induced protein 6
NM_178317	TRB2	TRB-2 protein
AW315959		13940 MARC 4BOV Bos taurus cDNA 5', mRNA sequence.
AW325368		16365 MARC 4BOV Bos taurus cDNA 5', mRNA sequence.
BE684800		186519 MARC 4BOV Bos taurus cDNA 5', mRNA sequence.
BI536463		393463 MARC 4BOV Bos taurus cDNA 5', mRNA sequence
BI536468		393469 MARC 4BOV Bos taurus cDNA 5', mRNA sequence.
BI537443		397313 MARC 4BOV Bos taurus cDNA 5', mRNA sequence.
BP102158		Transcribed locus
BP103904		BP103904 ORCS bovine liver cDNA Bos taurus cDNA clone ORCS25139 3', mRNA sequence.
BP104736		BP104736 ORCS bovine liver cDNA Bos taurus cDNA clone ORCS26135 3', mRNA sequence.
BP105513		BP105513 ORCS bovine liver cDNA Bos taurus cDNA clone ORCS27141 3', mRNA sequence.
BP107839		BP107839 ORCS bovine utero-placenta cDNA Bos taurus cDNA clone ORCS11248 3', mRNA sequence.
BP108716		Isolate UoG-BovSAGE-UK2 unknown mRNA
BP110155		Testis derived transcript (3 LIM domains)
BP110180		Transcribed locus
BP110819		BP110819 ORCS bovine utero-placenta cDNA Bos taurus cDNA clone ORCS11012 5', mRNA sequence.
BP111150		BP111150 ORCS bovine utero-placenta cDNA Bos taurus cDNA clone ORCS11443 5', mRNA sequence.
XM_614289		Similar to glucocorticoid induced transcript 1 (GLCCI1), mRNA.
XM_864694		Similar to tolloid-like 2, transcript variant 2 (TLL2), mRNA.

Quantitative PCR analysis of representative highly expressed genes in F1 and F2

Figure 2 shows the results of QPCR analysis of the eight representative genes that were highly expressed in F2 (cluster A) compared with F1 (cluster B) in microarray analysis. Messenger RNA expression for chemokine ligand 2 (CCL2), GADD45A, IGFBP5, plasminogen activator urokinase receptor (PLAUR), secreted phosphoprotein 1 (SPP1), selectin P (SELP), tissue inhibitor of matrix metalloprotease-1 (TIMP1) and TSP2 was greater in the F2 than in the F1 (P < 0.05). The results of QPCR analysis of the eight representative genes that were highly expressed in the F1 as compared with the F2 in microarray analysis are shown in Figure 3. The expression of AMH, CYP19, FSHR, glutathione peroxidase 3 (GPX3), PlGF, phospholipase A2 group 1B (PLA2G1B), stearoyl-CoA desaturase (SCD) and tribbles homolog 2 (TRB2) mRNA was greater in the F1 than in the F2 (P < 0.05).

Figure 2

QPCR analysis of representative eight genes (. These genes were highly expressed in F2 (cluster A) compared with F1 (cluster B) in microarray analysis. The expression of mRNA was normalized to the expression of GAPDH measured in the same RNA preparation. The black bar and the white bar indicate the F1 and the F2, respectively. Data are shown as the mean ± SEM. Different letters denote significant differences (P < 0.05).

Figure 3

QPCR analysis of representative eight genes (. These genes were highly expressed in F1 (cluster B) compared with F2 (cluster A) in microarray analysis. The expression of mRNA was normalized to the expression of GAPDH measured in the same RNA preparation. The black bar and the white bar indicate the F1 and the F2, respectively. Data are shown as the mean ± SEM. Different letters denote significant differences (P < 0.05).

Follicular fluid concentrations of E2 and P4 in follicles

We classified follicles into healthy or atretic based on the relative concentrations of E2 and P4 in FF (healthy: E2/P4 ratio ≥1, atretic: E2/P4 ratio <1). From a total of 11 follicles, eight were categorized into healthy while the other three were atretic. Table 4 shows the characteristics of the follicles used in experiment 2. Healthy follicles had higher E2 and lower P4 concentrations in FF than atretic follicles. The E2/P4 ratio in FF was significantly higher in healthy follicles than in atretic follicles.

Table 4

Follicular fluid concentrations of estradiol (E2) and progesterone (P4) in examined follicles used in experiment 2.

Follicle	E2 (ng/ml)	P4 (ng/ml)	E2/P4 ratio
Healthy	180.0 ± 44.9	15.9 ± 15.3	4.3 ± 0.8
Atretic	41.4 ± 5.3*	387.7 ± 121.7*	0.1 ± 0.1*

Values are mean ± SEM

* P < 0.05 vs. healthy follicle

In situ hybridization of representative genes identified in experiment 1

Figure 4 shows mRNA localization for GADD45A, IGFBP5 and TSP2 in healthy and atretic follicles by in situ hybridization. These genes were highly expressed in F2 than in F1 in microarray and QPCR analysis of experiment 1. IGFBP5 mRNA was localized in the GC and theca layer (TL) of atretic follicles but not in healthy follicles (Figure 4E, F, G and 4H). GADD45A (Figure 4A, B, C and 4D) and TSP2 (Figure 4I, J, K and 4L) mRNA were found in both GC and TL of atretic follicles but they were expressed in only GC of healthy follicles. No significant signals were detected with any sense probes (Figure 4B, D, F, H, J and 4L).

Figure 4

Localization of . These genes were expressed more in F2 than in F1 in QPCR analysis. (A, C, E, G, I and K) DIG-labeled anti-sense cRNA probes were used. (B, D, F, H, J and L) DIG-labeled sense cRNA probes were used. Seven-micrometer sections of bovine follicles were hybridized with each probe. GADD45A (A, B, C and D) and TSP2 (I, J, K and L) mRNA were found in both granulosa cells (GC) and theca layer (TL) of atretic follicle, whereas it was localized in only GC of healthy follicle. IGFBP5 mRNA (E, F, G and H) was localized in GC and TL of atretic follicle but not found in healthy follicle. Scale bar = 20 μm.

Localization of AMH, CYP19 and PlGF mRNA in healthy and atretic follicles are shown in Figure 5. These genes were expressed more in the F1 than in the F2 in experiment 1. AMH (Figure. 5A, B, C and 5D) and CYP19 (Figure. 5E, F, G and 5H) mRNA was localized in GC of healthy as well as atretic follicles. PlGF mRNA was found in GC and TL of only healthy follicles but not atretic follicles (Figure. 5I, J, K and 5L). No significant signals were detected with any sense probes (Figure. 5B, D, F, H, J and 5L).

Figure 5

Localization of . These genes were expressed more in F1 than in F2 in QPCR analysis. (A, C, E, G, I and K) DIG-labeled anti-sense cRNA probes were used. (B, D, F, H, J and L) DIG-labeled sense cRNA probes were used. Seven-micrometer sections of bovine follicles were hybridized with each probe. AMH (A, B, C and D) and CYP19 (E, F, G and H) mRNA was localized in granulosa cells (GC) of healthy as well as atretic follicles. PlGF mRNA (I, J, K and L) was found in GC and theca layer of healthy follicle but not atretic follicle. Scale bar = 20 μm.

Discussion

In this study, as expected, hierarchical cluster analysis of the microarray data classified F1 and F2 according to differences in gene expression profiles. In each follicular group, characteristic genes determining their developmental status were expressed. The F1 showed greater expression of genes responsible for enhancement of follicular E2 production than the F2. These genes were gonadotropin receptor (FSHR), steroidogenic enzymes (CYP17, CYP19 and HSD17B1) and inhibin-activin-follistatin system (INHA, INHBA and FST). It is well demonstrated that mRNA expression for FSHR, CYP17, CYP19, INHA and INHBA increases with the progress of bovine follicular development and is greater in DF than SF [8,12,29]. On the other hand, the F2 had greater expression of IGFBP5 mRNA than the F1. IGFBP5 mRNA expression dramatically increased in bovine atretic follicles compared with the DF [26]. Intrafollicular levels of IGFBP proteolytic activity and IGFBPs gene expression are important for bioavailability of free IGF within the follicle and play a crucial role for determining follicular dominance and fate [7,30]. Therefore, we evaluated the F1 were selected DF and the F2 were unselected SF.

Our evaluation of follicular status was confirmed to investigate CYP19 and IGFBP5 mRNA localization in healthy and atretic follicles in experiment 2 using in situ hybridization. CYP19 mRNA was abundantly expressed in healthy follicles but it was also expressed in atretic follicles while IGFBP5 mRNA was detected only in atretic follicles. Both CYP19 and IGFBP5 mRNA is hormonally regulated in bovine follicular cells [31-34], in addition, in situ hybridization is not quantitative and not be as sensitive as QPCR. Thus, small amounts of CYP19 and IGFBP5 mRNA may be detected or regulated in atretic and healthy follicles, respectively. Since we did not perform sample collection at a specific phase of follicular wave in experiment 1, detailed growth profiles of the follicles we used were unclear. However, our results demonstrate that randomly collected follicles can be divided into several groups by similarities of gene expression profiles among the follicles and suggest that gene expression profiles of examined follicles are closely associated with their development status.

Confirmation of microarray data by QPCR analysis successfully identified a set of genes differentially expressed between the F1 and F2. In addition, possible involvement of these genes in follicular development and/or atresia was further demonstrated to investigate mRNA localization in healthy and atretic follicles. The F1 showed greater expression of genes involved in follicular growth and survivability (AMH, PLA2G1B, SCD2 and TRB2) than the F2. High expression of these genes may be closely associated with the establishment and maintenance of follicular dominance. Although the functional role of AMH in antral follicle development is poorly understood, recent studies showed that both intrafollicular AMH concentration and AMH mRNA expression were highest in small antral follicles and then decreased with follicular growth, suggesting the involvement of AMH in bovine follicular recruitment and/or selection [16,35,36]. A recent study showed a significant decrease of AMH mRNA expression in late atretic follicles compared with healthy follicles [36], which is consistent with our present result. Furthermore, in agreement with previous studies [16,37], our in situ hybridization study showed that AMH mRNA was localized in only GC. High expression and clear localization of AMH mRNA in the GC of healthy large follicles implies that this growth factor has a plausible effect on the development of DF after follicular selection as well as recruitment.

PLA2 enzymes including PLA2G1B hydrolyze fatty acids from the sn-2 position of phospholipids with concomitant formation of lysophospholipids, which serve as precursor for lipid mediators such as lysophosphatidic acid (LPA) [38,39]. Released LPA has diverse biological activities including cell proliferation and differentiation, suppression of apoptosis and cytoskeleton modulation in reproductive tissues [39]. Because Diouf et al. reported that PLA2G1B mRNA expression in the GC of bovine preovulatory follicle decreased after hCG injection [40], PLA2G1B may mainly contribute to generation of LPA during DF growth before the LH surge.

SCD is a rate-limiting enzyme that catalyzes the synthesis of monounsaturated fatty acids, mainly palmitic and oleic acid [41]. Consistent with our result, SCD mRNA expression in bovine follicles was found to be highest in GC of DF than in cohort follicles before selection or SF [11,17]. Expression of SCD2 is hormonally regulated during follicular development because both SCD2 mRNA and protein expression in rat large follicles were stimulated by gonadotropin and IGF-I treatment [42]. Increase of monounsaturated acids synthesize by SCD2 activation during DF growth may be required to maintain membrane fluidity [43] and a major lipid reserve of oocytes [44].

Members of the TRB family including TRB2 interact and modulate the activity of mitogen-activated protein kinase (MAPK) which regulates cell proliferation, differentiation, apoptosis and survival [45]. These MAPK cascade protein levels were greater in DF than in SF [46]. In addition, it has been reported that TRB2 mRNA was constantly expressed between bovine small follicles and DF [11]. These studies and our present result suggest the potential role of TRB2 in the regulation of MAPK cascades in the growing DF.

The F2 are characterized by high expression of the genes involved in immune reaction (CCL2, SELP and SPP1). In bovine follicles, expression of CCL2 and SPP1 mRNAs and SELP protein was up-regulated in association with follicular development and ovulation [16,18,47]. Our results raise the possibility that these immune-related genes may be involved in bovine follicular atresia as well as follicular development and ovulation. Both CCL2 and SELP mediate induction of leukocyte emigration into extravascular inflammatory sites [48]. Although SPP1, also known as osteopontin, has diverse physiological functions, one of its potent actions is recruitment and retention of macrophages and T cells to inflamed sites [49]. Since number of leukocytes, lymphocytes and activated macrophages are increased in atretic follicles [50], CCL2, SELP and SPP1 participate in the regulation of inflammatory processes during follicular atresia to attract white blood cells.

The F2 are also characterized by high expression of genes regulating tissue remodeling (TIMP1 and PLAUR). Both plasminogen activator (PA)-plasmin and matrix metalloproteinase (MMP) systems play a crucial role in the degradation and remodeling of extracellular matrix associated with follicular development, ovulation and atresia [51]. Urokinase PA (uPA) receptor is a specific cell surface receptor for uPA and its principal role is to localize pericellular plasmin activity to induce extracellular matrix degradation [52]. A previous study showed that there was no difference in mRNA expression levels of uPA between non-atretic and atretic bovine follicles, whereas atretic follicles had lower FF protein level and mRNA expression of a PA inhibitor and higher FF plasmin activity than non-atretic follicles [53]. Therefore, the follicular PA-plasmin system may be primarily regulated by changes in their receptors and inhibitors' expressions. TIMP-1 is an intrinsic inhibitor of MMPs and preferentially binds to MMP-9 [54]. A previous study demonstrated that MMP-9 proenzyme (proMMP-9) protein in FF was detected only in atretic follicles but not in healthy follicles in cattle [55]. Atretic follicles may balance MMPs and TIMP-1 in response to an increase in proMMP-9 to control extracellular matrix degradation by MMP-9.

In the present study, we identified differential expression of two anti-apoptosis factors (GADD45A and GPX3) between the groups. Expression of GADD45A mRNA was greater in the F2 than in the F1. GADD45A controls cell cycle arrest, apoptosis induction and DNA damage repair in response to DNA damaging agents and growth arrest signals of genotoxic stress [56]. Our result implies that the atretic follicles suffered more severe DNA damage than healthy follicles. Indeed, we found in experiment 2 that the atretic follicles expressed GADD45A mRNA in both GC and TL whereas the healthy follicles expressed it only in GC. This result supports our microarray result and suggests an increase in the requirement of GADD45A activity for progression of apoptotic cell death in GC and TC during follicular atresia. On the other hand, GPX3 mRNA was found to be more greatly expressed in the F1 than in the F2. Glutathione peroxidase protects cells against oxidative damage to catalyze the reduction of free hydrogen peroxide and other hydroperoxides [57]. High oxidative stress can trigger apoptosis of follicular cells and induce atresia [58]. In cultured swine GC, GPX3 mRNA expression was upregulated by FSH treatment [59]. Thus GPX3 could prevent cell apoptosis from oxidative stress during growth of the healthy follicles. It is likely that follicular oxidative stress-response enzymes are expressed in a stage-dependent manner since mRNA expression of other anti-oxidative stress enzymes in bovine GC was increased in atretic DF than in healthy DF [60].

Providing a sufficient blood supply is essential for follicular growth [61,62]. A morphological study has demonstrated that bovine healthy DF has a high density and well developed capillaries in TL whereas atretic follicles has sparse and poorly developed capillaries [63]. In the present study, two genes regulating angiogenesis, PlGF and TSP2, were differentially expressed between the groups. PlGF was expressed most in F1 than in F2 and localized in both GC and TL of healthy follicles but not detected in atretic follicles. PlGF is a member of the vascular endothelial growth factor family and stimulates the proliferation of endothelial cells and supports angiogenesis [64,65]. Therefore, PlGF may contribute to follicular thecal angiogenesis via paracrine/autocrine action in healthy follicles as well as other angiogenic factors. In contrast to PlGF, TSP-2, a member of the TSP family, acts as a potent inhibitor of angiogenesis and induces endothelial cell apoptosis [66]. In experiment 1, TSP2 mRNA expression was greater in the F2 than in the F1. TSP2 mRNA level in the bovine follicles decreased in accordance with an increase in follicular diameter [67]. The same authors also showed that TSP protein was localized in both GC and TC of small follicles but in only in the GC of large follicles [67]. We demonstrated in experiment 2 that TSP2 mRNA was localized in both GC and TL of atretic follicles while it was expressed in only GC of healthy follicles. Recent studies have demonstrated that mRNA and protein expression of TSP-1, another antiangiogenic TSP, is upregulated in primate GC during progression of follicular atresia [68] and TSP1 mRNA abundance is decreased by IGF-I treatment in cultured porcine GC [69]. Thus, we speculate that TSP2 mRNA expression is maintained at high levels in follicular cells of atretic follicles whereas it decreases in healthy follicles. Highly expressed TSP2 mRNA in the follicles could negatively influence their angiogenesis. It may cause an insufficient supply of substrates essential for follicular growth, thereby affecting follicular hormone production and cell proliferation, and, as a result, inducing atresia.

Conclusion

Microarray and QPCR analysis enabled us to classify uncharacterized bovine follicles and to evaluate their representative follicular status according to differences in global gene expression profiles. Our present study demonstrates that the expression of stage-specific genes in F1 and F2 may be closely associated with follicular growth and atresia. Several genes identified in this study will provide information on the genomic actions of intriguing candidates for the determinant of bovine follicular development.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

KGH participated in the design of the study, collected the materials, carried out all experiments and drafted the manuscript. KU collected the materials, carried out the microarray experiments and analysis, and helped to carry out QPCR and in situ hybridization. MH was responsible for all animal care, collected the materials and carried out the microarray experiments. TT supervised the study, collected the materials and helped to draft the manuscript. All authors read and approved the final manuscript.