Introduction Together with more subtle genetic modifications such as gene expression changes and point substitutions, new genes with novel functions may have significantly contributed to the evolution of new phenotypes specific to humans and their closest evolutionary relatives. New duplicate genes may originate through (segmental) gene duplication by intra- or interchromosomal transposition of gene-containing segments [1,2]. Another mechanism, retroposition, generates new intronless gene copies (retrocopies) by reverse transcription of mRNAs derived from source genes (“parental” genes), followed by reintegration of the resulting cDNA in the genome [2,3]. Retroposition was commonly thought to generate nonfunctional gene copies (retropseudogenes) that accumulate disablements such as premature stop codons and frameshift mutations [4], because the copied mRNA is generally lacking regulatory elements. However, we and others have recently shown that retroposition has generated a significant number of new functional genes (retrogenes) in mammalian and invertebrate animal genomes [3,5,6]. Multiple studies have suggested a high rate of retroposition on the primate and rodent lineages [7−9], probably driven by the activity of L1 retrotransposable elements [10]. Thus, retroposition may also have provided abundant raw material for the formation of new genes on the primate lineage leading to humans, potentially generating many more retrogenes than the four primate-specific retrogenes (present in the human genome) with functional roles and/or expression in testis, brain, and lymphocytes previously described [11−14]. To assess the importance of retroposition for the creation of new genes on the primate lineage leading to humans, we systematically screened the human genome for retrogenes that emerged during the primate burst of retroposition. Our results suggest an important role of retroposition in the formation of new genes and phenotypes in the recent evolution of the human genome. Results Age Distribution of Human Retrocopies We identified 3,951 retrocopies (and their corresponding parental genes) in the human genome using a refinement of a previously published procedure [5] (see Materials and Methods). Among these, 705 retrocopies (∼18%) are found to be “intact,” i.e., they show no disablements such as premature stop codons or frameshift mutations when compared to the open reading frame (ORF) of their parental genes. To assess the age distribution of retrocopies, we calculated nucleotide divergence at silent sites (K S) between retrocopies and their parental genes (Figure 1). Assuming neutral mutation rates of 1−1.3 × 10−9 substitutions per site per year [15], the high number of retrocopies with K S ≍ 0.1 suggests that the burst of retroposition reached its peak approximately 38−50 million years ago (MYA) on the primate lineage, in agreement with previous estimates [7,8]. The vast majority of retrocopies (91%) also show a divergence at silent sites much lower than that observed between human and mouse genes (Figures 1 and S1), indicating that they arose after the human−mouse split. Therefore, our data are consistent with a high retroposition activity on the primate lineage. Rate of New Retrogene Formation To estimate the number of recent human retrogenes, we compared signatures of selective constraint between intact, potentially functional retrocopies, and retropseudogenes (assumed to be nonfunctional, i.e., evolving neutrally). To this end, we calculated the ratio of nonsynonymous to synonymous substitutions per site (K A/K S) for retrocopy/parental gene pairs with a synonymous divergence of less than 0.15. This value approximately reflects the deepest neutral divergence in the primate tree between humans and the most divergent extant primate lineage [16,17], the lemurs, and corresponds to around 63 million years of primate evolution [18]. This analysis reveals a difference in the K A/K S distributions between intact copies and retropseudogenes (which may show low K A/K S ratios by chance), with a highly significant excess of intact copies for K A/K S 600 bp) and age (>8 million years; i.e., presence in humans and African apes), characteristics estimated to provide sufficient statistical power for the simulation approach (see Materials and Methods). We sequenced these copies in all species carrying them. Sequence alignments show that eight of these 23 retrocopies are intact in all species, whereas the remaining copies carry one or more stop codons and/or frameshift mutations in one or more lineages (Table 1). Next, we used a simulation approach (see Materials and Methods), which is based on the basic assumption that under neutrality, an intact retrocopy will accumulate deleterious mutations (stop codons or frameshifts) over time that will disrupt its ORF and may eventually preclude gene function, whereas under functional constraint, natural selection will prevent the accumulation of deleterious mutations in the retrocopy sequence. Our simulation approach estimates the probability that a gene copy would have retained its ORF since the duplication event, in all or most species in which it is present, if it had evolved neutrally along all lineages of the species tree. In parallel, this approach tests whether the number of nonsynonymous substitutions that accumulated since the retroposition event along the different branches of the species phylogeny is consistent with neutral evolution. The simulations revealed that seven retrocopies are unlikely to have remained intact in all (or most) species if they had evolved neutrally throughout their evolutionary history, even after correcting for multiple tests (p 1, p 1) and/or a neutral phase of evolution after duplication [3,12]. Second, retrocopies with disablements in their ORFs (as defined by their parents) are treated as pseudogenes in this analysis, although new retrogenes may emerge from truncated coding regions [3,13]. It is also known that new splicing signals in a coding region that contains frameshifts or premature stop codons may evolve to define a new intron or to generate chimeric transcripts with nearby or “host” genes [3]. Finally, duplicate “pseudogene” copies may play functional roles by virtue of their RNAs regulating closely related paralogous genes [39,40]. At any rate, our results suggest that in addition to other types of duplications [1], retroposition significantly contributed to new gene formation in primates. Retrogenes and Male Functions It is remarkable that all seven retrogenes identified in this report are expressed predominantly or exclusively in testis, whereas their parents are all expressed ubiquitously. A preliminary survey of retrocopy transcription using expressed sequence tag databases suggests that this observation may reflect a general pattern (data not shown). Several factors may contribute to this effect. For example, chromatin remodeling [41] and abundance of RNA polymerase II complexes during late phases of male meiosis [42] lead to a state of “hypertranscription” [43], which may allow retrocopies to become initially transcribed in testis. This may also have facilitated transcription of new genes arising from pericentromeric segmental duplications [44,45]. Thus, there is a mechanistic bias that may favor testis expression of new genes. However, our results suggest that testis expression is often not merely a by-product of new retrogene formation but that natural selection may have favored the recruitment of testis-specific regulatory elements to enhance the beneficial effects of the initial mechanistically driven testis transcription. Consistently, we can infer a testis function for five of the seven primate retrogenes identified here and for two of the four previously identified retrogenes (TAF1L and UTP14C; [13,14]). Five retrogenes (eIF-2-gamma2, RBMXL1, KIF4b, TAF1L, and UTP14C) stem from the X chromosome and probably either substitute for their parental genes during male meiosis [30] or otherwise enhance male germline function [46]. For one retrogene (GMCL2), a function in sperm formation can be postulated based on studies of parental orthologs. Finally, PABP3 functionally adapted to late spermatogenesis both on the protein sequence level and by developing a highly specific expression pattern [20]. Sex- and reproduction-related genes are generally recognized as a class of rapidly evolving genes, particularly genes involved in male reproduction [47]. Possible causes include sperm competition, sexual conflict, and selection for reproductive isolation [48]. A comparison of the human and mouse genomes revealed an excess of lineage-specific expansions of genes related to reproduction as well as an accelerated protein evolution of such genes [49]. Together, these observations suggest that duplicate gene copies may have provided important raw material for rapid testis evolution in primates. Specifically, gene duplication may allow one copy of the duplicate pair to specialize in testis function, while the other is selectively preserved to sustain a role in somatic tissues [50−52]. Our data suggest that retroduplication may have provided a means to allow for such decoupling of functions in primates. Indeed, we show that selection to attain enhanced male germline function has progressively fixed and adapted retroposed gene copies on the primate lineage leading to humans. Materials and Methods Retrocopy screen We retrieved all peptide sequences (categories: known and novel) from the Ensembl ([53]; http://www.ensembl.org/index.html) database (version 29). To screen for retrocopies, these peptide sequences were used as queries in translated similarity searches against the complete human genome (NCBI genome release 35) sequence using tBLASTn [54]. Adjacent homology matches were merged in a series of parsing steps using Perl scripts, combining only nearby matches (distance 50%) and aligned to one another over more than 70% of the length of their sequence (minimum length: 50 amino acids). Next, we performed similarity searches of the merged sequences against all Ensembl genes (intron-containing and intronless) using FASTA. We kept only copies where the closest hit was an Ensembl peptide with multiple coding exons (putative parental gene). Merged sequences for which the closest match was an intronless gene were excluded from the data (e.g., to avoid intronless genes of other types such as olfactory receptor genes). We also confirmed the absence of introns in these retrocopies by mapping parental intron locations onto the alignments. We required that parental introns map within the alignments between parents and retrocopies and be larger than 80 bp. This threshold was chosen to ensure that real introns are missing in the retrocopies; 80 bp is larger than the gap size (40 bp) allowed in the merging step, it avoids mapping of small gaps in parental exons erroneously annotated as introns, and it takes into account that the majority of human introns are ∼80 bp or larger [55]. Samples Primate DNA samples were mainly obtained from the E CACC repository (Wiltshire, United Kingdom): chimpanzee (Pan troglodytes), gorilla (Gorilla gorilla), orangutan (Pongo pygmaeus), gibbon (Hylobates lar), Old World monkey (African green monkey, Cercopithecus aethiops sabaeus), and New World monkey (owl monkey, Aotus trivirgatus). Lemur (Lemur catta) and tupaia (Tupaia glis) DNA samples were obtained from Institut des Sciences de l'Evolution, Montpellier University 2. PCR and sequencing reactions PCR amplifications were performed in a Mastercycler gradient (Eppendorf, Hamburg, Germany) using either Taq DNA Polymerase or ProofStart DNA Polymerase from Qiagen (Valencia, California, United States). PCRs were performed according to the instructions of the manufacturer. For sequencing, amplified PCR products were reamplified using a pair of nested primers. The resulting PCR products were purified using the MinElute PCR Purification Kit or QIAquick Gel Extraction Kit from Qiagen. From these PCR products, both strands of the retrogene coding sequence were determined using the BigDye 3.1 cycle sequencing kit (PerkinElmer, Wellesley, California, United States). The sequencing reactions were run on an ABI 3730 automated sequencer (Applied Biosystems, Foster City, California, United States). Parental and retrogene expression patterns were analyzed using PCR and a cDNA panel of 20 different human tissues. Experiments were repeated twice to confirm the expression pattern. Unique primer pairs were designed for both parental gene and retrogene, based on ClustalX alignments of parental and retrogene cDNA sequences. The cDNA panel was synthesized using the FirstChoice Human Total RNA Survey panel from Ambion (Austin, Texas, United States) and a SuperScript II First-Strand Synthesis System RT-PCR (Invitrogen, Carlsbad, California, United States). Reactions without reverse transcriptase were done in parallel as negative controls for all 20 tissues. RT-PCR amplifications were performed in a Mastercycler gradient (Eppendorf) using JumpsStart DNA Polymerase (Sigma-Aldrich, St. Louis, Missouri, United States) using standard conditions as recommended by the supplier. Products were purified using the MinElute PCR Purification Kit from Qiagen and sequenced using the same pair of primers. Obtained sequences for each retrogene were then aligned with both retrogene and parental gene sequences using ClustalX. To ensure that RT-PCR products were derived from the retrogene, nucleotides at diagnostic sites that discriminate between retrogene and parental gene were manually confirmed. All oligonucleotide sequences used for PCR and sequencing are available upon request. Age of retrocopies We estimated the age of retroposition events by calculating coding sequence divergence at synonymous sites (K S) between each retrocopy and the corresponding parental gene. The same analysis was performed for parental genes and their mouse orthologs. Codon sequences were aligned on the basis of the translated sequence alignment using the EMBOSS package [56]. In all alignments, the coding sequence of the parental gene was used as a reference. Pairwise K S statistics were estimated using the YN00 program of PAML [57] version 3.14. We note that the ages of retrocopies may be slightly underestimated by this approach, because silent sites are not always completely neutral ([58] and references therein). Using a phylogenetic dating approach, we determined the age of individual retrocopies by screening for their presence or absence in primate genomes using PCR with primers flanking the insertion site. We confirmed that the insertion site in species not carrying the copy reflects the expected size of the ancestral state (before retrocopy insertion [12]). For five of the retrogenes analyzed in detail, the ancestral state of the insertion site was further confirmed by sequencing. For the two retrogenes (NACA2 and PABP3) present in all anthropoid primates (hominoids, Old World monkey, and New World monkey), we confirmed their absence in lemur and tupaia using several different primer pairs located in their coding regions, as the insertion site could not be amplified using primers in the flanking region. Rate of retrogene formation Pairwise K A and K S statistics for all retrocopies were estimated using the YN00 program of PAML [57] version 3.14. To estimate K A/K S on the retrocopy lineage itself, we performed the same analysis but compared the retrocopy and the ancestral sequence of the retrocopy at the time point of retroposition (estimated by a maximum likelihood procedure; using the codeml program of PAML [57] and the mouse ortholog of the parent as outgroup). K A/K S is influenced by the GC content at synonymous sites of the parent as well as by the GC content of the genomic region surrounding the retrogene [59]. In particular, retrocopies derived from parental genes with high GC that insert into regions of low GC may show low K A/K S driven by local adaptation to local GC. To test whether GC differences between intact and retropseudogene copies with low K A/K S ( 60% at 4-fold degenerate sites) parents that inserted into low-GC (lower than median value of GC) regions (52 of 130 intact retrocopies versus 60 of 172 retropseudogene copies, p = 0.4). Thus, the difference in the distributions for K A/K S 1) on the retrogene lineages, we compared model M1 and model A as implemented in codeml from the PAML package [57] using likelihood ratio tests [62]. Model M1 assumes two classes of sites for the sequences in the whole phylogeny: sites under purifying selection (K A/K S 1. We also compared this model A to a modified model where K A/K S is fixed at one. Sites under positive selection in the retrogene lineages were identified using the Bayesian approach as implemented in codeml [63]. Note that with respect to CDC14B2, the human and chimpanzee sequences have lost the original translation initiation codon (methionine) used by the parental gene (which may have led to the annotation of this gene as a VEGA pseudogene, OTTHUMG00000033880) and gained a putatively new methionine start codon at position 31. The selection tests show similar (statistically significant) results when either the original full-length sequence alignment or a shorter alignment starting from position 31 is used (data not shown). Supporting Information Figure S1 Distribution of K S between Parental Genes and Their Orthologs in Mouse (644 KB EPS). Click here for additional data file. Figure S2 K A /K S on the Retrocopy Lineage: Comparison of the K A /K S Distributions for Intact Retrocopies and Retropseudogenes The mode of the K A /K S distributions is smaller than one (usually expected under neutrality), owing to the effect previously described [59]. White bars correspond to intact retrocopies, and dark bars to retropseudogene copies. (636 KB EPS). Click here for additional data file. Figure S3 Simulation Results for CDC14B2 See legend of Figure 3. (7.9 MB PDF). Click here for additional data file. Figure S4 Simulation Results for eIF-2-gamma2 See legend of Figure 3. (6.6 MB PDF). Click here for additional data file. Figure S5 Simulation Results for GMCL2 See legend of Figure 3. (6.6 MB PDF). Click here for additional data file. Figure S6 Simulation Results for KIF4B See legend of Figure 3. (6.6 MB PDF). Click here for additional data file. Figure S7 Simulation Results for NACA2 See legend of Figure 3. (6.7 MB PDF). Click here for additional data file. Figure S8 Simulation Results for PABP3 See legend of Figure 3. (6.6 MB PDF). Click here for additional data file. Figure S9 Phylogenetic Trees for the Seven Retrogenes Identified in This Study Maximum likelihood K A /K S values and the estimated number of nonsynonymous versus synonymous substitutions (in parentheses) for each branch are indicated. (791 KB EPS). Click here for additional data file. Table S1 Dated Retrocopies (107 KB DOC). Click here for additional data file. Accession Numbers The GenBank (http://www.ncbi.nlm.nih.gov/Genbank/) accession numbers for the primate sequences generated for this paper are DQ120612−DQ120720. They are detailed in Table S1. The Ensembl (http://www.ensembl.org/) accession numbers for other genes discussed in this paper are GMCL1 (ENSG00000087338) and PABP1 (ENSG00000152520).
