Abstract

Among closely related taxa, proteins involved in reproduction generally evolve more rapidly than other proteins. Here, we apply a functional and comparative genomics approach to compare functional divergence across a deep phylogenetic array of egg-laying and live-bearing vertebrate taxa. We aligned and annotated a set of 4,986 1 : 1 : 1 : 1 : 1 orthologs in Anolis carolinensis (green lizard), Danio rerio (zebrafish), Xenopus tropicalis (frog), Gallus gallus (chicken), and Mus musculus (mouse) according to function using ESTs from available reproductive (including testis and ovary) and non-reproductive tissues as well as Gene Ontology. For each species lineage, genes were further classified as tissue-specific (found in a single tissue) or tissue-expressed (found in multiple tissues). Within independent vertebrate lineages, we generally find that gonadal-specific genes evolve at a faster rate than gonadal-expressed genes and significantly faster than non-reproductive genes. Among the gonadal set, testis genes are generally more diverged than ovary genes. Surprisingly, an opposite but nonsignificant pattern is found among the subset of orthologs that remained functionally conserved across all five lineages. These contrasting evolutionary patterns found between functionally diverged and functionally conserved reproductive orthologs provide evidence for pervasive and potentially cryptic lineage-specific selective processes on ancestral reproductive systems in vertebrates.

1. Introduction

Over the past 550 million years, evolutionary processes have generated a diverse array of vertebrate species. Taxa that include fishes, birds, reptiles, and mammals evolved unique suites of adaptations allowing them to prosper in the most extreme sea, air, and land environments. Vertebrate diversity spans morphological innovations, developmental programming, cellular responses, as well as behaviors and life histories, and such differences become increasingly evident when taxa are compared across deep phylogenies. Studying the evolutionary patterns of functional change across this subphylum provides an opportunity to understand the evolutionary processes that have been important throughout vertebrate evolution. Yet, to date, clear common functional signatures that are in rapid flux across all vertebrate taxa have not been identified indicating the historical presence of a variety of niche- and lineage-dependent selective processes.

While functional evolutionary signals are not apparent across diverse phylogenetic lineages, when more closely related species such as sister species or multiple species within a single genus are compared, reproductive traits consistently reveal high diversity among species. This reproductive signature has been known for centuries, beginning with Linnaeus' binomial classification system [1]. Charles Darwin, in his 1871 treatise on sexual selection, also catalogued highly differentiated secondary sexual organs between closely related bird and mammal species [2]. Over a century later, William Eberhard described the diversity of morphological differences found in male secondary sexual traits, including vertebrate genitalia [3]. Both Darwin and Eberhard explain this higher male variance as the result of female mate choice or male-male competition on sexually selected traits within populations. The last three decades have amassed more vertebrate examples including cichlids [4], frogs [5], and primates [6] indicating that selection on reproductive traits may be a common underlying evolutionary process in vertebrates.

Studying rates of morphological character change demonstrates how certain functional classes evolve relative to others and provides a lens into evolutionary processes of the past. While this framework works well on closely related species, signatures diminish when applied to distantly related taxa due to the presence of lineage-specific rates of development, selective constraints, and genetic architectural differences [7]. In addition, there are many processes in which the selected phenotype may be hidden or cryptic to human observers. Such phenotypes often occur at the molecular level and include immune response [8], gametic interactions [9], and pheromonal exchange [10]. To systematically understand the relative roles of different functional classes in the evolutionary history of vertebrates, and hence the role of certain selective processes, it would be instructive to employ a common and unbiased framework on a representative sample of taxa.

With the availability of annotated genomic sequences across an ever-expanding number of taxa in addition to associated functional data (e.g., ESTs, GO) that can link genes to function, an operational framework is emerging that compares rates of functional change across varying degrees of phylogenetic relatedness [11, 12]. By applying this functional and comparative genomics approach, we now can use normalized information from sequences to infer how functional categories of genes have changed in the past. Combining the two domains of time and function can provide valuable information about the history of these lineages, in particular, how certain selective forces act upon certain reproductive processes such as gamete recognition, oogenesis, spermatogenesis, and adult behavior.

In this paper, we quantify the rates of change among reproductive and non-reproductive genes in five distantly related vertebrate lineages. We functionally categorize ~5,000 orthologs using available testis, ovary, and non-reproductive EST libraries in each species and find that individual vertebrate lineages generally follow a pattern of greater divergence in genes solely expressed in the gonads compared to genes expressed in non-reproductive tissue. In most cases, the testis appears to be driving gonadal divergence. However, an opposing pattern emerges when we compare evolutionary rates among the much smaller subset of tissue-expressed genes that have remained functionally conserved across vertebrates ( ). Using this framework, we are beginning to unmask a pattern of rapid and cryptic molecular evolution on lineage-specific reproductive features that are part of conserved developmental processes, thus, providing a common underlying genetic basis of functional evolutionary change in the vertebrate subphylum.

2. Materials and Methods

2.1. Orthology and Estimates of Divergence

Protein coding genes from A. carolinensis, G. gallus, D. rerio, M. musculus, and X. tropicalis were used in this analysis. Orthologs for each species pair were obtained from BioMart (http://uswest.ensembl.org/biomart/index.html). Orthologs were filtered so that only transitive sets of 1 : 1 : 1 : 1 : 1 orthologs remained, producing 4,986 sets of 5-species orthologs. We excluded all paralogous relationships (including 10,122 1 : 1 :  2 : 1 : 1 relationships, where “2” denotes paralogous sequences from the zebrafish lineage) in order to maintain a relatively ambiguous ortholog set. The protein coding CDS and amino acid sequence of each gene’s longest transcript were also obtained from BioMart: in the case of transcript length ties, the transcript with the lower incremental Ensemble ID number was used. Multiple sequence alignments for each orthologous set of proteins were generated using MUSCLE (version 3.8; [13]) and then back-translated using corresponding CDS and a custom Perl script (available from CJG on request). All 1 : 1 : 1 : 1 : 1 alignments in addition to their associated functional assignments will be made available via lizardbase (http://www.lizardbase.org/) as an active link to current A. carolinensis annotations in lizardbase’s genome browser, JBrowse, and lizardbase’s Resources Page. All alignments will also be made available on the Resources page in lizardbase.

A protein distance matrix was calculated for each protein alignment using the Jones-Taylor-Thornton (JTT) model in the prodist program from the Phylip suite of phylogenetic programs (version 3.69; [14]). Consensus phylogenetic trees were generated using concatenated sequences from both CDS and its associated protein sequences (See in Supplementary Material available online at http://dx.doi.org/10.4061/2011/274975 Supplementary Figure  1). For a given gene from each species, the mean of its four orthologous protein distances was used as one of two estimates of sequence divergence. A matrix of nonsynonymous substitutions per nonsynonymous site, dN, was also estimated for each codon alignment using Nei and Gojobori’s method [15] using the SNAP Perl program [16], and its mean dN across four orthologs was used as an estimate of sequence divergence.

2.2. Functional Annotation Using EST Libraries and Gene Ontologies

ESTs from each of the five species were filtered as “normal adult” tissue from NCBI’s dbEST (downloaded in October 2009) and assigned to species-specific tissue libraries (see Supplementary Table  1) based on either organ or tissue fields in the Genbank record. EST sequences were locally indexed and aligned to genes from the same species using a standalone version of blastn (version 2.22; [17]). EST-to-gene alignments of at least 100 nucleotides, 90% identity, and an -value of were used as alignment criteria. For each of the five species lineages, genes with at least three ESTs (i.e., hits) meeting the above alignment criteria were assigned to seven non mutually exclusive functional classes: (1) genes with hits in only the testis were classified as testis-specific; (2) genes with hits in the testis and another tissue(s) were classified as testis-expressed; (3) genes with hits in only the ovary were classified as ovary-specific; (4) genes with hits in the ovary and another tissue(s) were classified as ovary-expressed; (5) genes with hits in only the testis and/or ovary were classified as gonadal-specific; (6) genes with hits in the testis and/or ovary, in addition to non-reproductive tissue(s), were classified as gonadal-expressed; (7) genes with hits from an assortment of non-reproductive tissues (see Supplementary Table  1) that were neither testis nor ovary were classified as non-reproductive. Thus, for each of the five species, genes with sufficient EST coverage fell into at least one functional class (Table 1). The difference between the mean dN of each reproductive class and the mean dN of the non-reproductive class was tested using an unpaired two-sample two-sided Wilcoxon rank sum test ([18]; Figure 1).

We also compared evolutionary rates in functionally conserved genes, that is, those orthologs that do not change functional class across all five lineages, according to our EST annotations. Interestingly, we were not able to identify a single gonadal-specific gene, but were able to identify functionally conserved subsets of testis-expressed ( ), ovary-expressed ( ), and non-reproductive ( ) orthologs. Figures 2(a)2(g) provides Venn diagrams for all species combinations in each functional class. Figure 3 compares dN across four (nonzero) functional classes.

To complement the functional annotations generated by ESTs, we linked the 10% most diverged orthologs to the GO categories, Biological Process (BP) and Cellular Component (CC) in each species. GeneMerge [19] was used to test for statistically significant over-represented functional terms. A “word cloud” that relates the frequency of each GO term to its font size was generated for four of the five species (Figure 4, Supplementary Figure  2). X. tropicalis was excluded from this analysis due to its sparse GO term set.

3. Results and Discussion

In this study, we chose five distantly related vertebrate species that fit the following criteria: (1) the presence of a well-assembled and freely available genome sequence, (2) the existence of well-curated gene models, (3) the availability of appreciable numbers of testis, ovary, and non-reproductive ESTs at dbEST, and (4) the condition that all five species, together, represent divergent clades thus presenting a deep vertebrate phylogeny with a diverse breadth of functional differences. After filtering out alignments that were of poor quality or had ambiguous orthologous relationships, a consensus tree-based off-concatenated CDS sequences from 4,986 orthologs was generated using the five vertebrate species. The tree’s topology was well supported in 100% of 1000 bootstrap replicates (Supplementary Figure  1). A concatenated protein tree-demonstrated the same topology and support (not shown) and mirrored published vertebrate phylogenies (e.g., [20]). We note that these ~5,000 orthologs represent a relatively “well-behaved” and conserved gene set that do not possess paralogs in any of the five lineages. This study focuses on 1 : 1 : 1 : 1 : 1 orthologs and ignores complications arising from neo-/subfunctionalization caused by gene duplication events [21, 22], particularly those found in the zebrafish lineage after an ancient duplication event [23].

We used the extensive EST libraries publically available for each species in order to categorize genes into functional classes. Our objective was to generate a standardized sample of genes in each of the reproductive and non-reproductive functional classes, for each species. Historically, EST libraries were originally developed to assist in the genome annotation process (e.g., [24]). The quantity, quality, specificity, and tissue-diversity of EST libraries vary considerably across species (see Supplementary Table  1) and are largely a function of each research community’s priorities and preferences for each of the five sequenced genomes. Since our principal objective is to compare reproductive versus. non-reproductive levels of molecular divergence in vertebrates, we sought to generate pooled gene samples derived from testis and ovary (i.e., reproductive) and non-reproductive tissue (any adult tissue that does not contain a sex-specific organ or tissue). In addition, genes from tissue-specific (or tissue-limited) classes were differentiated from “tissue-expressed” genes that are expressed more ubiquitously. This approach enables us to compare functional gene classes using relatively large sample sizes and ample statistical power.

A total of seven functional classes were assigned to genes in each of the five species (see Section 2). Table 1 summarizes the number of genes that are contained in each functional class for each species. It is important to note that the proportion of reproductive (e.g., testis, ovary) to non-reproductive genes in each species is not necessarily indicative of the total fraction of reproductive genes found in each genome but, again, reflects each community’s specialized interests in generating certain libraries. In addition, overall EST library coverage can be different by an order of magnitude. For example, at last count, the mouse has nearly 5 million ESTs deposited in dbEST, while the green lizard has only 150,000 ESTs. The broader EST coverage in mouse may explain why our screen failed to identify any ovary-specific genes in this taxon. In contrast, since the anoles EST set includes only three non-reproductive tissues at a lower coverage than other species, this may also explain the relatively high number of ovary-specific genes in this species. With such large differences in EST coverage in each of the five species, it is important to understand the limits of these analyses.

Overall, our results provide evidence of a general pattern of rapid reproductive change over deep vertebrate lineages. Each of the five vertebrates demonstrate significantly higher protein divergence in gonadal genes compared to non-reproductive genes (Figure 1). Rapidly evolving testis genes appear to be driving much of the pattern of higher gonadal-specific gene divergence in these lineages: four of the five taxa—zebrafish, Xenopus, chicken, and mouse—all share significantly higher testis-specific divergence. Interestingly, these three taxa include two of the more basal taxa, Xenopus and zebrafish (Supplementary Figure  1), supporting that this pattern spans broad phylogenetic groups across the vertebrate subphylum. In green lizards, we observe a contrasting pattern of gonadal divergence as ovary-specific genes appear to be driving the significantly higher divergence of gonadal genes (but see caveat above). Thus, while we see a general pattern of significantly higher divergence among reproduction-specific genes across all vertebrate lineages, there may be large differences in the subset of reproductive genes that are diverging.

In Drosophila, we also see a similar pattern of rapidly evolving gonadal genes from EST libraries. Reproductive genes from the testis and ovary and non-reproductive genes from the brain have been used to characterize sexually dimorphic expression patterns [2527] as well as to compare the evolution of reproductive genes relative to non-reproductive genes [2830]. A recent study using 12 genomes in Drosophila and an extensive EST set from D. melanogaster also found that rates of evolution among testis-expressed genes are significantly higher than genes expressed in the ovary or head [12]. A number of studies in mammals have also demonstrated a similar pattern of higher divergence rates in male reproductive genes [11, 3133].

This higher divergence of reproductive genes, and in particular, male-specific proteins, supports the hypothesis that sexual selection may be an important driver of evolutionary change and extends sexual selection theory to the level of molecules such as those found in gametogenesis and fertilization [3436]. The strength of this molecular signature indicates the pervasive and cryptic nature of this process: much of this pattern would remain hidden without a comparative and systematic treatment of genome-wide sequence data. We also note that reproductive proteins, particularly those regulating sperm development, are of particular interest to researchers studying mechanisms of reproductive isolation because hybrid male sterility may be the product of the rapid evolution of male reproductive genes: spermatogenesis appears to be a selected target of hybrid male fertility breakdown [3741]. In addition, there is mounting evidence that positive selection drives the evolution of genes controlling key transitions during both spermatogenesis and oogenesis [42, 43].

Other functional classes of testis-associated genes have also been found in Drosophila. Genes encoding proteins secreted by male accessory glands (Acps), the ejaculatory duct, and the ejaculatory bulb, as well as many components of D. melanogaster seminal fluid, are known to be rapidly evolving. These proteins are transferred from the male to the female along with sperm during mating and mediate a series of postmating events [4446]. Furthermore, there is ample evidence of adaptive evolution at several loci that encode D. melanogaster seminal fluid proteins [4752]. Whether a similar signal among secretory reproductive classes is found in vertebrate lineages is an intriguing question.

While a clear pattern of rapid testis-specific divergence emerges from our lineage-specific annotations, we then asked whether the same evolutionary pattern holds across genes that have maintained a similar function across all five vertebrate species. In other words, what are the relative rates of evolutionary change across functionally conserved classes? Surprisingly, the numbers of genes per class were drastically reduced to the point that only four classes—testis-expressed, ovary-expressed, gonadal-expressed, and non-reproductive genes—share genes in common across all five species (Figure 2). Furthermore, a decreasing but nonsignificant trend of evolutionary rates was found among these four functional classes: . Overall, this functionally conserved group describes a subset of the data with a contrasting evolutionary pattern, thereby demonstrating that testis-specific genes are affected by a variety of evolutionary forces. In a recent study, Dean et al. [33] performed a genomic and proteomic study on six tissue types from the male reproductive tract of mouse (excluding testis) and found that one tract, the seminal vesicle, had significantly higher rates of divergence while the other five tracts showed significantly lower rates of divergence when compared to other proteins. Our results demonstrate a similar high variance of evolutionary rates within the testis.

A. carolinensis was the outlier of the five vertebrate taxa with a significantly higher divergence among ovary-specific genes. Ovaries have also been shown to be sites of rapid divergence in D. melanogaster as part of a molecular coevolutionary process between sperm and egg. A number of rapidly evolving genes have been found expressed in the female reproductive tract and potentially secreted [53] or induced in the female reproductive tract by mating [5456]. Further characterization of the green anole genome, in addition to other Lepidosauria genomes and genomic resources that will soon be available, will allow us to address whether female lizards are indeed driving sexual selective processes and whether this is a common lineage-specific process among squamate reptiles or simply an artifact of EST functional annotation.

While aligning orthologs to ESTs offers a powerful approach for functional annotation, it is important to procure a more granular understanding of process, function and localization. Therefore, we took the 10% most diverged orthologs and associated each species’ corresponding gene to its Gene Ontology (GO), namely, Biological Process (BP) and Cellular Component (CC). A word cloud in which the font size is a function of the frequency of statistically over-represented functional phrases in diverged orthologs is shown for A. carolinensis in Figure 4. GO-associated word frequencies for fish, mouse, and chicken are found in Supplementary Figure  2. The density of each word cloud for each species reflects the amount of curation effort in Gene Ontology within these species’ communities. As expected, we don’t see much overlap between the GO and EST approaches to functional annotation. Reproductive function is a poorly annotated ontological class harboring a level of characterization that will not substantially improve until more geneticists and molecular biologists study reproductive loci in greater detail.

4. Concluding Remarks

Patterns parsed from extant genomes can inform us about the underlying evolutionary processes that have acted upon lineages in the past. As a functional class, reproductive-specific genes are more rapidly evolving than other functional gene classes, and it appears that testis genes are driving this pattern of divergence in the majority of vertebrate lineages. This work sets the stage for a more nuanced analysis of divergence leveraged against function across diverse taxa. With more genomes and ESTs generated, greater effort can be afforded to better estimate the probability that a gene is a member of a particular functional class, even when the number of ESTs and libraries are quite different between species. Newer data types such as RNAseq will certainly help solve the sampling bias problem with better coverage and more tissues sampled. Future studies that include paralogous sequences to evaluate birth/death processes and de novo gene functionalization models (including incorporating the large number of paralogs from zebrafish) in the context of functional class will also be useful in addressing the role of reproductive genes in vertebrate evolution.

It is remarkable that across very distant phylogenetic lineages, we detect the same evolutionary patterns found among closely related species: high lineage-specific reproductive diversity and, in particular, a high variance in male reproductive characters. These parallel patterns support the contention that sexual selection on both morphological and molecular characters may be an important, common, and pervasive feature of vertebrate evolution.

Abbreviations

GO:Gene ontology
EST:Expressed sequence tag
CDS:Coding sequence
BP:Biological process
CC:Cellular component
dN:Nonsynonymous substitutions per nonsynonymous site.

Acknowledgments

The authors would like to thank Dr. Tonia Hsieh (Temple University) for guidance and support during all stages of this project. They also thank Dr. Ed Braun (University of Florida) for technical advice and the Beaty Biodiversity Research Centre at the University of British Columbia for use of their SciBorg cluster. C. J. Grassa dedicates his contributions to Thomas and Sarah Grassa.

Supplementary Materials

A catalog of the EST libraries used in this study are listed in Supplementary Table 1.

Supplementary Figure 1 displays the five-species vertebrate phylogeny based off concatenated CDS sequences.

The word-size frequency distribution of Gene Ontology (GO) terms for the most diverged orthologs in M. musculus, G. gallus, and D. rerio are found in Supplementary Figure 2.

  1. Supplementary Materials