Profiling of Luteal Transcriptome during Prostaglandin F2-Alpha Treatment in Buffalo Cows: Analysis of Signaling Pathways Associated with Luteolysis

In several species including the buffalo cow, prostaglandin (PG) F2α is the key molecule responsible for regression of corpus luteum (CL). Experiments were carried out to characterize gene expression changes in the CL tissue at various time points after administration of luteolytic dose of PGF2α in buffalo cows. Circulating progesterone levels decreased within 1 h of PGF2α treatment and evidence of apoptosis was demonstrable at 18 h post treatment. Microarray analysis indicated expression changes in several of immediate early genes and transcription factors within 3 h of treatment. Also, changes in expression of genes associated with cell to cell signaling, cytokine signaling, steroidogenesis, PG synthesis and apoptosis were observed. Analysis of various components of LH/CGR signaling in CL tissues indicated decreased LH/CGR protein expression, pCREB levels and PKA activity post PGF2α treatment. The novel finding of this study is the down regulation of CYP19A1 gene expression accompanied by decrease in expression of E2 receptors and circulating and intra luteal E2 post PGF2α treatment. Mining of microarray data revealed several differentially expressed E2 responsive genes. Since CYP19A1 gene expression is low in the bovine CL, mining of microarray data of PGF2α-treated macaques, the species with high luteal CYP19A1 expression, showed good correlation between differentially expressed E2 responsive genes between both the species. Taken together, the results of this study suggest that PGF2α interferes with luteotrophic signaling, impairs intra-luteal E2 levels and regulates various signaling pathways before the effects on structural luteolysis are manifest.


Introduction
Corpus luteum (CL) is a transient endocrine structure formed from the remnants of the ovulating follicle. Through secretion of progesterone (P 4 ), it plays a pivotal role in the control of reproduction in mammals. During non-pregnant cycles, regression of CL precedes initiation of new reproductive cycle to allow reovulation and another chance for conception to occur [1]. The structure and function of CL are controlled by intricate interplay between luteotrophic and luteolytic factors. In several species including bovines, the influence of lytic factors dominate over luteotropic factors in the control of luteal function. Of the many factors, prostaglandin (PG) F 2a molecule is recognized as the physiological luteolysin responsible for luteal regression in mammals [2][3][4].
PGF 2a or its synthetic analogues initiate a complex cascade of events within the CL leading to inhibition of steroidogenesis resulting in rapid fall in circulating P 4 levels (functional luteolysis) [5][6][7] and initiation of apoptosis (structural regression) in the bovine species [6,[8][9][10]. In several species, many of the endocrine, biochemical and morphological events observed during PGF 2ainduced luteolysis have been reported to be broadly similar [11][12][13][14][15][16][17][18][19][20]. Disruption of transport of cholesterol due to inhibition of StAR expression has been implicated as the cause of inhibition of steroidogenesis after PGF 2a treatment under both in vivo and in vitro conditions [21,22]. However, in a more recent study, an initial increase in StAR expression post intrauterine PGF 2a treatment has been observed [23], suggesting factors other than StAR are responsible for inhibition of steroidogenesis. Studies on molecular analysis of structural regression of bovine CL leading to apoptosis has previously been reported by others [24,25] and us [6,9], but the actual mechanisms responsible for functional luteolysis has not been systematically examined. Recent studies on transcriptome analysis of CL tissues have expanded our understanding of effects of PGF 2a during luteolysis or refractory state of CL to PGF 2a effects [26][27][28].
In contrast to a well-established role for pituitary LH in the control of CL function in primates and rodents [29,30], the role of LH in the regulation of CL function in bovines has not been clearly established. In cattle, an earlier study had suggested a role for LH in CL function [31], while other studies have suggested a minor role for LH especially during pregnancy [32][33][34]. To understand mechanisms involved in the initial loss of CL function in response to PGF 2a treatment, it is essential to examine molecules downstream of LH/CG receptor (LH/CGR) activation during PGF 2a -induced luteolysis. The role of intraluteal molecules during luteolysis forms another crucial area for examining the PGF 2a actions during the initial period of PGF 2a treatment [35]. In bovines, even though expression of CYP19A1 that codes for aromatase enzyme responsible for converting androgens to estrogens, gets down regulated after ovulation, its expression is restored albeit at low levels in the developed CL and is never the less responsible for low levels of intraluteal estradiol-17b (E 2 ) [36] in the buffalo cow. Since luteal tissue expresses both a and b receptors [37], it remains to be determined whether biosynthesis of E 2 and its receptor signaling is disrupted following PGF 2a treatment contributes to the loss of luteal function. Interestingly, E 2 is the primary luteotrophic hormone in rabbits and experimental withdrawal of E 2 has been shown to inhibit luteal function and activate apoptotic pathways [38,39].
In view of the above, the following objectives were addressed in the present study: 1) analysis of global gene expression changes in CL of buffalo cows administered with luteolytic dose of PGF 2a treatment, 2) examination of effects of PGF 2a on molecules downstream of LH receptor signaling and 3) examination of the effects of PGF 2a on CYP19A1 expression, E 2 biosynthesis and expression of E 2 responsive genes in the CL.

Ethics Statement
All procedures in animals were approved by the Institutional Animal Ethics Committee, Indian Institute of Science, Bangalore, India.

Experimental protocol, blood and CL collection schedule
Non-lactating buffalo cows (Bubalus bubalis, Surthi breed) aged 5-6 years, with a history of normal cyclicity were monitored daily for onset of behavioral signs of estrus such as bellowing, restlessness and mucous discharge from the vulva. The estrous cycles were synchronized with two injections of PGF 2a and the day of onset of estrus was designated as day 0 of estrous cycle. To verify the presence of functional CL, blood samples were collected on days 3 to 7 of the cycle for monitoring circulating P 4 concentration. In this experiment, Juramate (500 mg i.m.) was administered on day 11 of estrous cycle and blood samples (n = 5-6 animals) were collected immediately before and at various time points post PGF 2a treatment. CL was collected at slaughter from untreated control animals (0 h), as well as from animals at 1, 2, 3, 6 and 18 h post PGF 2a injection (n = 3-4 animals/time point). For serum E 2 estimation, blood sampling has also been done at earliest time points, i.e. at 1 and 2 h post PGF 2a treatment. CL tissues were also collected at varied points throughout an estrous cycle with day 5, 11, 14, 17 and 20 CL being designated as early (E), mid (M), mid-late (ML), late (L) and very late (VL) CL, respectively. Ovaries containing CL were collected and washed in sterile ice cold PBS and transferred into Dulbecco's Modified Eagles Medium supplemented with penicillin (500 U/ml) and streptomycin (50 mg/ml) and transported to the laboratory on ice within 30 min of collection. Following retrieval of CL from the ovary, CL was cut into ,50 mg pieces and 3-4 pieces were placed in sterile cryovial and flash frozen in liquid nitrogen before storing at 270 uC for various analyses. The progesterone analysis data and a portion of CL tissue used for different analysis have been reported recently [40].

Hormone assays
For the measurement of luteal P 4 and E 2 , wet weight of one of the CL pieces was determined before homogenization in a known volume of 1X PBS solution. Undiluted CL tissue lysate (E 2 ) or at a dilution of 1:500 (P 4 ) was used in the assay. Hormone extraction from lysate was carried out using diethyl ether and reconstituted in GPBS at 37uC for 1 hour in an incubator shaker. Luteal and serum P 4 and E 2 concentrations were determined by radioimmunoassay method as reported previously [41]. The sensitivity of the assay was 10 pg/tube and 3.9 pg/tube for P 4 and E 2 , respectively and the inter-and intra-assay coefficients of variation were , 10%.

Isolation of genomic DNA and fragmentation analysis
Genomic DNA was extracted from individual CL using standard phenol-chloroform extraction method. The isolated genomic DNA was precipitated, dissolved and quantitated spectrophotometrically for fragmentation analysis as described elsewhere [6]. Briefly, genomic DNA from different CL tissue samples were analyzed for quantitation of low molecular weight DNA fragmentation using protocol described previously [42,43] with few modifications. Genomic DNA (4 mg) was labeled with 50 mCi of [a 32 P] dCTP by incubation at 37uC for 1 h with 10 units of terminal deoxynucleotidyl transferase enzyme. The DNA samples were then resolved on a 2% agarose gel and transferred onto a nylon transfer membrane using the alkaline (0.4 M NaOH and 1 M NaCl) buffer with capillary blotting transfer method for 16 h. The membranes post transfer were cross linked using 1200 mJoules of UV for 1 min in Stratalinker. The membranes were covered with plastic wrap and were exposed to an intensifying screen at room temperature and the low molecular weight labeled DNA signals were quantitated densitometrically using a PhosphorImager (Typhoon 9210; Amersham Biosciences) with exposure of 6 h or longer.

RNA isolation
Total RNA was extracted from control and PGF 2a treated CL tissues using TRI Reagent according to the procedure as reported previously [41]. RNA was quantitated spectrophotometrically using ND-1000 (NanoDrop, Thermo Scientific, Wilmington, DE). The quality and quantity of RNA were determined by electrophoresis on a 2% (w/v) formaldehyde agarose gel along with RNA samples of known concentration, and A 260 : A 280 ratio was .1.8.
Additionally, the quality of RNA samples was assessed using Agilent Bioanalyzer 2100 and RNA samples with RNA integrity number $9 were used for microarray analysis.

Quantitative real time PCR (qPCR)
The qPCR analysis was carried out as described previously from the laboratory [44]. The cDNA samples equivalent to 10 ng of total RNA were subjected to validation analysis on Applied Biosystems 7500 Fast Real Time PCR system with SDS v 1.4 program employing Power SYBR green 2X PCR master mix. The details of primers employed along with the annealing temperature and expected product size are provided in Table S1. Analysis of expression of each gene included a no template control (NTC) and generation of a dissociation curve. Expression level of individual gene was normalized by using L19 expression level as calibrator (internal control) for each cDNA sample. The relative expression of each gene was determined using the DDC t method which calculates the fold change in gene expression using the formula: Fold change = 2 2DDCt , where C t = Threshold cycle i.e. the cycle number at which the relative fluorescence of test samples increases above the background fluorescence and DDC t = [C t gene of interest (unknown sample) -C t of L19 (unknown sample)] -[C t gene of interest (calibrator sample) -C t of L19 (calibrator sample)]. PCR for each sample was set up in duplicates and the average C t value was used in the DDC t equation.

Immunoblot analysis
Immunoblot analysis of the total protein lysates from CL tissues was carried out as per the procedure reported previously [6].

Protein kinase A assays
Protein kinase A assays were performed according to the manufacturer's instructions provided with the kit and as described in detail previously [45].

Microarray target preparation and hybridization
Transcriptome analysis of CL tissues was performed on GeneChip Bovine Genome Array (Affymetrix, Santa Clara, CA) that contains 24,072 probe sets representing ,23,000 bovine transcripts and few housekeeping/control genes. Use of heterologous array hybridization screening for gene expression changes of closely related species has been validated for many species [46]. Affymetrix Bovine GeneChip platform is very efficient and can be easily extended to other species for which genetic sequences are available [47,48]. The water buffalo and domestic cattle, both belonging to the Bovidae family, are closely related. Moreover, use of robust multiarray average (RMA) algorithm for background correction, normalization across arrays and summarization is well suited for cross species microarray analysis, since the software considers hybridization from perfect match probes and not the mismatch probes [47]. More recently, microarray analysis has been performed on buffalo granulosa cells employing GeneChip Bovine Genome Array and analysis data indicated excellent hybridization [49]. The entire procedure of microarray target preparation, hybridization and scanning were performed at Center for Integrated Biosystems, Utah State University, Utah by employing standard Affymetrix protocols (Chip Expression Analysis Technical Manual rev. 3,2001). Further, scanning was performed using Affymetrix Genearray Scanner 3000 series.

Analysis of microarray data
Total RNA isolated from control and PGF 2a -treated CL tissue was labeled and hybridized to Affymetrix GeneChip Bovine Genome Arrays as per the manufacturer's recommendations. 3 Affymetrix GeneChip Bovine Genome Arrays were used for RNA from individual CL as biological replicates (n = 3/time point; 0, 3, 6 and 18 h post PGF 2a treatment). A total of 12 arrays were used. A detailed description of procedures and subsequent generation of processed image files of microarray analysis have been reported previously [44,49]. The microarray procedure and the data analysis were performed as per Minimum Information About Microarray Experiments (MIAME) compliance. The raw data and the completed analysis of microarray data files have been deposited at NCBI's Gene Expression Omnibus, http://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE27961 and are accessible through GEO series accession number GSE27961. The processed image files (.cel) were normalized across arrays using RMA software [50] and log-transformed (base 2), which allowed direct comparison of probe set values between all samples used in the experiment normalization. GeneSifter (VizX Labs; Seattle, WA, data not shown), GeneSpring GX (Agilent Technologies, Inc. Santa Clara, CA) and R/Bioconductor (FHCRC labs; Seattle, WA) microarray expression analysis softwares were used to identify differentially expressed transcripts. The differential expression of genes was calculated by averaging the normalized samples and performing a pairwise analysis. Statistical significance for differentially expressed genes among control and treatment groups was determined using Student's t-test (two tail, unpaired) employing Benjamini and Hochberg correction factor for false discovery rate [51]. As many of the original annotations for the Affymetrix Bovine Genome Chip have been reported to be erroneous [52], a computational tool, AffyProbeMiner, was employed to redefine the chip definition files taking into account the most recent genome sequence information [53]. AffyProbe-Miner analysis results suggested that the original Affymetrix gene annotation were not compromised and the identified differentially expressed genes were transcript consistent and did not hybridize to multiple transcripts. Analysis of the data by Bioconductor analysis tool employing $2 fold change cut-off and statistical filters provided a number of differentially expressed genes. The linear regression analysis was performed on fold change expression values for selected differentially expressed genes obtained from qPCR and microarray analyses and a statistically significant (p, 0.05) correlation between the two analyses was determined as reported previously [44,49]. The list of top 15 differentially expressed (both up-and down-regulated) genes at different time points post PGF 2a administration are represented in Table S2, S3, S4, S5, S6, S7. The differentially expressed genes were clustered by hierarchy analysis by GeneSpring analysis for all the probe sets at each time point post PGF 2a administration and represented as dendrograms or heat map in Figure S1. Further, to examine the molecular function and genetic networks, the microarray data was analyzed using Ingenuity Pathways Analysis tool (IPA version 8.7, Ingenuity Systems Inc., Redwood City, CA, USA; http://www. ingenuity.com), a web-based software application that enables identification of biological mechanisms, pathways, and functions from the differentially expressed genes.

Statistical analysis
The data were expressed as mean 6 SEM. Statistical evaluation of mean differences of E 2 and P 4 , protein blots and qPCR fold change in expression of genes among different treatment groups were analyzed using one-way ANOVA test followed by the Newman-Keuls multiple comparison tests (PRISM GraphPad; GraphPad Software Inc., San Diego, CA). A p-value of ,0.05 was considered to be significant.

Results
Effects of PGF 2a on circulating P 4 , StAR expression and biochemical integrity of DNA (DNA laddering) in the CL tissue Following administration of PGF 2a , circulating P 4 levels declined significantly within 1 h and continued to be lower at 2, 3, 6 and 18 h time points examined (p,0.05; Figure S2A). The luteal P 4 level too declined significantly by 3 h and further declined by 6 and 18 h post PGF 2a treatment (p,0.05; Figure   S2A). Analysis of StAR protein pattern paralleled P 4 decline at alltime points examined ( Figure S2B). Analysis of low molecular weight DNA fragments indicated that although a slight increase in low molecular weight fragments could be visualized in CL tissues at 6 h post PGF 2a treatment compared to untreated CL tissue, but a clear increase of laddering pattern of DNA was observed in CL tissues collected from animals, 18 h post PGF 2a treatment ( Figure  S2C). These results verified that CL tissue underwent a time related rapid functional loss followed by evidence of structural loss after initiation of PGF 2a treatment.
Genome wide analysis of differentially expressed genes in the CL of buffalo cow post PGF 2a treatment The mechanism of PGF 2a -induced luteolysis is a complex process involving expression changes in genes associated with regulation of steroidogenesis, apoptosis, expression changes of gene products of several factors such as enzymes involved in PGF 2a biosynthesis and cross talk with other signaling pathways that include luteotrophic factors. In the present study, efforts were made to identify temporal changes in the profile of transcriptome of CL tissue in response to PGF 2a treatment. The microarray analysis provided a set of differentially expressed genes based on various high stringency statistical filters such as a t-test with p, 0.05 and multiple hypothesis testing (Benjamini and Hochberg comparison test) to eliminate the false positives. Further, the differentially expressed genes that passed the statistical filters were classified based on whether the identified differentially expressed gene was present in all the samples examined. Another parameter employed was a sliding scale fold change cut-off filter of $2 (except in identification of E 2 target genes, in which changes ,1.5 fold was considered for analysis) with some degree of confidence in our candidate list to narrow down the list.
The whole transcriptome analysis data presented as Venn diagrams in Figure 1A and B revealed a total of 127, 774 and 939 genes differentially expressed in CL tissues collected at 3, 6 and 18 h post PGF 2a treatment, respectively. Also, the distribution of nearly 40 differentially expressed genes which were common across different time points post PGF 2a treatment were identified as up regulated (20 genes), down regulated (16 genes) or up regulated at one time point, but down regulated at another time points (4 genes) are shown in the Venn diagrams ( Figure 1A&B). It was observed that most cluster of genes appeared to be modulated in a synchronous or synergistic way following PGF 2a treatment, i.e. a higher percentage of the genes regulated at the 3 h time point had a similar regulation at 6 and 18 h time points. qPCR validation of differentially expressed genes from the microarray analysis data Figure 2 shows fold changes in mRNA expression of 20 differentially expressed genes of both microarray and qPCR analyses. Of the 20 differentially genes analyzed by qPCR analysis, the results of 15 genes (CYP19A1, TGFBR3, AGTR1, VEGFR2/ KDR, VEGFR3/FLT4, AREG, PTGS2, TIMP3, BMP2, TGFBR1, PTGFR/FPR, GHR, IGF1R, IGF2R and StAR) were also utilized for validation of microarray data (i.e. for validation and correlation between microarray and qPCR analyses was performed and the results presented in Figure 1C). The goodness of fit analysis showed that the fold change in qPCR analysis of mRNA expression data corroborated well with the microarray analysis data ( Figure 1C). At 3 h post PGF 2a administration the correlation coefficient 'r' value was 0.69 (p,0.003, Figure 1C), while the 'r' value at both 6 and 18 h time points was 0.89 (p, 0.0001, Figure 1C). Many of the differentially expressed genes were selected for qPCR analysis to further characterize their expression in relation to their function such as signaling pathways (for pathway activity analysis), key enzymes/proteins associated with the steroidogenesis, growth factor activation status, prostaglandin biosynthesis, activation of cytokine signaling and genes associated with extracellular matrix (Figure 2). Expression analysis data of select genes presented in Figure 2 also provide information as to the time course changes in expression post PGF 2a treatment. Comparison of expression data of select genes between microarray and qPCR analyses at 3 h time point revealed that the fold expression changes by microarray analysis was lower compared to other time points, but fold change in expression of same genes by qPCR analysis was higher (  (Tables S2, S3, S4, S5, S6, S7).

Pathway analysis of differentially expressed genes employing Ingenuity pathway analysis (IPA) software
The differentially expressed genes identified in CL tissues post PGF 2a treatment were classified into different functions. The IPA analysis revealed that 85-90% of the differentially expressed genes were observed to be function or pathway eligible. Gene classification according to the canonical signaling pathways indicated that many of the differentially expressed genes were associated with process and/or function of IGF-1 signaling, steroidogenesis, chemokine signaling, prolactin signaling, cellular growth and proliferation, extracellular matrix modulation and apoptosis. Further, performing IPA on the differentially expressed genes, 42 networks were identified, out of which 16 networks had 20 or more focus genes in each network. Figure 3A&B shows different canonical signaling pathways and genes identified in Network for 0 vs. 3 h time point (score 48, 21 focus molecules, Figure 3B). The network of genes are associated with growth factor signaling, cytokine signaling, steroidogenesis and extracellular matrix The classification of canonical pathways and network of genes for 6 and 18 h time points post PGF 2a treatment are provided in Figure S3 and S4, respectively.

Effect of PGF 2a administration on expression of orphan nuclear transcription factors associated with regulation of steroidogenesis in the CL
Since the effects of PGF 2a on PKC activation pathway are well recognized, its effects on steroidogenesis at the molecular and biochemical levels were examined. Both microarray and qPCR analyses of mRNA expression of orphan nuclear transcription factors NR5A1/SF-1 and NR5A2/LRH-1 involved in the regulation of expression of steroidogenic genes showed down regulation, excepting at 3 h time point for SF-1, at all-time points post PGF 2a treatment (Figure 2). Immunoblot analyses for protein of both these genes also showed down regulation similar to their mRNA expression patterns ( Figure S5A&B).

Effect of PGF 2a treatment on LH/CGR activation
Although LH/CGR mRNA expression by microarray and qPCR analyses indicated marginal decrease at 6 and 18 h post PGF 2a treatment (Figure 2), but the protein expression was significantly (p,0.05) lower post PGF 2a treatment ( Figure 4A). To further examine whether PGF 2a treatment affected signaling molecules downstream of LH/CGR activation, PKA activity and CREB phosphorylation levels were determined and the results indicated decreased PKA activity ( Figure 4B) and lowered pCREB levels ( Figure 4C) at 3 and 18 h time points examined. Also, a decreased phosphorylated Akt ( Figure 4D) and FKHR ( Figure 4E) levels post PGF 2a treatment suggesting inhibitory effects on survival signaling pathway.
Effect of PGF 2a treatment on circulating and luteal E 2 levels and expression of E 2 responsive genes in the CL Of the many steroidogenic genes which were affected by PGF 2a treatment, mRNA expression of CYP19A1 gene that codes for aromatase enzyme responsible for conversion of androgens to estrogens was down regulated (see Figure 2) at all the time points examined. Figure 5A shows circulating and luteal E 2 levels before and after PGF 2a treatment. Due to variations in circulating E 2 concentrations between animals, the pre-treatment concentration of E 2 in each animal was set as 100% and values post PGF 2a treatment were expressed in relation to 100%. As can be seen from Figure 5A, E 2 concentrations decreased significantly within 1 h after PGF 2a treatment and continued to be lower until 6 h. Since Transcriptome Profiling in Buffalo Cows CL PLOS ONE | www.plosone.org % decrease in E 2 was not seen after 6 h post treatment this was perhaps due to contribution of E 2 from increased growth and development of ovarian follicles. To further validate the effects of PGF 2a on E 2 synthesis and secretion, luteal E 2 levels were monitored and the results are presented in Figure 5A. As can be seen from the figure, luteal E 2 levels continued to decline post PGF 2a treatment.
Luteal expression of E 2 receptors, a and b, during the luteal phase of estrous cycle and their expression patterns post PGF 2a treatment are presented in Figure S6. Since estrogen has both local and systemic effects, expression of estrogen responsive genes in the CL tissue was examined. For examining the effects of PGF 2a treatment on expression of E 2 responsive genes in the CL tissue, in addition to the classical E 2 responsive genes [54], genes identified as target of E 2 action by PCR array human estrogen signaling (Qiagen, SABiosciences, PHAS-005A) as well as contents of estrogen responsive gene data base (ERGDB) were included. In all, a total of 89 genes identified as potential E 2 responsive genes were considered for analysis. From the list of differentially expressed genes obtained by microarray analysis in CL tissue at different time points post PGF 2a treatment, as many as 57 out of 89 E 2 responsive genes (both up and down regulated) were identified to be differentially expressed at least 1 fold or more. The data is represented as Venn diagram in Figure 5B. The list of E 2 responsive genes for 3 h time point is shown in Tables S8, S9. The list for other time points is not shown since expression of large number of E 2 responsive genes were commonly regulated at alltime points. The list of E 2 responsive genes commonly up and down regulated at all-time points post PGF 2a is shown in Tables S10, S11. Twenty one of the 43 up regulated genes were observed to be commonly regulated at all-time points examined, while 37 genes identified as down regulated, 14 were found to be commonly regulated at all-time points post PGF 2a treatment ( Figure 5B).
The buffalo CL much like the cattle is not regarded as a site of high E 2 production since much of the circulating E 2 in bovines is contributed by waves of follicular growth occurring throughout the estrous cycle [55,56]. On the other hand, in species such as primates, CL is the predominant source of E 2 production throughout the luteal phase period and circulating E 2 contributed by the functional CL is considerable [57][58][59]. To further validate our hypothesis that decreased estrogen production following inhibition of CYP19A1 expression might adversely affect CL function and its eventual survival, the previously published microarray data of the differentially expressed genes from the CL tissues of macaques receiving PGF 2a treatment for 24 h ( [44], GEO accession number GSE8371) was mined for E 2 responsive genes (89 genes identified as potential E 2 responsive genes, see above) for purposes of comparing the number of E 2 responsive genes that were differentially expressed in CL of macaques (in which E 2 is secreted in higher amounts) to that of the buffalo cow CL (in which E 2 secretion is low). The mined data comprising up and down regulated genes of CL from macaques at 24 h time point (the single time point for which microarray analysis available) vs. E 2 target genes identified at 3, 6 and 18 h time points of buffalo CL tissue post PGF 2a treatment are shown in Figure 6 (A-C). Forty seven out of the 89 genes regarded as E 2 responsive genes were identified as up regulated in the monkey CL, 24 h post PGF 2a treatment and 19 of the 47 genes were also observed to be commonly regulated in CL tissues of buffalo cows post 3 h PGF 2a treatment ( Figure 6A). Remarkably, the number of down regulated genes identified in CL tissues at 24 h post PGF 2a treatment in macaques and at 6 and 18 h post PGF 2a treatment in buffalo cows were similar and 14 of the genes were found to be commonly regulated between the two species ( Figure 6B&C). The list of E 2 responsive genes commonly up and down regulated from macaques at 24 h time point vs. genes identified at 3, 6 and 18 h time points of buffalo CL tissue post PGF 2a treatment is shown in Tables S12, S13, S14, S15, S16, S17.

Effect of PGF 2a administration on expression of different regulatory molecules associated with cell survival and apoptosis
To examine whether down regulation of E 2 receptor expression results in activation of apoptotic changes in the luteal tissue, expression of Bax and Bcl-2 genes were determined in the luteal tissue post PGF 2a treatment. Figure 7A&B illustrates expression of (both mRNA and protein) Bcl-2 family members, Bax (proapoptotic) and Bcl-2 (anti-apoptotic) during PGF 2a -induced luteolysis. The Bcl-2 mRNA expression was unchanged post PGF 2a treatment, but a significant increase in Bax mRNA expression was seen at 18 h post PGF 2a treatment ( Figure 7A). Similarly, at the protein levels, Bcl-2 protein levels did not change, but a significant increase was seen for Bax protein at 18 h time point post PGF 2a treatment ( Figure 7B). As indices of activation or inhibition of mitochondrial permeability to apoptogenic molecules, changes in Bax and Bcl-2 mRNA expression and protein levels were expressed as Bax/Bcl-2 ratios. A significant increase in both Bax mRNA and protein levels post PGF 2a treatment was observed ( Figure 7A&B). As can be seen from Figure 7C, Bax/ Bcl-2 ratio for both mRNA and protein levels, increased significantly by 18 h post PGF 2a treatment. Also, a significant increase in Bax/Bcl-2 protein was observed by 6 h post PGF 2a treatment suggesting increased mitochondrial permeability to apoptogenic molecules ( Figure 7C).
It has been reported that PI3K activity is essential for the activation of Akt for E 2 actions [60,61]. In present study, although a significant change in the expression of PI3K p85 mRNA by qPCR analysis was not observed ( Figure 7D), but a significant decrease in phospho levels of PI3K p85 was observed at all-time points post PGF 2a treatment ( Figure 7D).

Schematic diagram illustrating a model for interaction amongst different signaling pathways during PGF 2ainduced luteolysis in buffalo cows
Extensive studies carried out previously by several groups have identified different molecules of the signaling pathways for LH and PGF 2a actions in the luteal tissue of various species [62][63][64][65][66]. In the present study, some of the signaling molecules of both these pathways were analysed as also genes associated with steroid biosynthesis, cell survival and apoptosis and the data for expression of various genes and their protein levels in the luteal tissue at different time points post PGF 2a treatment are presented in Figure 8. A schematic diagram of LH and PGF 2a signaling pathways and E 2 receptor signaling pathways and their likely interactions is depicted in Figure 8. In the proposed model, an association between cell survival and E 2 action could be suggested for regulating luteal function in the bovine luteal tissue. The down regulation of E 2 signaling by PGF 2a as suggested by the data indicates activation of apoptotic pathway (Figure 8).

Discussion
In bovines, during non-conception reproductive cycles, the onset of spontaneous luteolysis is initiated by increased pulsatile secretion of PGF 2a from the uterus [62,67]. Several recent studies have compared expression changes during periods of refractoriness and responsiveness of CL to PGF 2a treatment [23,28,68] and few studies have previously reported global gene expression changes during PGF 2a induced luteolysis [23,26]. In the present study, employing high throughput screening technique, genomewide gene expression changes in the CL of buffalo cow during PGF 2a -induced luteolysis were determined. The primary objective of this study was to examine temporal expression changes during early periods of PGF 2a treatment in which functional loss in CL will be manifest. The global transcriptome changes post PGF 2a treatment revealed temporal increases in the expression of a number of differentially expressed genes at 3, 6 and 18 h time points examined, and of these several of them were immediate early genes. The analysis of differentially expressed genes during induced luteolysis also revealed cluster of genes associated with steroidogenesis, angiogenesis, cell survival and apoptosis and these have major bearing on initiation of functional luteolysis and apoptotic processes. Further, comparison of differentially expressed genes data from the present study with other studies following administration of PGF 2a treatment in bovine and primate species indicates gene expression changes related to angiogenesis, immune system and participation of functional changes in non-steroidogenic cells [44,23].
Analysis of differentially expressed genes at 3 h post PGF 2a treatment indicated that many of the genes were associated with canonical cell to cell signaling pathways (chemokine, IGF-I and prolactin signaling), cellular growth and proliferation. Several of the products of the up-regulated genes (TNFRSF12A, DAPL1, GEM, PPP4R4 and OSR2) were associated with activation of cellular signaling resulting in inhibition of aggregation of denatured proteins as well as processes associated with protection against metal toxicity and oxidative stress (FABP4, HSPH1, MT1A, MT2A, HP and CBR1). Many of the down-regulated genes belonged to steroidogenesis process (CYP19A1, NR5A1/SF-1, NR5A2/LRH-1, PLA2G1B and PTGR1), cell cycle regulation and tissue development (INSL3, AHSG, BEX2, CDKN1C and HSPB3). Earlier studies involving intra-uterine infusions of low [23] and pulsatile [3] administration of PGF 2a also observed similar changes in many of the genes during functional luteolysis. It should be pointed out that expression of genes regulating steroidogenesis i.e. P 4 biosynthesis (LH/CGR, P450scc and HMGCR) or expression of genes involved in the uptake and trafficking of cholesterol (SR-B1 and StAR) did not change significantly at 3 h post PGF 2a treatment. However, changes in expression of many of these genes were observed at later time points. It has been reported that P 4 decline occurs as early as 30 min after PGF 2a treatment by others and us [69,70] suggest To examine the mechanism of PGF 2a -induced luteolysis, the differentially expressed genes belonging to some processes and functions of luteal cells and few top early differentially expressed genes were subjected to additional analysis. More specifically, as per the objectives of the present study stated in the introduction section, effects of PGF 2a on luteotrophic signaling and role of intra luteal factors were further examined from the differentially expressed genes post PGF 2a treatment. Several studies have reported inhibition of LH-stimulated P 4 secretion by PGF 2a [71,72]. One of the reported actions of PGF 2a was inhibition of adenylate cyclase activity [73,74] that lie downstream of LH receptor activation, but this effect has been reported to be a transient one and moreover, this effect appears to be limited to small luteal cells [75,76]. It is not clear to what extent PGF 2ainhibited adenylate cyclase activity will adversely affect overall P 4 production in vivo, since large luteal cells, besides having large volume and being most abundant, have constitutive steroid synthesis and therefore high P 4 production capacity [63]. On the other hand, the LH-responsive small luteal cells contribute very small fraction of total P 4 produced from the luteal tissue. It should be pointed out that fewer pulses of LH are secreted during the luteal phase that may have little or no impact on the overall P 4 production capacity of the CL [77][78][79][80]. It has been reported that LH/CGR expression gets inhibited post PGF 2a administration in several species [71,72]. In the present study, LH/CGR mRNA expression did not change within 3 h, but decreased at 6 h and thereafter. Surprisingly, LH/CGR protein levels decreased within 3 h post PGF 2a treatment and this may be the reason for immediate inhibitory effect observed in luteal cells as reported in bovines [81] and other species [11,72]. The decrease in PKA activity and pCREB levels at later time points post PGF 2a treatment suggest that aspects of LH/CGR signaling appears to be modulated by PGF 2a directly at the cellular level through inhibition of LH/CGR protein or cAMP levels due to inhibition of cAMP phosphodiesterase enzyme [82,83]. A decrease in levels of SF1 and LRH1 associated with regulation of steroidogenesis indicates regulation of expression of early genes involved in modulation of steroidogenesis [84].
One of the novel findings of this study is the early and consistent down regulation of CYP19A1 expression post PGF 2a treatment. It has been established that in response to LH surge, CYP19A1 expression gets down regulated and may even be absent post ovulation [85][86][87]. However, others [36] and we (in the present study) have observed CYP19A1 expression in the non-pregnant CL. It was of interest to examine the importance of PGF 2ainduced down regulation of CYP19A1 expression. The aromatase enzyme, encoded by CYP19A1 gene, a complex composed of two proteins, an ubiquitous NADPH cytochrome P450 reductase and the cytochrome P450 aromatase, is essential for conversion of androstenedione to E 2 . It has been determined that the regulatory sequence, cAMP-response element (CRE) on the promoter region of CYP19A1 gene is critical for control of expression of aromatase in mammals [88]. For driving expression of CYP19 gene, many tissue specific promoters have been described for several species and in the ovary, promoter II sequence has been reported to be responsible for driving aromatase expression [89][90][91]. Although the promoter II sequence in several species including the cow presents strong homologies, differences in the regulatory sequence have been observed amongst species [92,93]. In the cow, one nucleotide deletion and two substitutions in the CRE-like sequence (CLS) leads to inactivation of the promoter II site and this has been suggested to be responsible for the presence of very low aromatase activity in luteinized granulosa cells or in the CL tissue [92,94]. Interestingly, CYP19 gene expression is most abundant in CL of pregnant rats [92], and this has important role since large quantities of androstenedione, the substrate for the aromatase enzyme, is synthesized in the placenta but gets converted to E 2 in the luteal tissue [95,96]. In bovines, even though CYP19 expression may be lower in non-pregnant CL, but it may be higher during pregnancy [97].
Intra luteal E 2 plays an important role in the structure and function of CL ranging from hypertrophy of luteal cells to increased transport of cholesterol for steroidogenesis [98][99][100]. Estrogens have been reported to have both luteotrophic and luteolytic function across different species [15,36,38,39,64,101], which suggest autocrine/paracrine role within the CL. The biological effects of E 2 are mediated by its binding to the structurally and functionally distinct estrogen receptor (ER), a and b. It has been documented that besides classical E 2 signaling, ligand activated ERa regulates a series of non-genomic events in the cytoplasm including PI3-kinase (PI3K) activation [61,102,103]. The present observation of decreased ER expression and inhibition of Akt phosphorylation post PGF 2a treatment point to loss of PI3K activity during luteal regression as reported previously by others [25,37]. Further, in order to gain more information on the effects of E 2 on the CL tissue, the microarray data of PGF 2a treatment study was mined for E 2 -responsive genes. The observation that a large number of E 2 -responsive genes was observed to be differentially expressed post PGF 2a treatment was surprising since intra luteal E 2 levels are considerably lower in the bovine CL tissue (and also in the circulation) and therefore, it was reasoned that influence of E 2 on luteal function, if any, will be limited. However, to further explore the possible E 2 effects on luteal tissue, effects of PGF 2a treatment on expression of genes in CL tissues of macaques, the species considered to secrete large quantities of E 2 [104] were determined. Mining of microarray data for E 2 -responsive genes from CL tissues of macaques receiving PGF 2a treatment [44] revealed that many of the E 2 responsive genes were also observed to be differentially expressed post PGF 2a treatment in the macaque CL. The observation of many of the differentially expressed E 2 -responsive genes in CL of both macaque and bovine species post PGF 2a treatment provides compelling evidence for regulation of CL function by E 2 . Furthermore, PGF 2a appears to interfere with luteal E 2 secretion and perhaps E 2 actions.
Earlier studies by others and us have provided evidence supporting involvement of cohort of conserved cell death regulatory factors, required for maintenance of growth and development of CL across species [8,9,105,106,107]. Apoptosis at cellular levels is regulated by the intricate changes in the members of the Bcl-2 family proteins that are classified either as anti-apoptotic (Bcl-2 and Bcl-xL) or pro-apoptotic (Bax, Bad, Bik, Bid) molecules. Bcl-2 has been designated as estrogen responsive gene [108,109] and up regulation of Bcl-2 expression has been identified as a critical mechanism for promoting cell survival. Akt mediated phosphorylation of cytosolic proteins is known to play a critical role in the regulation of different metabolic pathways. The promoter region of Bcl-2 contains a cAMP responsive element (CRE) site which is flanked by two estrogen response element (ERE) sites [103,110,111,112]. Based on the previous studies, the transcription factor, CREB can be considered as a positive regulator of Bcl-2 expression [111,113]. Further, Akt, a target of E 2 signaling through PI3K pathway has been shown to activate CREB [103,114] thereby regulating Bcl-2 expression via differential activation of ERs. In our present study, an increase in the Bax/Bcl-2 ratio at mRNA and protein levels further suggests an important role played by permeability of these regulatory proteins across mitochondria as reported during induced and spontaneous luteolysis in bovines [8,9,63]. Another regulatory molecule FKHR which is target of Akt has been proposed to participate in both metabolic and cell survival pathway [115,116], has also been observed to decline post PGF 2a treatment, activating the apoptotic pathway.
The other effects of PGF 2a involve interaction with intraovarian paracrine and endocrine factors [117,118] to mediate intracellular communication within the CL [119]. PGF 2a -induced structural luteolysis depends on cell composition (immune or endothelial cell) [120] and contact [121] as the vasoconstrictor property of PGF 2a may decrease ovarian blood flow and cause apoptosis of endothelial cells [65,[122][123][124] restricting access of gonadotropins and oxygen to steroidogenic cells. The possible role of apoptosis as a mechanism for structural luteolysis has received considerable attention in several species [125,126]. The PGF 2a -induced structural luteolysis mediated by apoptosis was confirmed by the presence of DNA laddering in 18 h PGF 2a -treated CL, which has been demonstrated by us and others through studies involving the role of Bax, Bcl-2, Fas/FasL and caspases during spontaneous or PGF 2a -induced apoptosis in ovine, bovine, rodent CL [9,107,127]. Thus, pathway analysis of microarray data at 6 and 18 h post PGF 2a administration using IPA revealed down-regulated genes to be associated with regulation of steroidogenesis (transcription factors belonging to steroidogenesis, biosynthesis of steroids, intracellular trafficking and lipid metabolism) and angiogenesis (VEGF signaling and coagulation system). Whereas, up-regulated genes to be associated with cell-cell interaction and signaling (EGF, p38 MAP kinase, NF-kB, TGF-b and apoptotic signaling), prostaglandin metabolism and tissue remodelling. Analysis of top differentially regulated genes indicated that most genes belonged to steroidogenesis, prostaglandin metabolism, tissue remodeling or ECM modulation, angiogenesis, TGF-b signaling, cellular stress activated binding proteins, intracellular trafficking, cell survival and apoptotic signaling.
In summary, the results of the present study describe gene expression changes at early time points post PGF 2a treatment. It was observed that PGF 2a treatment caused down regulation of various components of LH receptor signaling, decreased CYP19A1 expression and inhibited intra luteal E 2 levels. A number of E 2 responsive genes were identified to be differentially expressed post PGF 2a treatment. In conclusion, based on the changes of key genes encoding proteins involved in regulating CL structure and function, we propose a model depicting a cross talk between PGF 2a , LH/CGR and E 2 signaling during luteolysis in bovines. PGF 2a -induced luteolysis involves down and up regulation of genes involved in luteal steroidognesis and susceptibility of cell system to apoptotic signals, respectively. The distinct changes that follow post PGF 2a treatment is triggered by down regulation of LH/CGR and ER signaling. E 2 -ER/PI3K-Akt signaling regulates many transcriptional regulatory molecules such as CREB and FKHR, and members of Bcl-2 family proteins such as Bax and Bcl-2. Thus, the present study suggests key role of inhibition of E 2 signaling in the regulation of various changes observed during PGF 2a -induced luteolysis in bovines.   Microarray data analysis was carried out to obtain a set of differentially expressed genes based on statistics, a Student's t-test (two tail, unpaired) with p,0.05 and multiple hypothesis testing (Benjamini and Hochberg comparison test) to reduce the false positives. The identified differentially expressed genes were transcript consistent and did not hybridize to multiple transcripts, as suggested by the AffyProbeMiner analysis. A Bioconductor analysis was performed with $2 fold change as cut-off with statistical filters for identification of differentially expressed genes. Whereas, the top 15 differentially UP regulated genes at 3 h post PGF 2a treatment are represented in this table. Probe Set ID: The identifier that refers to a set of probe pairs selected to represent expressed sequences on an array; Fold Change: It is a number describing changes in expression level of a gene compared between control and treatment; Gene ID: Gene symbols extracted from Entrez Gene or UniGene; Gene Title: Gene name extracted from Entrez Gene or UniGene. (TIF)    responsive genes in bovine CL were identified based on the available list of classical E 2 responsive genes, genes employed in PCR array human estrogen signaling and the data base, ERGDB. Microarray data analysis was carried out to obtain a set of differentially expressed genes based on statistics, a Student's ttest (two tail, unpaired) with p,0.05 and multiple hypothesis testing (Benjamini and Hochberg comparison test) to reduce the false positives. The identified differentially expressed E 2 responsive genes were transcript consistent and did not hybridize to multiple transcripts, as suggested by the AffyPro-beMiner analysis. A Bioconductor analysis was performed with $1 fold change as cut-off and statistical filters for identification of differentially expressed E 2 responsive genes. Whereas, the top 15 differentially UP regulated genes at 3 h post PGF 2a treatment are represented in this table. Probe Set ID: The identifier that refers to a set of probe pairs selected to represent expressed sequences on an array; Fold Change: It is a number describing changes in expression level of a gene compared between control and treatment; Gene ID: Gene symbols extracted from Entrez Gene or UniGene; Gene Title: Gene name extracted from Entrez Gene or UniGene. (TIF)  Table S12 List of common up regulated E 2 responsive genes in monkey (24 h) and bovine (3 h) CL. The previously published microarray data of the differentially expressed genes from the CL tissues of macaques receiving PGF 2a treatment for 24 h [GEO accession number GSE8371] was mined for E 2 responsive genes for purposes of comparing the number of E 2 responsive genes that were differentially expressed in CL of macaques to that of the buffalo cow CL [GEO accession number GSE27961]. The mined data comprising common UP regulated E 2 responsive genes (19 genes) of macaques CL at 24 h vs. E 2 responsive genes of bovine CL at 3 h post PGF 2a treatment are represented in this Table. (TIF) Table S13 List of common down regulated E 2 responsive genes in monkey (24 h) and bovine (3 h) CL. The mined data comprising common DOWN regulated E 2 responsive genes (8 genes) of macaques CL at 24 h vs. E 2 responsive genes of bovine CL at 3 h post PGF 2a treatment are represented in this Table. (TIF) Table S14 List of common up regulated E 2 responsive genes in monkey (24 h) and bovine (6 h) CL. The mined data comprising common UP regulated E 2 responsive genes (17 genes) of macaques CL at 24 h vs. E 2 responsive genes of bovine CL at 6 h post PGF 2a treatment are represented in this Table. (TIF) Table S15 List of common down regulated E 2 responsive genes in monkey (24 h) and bovine (6 h) CL. The mined data comprising common DOWN regulated E 2 responsive genes (13 genes) of macaques CL at 24 h vs. E 2 responsive genes of bovine CL at 6 h post PGF 2a treatment are represented in this Table. (TIF) Table S17 List of common down regulated E 2 responsive genes in monkey (24 h) and bovine (18 h) CL. The mined data comprising common DOWN regulated E 2 responsive genes (13 genes) of macaques CL at 24 h vs. E 2 responsive genes of bovine CL at 6 h post PGF 2a treatment are represented in this Table. (TIF)