Gene expression patterns in four brain areas associate with quantitative measure of estrous behavior in dairy cows.

BACKGROUND: The decline noticed in several fertility traits of dairy cattle over the past few decades is of major concern. Understanding of the genomic factors underlying fertility, which could have potential applications to improve fertility, is very limited. Here, we aimed to identify and study those genes that associated with a key fertility trait namely estrous behavior, among genes expressed in four bovine brain areas (hippocampus, amygdala, dorsal hypothalamus and ventral hypothalamus), either at the start of estrous cycle, or at mid cycle, or regardless of the phase of cycle.
RESULTS: An average heat score was calculated for each of 28 primiparous cows in which estrous behavior was recorded for at least two consecutive estrous cycles starting from 30 days post-partum. Gene expression was then measured in brain tissue samples collected from these cows, 14 of which were sacrificed at the start of estrus and 14 around mid cycle. For each brain area, gene expression was modeled as a function of the orthogonally transformed average heat score values using a Bayesian hierarchical mixed model. Genes whose expression patterns showed significant linear or quadratic relationships with heat scores were identified. These included genes expected to be related to estrous behavior as they influence states like socio-sexual behavior, anxiety, stress and feeding motivation (OXT, AVP, POMC, MCHR1), but also genes whose association with estrous behavior is novel and warrants further investigation.
CONCLUSIONS: Several genes were identified whose expression levels in the bovine brain associated with the level of expression of estrous behavior. The genes OXT and AVP play major roles in regulating estrous behavior in dairy cows. Genes related to neurotransmission and neuronal plasticity are also involved in estrous regulation, with several genes and processes expressed in mid-cycle probably contributing to proper expression of estrous behavior in the next estrus. Studying these genes and the processes they control improves our understanding of the genomic regulation of estrous behavior expression.

Background

Maintaining good fertility and thereby optimum reproductive performance in dairy cows is of great economic importance for the dairy industry. Knowledge on factors influencing fertility is already being applied to improve or regulate fertility. For example, the importance of limiting negative energy balance in early lactation cows for proper reproductive performance is well recognized [1,2]. Insight into the hormonal regulation of estrous cycle has found practical application to artificially regulate the cycle in farm animals and to manage or treat fertility related problems. However, current understanding of genomic factors underlying fertility is limited and this obstructs the development of novel genomic tools and managemental strategies for improving and optimizing reproductive performance, such as biomarkers to monitor the fertility status of cows. Studying the genomic factors underlying fertility may help to optimize nutritional or management systems that improve reproductive performance [3] and also to explain the genetic basis for the decline in several fertility traits of high producing dairy cows. Currently it is known that this decline may be partly attributed to physiological adaptations by the cow to high milk production [4].

Among the fertility traits, the expression of estrous behavior (heat), a key fertility trait that marks the fertile period in cows, has decreased both in duration and intensity over generations of cows selected for high milk yield [5]. Short heat periods and the absence of clear behavioral signs of heat cause farmers to fail to detect heat or to misjudge the optimum time of insemination of their cows, resulting in financial losses due to prolonged interval from calving to first insemination, reduced conception rates and increased calving intervals.

In an effort to understand the genomic regulation of estrous behavior expression in dairy cows, a microarray experiment was set up to study gene expression levels in 4 different brain areas and the anterior pituitary of cows sacrificed at either the start of estrus (day0 of estrous cycle) or at mid-cycle (day12). Differential gene expression analysis between day0 and day12 cows for each of these tissues revealed a limited number of significant genes in the anterior pituitary alone and none in the brain areas (detailed results not reported here). When the trait of interest is quantitative, as in this case (estrous behavior quantified as heat score), the grouping of individuals into qualitative classes dilutes the available information. Therefore, the association between gene expression and phenotypic trait may be better analyzed using the individual quantitative measurements. Using this approach in an earlier study, we identified a set of a few hundred probes out of approximately 24,000 probes on a bovine microarray corresponding to genes whose level of expression in the anterior pituitary of experimental cows associated with the degree to which the individual cows expressed estrous behavior [6]. Among these probes were genes encoding hormones like FSH and prolactin, whose roles in estrous regulation are well-known. Further, biological processes relevant to estrous behavior were over-represented in this set of genes. These results give confidence in the association analysis methodology followed, though experimental validation is needed to determine to what extent the associated genes regulate estrous behavior. In addition to the anterior pituitary, it is likely that a number of areas of the cow brain also have genes whose levels of expression associate with the degree to which cows express estrous behavior. Studies in rodents have revealed that estrogen dependent female reproductive behavior happens via well orchestrated genomic responses in the forebrain with the hypothalamus playing a major role [7,8]. Areas in the limbic region of the forebrain like the amygdala and hippocampus were found to have functions related to sexual behavior and associated emotional responses [9,10]. As yet, there have been no studies in cows linking gene expression in the brain to estrous behavior. Identifying and studying genes whose level of expression in different brain areas of cows associate with the degree to which these cows express estrous behavior will help improve our understanding of genomic factors underlying fertility.

The objective here is to identify and study those genes that associated with estrous behavior, among genes expressed in four bovine brain areas (hippocampus, amygdala, dorsal hypothalamus and ventral hypothalamus), either at the start of estrous cycle, or at mid cycle, or regardless of the phase of cycle.

Results

The degree of estrous behavior expression was quantified as a cow's average heat score using heat scores recorded from at least two consecutive cycles (Table 1). The data from one of the day0 cows was excluded from further analysis because of its high outlier value of 1750 for average heat score, which we attributed to that cow's several consecutive attempts to mount other cows during one observation period. The average heat scores for the remaining 13 day0 cows ranged from 0 to 405, with a mean value of 178.4 (SD 125.7), and the average heat scores for the 14 day12 cows ranged from 2 to 505, with a mean value of 244.7 (SD 175.4). The average heat scores were used with the corresponding gene expression data to run the three analyses per brain area as summarized in Table 2.

Table 1

Average heat scores of the experimental cows sacrificed at day0 or day12 of estrous cycle

day0 Cow Nr	Heat score	day12 Cow Nr	Heat score
d0_5006	0	d12_1528	2
d0_9284	37	d12_9303	4
d0_3739	43	d12_8860	5
d0_8855	53	d12_1773	75
d0_1194	137	d12_8873	157
d0_1821	175	d12_7724	198
d0_8870	191	d12_1520	257
d0_5507	200	d12_7942	275
d0_1786	206	d12_1607	318
d0_3472	246	d12_1822	368
d0_6487	248	d12_1638	383
d0_3747	378	d12_6956	404
d0_7008	405	d12_8857	475
d0_5125 *	1750	d12_2540	505

* Data from this cow was not used in the analysis due to outlier heat score value

Table 2

Description of the three analyses and their objectives

Analysis	Data	Objective
day0	Gene expression data from the tested brain area of Day0 cows and their average heat scores	To identify genes of which the expression in the tested brain area at the start of estrus was associated with estrous behavior
day12	Gene expression data from the tested brain area of Day12 cows and their average heat scores	To identify genes of which the expression in the tested brain area around mid cycle (diestrus) was associated with estrous behavior
day0 + day12	Gene expression data from the tested brain area of Day0 and Day12 cows and their average heat scores	To identify genes of which the expression in the tested brain area was associated with estrous behavior regardless of the phase of estrous cycle

Additional file 1 lists the associated genes found in each analysis and describes the pattern of association between gene expression and heat score for each of the genes in the list. The patterns noticed were: linear (positive or negative slopes) or quadratic (positive/convex shaped or negative/concave shaped curves). The total number of heat score associated probes found common to all 5 Gibbs sampling chains per analysis in each brain area is provided in Figure 1. Additional file 2 depicts Venn diagrams that show the number of overlapping probes between the different brain areas per analysis. The overlap was highest between the DH and VH and then between the AM and HC. The figure 1 and the additional file 2 includes the results of the re-analysis on the AP as well, for which the number of associated probes found for the three analyses (82, 63 and 23 for day0, day12 and day0+day12 respectively) were now considerably lower than the numbers (177, 142 and 118 for day0, day12 and day0+day12 respectively) from the association analysis with third order polynomials and one Gibbs sampling chain per analysis as done earlier [6]. For the different brain areas, the percentage of associated probes found common to all chains per analysis varied from 50-80%. For day0, a relatively high number of heat score associated probes were detected for AM and AP (146 and 82 respectively) whereas for day12 this was true for VH (104). For DH, the numbers were approximately equal at day0 and day12 (51 and 50 respectively). Figure 2 provides the association patterns for three genes whose expression values were found to be associated with heat score at day0. For AVP and MCHR1, a linear trend was observed in HC and AM respectively whereas for OXT, a quadratic trend was observed in DH.

Figure 1

Number of estrous behavior associated probes found in 4 brain areas and anterior pituitary. AM - Amygdala; HC - Hippocampus; DH - Dorsal Hypothalamus; VH - Ventral Hypothalamus; AP - Anterior Pituitary.

Figure 2

Association patterns of gene expression profile with heat scores for three estrous behavior associated genes at day0.

Of the total 640 unique probes found associated with heat score in all analyses and brain areas, 372 had gene annotations. Sets of genes associated with functional categories which we group as: "transcription and regulation of gene expression", "detection of and responses to stimuli" and "signaling pathways" each made up about 15% of the genes that were identified. In the category "detection of and responses to stimuli" we included the genes involved in neurotransmission through encoding for neurotransmitters/hormones (OXT, POMC, AVP, CCK, CGA) or neurotransmitter receptors (GABRA6, HTR2A, MCHR1, CHRM1, CHRM3, DRD2, CHRNA5) and metabolizing enzymes (PTGDS, PTGIS, PTGR1, ACHE, SULT4A1). We could also identify 5 more sets of genes, each consisting of between 5 to 10% of the total genes. These included genes associated with: "transport and localization" (positioning of a substance or cellular entity, like sorting nexin family member genes), "transporter activity" (e.g. calcium channels and Na+/K+ transport), "metabolism" (e.g. glucose and amino acid metabolism or steroidogenesis), "cell cycle" (here including processes linked to DNA folding and repair) and "multi-organism processes". In the last category we included immune system related genes like CTLA4, IL1RL1, MARCO, FCRLA and IL33.

A detailed list of all significant enriched GO and KEGG pathway terms (p < 0.10) found in the different brain areas and analyses is provided in Additional file 3. For illustration, a summary of only those GO biological processes that cleared a FDR cutoff of 20% is provided in Table 3. Several relevant processes related to behavior, neurotransmission and signaling, and ion regulation were found especially in day0. The presence of enriched processes related to behavior like 'grooming behavior', 'regulation of female receptivity' and 'female mating behavior' within the associated genes found in HC indicates the key role of the genes AVP and OXT that contribute to these processes. Terms related to neurotransmission and signaling and the associated ion regulation terms too have biological implications related to estrous behavior as evident from the discussion below on some of the genes that contribute to these processes.

Table 3

Significantly enriched Gene Ontology Biological Process terms in the estrous behavior associated gene lists*

Gene Ontology Biological Process term	day0	day12	day0 + day12
Behavior related terms
regulation of female receptivity; reproductive behavior in a multicellular organism; female mating behavior; maternal aggressive behavior; maternal behavior; parental behavior	HC	DH	HC
grooming behavior			HC
locomotory behavior	AP+
behavior	AP+
Neurotransmission and signaling related terms
neurological system process	AP+
regulation of neurotransmitter levels	AP
transmission of nerve impulse	AP+
regulation of synaptic transmission	DH
cell surface receptor linked signal transduction		DH+
cell-cell signaling	HC+
second-messenger-mediated signaling			VH+
response to cAMP	DH
Wnt receptor signaling pathway	DH
phosphoinositide-mediated signaling			VH
Ion regulation related terms
ion transport	AP+
elevation of cytosolic calcium ion concentration			DH
negative regulation of ion transmembrane transporter activity	AM
cytosolic calcium ion homeostasis	AM

* Only terms with p < 0.10 and fdr <0.20 are included. Terms typical of other organs (these terms appeared here due to genes with different functions in different organs) have been removed from this summary table. Where the term is supported by 5 or more genes in the gene list, the corresponding brain area is followed by a sign whereas in cases where the term is supported by only 2-4 genes, the sign is absent.

Discussion

Variation in the behavioral trait to express estrous behavior was found to be associated with variation in the expression of genes in the cow's brain areas: DH, VH, HC and AM and also the AP, an endocrine gland pivotal in synchronizing estrous behavior with hormonal and ovarian events preceding ovulation. The choice of these tissues for this study was based on their reported involvement in regulating female sexual behavior. Although genes differ in their influence on specific traits, it is tempting to consider the AM to be of relatively higher importance for regulating estrous behavior as the largest number of associated genes was found here at day0 and because AM is known for its central role in regulating emotions. Key genes and biological processes as identified from the lists of heat score associated genes are discussed next and linked to estrous behavioral expression, though the links are not always as expected.

Genes and biological processes associated with estrous behavior in line with expectations

An association with heat scores was detected for the gene, which encodes for the alpha subunit of glycoprotein hormones (FSH, LH, TSH), with associations being time and brain area specific. It seems that cows with clear expression of estrous behavior have relatively high expression of this gene in the hypothalamus around mid-cycle and low expression around start of estrus.

The gene was found associated in the HC and DH whereas the gene was found associated in the HC and AM. The known functional properties of these genes contributed to several enriched GO terms in the DAVID analysis, especially those related to female mating behavior. Oxytocin is released within the brain where it acts on specific oxytocin receptors to elicit effects like female sexual receptivity, grooming behavior and partner bonding [11]. In the presence of estrogen, oxytocin exerts an anxiolytic effect thereby favoring courtship and mating [12,13]. Similar to oxytocin, AVP is associated with sexual behavior and bonding [14] and its expression is under control of the steroids progesterone and estrogen [15,16]. In mice, the absence of estrogen receptors (ERα and ERβ) impairs social recognition similar to the effect of OXT gene deletion [17]. Studies with rodents and humans [18] demonstrate that oxytocin and vasopressin modulate complex socio-sexual behavior, typically under the influence of reproductive steroid hormones. The present association of OXT and AVP in several brain areas is in line with the above mentioned findings and suggests a major influence of these genes on estrous behavior expression in dairy cows.

The genes and , which associated with heat score in at least one brain area, are known to modulate emotional states like anxiety and satiety [19-21] or even sexual motivation. In sheep, dopamine-mediated D2 receptor (DRD2) signaling in the mediobasal hypothalamus affects female sexual motivation and receptivity [22]. Interactions between monoamines (dopamine, serotonin, noradrenaline) and steroid hormones play a major role in the integration of reproductive behavior and gonadal function [23]. The perception and awareness of male-related cues differs with the stage of estrous cycle, with releases of monoamines (linked to HTR2A and DRD2) and gamma-aminobutyric acid (GABA) (linked to GABRA6) in the mediobasal hypothalamus being triggered by such cues only when ewes are in estrus [24]. We found serotonin receptor 2A (HTR2A) associated with heat score at day12 in VH. Studies in female rats and hamsters have shown the inhibitory and facilitatory effects of serotonin receptor agonists and antagonists on the hypothalamic regulation of sexual receptivity [25,26], and this regulation is also mediated by GABAergic neurons interacting with serotonin containing neurons [27].

Noteworthy, are the heat score associated genes (that regulate prostaglandin functioning. In the central nervous system, prostaglandins are involved in functions like thermoregulation and influencing neuronal morphology. Prostaglandins are also known to be under the influence of estradiol [28] and are capable of directly affecting neurons that synthesize and secrete gonadotropin-releasing hormone [29].

Another heat score associated gene, , a carrier of thyroid hormone, is known to influence anxiety [30], behavioral activity [31] and mental functions [32]. The melanin-concentrating hormone receptor, MCHR1, plays a role in metabolic rate and feed intake [33]. Changes in anxiety behavior and feeding motivation are likely to facilitate mating.

The present association between heat scores and the expression of and several cholinergic receptors (and may be explained by the effect of the neurotransmitter, acetylcholine on arousal, plasticity and reward. The products of the muscarinic cholinergic receptor genes, CHRM1 and CHRM3, are Gq-protein coupled receptors whose activation releases intracellular Ca2+ via the phospholipase C - inositol 1,4,5-trisphosphate signaling pathway [34]. The genes for phospholipase C and inositol triphosphate kinase () and several protein kinases were also found associated to heat score. These findings may be explained based on the hypothesis put forward by Kow and Pfaff [35] that the membrane actions of estrogen can modulate the genomic actions of estrogen and that this transcriptional potentiation was mediated via signaling pathways requiring the activation of certain protein kinases and increased intracellular Ca2+.

Genes and biological processes unexpectedly associated with estrous behavior

The finding of expected processes related to estrus as described in the previous section supports assumed neurophysiological mechanisms underlying female sexual behavior in dairy cows. Here, examples are given of more novel candidate genes and mechanisms. The gene encodes the protein tachykinin (or neurokinin B), which in humans has been considered a critical regulator of gonadotropins (LH, FSH) via regulating GnRH secretion [36]. The present findings encourage further investigation towards the importance of tachykinin associated mechanisms in dairy cow fertility. Behavioral changes during estrus represent changes in central perception and processing of information, i.e. cognition, and some of the genes that associated with heat scores have been linked to cognition, e.g. . PEBP1 has been suggested to be involved in chronic stress-induced memory impairment [37], MOBP has been linked to mood disorders [38], LTA4H to depression in women [39], and KCNN2 to anxiety and stress responses [40].

The gene in mice seems to establish a cross-talk between reproduction and metabolic events [41] and was here associated with heat scores. Also the gene, which encodes transcription factors involved in activation of growth hormone and prolactin, was found associated with heat score in the AP on day0. This gene may contribute in part to the generally observed reduction in estrous behavior in high producing cows. The heat score associated gene is noted here because it has been linked to egg productivity in laying hens [42] and, thus, may have a fertility related function.

A number of heat score associated genes have been linked to the immune system, e.g. and , indicating the importance of cell-cell interactions. It has been shown that immunoglobulin superfamily proteins may play important roles in brain developmental processes and the functioning of neuronal networks in adults because they provide the ideal structure for protein-protein interactions and, thus, cell-cell interactions [43]. Remodeling of synaptic networks, which occurs during estrogen promoted female sexual behavior [7], may also be facilitated by immunoglobulins. Together with the many associated genes that are known to regulate cell fate, activity and morphology, this seems to underline the importance of neural tissue plasticity in the appropriate expression of estrous behavior.

The gene product Sodium/potassium-transporting ATPase subunit alpha-3 (), which was found to be associated to heat score at estrus in AP in our earlier study was also found in VH. ATP1A3 has been implicated in rapid-onset dystonia parkinsonism (RDP), characterized by sudden onset of neurological symptoms over hours to a few days [44], suggesting a role in the sudden onset of behavioral changes like during estrus.

Estrous behavior associated genes expressed at estrus and mid-cycle

Including data from both day0 and day12 cows in a combined analysis not only revealed genes that were associated with estrous behavior regardless of phase of cycle but also resulted in greater statistical power and helped reveal some associated genes which could not be found in the separate analysis. For example, DRD2 (dopamine receptor D2) gene has been found associated in DH. However, care needs to be taken to interpret the results of the combined analysis where an associated gene was also found in one or both of the separate analyses. There could be cases where the day0+day12 combined analysis found an associated gene because of the effect of a strong association found in one of the separate analyses. An example of this case is HTR2A which did not associate with heat scores on day0, but did on day12 (R2 = 0.27) and also weakly (R2 = 0.1) in the day0+day12 combined analysis, probably as a carryover effect.

Interesting observations were also made by investigating associated genes that appeared in several brain areas and analyses. For example, the gene found in VH at day0 and in the combined analysis and also in the AP, is known to mediate the uptake of glutamate into synaptic vesicles at presynaptic nerve terminals of excitatory neural cells and may also mediate the transport of inorganic phosphate [45]. The gene could be a contributor to neurotransmission associated with estrous behavior expression and was one of the genes contributing to several ion transport related GO terms found enriched in the DAVID analysis. The gene , found in several brain areas across all analyses, was seen to contribute to several GO terms related to binding and signal transduction. It has been reported to have a role as an antagonist of Wnt signaling pathway which may influence the development of neurons in dorsal midbrain [46], suggesting again a link between neuronal plasticity and estrous behavior. The gene , also found associated in several brain areas and analyses, is a member of the anoctamin family which has been implicated in calcium ion-activated chloride channels that perform several important functions including neuronal excitability.

To summarize, several genes were found here whose expression levels in the brain areas were associated with the degree to which cows express estrous behavior. For some of these genes, there is a known function relating them to processes regulating estrous behavior, while for others, the current association suggests such a relation. We propose that the genes OXT and AVP play major roles in regulating estrous behavior along with genes affecting neurotransmission and neuronal plasticity. Genes, whose expression in mid-cycle associated with estrous behavior, may contribute to preparing the cow for the next estrus. Figure 3 depicts the key findings. This study may assist in the search for biomarkers for estrus detection or in screening for the most likely genes within QTLs associated with fertility. Further research with more animals sampled at multiple time points in the estrous cycle may improve our understanding of the dynamic regulation of estrous behavior over the cycle. Nevertheless, this first study provides an understanding of some of the genes and processes related to expression of estrous behavior in dairy cows.

Figure 3

Schematic representation of key genes found associated with estrous behavior and their relationship with estrus. Some of the key genes whose expressions in the brain areas were found associated with estrous behavior are depicted here along with mechanisms affected during estrus that are important for estrous behavior. GnRH released from the DH stimulates release of FSH and LH from the AP which in turn influences the ovary to release estrogen and progesterone (Hypothalamic-pituitary-gonadal axis). In the brain areas, estrogen influences expression of OXT and AVP, which play major roles in the regulation of sexual behavior. The neurotransmitter receptors, GABRA6, HTR2A and DRD2 affect sexual motivation, whereas CHRM1, CHRM3 and CHRNA5 along with "immune related" genes like CTLA4, IL1RL1 and MARCO affect neuronal plasticity. The genes POMC, MCHR1 and TTR are involved in altering anxiety, stress and feeding motivation.

Conclusions

The study predicted estrous behavior associated genes in the 4 brain areas at two time points of the estrous cycle. These included genes expected to be related to estrous behavior as they influence states like socio-sexual behavior, anxiety, stress and feeding motivation (OXT, AVP, POMC, MCHR1), but also genes whose association with estrous behavior is novel and warrants further investigation. Studying these genes and the processes they control improves our understanding of the genomic regulation of estrous behavior expression, ultimately leading to better management strategies or tools to improve or monitor reproductive performance.

Methods

Data recording, sample isolation and microarray hybridization

Estrous behavior was recorded in 28 healthy Holstein Friesian cows from 30 days in milk (DIM) onwards till their time of sacrifice which varied between 77 and 139 DIM. The estrous behavior recorded in these cows were quantified as heat scores and the scores from multiple consecutive cycles per cow were averaged to obtain the average heat score per cow [6]. Of the 28 cows, 14 were sacrificed at the start of their estrous cycle (day0) and 14 around mid cycle (day12). Following sacrifice, brain tissue samples were collected from 4 brain areas: amygdala (AM), hippocampus (HC), dorsal hypothalamus (DH) and ventral hypothalamus (VH). Since the brain areas are not clearly demarcated, a standardized system was used for positioning and dissecting the brain which ensured that the collected samples representing each brain area in all cows were identical. The procedure followed for collecting brain samples is described in detail in Additional file 4. The study was approved by the Animal Care and Ethics Committee of the Animal Sciences Group of Wageningen University and Research Centre, Lelystad.

RNA extracted from brain tissue samples were hybridized on Bovine 24K oligonucleotide (70-mer) microarrays designed and produced by the Bovine Oligonucleotide Microarray Consortium (BOMC), USA http://www.bovineoligo.org/. This array was amongst the few whole genome bovine spotted microarrays available at the time this experiment was performed in the year 2007-08. The choice of microarray technology over RNA sequencing, which was not a well established technology then, was appropriate for fulfilling the objectives of this study in a precise and cost effective manner. The procedures followed for RNA preparation and microarray hybridization were as described in our earlier study [6]. The purity of the total RNA was assessed using the A260/280 ratios given by NanoDrop spectrophotometer and found to have a ratio above 1.8 indicating good quality RNA. In addition, we also performed agarose gel electrophoresis for a visual inspection of the RNA integrity. There was no indication of degraded RNA as deduced from the intensities of the 28S and 18S bands. Based on the combination of results of NanoDrop and gel electrophoresis, we were satisfied with the quality of RNA extracted by our extraction procedure. Five μg RNA was used per labelling using the RNA MICROMAX TSA labelling and detection kit (Perkin-Elmer). A total of 224 arrays (i.e. 28 cows × 4 brain areas × 2 dye swaps) were prepared in a common reference design with the dye labels swapped between individual samples from each brain area and a reference sample consisting of equal proportions of RNA from all four brain areas as well as AP of all cows. Processed slides were scanned using GenePix 4200A (Molecular Devices), with identical settings and the images processed using GenePix Pro 6 software (Molecular Devices). Regarding microarray quality and validation, a recent re-annotation of the probes in this array following the drafting of the bovine genome ensured that only good quality probes were considered in further analysis (details in the second following section). All array hybridizations and processing were done in a standard lab with experienced personnel using a standard microarray platform. Moreover, the experimental design used 14 biological replicates per time point which is relatively large and sufficient for a sound statistical analysis and the use of dye swap technical replication also improves the analysis. Given the above facts, we believe that a PCR validation of the individual probes was not required, the more so because we analysed clusters of genes instead of individual genes.

The array experiment data is publicly accessible on ArrayExpress (accession: E-TABM-916 on http://www.ebi.ac.uk/arrayexpress).

Microarray data analysis and identification of estrous behavior associated genes

Microarray data pre-processing and analysis was done using the LIMMA (Linear models for microarray data) package [47] within Bioconductor project [48] of R statistical programming language http://www.r-project.org, identical to the procedures described and used in our earlier study [6]. Gene expression levels expressed as M-values (log differential-expression ratio of sample versus reference) were thus obtained for each brain area per cow. For each gene, we then modeled its expression level across all cows as a function of the cows' average heat score. For this we used the Bayesian hierarchical mixed model developed by Jia et al. [49] that employs orthogonal polynomials to quantify linear and non-linear associations between quantitative phenotypes and gene expression data. We used the algorithm of Jia et al. [49] coded in SAS® language which the authors kindly provided. The program was run on SAS® software, Version 9.1 of the SAS® System for Windows. Similar to the approach in our previous study on the AP [6], here we used the SAS program for three separate analyses on data from each of the 4 brain areas to achieve the objectives of this study (summarized in Table 2). However, in contrast to the procedure followed in our earlier study, where gene expression was modeled as a function of average heat score with third order polynomials, here we used second order polynomials. The Bayesian hierarchical mixed model used here remains the same as was explained in our earlier study [6]. However, we limited our analyses to two orders, linear and quadratic, because the resulting relationships for third order polynomials were observed to be less reliable considering the low number of animals used in this study. For each gene, posterior probabilities were obtained for the two regression coefficients to be different from zero. Genes were considered to have a significant association when the posterior probability of at least one of the regression coefficients was larger than 0.80, thereby limiting the FDR to below 1% [50]. We performed the Gibbs sampling chains multiple times (arbitrarily set to 5 chains with each chain running 10,000 iterations of which 5000 were burn-in) for each analysis of a particular brain area and only those genes that appeared significant in all chains were considered in the functional analysis (see next section). By doing so, we limited the variations due to Gibbs sampling.

Functional analysis of estrous behavior associated genes

The original annotation of the bovine 24 K oligonucleotide microarray provided by BOMC dates back to June 2007. For our analysis, we used the bovine oligonucleotide array probe re-annotation (Version 5) based on Ensembl http://www.ensembl.org release 56 (October 2009) provided on the EADGENE website by the authors of the oligo-set re-annotation pipeline, sigReannot [51]. For the re-annotation, out of the 23,496 probes (all control probes removed) on the bovine oligonucleotide array, only 16,620 probes that were assigned a quality score between 1 and 4 based on their specificity to hits on the bovine genome were considered. Probes with quality scores between 5 and 7 had either no hits or multiple hits and were not annotated as they were not specific.

For gaining insight into biological processes underlying estrous behavior, we performed functional analyses of the sets of estrous behavior associated genes (study sets) identified in the 3 analyses for the 4 brain areas and the anterior pituitary. To increase the accuracy of functional analysis, we used the re-annotated probes information. The probe re-annotation also provides the Ensembl gene ID of the orthologous human genes for 14,585 probes, which we used in the functional analysis instead of bovine genes so as to benefit from the greater gene annotation information available for human species.

Functional analysis was done using DAVID bioinformatics resources [52,53], a freely available web-based tool http://david.abcc.ncifcrf.gov/ that integrates biological data from several sources including GO [54] and biological pathway databases. In order to get an indication of the major processes over-represented among the annotated probes (genes) in the study sets (all linear/quadratic associated genes in each of the 3 analyses for the 4 brain areas), we interpreted the significance of the DAVID results based on EASE score (p value derived from a modified Fisher's exact test) [55] threshold set at 0.10. We did not consider multiple testing correction as it was too conservative for our purpose of discovery of biological themes in our study sets. Genes in the study sets were tested for enriched GO terms and KEGG [56] pathway terms using DAVID functional annotation tool. The population set against which the study set genes were tested consisted of the Ensembl IDs of the orthologous human genes, of which 11,589 remained after removing duplicates..

List of abbreviations

AM: Amygdala; AP: Anterior Pituitary; BOMC: Bovine Oligonucleotide Microarray Consortium; DAVID: Database for Annotation, Visualization and Integrated Discovery; DH: Dorsal Hypothalamus, EADGENE: European Animal Disease Genomics Network of Excellence for Animal Health and Food Safety; FDR: False Discovery Rate; FSH: Follicle stimulating hormone; GnRH: Gonadotropin-releasing hormone; GO: Gene Ontology; HC: Hippocampus; KEGG: Kyoto Encyclopedia of Genes and Genomes; LH: Luteinizing hormone; QTL: Quantitative trait loci; SD: Standard Deviation; TSH: Thyroid stimulating hormone; VH: Ventral Hypothalamus;

Authors' contributions

AACW performed the microarray experiments. AK performed the data analysis and drafted the manuscript. HAM and RFV helped with the statistical analysis. BB, HW, MFWP and MAS helped with the biological interpretation of the results. All authors read and helped in improving the manuscript.

Additional file 1

Association patterns between gene expression and heat score of estrous behavior associated genes identified in the four brain areas and anterior pituitary. See legend at end of table.

Click here for file

Additional file 2

Venn diagrams showing the number of overlapping probes between the different brain areas per analysis. The figures in brackets represent the total number of estrous behavior associated probes found in each brain area per analysis. Here, AM - Amygdala; HC - Hippocampus; DH - Dorsal Hypothalamus; VH - Ventral Hypothalamus; AP - Anterior Pituitary.

Click here for file

Additional file 3

Significantly enriched Gene Ontology and KEGG pathway terms in the estrous behavior associated gene lists. Table shows all terms with p < 0.10 and includes gene ontology terms in the 3 categories: Biological Process, Molecular function and Cellular component.

Click here for file

Additional file 4

Collection of brain samples and pituitary from the experimental cows. Procedure followed for collection of anterior pituitary and brain samples: amygdala, hippocampus, dorsal hypothalamus and ventral hypothalamus (with pictures).

Click here for file