Changes in gene expression associated with reproductive maturation in wild female baboons.

Changes in gene expression during development play an important role in shaping morphological and behavioral differences, including between humans and nonhuman primates. Although many of the most striking developmental changes occur during early development, reproductive maturation represents another critical window in primate life history. However, this process is difficult to study at the molecular level in natural primate populations. Here, we took advantage of ovarian samples made available through an unusual episode of human-wildlife conflict to identify genes that are important in this process. Specifically, we used RNA sequencing (RNA-Seq) to compare genome-wide gene expression patterns in the ovarian tissue of juvenile and adult female baboons from Amboseli National Park, Kenya. We combined this information with prior evidence of selection occurring on two primate lineages (human and chimpanzee). We found that in cases in which genes were both differentially expressed over the course of ovarian maturation and also linked to lineage-specific selection this selective signature was much more likely to occur in regulatory regions than in coding regions. These results suggest that adaptive change in the development of the primate ovary may be largely driven at the mechanistic level by selection on gene regulation, potentially in relationship to the physiology or timing of female reproductive maturation.

Transcriptome Analysis and Life History in a Wild Primate Population

Nonhuman primates are valuable sources of insight into human evolution. Until recently, however, such insight was limited by the dearth of genetic resources for most primate species. In addition, studies of primates in their natural habitats, while rich in behavioral and ecological detail, have rarely included extensive genetic or genomic components. This situation is changing now that genomic resources are increasingly available, and gene regulatory studies of captive primates have set the stage (reviewed in Tung et al. 2010). However, we still know relatively little about variation in gene expression in wild primates.

Collecting functional genomic data on such systems could provide important context for the evolution of gene regulation in humans. Specifically, studying changes in gene expression during maturational milestones in nonhuman primates may provide insight into the loci that contributed to shifts in developmental timing and physiology during human evolution (Uddin et al. 2008; Somel et al. 2009; Gunz et al. 2010). Some examples of these shifts include relatively late menarche in human hunter-gatherers compared with nonhuman primates (reviewed in Blurton Jones et al. 1999); a skeletal growth spurt that accompanies female maturation in humans that appears to be absent in nonhuman primates (Bogin and Smith 1996); and short interbirth intervals in humans relative to body size (reviewed in Mace 2000). Circumstantial evidence suggests a role for gene regulation in these changes. Indeed, sequence-based analyses have revealed that the regulatory regions of many development-related genes have undergone positive selection within primates (Haygood et al. 2010) and that rapidly evolving regulatory regions near duplicated genes in humans are enriched for genes related to pregnancy and reproduction (Kostka et al. 2010).

Yellow baboons (Papio cynocephalus) are close human relatives (∼94% sequence similarity: see Silva and Kondrashov 2002) that, like humans, are large-bodied terrestrial savanna primates with long life histories and nonseasonal reproduction. They also inhabit African savanna environments similar to those relevant for early humans (Potts 1998; Behrensmeyer 2006). Yellow baboons have been the subjects of extensive study in the wild (Altmann SA and Altmann J 1970; Jolly 1993; Rhine et al. 2000; Buchan et al. 2003; Wasser et al. 2004; Alberts et al. 2006), including in the Amboseli basin of Kenya where individually recognized baboons have been monitored since 1971 (Altmann SA and Altmann J 1970; Buchan et al. 2003; Alberts et al. 2006). This system therefore presents an exceptional opportunity to integrate functional genomic data sets with detailed life history information about the same animals.

Here, we take advantage of life history and behavioral data from the Amboseli baboon population, combined with an unusual circumstance in which we were able to collect fresh tissue from seven known females (four premenarcheal juveniles and three multiparous adults). Six of these seven females died in an episode of conflict with the local human population (the Maasai community in Amboseli) that perceived the baboons as a threat to their livestock; the seventh died of natural causes a few days later. The bodies of all seven females were collected within a few hours of their death, with the help of the Maasai community. We used these data and samples to investigate gene expression changes related to the onset of sexual maturity in females and to examine differential expression in maturity-related genes among genes inferred to have evolved under lineage-specific selection in primates. We focused specifically on expression differences in the ovary, an organ that plays a central role in reproductive maturation. We present a genome-wide analysis of ovarian gene expression changes in these seven female baboons from this natural population using RNA-Seq.

Expression Differences by Life History Stage

RNA-Seq libraries were made using ovarian RNA from three adult and four juvenile females (supplementary fig. S1 and table S1, Supplementary Material online). We obtained ∼15 million reads per individual (supplementary table S1, Supplementary Material online), and we measured the expression of a total of 9,770 genes in the baboon ovary. Ninety-seven genes (∼1% of genes in the data set) were differentially expressed between the juveniles and the adults (False discovery rate [FDR]–adjusted P value < 0.05) (fig. 1). This result is consistent with the expectation that intraspecific differential expression, particularly within a population and within sex, is likely to be less common than interspecific differential expression between different primate species (Babbitt et al. 2010; Blekhman et al. 2010; Xu et al. 2010). Of the differentially expressed genes, 79 were upregulated in the adults and 18 were upregulated in the juveniles. This imbalance in upregulated expression toward the adult females was expected, as the adult ovary is much more metabolically active than the premenarcheal ovary (reviewed in McGee and Hsueh 2000).

F

MA plot of the normalized data. Each dot represents a single gene, and significantly differentially expressed genes are colored by higher expression levels in adults (red) or juveniles (blue).

To evaluate the global effect of maturation stage on gene expression variation, we performed a principal components (PCs) analysis. The first three PCs in this analysis explained ∼67% of variation in the gene expression data (supplementary fig. S3, Supplementary Material online). None of these PCs clearly differentiated adult and juvenile tissues, although PC2 exhibited the strongest (albeit nonsignificant) relationship with life history stage (Mann–Whitney test, W = 11, P value = 0.1143). In contrast, when we examined only those genes that were significantly differentially expressed (n = 97), PC1 alone explained 70% of the variation in the gene expression data. PC1 also exhibited a trend toward higher values for juveniles than for adults (Mann–Whitney test, W = 0, P value = 0.05714 and supplementary fig. S3, Supplementary Material online).

Little is known about ovarian gene expression in human or mouse models during either the premenarcheal stage or in noncycling adult tissue, as most studies concerning ovarian gene expression have focused on embryonic sex specification (Nef et al. 2005; Small et al. 2005), fertility disorders (reviewed in Matzuk and Lamb 2008), or cancer states (e.g., Wang et al. 1999; Welsh et al. 2001; Haviv and Campbell 2002; Adib et al. 2004). To explore patterns in expression differences between these life history stages, we performed categorical enrichment analyses using the GO (Gene Ontology Consortium 2000) and PANTHER (Mi et al. 2005) ontology databases. The enrichments were performed in two ways: first, using the absolute rankings of gene expression differences between adults and juveniles, regardless of the direction of the difference; and second, using only genes that were more highly expressed in the more metabolically active adult tissue (table 1, supplementary tables S2 and S3, Supplementary Material online).

Table 1

PANTHER Biological Process Categorical Enrichments

NOTE.—Categorical enrichments were performed using a Wilcoxon rank test. The right-hand column shows the total number of genes evaluated. Categories that evaluated fewer than ten genes are not shown. Categories in white have a B-H corrected P value < 0.05 (Benjamini and Hochberg 1995), and categories in gray have a nominal P value < 0.05.

Several patterns emerged from these analyses. First, we identified a number of ontology categories generally associated with blood, including “immunity and defense” and “angiogenesis” (table 1). The cortex of the ovary becomes highly vascularized after the onset of maturity (Redmer and Reynolds 1996; Abulafia and Sherer 2000), a maturational process that could account for some of the observed enrichments. In addition, follicular development in the mature ovary is correlated with increased inflammation (reviewed in Bukovsky and Caudle 2008). In keeping with this change, we identified cytokine, chemokine, and macrophage-related immune activities among the significant categories of genes that show increased expression in the adults (supplementary table S3, Supplementary Material online). Second, and perhaps unsurprisingly, genes involved in developmental processes (i.e., developmental processes and mesoderm development) tended to be enriched for differential expression (table 1 and supplementary table S2, Supplementary Material online). These enrichments emphasize that the physiological distinctions between the mainly quiescent juvenile ovary and the mature ovary are likely related, at least in part, to differences in gene regulation.

At the level of individual genes, we found a significant upregulation in the adult ovary of genes essential for ovarian function and folliculogenesis (table 2), including genes such as VGF (VGF nerve growth factor inducible), MMP19 (matrix metalloproteinase-19), and ADAMTS1 (a disintegrin and metalloproteinase motif 1) (fig. 2). MMP19 and ADAMTS1 function to remodel the extracellular matrix as follicles develop (Jo and Curry 2004; Brown et al. 2010). The role of VGF is less clear, but its essential role has been demonstrated in VGF−/− mice, which produce many primary follicles but few mature follicles (Hahm et al. 1999; Jethwa and Ebling 2008). Fewer genes are upregulated in the juveniles; however, one intriguing example is RSPO1 (R-spondin1), which is known to be critical for early human ovary development and specification (Tomaselli et al. 2011). Our data indicate that it continues to be expressed until the stages right before puberty (fig. 2).

Table 2

Differentially Expressed Genes (FDR-adjusted P value < 0.05) in the Adult and Juvenile Baboon Ovary

Gene ID	Log FC	P Value	P Value FDR	Gene ID	Log FC	P Value	P Value FDR
serpina3	5.56053945	2.44 × 10−16	2.39 × 10−12	EPO	−2.3597099	6.37 × 10−05	0.01219865
ADAMTS4	5.09716829	5.79 × 10−15	2.83 × 10−11	SLC7A8	2.85383446	6.49 × 10−05	0.01219865
REN	3.98944937	5.39 × 10−12	1.76 × 10−08	ERRFI1	2.13531073	7.99 × 10−05	0.01472884
TFPI2	5.09331391	1.82 × 10−11	3.58 × 10−08	MYOC	−4.5033477	8.42 × 10−05	0.0152368
ADAMTS6	4.69173617	1.83 × 10−11	3.58 × 10−08	Mmp1	8.22694892	8.90 × 10−05	0.01581064
Melk	3.76852295	2.12 × 10−10	3.44 × 10−07	rspo1	−2.7082948	0.00010143	0.01769646
LRG1	4.68207316	4.18 × 10−10	5.83 × 10−07	HIF1A	2.20582485	0.00010339	0.01772116
FABP4	9.28927774	3.47 × 10−09	4.24 × 10−06	rpl21	3.86520992	0.0001073	0.01807402
CH25H	3.20323737	1.81 × 10−07	0.00018065	OSMR	2.19076016	0.00011098	0.01837741
SOST	−4.8295275	1.85 × 10−07	0.00018065	TIMP1	2.00766916	0.00011652	0.01897307
IL1RL1	9.32559106	3.14 × 10−07	0.0002791	CYP21A2	2.59714864	0.00012208	0.0193343
F3	2.72902443	4.14 × 10−07	0.00033246	PAPPA	3.66989958	0.00012459	0.0193343
RGS2	2.94854719	4.42 × 10−07	0.00033246	Adamts1	−3.7455767	0.0001265	0.0193343
GALNT9	−3.4824072	6.16 × 10−07	0.00042976	Trem1	1.97618174	0.0001301	0.0193343
TNFAIP6	9.06613349	6.82 × 10−07	0.00044438	RAB38	−2.3456431	0.00013012	0.0193343
CRTAC1	−3.1329928	7.44 × 10−07	0.00045431	Sgip1	2.17468855	0.00013061	0.0193343
C19orf26	2.79022659	1.08 × 10−06	0.00062056	DDX21	2.13304855	0.00013343	0.01945687
LdhA	2.80233251	1.91 × 10−06	0.00102984	f2rl1	3.0581654	0.00014459	0.0207748
stc1	2.95769089	2.00 × 10−06	0.00102984	SBNO2	2.8690604	0.00015474	0.02160406
Gdf15	2.63705264	2.61 × 10−06	0.00127693	S100A8	3.74690504	0.00015479	0.02160406
FCER1G	3.18705038	2.77 × 10−06	0.00129071	DST	3.67689535	0.00016231	0.02233414
VGF	3.53951033	3.72 × 10−06	0.00165083	AADAC	7.97969292	0.00016903	0.02293631
MMP19	2.49571312	5.03 × 10−06	0.00213603	CHI3L1	2.26915607	0.00017834	0.0238048
Fosl2	2.56131942	5.46 × 10−06	0.00222307	H6pd	2.1574967	0.0001803	0.0238048
SERPINE1	3.00295543	6.06 × 10−06	0.00226587	ADAMTS16	−2.5741504	0.00018334	0.02388255
S100A9	2.80444185	6.24 × 10−06	0.00226587	HLA-DQB1	3.33877719	0.00019359	0.02488701
ADPRHL1	3.46982396	6.26 × 10−06	0.00226587	CTSG	2.8539012	0.00021891	0.02777611
Cd163l1	2.93708735	7.60 × 10−06	0.00265144	RPF2	2.32131206	0.00022495	0.02817639
socs3	2.49813522	9.29 × 10−06	0.00306778	Cd48	7.59418959	0.00023965	0.02877211
ifi30	2.43485712	9.42 × 10−06	0.00306778	tnfrsf11b	2.68764743	0.0002414	0.02877211
CHGB	2.68694397	1.04 × 10−05	0.00327812	C10orf10	1.97584087	0.00024149	0.02877211
Cntn4	6.69229158	1.41 × 10−05	0.00430415	KCNN4	−8.0092773	0.00024933	0.02934901
Il1r1	2.38670209	1.64 × 10−05	0.00484285	IL8	4.44362712	0.00025339	0.02947116
AG2	7.88467942	1.87 × 10−05	0.00507962	ZFP36	1.96860745	0.00026324	0.02998838
GADD45A	2.30721084	1.89 × 10−05	0.00507962	DLK1	1.87823533	0.00026397	0.02998838
LMO1	2.31100372	1.97 × 10−05	0.00507962	GALR3	−3.0699664	0.00027127	0.0304633
TNFAIP3	−3.1164769	2.00 × 10−05	0.00507962	ANKRD31	7.59845327	0.00027525	0.03055914
ANKRD1	7.90497765	2.01 × 10−05	0.00507962	TRIB1	2.3659021	0.00029898	0.03274912
gpr84	2.20041148	2.08 × 10−05	0.00507962	apol3	2.42154371	0.00030168	0.03274912
NUP35	2.25111635	2.11 × 10−05	0.00507962	PPARGC1A	−7.848901	0.00031593	0.03366667
LCNL1	−8.1713233	2.13 × 10−05	0.00507962	PTCHD1	2.54755259	0.00031702	0.03366667
cebpd	2.21602916	2.28 × 10−05	0.00530536	EFNA5	−2.683381	0.00037376	0.03926515
NR5A2	3.01054291	2.37 × 10−05	0.00537637	LGALS3	1.95933307	0.00038655	0.04017603
TMEM49	2.23411164	3.46 × 10−05	0.00767944	c2cd4c	1.86851938	0.00039556	0.04042629
GZMB	4.67164885	3.73 × 10−05	0.00810511	NFIL3	−3.0702127	0.00039723	0.04042629
SLC16A10	6.39312229	4.21 × 10−05	0.00894971	WNT6	−1.9481788	0.0004573	0.04606028
ptgds	−2.2794722	4.40 × 10−05	0.00913917	CAMP	7.71716682	0.00047286	0.04714096
HPGDS	8.25515497	5.22 × 10−05	0.01062847	ism1	1.96334052	0.00050099	0.04944076
CD163	2.06718683	5.90 × 10−05	0.01153785

NOTE.—FC means fold change.

F

Boxplot diagrams of four representative differentially expressed genes involved in ovarian function and folliculogenesis. Juvenile expression data are in light blue, and adult expression data are in dark blue.

Changes in gene regulation could reflect differences in alternative splicing and exon usage between juveniles and adults in addition to changes in transcript abundance (e.g., Barberan-Soler and Zahler 2008; Revil et al. 2010). To investigate this possibility, we looked for differential exon expression (FDR-adjusted P value < 0.05)—a proxy for alternative splicing in a transcriptome without alternative splicing gene models—in genes with more than one exon. Specifically, we identified cases in which at least one exon, but not all exons, were differentially expressed (supplementary table S4, Supplementary Material online). To avoid false positives due to limited power, if one exon was differentially expressed, we relaxed the FDR-adjusted P value for differential expression to 0.15. Thus, evidence for exon-specific differential expression required relatively strong evidence for differential expression in at least one exon and a relative absence of evidence for differential expression in at least one other exon. Twenty-four genes exhibited this pattern, including STC (stanniocalcin) and GCLC (gamma-glutamylcysteine synthetase, catalytic subunit), both of which are thought to be important in ovarian development and function (Paciga et al. 2002; Luderer et al. 2003; Luo et al. 2004; Hoang et al. 2009).

Differential Expression in the Ovary and Lineage-Specific Selection in Primate Noncoding Regions

Many of the genes expressed in juvenile and adult baboon ovaries are also likely to be expressed in juvenile and adult ovaries of other primates, including humans. Thus, genes that we identified as differentially expressed across life history stages in baboons might be informative for identifying genes important in female life history evolution in humans or in primates more generally. To gain insight into the patterns of natural selection that may have acted on such genes, we therefore integrated the novel functional data from this study with evidence for selection in primates from previous studies.

We obtained estimates of positive selection on the lineage leading to humans for protein-coding regions from Nielsen et al. (2005) and for putative regulatory regions 5 kb upstream of genes from Haygood et al. (2007). Both studies took a similar approach to identify selective targets: Specifically, they compared the rate of nucleotide evolution in the focal region (protein-coding regions in Nielsen et al. 2005 and upstream regulatory regions in Haygood et al. 2007) with the rate of nucleotide evolution in a nearby region thought to be evolving neutrally (the general approach is reviewed in Yang and Bielawski 2000). An elevated rate of nucleotide evolution in the focal region relative to the nearby neutral region was interpreted as a signature of adaptive change. Likelihood ratio tests were then used to identify cases in which these rates differed across different branches of a species tree; we identified possible targets of lineage-specific selection by locating elevated rates of evolution in protein-coding or regulatory regions that occurred only on specific branches of the tree.

Combining our data with results from these studies (Nielsen et al. 2005; Haygood et al. 2007), we identified 225 genes that were both included in Haygood et al. (2007) and were differentially expressed in this study (P < 0.05 for differential expression; we relaxed this threshold to increase the sensitivity of this analysis). Of these 225 genes, we found 19 differentially expressed genes that were associated with signatures of selection in noncoding regions on the human lineage (P < 0.05 for the test for selection). In contrast, we found that none of our differentially expressed genes overlapped with signatures of positive selection in coding regions (of a total of 35 genes that were differentially expressed in this study and were included in Nielsen et al. (2005)). We did not observe a significant enrichment of ovarian differentially expressed genes among genes with a history of positive selection on the human branch. However, the target of selection in genes that were both differentially expressed between reproductively mature and immature ovarian tissue and also exhibited evidence for selection in the lineage leading to humans, was much more likely to have been a gene regulatory region than a coding region (Fisher's exact test, P = 2.367 × 10−08). If historical selection pressures on these loci were related to female maturation, changes in gene regulation may therefore have played an important role in the evolution of these traits in humans.

Genes expressed in reproductive tissue tend to be rapidly evolving exhibiting signatures of selection in multiple lineages (reviewed in Swanson and Vacquier 2002). We therefore examined whether differentially expressed genes were likely to be members of this rapidly evolving class or if they were specific to selection on the human branch. We asked whether noncoding regions that appear to have been positively selected on the chimpanzee (Pan troglodytes) lineage (Haygood et al. 2007) were similarly enriched for differential expression. We observed a similar number of differentially expressed genes by life history stage that correspond to positively selected regulatory regions in chimpanzees (21 in chimpanzees vs. the 19 seen in humans). Interestingly, 10 of these regions are shared between the two species, significantly more than expected by chance (hypergeometric test, P = 7.595 × 10−18; table 3). These results suggest that positive selection on the specific aspects of ovarian maturation controlled by these genes may be a general characteristic of primate evolution. Indeed, genes involved in reproductive and immune pathways that evolved under selection in humans are often also under selection in other primates (reviewed in Vallender and Lahn 2004) and in mammals more generally (Kosiol et al. 2008). Our data suggest that this shared pattern of positive selection may apply to regulatory regions of reproductively important genes as well.

Table 3

The Overlap Get of Genes That 1) Show Significant P Values for Selection in Noncoding Regions in Both Humans and Chimpanzees and 2) Also Show Evidence for Significant Differential Expression by Life History Stage in the Baboon Ovary Gene Expression Data

Gene ID	Selection P Value Human Noncoding	Selection P Value Chimpanzee Noncoding	P Value Differential Expression Baboons
serpina3	0.01	0.04	2.44 × 10−16
CHGA	0.01	0	0.01
LMO1	0.01	0.03	2.00 × 10−05
OSMR	0	0.03	0
DRG1	0.03	0.05	0.02
CAMP	0	0	0
dusp5	0	0.03	0.03
pfkfb3	0.03	0	0
SCUBE3	0.01	0.04	0.04
vwa2	0.02	0.01	0.04

Conclusion

The timing of female sexual maturity is one of many life-history traits that have shifted during primate and human evolution, probably in response to selection. Our results suggest that there has been repeated selection on the cis-regulatory regions of some sexual maturity-related genes in multiple primate lineages. These loci are therefore of special interest in relationship to phenotypic evolution during reproductive maturation. Thus, examining the overlap of signatures of selection and differential gene expression from samples obtained from natural populations may serve as a useful filter for identifying loci of particular evolutionary or phenotypic interest. Although such opportunities will be uncommon, they promise to enrich our ability to interpret the phenotypic relevance of sequence-based signatures of selection.

Materials and Methods

Study Subjects

Samples used in this study were obtained from seven healthy females from the Amboseli baboon population in Kenya (supplementary fig. S1 and table S1, Supplementary Material online) and retrieved within 5–8 h of death. Tissue was stored in RNAlater (Ambion, Austin, TX) and transported to −20 °C storage in Nairobi within 24 h. Upon transport to the United States, samples were stored at −80 °C. To minimize the effects of cell type heterogeneity in the ovary, we sampled from the lateral ovarian cortex.

Sample Preparation and Sequencing

Four micrograms of total RNA were isolated for each sample using an RNeasy kit (Qiagen, Valencia, CA) (supplementary table S1, Supplementary Material online) and used as input for the mRNA-Seq 8-Sample Prep Kit (Illumina, San Diego, CA). Libraries were sequenced on an Illumina GAIIx (one lane per sample) at the Yale University Keck Sequencing Core Facility. Approximately 15 million of 75 bp sequences resulted from each lane of sequencing.

Baboon Gene Models

The current publicly available baboon genome assembly (Pham_1.0, 20 November 2008) contains 387,373 linear scaffolds with approximately 5.3× coverage of the genome but has not yet been assembled into chromosomes (http://www.hgsc.bcm.tmc.edu/project-species-p-Papio%20hamadryas.hgsc). We mapped the RNA-Seq reads to the subset of these scaffolds (134,448 scaffolds with mean length of 20,246 bp) that mapped unambiguously to the macaque genome (Mmul_051212, rhemac2) using lastz (Harris 2007). Overall, the subset covered 94.9% of the current rhesus macaque assembly. Gene models were obtained by mapping human RefSeq exons to the baboon genome with lastz in Galaxy (Taylor et al. 2007) with a 90% similarity cutoff based on previous estimates of human–baboon sequence conservation (Silva and Kondrashov 2002).

Mapping Reads, Data Normalization, and Patterns of Differential Gene Expression

The RNA-Seq reads were mapped to the baboon scaffolds using “bowtie” (Langmead et al. 2009). Reads were defined as being within exon models using HTSeq (http://www.huber.embl.de/users/anders/HTSeq/doc/overview.html). Gene counts are the sum of the exon expression counts. The overall distributions of read counts were similar across all individuals and, more importantly, were not different between juveniles and adults, our primary axis of comparison (supplementary fig. S2, Supplementary Material online). Both exon counts and gene counts were normalized using the edgeR package (Robinson et al. 2010) in R (R Development Core Team 2008).

To evaluate the effect of maturation stage on specific genes, we used a generalized linear model with a negative binomial error structure to model variation in gene expression for each gene. Gene expression counts represented the response variable, and life history stage was modeled as a binary explanatory variable (juvenile or adult). We eliminated seven genes from this analysis that exhibited a significant relationship between gene expression and admixture-related genetic background as well as a relationship with life history stage (admixture between P. cynocephalus and a sister taxon, P. anubis, has previously been documented in this population and presented a possible confounder: Alberts and Altmann 2001; Tung et al. 2008). FDR corrections for multiple comparisons were performed using the Benjamini–Hochberg method (Benjamini and Hochberg 1995) at an FDR of 5% (fig. 1).

Categorical Enrichment Analyses and Alternative Exon Usage

To determine functional categorical enrichment for the differentially expressed genes, we employed the PANTHER (HMM Library Version 6.0) (Mi et al. 2005) and GO (Gene Ontology Consortium 2000) databases and computed enrichment scores using Wilcoxon rank tests. Our background set of genes was composed only of genes measured in this study.

Supplementary Material

Supplementary figures S1–S3 and tables S1–S4 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).