1. Genes and Chromosomes
  2. Genomics and Evolutionary Biology
Download icon

The life cycle of Drosophila orphan genes

Research Article
Cited
34
Views
2,870
Comments
0
Cite as: eLife 2014;3:e01311 doi: 10.7554/eLife.01311

Abstract

Orphans are genes restricted to a single phylogenetic lineage and emerge at high rates. While this predicts an accumulation of genes, the gene number has remained remarkably constant through evolution. This paradox has not yet been resolved. Because orphan genes have been mainly analyzed over long evolutionary time scales, orphan loss has remained unexplored. Here we study the patterns of orphan turnover among close relatives in the Drosophila obscura group. We show that orphans are not only emerging at a high rate, but that they are also rapidly lost. Interestingly, recently emerged orphans are more likely to be lost than older ones. Furthermore, highly expressed orphans with a strong male-bias are more likely to be retained. Since both lost and retained orphans show similar evolutionary signatures of functional conservation, we propose that orphan loss is not driven by high rates of sequence evolution, but reflects lineage-specific functional requirements.

https://doi.org/10.7554/eLife.01311.001

eLife digest

New genes are added to most genomes on a steady basis. A new gene can either begin as a copy of an existing gene from elsewhere in the genome, or is created entirely ‘from scratch’ from a DNA sequence that had not previously encoded for a protein. New genes that are not found in other related species are called orphan genes—and these genes can account for up to 30% of all the genes in the well-studied genomes. However, for reasons that are not fully understood, the total number of genes in most genomes remains fairly constant despite these regular additions. Now, Palmieri et al. have investigated this paradox by following the evolutionary fate of orphan genes in a small group of related species of fruit fly.

Palmieri et al. discovered that most orphan genes are very short-lived, even though they showed clear signals of carrying out important functions. Most orphan genes died out quickly due to mutations that made them unable to be expressed as functional proteins, and a small number were deleted entirely from the genome. Unexpectedly, new orphan genes were more likely to die out than those that had been around for a while.

Palmieri et al. also found that the expression levels of orphan genes determined their probability of dying with those genes that were expressed to the highest levels being most likely to persist longer. Furthermore, genes that were expressed more in males than in females were also less likely to die. The next challenge will be to identify the mechanisms that determine which orphan genes survive and which do not.

https://doi.org/10.7554/eLife.01311.002

Introduction

Orphans are genes with limited phylogenetic distribution and represent a considerable fraction (up to 30%) of the gene catalog in all sequenced genomes (Khalturin et al., 2009). Studies conducted in different eukaryotes showed that orphans emerge at high rates (Domazet-Loso et al., 2007; Wissler et al., 2013). While gene duplication and exaptation from transposable elements often result in orphan genes (Toll-Riera et al., 2009), they also originate frequently de novo from non-coding DNA (Cai et al., 2008; Heinen et al., 2009; Knowles and McLysaght, 2009; Wu et al., 2011; Yang and Huang, 2011; Xie et al., 2012; Neme and Tautz, 2013; Wu and Zhang, 2013), probably through intermediate proto-genes (Carvunis et al., 2012). Compared to evolutionary conserved genes, orphans are overall shorter (Lipman et al., 2002), fast evolving (Domazet-Loso and Tautz, 2003), have lower and more tissue-restricted expression (Lemos et al., 2005). Moreover, they often show testis-biased expression (Levine et al., 2006; Begun et al., 2007), probably due to frequent origination in testis (Kaessmann, 2010).

In Drosophila the rate of orphan emergence is particulary high (Domazet-Loso et al., 2007) and many orphans become quickly essential (Chen et al., 2010). Although the function of only a few orphan genes has been studied, it has been proposed that orphans might serve an important role in speciation and adaptation to different environments (Khalturin et al., 2008; Khalturin et al., 2009; Colbourne et al., 2011). The high rate of orphan origination would predict an increase in gene content over time. However, gene content in eukaryotes is remarkably stable compared to genome size, as highlighted by Tautz and Domazet-Loso (2011). To solve this paradox Tautz and Domazet-Loso proposed that orphans have only a short lifetime (‘rapid-turnover’ hypothesis) (Tautz and Domazet-Loso, 2011). Thus, although orphans are continuously created, most of them might be lost in a relatively short evolutionary time. Relaxed selective costraints in orphans (Cai and Petrov, 2010) might also contribute to the high rate of orphan loss.

Moreover, since orphans are typically identifed by the comparison of distantly related species, their evolutionary stability has been so far neglected. This contrasts the comprehensive analysis of evolutionary patterns of gains and losses of non-orphan genes (Hahn et al., 2007). In this study, several partially interrelated factors affect gene loss, including gene expression levels, number of protein–protein interactions, gene dispensability, and rate of sequence evolution (Krylov et al., 2003; Cai and Petrov, 2010).

This study focuses for the first time on the evolutionary stability of orphan genes. We investigate the factors contributing to orphan loss and find that orphan age, male-biased gene expression, and microsatellite content are correlated with orphan stability. Surprisingly, differences in evolutionary rates cannot explain orphan loss and we propose that orphan loss is driven by lineage-specific evolutionary constraints. Overall, orphan genes are lost at a significantly higher rate than non-orphan genes, supporting the ‘rapid-turnover’ hypothesis.

Results

Orphans are commonly detected by BLASTing the genes of a given organism against a set of outgroup species (Domazet-Loso and Tautz, 2003; Toll-Riera et al., 2009). A BLASTP cutoff of 10−3–10−4 was found to be optimal to maximize sensitivity and specificity in Drosophila (Domazet-Loso and Tautz, 2003). To identify orphans we used a BLASTP cutoff of 10−4 combined with a TBLASTN cutoff of 10−4, to exclude genes with unannotated orthologs in other species. Following these criteria, we searched in Drosophila pseudoobscura for genes with no sequence conservation in 10 Drosophila species outside the Drosophila obscura group (Figure 6—figure supplement 1). In total, we identified 1152 orphans, corresponding to 7% of all the D. pseudoobscura genes. Our estimate is slightly lower than a previous one (Zhang et al., 2010), due to our different filtering procedure, but still consistent with a high rate of orphan gain in Drosophila (Domazet-Loso and Tautz, 2003; Domazet-Loso et al., 2007; Zhou et al., 2008; Wissler et al., 2013). Our data clearly indicate that orphan genes are subject to purifying selection, as they show several hallmarks of functional protein-coding sequences (Figures 1, 2). A comparison of orphan genes preserved between D. pseudoobscura and D. affinis resulted in a distribution of dN/dS significantly lower than 1 with a median of 0.44 (Figure 1—figure supplement 1, one-sided Wilcoxon signed-rank test, p<1.0 × 10−15), as expected for protein-coding sequences. Moreover, dN/dS for orphans is significantly lower (Mann–Whitney test, p=2.7 × 10−14) than dN/dS calculated on a random set of intergenic regions with the same length distribution of orphans (see ‘Materials and methods’, section ‘Evolutionary rates’) (Figure 1A). Consistent with this, we also found orphans to be more conserved than intergenic regions (Figure 1B, Figure 1—figure supplement 2). The codon usage bias of orphans is intermediate to that of old genes and intergenic regions (Figure 1C).

Figure 1 with 2 supplements see all
Orphans are subject to purifying selection.

(A) dN/dS of D. pseudoobscura and D. affinis orthologs. dN/dS is lowest for old genes, but also orphan genes have dN/dS smaller than one. A comparison of orphans and intergenic regions shows that dN/dS for orphans is significantly smaller (Mann–Whitney test, p=9.5 × 10−10), indicating purifying selection on orphan genes. Intergenic regions were of similar length and chromosomal position as the orphan genes. (B) Sequence similarity in HSPs obtained from BLASTing D. pseudoobscura genes against the D. affinis genome. Orphans are more conserved than intergenic regions (Mann–Whitney test, p=0.00238) and less conserved than old genes (Mann–Whitney test, p<1.0 × 10−15). (C) Codon usage was measured by the Codon Adaptation Index (Sharp and Li, 1987). The codon usage of orphans is significantly higher than that of intergenic regions (Mann–Whitney test, p<1.0 × 10−15) indicating that orphans are subject to purifying selection. In comparison to old genes, orphans have a significantly lower codon usage bias (Mann–Whitney test–p<1.0 × 10−15). Overall, all three analyses demonstrate that orphans are not annotation artifacts, but evolutionary conserved genes.

https://doi.org/10.7554/eLife.01311.003
pN/pS for old genes, orphans, and intergenic regions.

Orphans show a pN/pS intermediate between old genes and intergenic regions. Nevertheless, pN/pS is significantly smaller for orphans compared to intergenic regions (Mann–Whitney test, p<1.0 × 10−15), indicating coding purifying selection acting on orphans.

https://doi.org/10.7554/eLife.01311.006

To further test for purifying selection acting on orphans, we used a polymorphism dataset of 45 strains being re-sequenced for the third chromosome of D. pseudoobscura (‘Materials and methods’). We calculated the ratio of synonymous to non-synonymous polymorphism (pN/pS), since it provides an indication of purifying selection. We found that pN/pS for orphans is significantly lower compared to intergenic regions (Mann–Whitney test, p=0.02182) (Figure 2), and significantly greater for old genes (Mann–Whitney test, p<1.0 × 10−15), consistent with purifying selection operating on orphans.

In agreement with studies in other species (Domazet-Loso and Tautz, 2003; Toll-Riera et al., 2009; Wolf et al., 2009; Cai and Petrov, 2010; Capra et al., 2010; Carvunis et al., 2012), we also find that orphan genes are shorter (median length for orphans = 344 bp, median length for old genes = 1470 bp), have a lower GC content (median GC content for orphans = 0.54, median GC content for old genes = 0.55), are expressed at lower levels (expression in D. pseudoobscura males: mean expression for orphans = 29 FPKM, mean expression for old genes = 41 FPKM) than old genes (Figure 3). Using CD-Hit (Li and Godzik, 2006), we found the fraction of genes with a paralog (>90% protein similarity) to be similar for orphans (6.9%) and old genes (6.4%). Orphans are more enriched in microsatellites, also consistent with previous findings in vertebrates (Toll-Riera et al., 2012) and rice (Guo et al., 2007). Furthermore, unlike mammals (Toll-Riera et al., 2009), none of the D. pseudoobscura orphans was found to be associated with transposable elements (see ‘Transposons detection’).

Comparison of orphans and genes conserved among 10 Drosophila species outside of the obscura group.

Orphans differ from old genes in various features: (A) gene length (B) GC content, (C) dN/dS (D) percentage of microsatellites in coding sequence (E) Codon Adaptation Index (F) gene expression level in D. pseudoobscura males (G) gene expression level in D. pseudoobscura females (H) sex-biased expression. Orphans are shorter (Mann–Whitney test, p<1.0 × 10−15), have lower GC content (Mann–Whitney test, p=3.9 × 10−7), lower codon usage bias (Mann–Whitney test, p<1.0 × 10−15), lower expression (Mann–Whitney test, p<1.0 × 10−15), higher proportion of microsatellites (Mann–Whitney test, p=1.8 × 10−4) and higher dN/dS (Mann–Whitney test, p<1.0 × 10−15) compared to old genes. Moreover, orphans are more enriched in male-biased genes compared to old genes (χ2-test, p<1.0 × 10−15).

https://doi.org/10.7554/eLife.01311.007

The distribution of orphans is heterogeneous across chromosomes (χ2-test, p<1.0 × 10−15), with the X chromosome having the highest fraction of orphans. In the obscura group, the two X-chromosome arms have a different evolutionary history. XL corresponds to Muller’s element A and is homologous to the X chromosome in D. melanogaster. XR, however, has been recently derived from an autosome (Muller’s element D, 3L in D. melanogaster). Analyzing the old-X and neo-X chromosomes separately, we observed a striking difference in the number of orphans despite similar chromosome sizes, with the old-X responsible for the excess of X-linked orphan genes, and the neo-X showing a similar number of orphans as the autosomes (Figure 4). For each chromosomal arm, we computed genomic features in 100 kb windows to correlate them with the difference in orphan content between old-X and neo-X. We found that average GC content, microsatellite density, transposon density, and length of intergenic regions differ between the two chromosomal arms (Figure 5).

Chromosomal distribution of old genes and orphan genes.

Orphans are overrepresented on the old-X. The number of orphan genes on the neo-X (XR) is significantly lower than on the old-X (XL) (χ2-test, p<1.0 × 10−15).

https://doi.org/10.7554/eLife.01311.008
Comparison of genomic features among autosomes, old-X and neo-X.

(A) GC content in 100 kb windows, (B) Microsatellite density in 100 kb windows, (C) Transposon density in 100-kb windows, (D) Length of intergenic regions, (E) Recombination rate. GC content is significantly greater on the neo-X compared to old-X for 10 kb windows (Mann–Whitney test, p=0.00020), but not for 100 kb windows (Mann–Whitney test, p=0.1092). Microsatellite density is significantly higher on the neo-X for both windows of 10 kb (Mann–Whitney test, p=1.9 × 10−12) and 100 kb (Mann–Whitney test, p=0.00025). Transposon density is significantly lower on the neo-X for both windows of 10 kb (Mann–Whitney test, p<1.0 × 10−15) and 100 kb (Mann–Whitney test, p=4.6 × 10−12). Intergenic regions are significantly shorter on the neo-X compared to the old-X (Mann–Whitney test, p=7.4 × 10−9). Recombination rate does not differ significantly between old-X and neo-X (Mann–Whitney test, p=0.629).

https://doi.org/10.7554/eLife.01311.009

We hypothesized that this pronounced difference between the two chromosome arms might reflect a different history of X-linkage. If orphan genes emerge at a higher rate on the X-chromosome (Levine et al., 2006), the shorter history of X-linkage on the neo-X could explain the paucity of orphans on the neo-X compared to old-X. In this case, the difference in orphan number between old-X and neo-X chromosomes should date back to the time before the origin of the neo-X, with a similar number of orphans originating after the creation of the neo-X. We therefore used the genomic sequences of five members of the D. obscura group (D. pseudoobscura [Richards et al., 2005], D. miranda [Zhou and Bachtrog, 2012], and the de novo assembled D. persimilis, D. lowei, and D. affinis) to date the origin of the orphan genes to different ancestral nodes in the phylogenetic tree of these species (Beckenbach et al., 1993). We distinguished five groups of genes: old genes (non orphans) and four different orphan age classes (Figure 6). Surprisingly, we observed a consistent paucity of orphans on XR relative to XL across all age classes (Figure 7). This persistent difference in orphan number between XL and XR in all age classes suggests that X-linkage is not sufficient to explain the enrichment of orphans on XL. We conclude that the former autosome differs from the ancestral X chromosomal arm by a yet unidentified feature that affects the emergence of new orphans.

Figure 6 with 1 supplement see all
Orphan gain and losses in the Drosophila obscura group.

Schematic phylogenetic tree of the Drosophila obscura group species according to Beckenbach et al. (1993) with D. melanogaster as outgroup. Genes conserved between D. pseudoobscura and 10 non-obscura Drosophila species correspond to age class 5 (old genes). For each age class the number of gene gains is shown in black. Orphans lost at a given branch are indicated in red. Note that losses at internal branches cannot be calculated, since all the orphans are present in D. pseudoobscura. Losses in D. affinis cannot be unambiguously assigned due to the absence of an additional obscura outgroup.

https://doi.org/10.7554/eLife.01311.010
Chromosomal distribution of orphans of different age classes.

In each age class orphans are underrepresented on the neo-X (XR) compared to old-X (XL) (Age class 4: χ2-test, p=6.3 × 10−9; age class 3: χ2-test, p=4.4 × 10−5; age class 2: χ2-test, p=0.00590; age class 1: χ2-test, p=0.00876; D. pseudoobscura specific: χ2-test, p=0.00030).

https://doi.org/10.7554/eLife.01311.012

The analysis of orphans that have putatively lost their function via the acquisition of a stop codon or a frame shift causing insertion/deletion (pseudogenized/lost orphans) reveals another interesting feature of the XL–XR fusion. The oldest orphans in our dataset (age class 4) show a pronounced excess of pseudogenized orphans on XR in D. affinis and D. miranda (Figure 8A). This trend was not observed for orphans that emerged on XR after the XL–XR fusion (Figure 8B,C), nor for old genes (Figure 8D) and is not due to an increased rate of orphan gain on XR (Figure 9). Since the oldest orphans (age class 4) on XR are a mixture of autosomal (i.e., before the fusion) and sex-chromosomal (i.e., after the fusion) orphans, we speculate that the high rate of pseudogenization of orphans on the XR may reflect the new X-linkage of previously autosomal orphans. A previous study (Meisel et al., 2009) found that the XR chromosome has experienced a burst of gene duplications to autosomes after its creation. It is plausible that after the conversion of the XR from autosome to sex-chromosome, orphans might have been duplicated to autosomes, whereas the XR ancestral copy would have become pseudogenized. To test this hypothesis, we looked for evidence of gene duplications for the orphans lost on the XR at node 4 (Figure 6). We aligned the sequences of these genes in D. lowei and D. miranda to the respective genomes using BLASTN (cutoff 10−5). Upon manual inspection of the alignments, we found that only 1 out of 21 genes in D. miranda (gene ID: GA23486) and 1 out of 14 genes in D. lowei (gene ID: GA23807) had a second hit on an autosome covering at least 50% of the length of the query gene. Other genes either produced a single best hit on the XR chromosome or spurious short hits on other chromosomes (data not shown). Thus, we conclude that duplication of orphans cannot explain the excess of pseudogenized orphans on XR. Nevertheless, our analysis clearly indicates that the emergence of the neo-X chromosome influenced the orphan dynamics on XR, affecting rates of both gain and loss, thus we excluded this chromosome arm for our analyses of the rate of orphan turnover.

Orphans predating the XL-XR fusion are preferentially lost on the neo-X.

For three terminal branches (D. lowei, D. miranda, and D. persimilis) the fraction of lost genes for each age class is shown. Each autosome and both X-chromosome arms are shown in different color. At node 4, where the neo-X originated, we observed the highest rate of orphan pseudogenization on the neo-X (A). Notably, this effect is not seen for younger orphans (B and C) neither for old genes (D).

https://doi.org/10.7554/eLife.01311.013
No change in orphan gain on the neo-X chromosome.

The percentage of orphan genes on the neo-X chromosome remains constant through time (indicated by age classes).

https://doi.org/10.7554/eLife.01311.014

For each age class, we determined the number of pseudogenized orphans (Tautz and Domazet-Loso, 2011). In the D. persimilis lineage, orphan pseudogenization can be studied for three different age classes. If orphans of all age classes were functionally equivalent, no difference in the rate of orphan pseudogenization would be expected. We observe, however, that the fraction of orphan pseudogenes decreases with age (Figure 10). The D. miranda lineage also shows a higher loss of young orphan genes. The relatively small number of observations, however, precludes statistical testing of this trend. Overall, orphan genes are lost significantly more often than old genes (Fisher’s exact test, p=3.3 × 10−8), consistent with the rapid turnover hypothesis. The unequal conservation of orphans of different age classes is also apparent after normalizing by coding sequence length (Figure 11), to account for the fact that longer coding sequences (CDS) have a greater chance of acquiring ORF-disrupting mutations. When looking at the distribution of premature termination codons (PTC) along the open reading frame (ORF) of all genes, we observed that PTCs are enriched at the beginning and at the end of the ORF (Figure 12), consistent with previous results in D. melanogaster (Lee and Reinhardt, 2012) and D. pseudoobscura (Hoehn et al., 2012). Since ORF-disrupting mutations occuring at the end of the ORF might have little impact on gene function, we redefined pseudogenes by considering only ORF-disrupting mutations localized in the first half of the ORF and confirmed that orphans of age class 3 are lost more often than those of age class 4 (Figure 13). Age class 2 was intermediate, most likely not reflecting a biological phenomenon, but due to a high sampling variance associated with the small number of observations (9 orphans). Finally, the pattern is also robust to a more conservative criterion for ortholog assignment (see ‘Annotation of the obscura species’, Figure 14).

Young orphan genes are more likely to be lost.

The barplot shows the fraction of orphans that has acquired a frameshift or premature stop codon (i.e., lost function). For D. lowei, D. miranda, and D. persimilis, the fraction of lost orphans is shown for different age classes. Orphans are more likely to be lost than old genes. Both the D. miranda and D. persimilis lineage show that younger orphans are more likely to lose function than older ones.

https://doi.org/10.7554/eLife.01311.015
Young orphan genes are more likely to be lost: accounting for CDS length.

To test if the short CDS of orphans affects the pattern that young orphans are more likely to lose function, we normalized the percentage of losses by the median CDS length of genes at that node.

https://doi.org/10.7554/eLife.01311.016
Distribution of premature stop codons (PTCs) along the ORF for all genes containing PTCs.

PTCs are enriched at the beginning and at the end of the ORF in each species.

https://doi.org/10.7554/eLife.01311.017
Young orphan genes are more likely to be lost: considering only frameshifts and premature stop codons occurring in the first half of the ORF.

We repeated the analysis shown in Figure 10 by considering only frameshifts and premature stop codons occurring in the first half of the ORF to define a conservative set of pseudogenes, since disrupting mutations occurring at the end of the ORF are likely to have little impact on the gene function.

https://doi.org/10.7554/eLife.01311.018
Young orphan genes are more likely to be lost: the conservative set of orthologs.

We repeated the analysis shown in Figure 10 by restricting it to orthologs for which at least one flanking gene is identified in the same contig (see ‘Annotation of the obscura species’). Due to the substantially reduced number of orphans in the older age classes, we combined age class 3 and 4.

https://doi.org/10.7554/eLife.01311.019

To determine features associated with the differences in disabling mutations among orphans from different age classes, we contrasted orphans lost in D. lowei and/or D. persimilis (lost orphans) vs orphans conserved in all the obscura species (conserved orphans). Genes in both classes evolve at the same rate, are of similar length, and have similar codon usage bias (Figure 15A–E). Conserved orphans have a higher GC content, contain fewer microsatellites, are expressed at a higher level and are more male-biased (Figure 15B,D–F,G,H) compared to lost orphans. Conserved orphans tend to increase their expression level as they become older (Figure 16A), whereas the opposite pattern is true for lost orphans (Figure 16B).

Features of conserved orphans vs lost orphans measured in D. pseudoobscura.

(A) Gene length (B) GC content, (C) dN/dS (D) percentage of microsatellites in coding sequence (E) Codon Adaptation Index (F) gene expression levels in D. pseudoobscura males (G) gene expression levels in D. pseudoobscura females (H) sex-biased expression. Gene length (Mann–Whitney test, p=0.7235) and evolutionary rates (Mann–Whitney test, p=0.5835) are not significantly different between conserved and lost orphans. Lost orphans have higher GC content (Mann–Whitney test, p=0.00325), lower expression in D. pseudoobscura males (Mann–Whitney test, p=0.00012) and females (Mann–Whitney test, p=0.00230) and a higher microsatellite content (Mann–Whitney test, p=0.00049) compared to conserved orphans. Lost orphans are enriched in unbiased genes compared to conserved orphans (χ2-test, p=0.02611).

https://doi.org/10.7554/eLife.01311.020
Conserved and lost orphans differ in their gene expression pattern.

Expression intensity and sex bias in D. miranda for orphans conserved in all the obscura species (conserved orphans) vs orphans that pseudogenized in D. lowei and/or D. persimilis (lost orphans). Expression is calculated in males for orphans of age classes 3 and 4. Expression level increases with age for conserved orphans (A), while it decreases for lost orphans (B).

https://doi.org/10.7554/eLife.01311.021

Orphan genes are frequently expressed in the testis (Levine et al., 2006; Begun et al., 2007) and have a male-biased gene expression pattern (Metta and Schlötterer, 2008). This pattern could be generated by pervasive gene expression in testis, which facilitates the functional recruitment of non-specific expression (Kaessmann, 2010). Another explanation is that expression in testis does not require a complex architecture of regulatory modules (Sassone-Corsi, 2002; Kleene, 2005; Kaessmann, 2010), so that fewer substitutions are required to obtain a functional regulatory module for expressing a novel gene in testis compared to other tissues. We scrutinized these explanations by comparing the fraction of male-biased genes among orphan genes from different age classes. Unexpectedly, the fraction of male-biased genes increases with the age of the orphan genes (Figure 17). This increase of male-biased orphans among the older age classes is the result of a preferential loss of orphans with an unbiased gene expression (Figure 18). To confirm that male-biased gene expression is associated with orphan retention rather than emergence, we analyzed the sex-bias in D. miranda for orphans with and without an open reading frame. Consistent with the gene expression pattern in D. pseudoobscura, we found that lost orphans have a significantly lower male-bias in D. miranda (Figure 19). We conclude that the previously reported male-biased gene expression of orphan genes is not the result of a preferential recruitment of male-biased transcripts, nor do orphans gradually acquire male-biased gene expression. Rather, male-biased orphans are more likely to be retained.

The proportion of male-biased orphans increases with age.

Sex-biased expression was measured in D. pseudoobscura for orphans belonging to different age classes and for old genes (age class 5).

https://doi.org/10.7554/eLife.01311.022
Conservation of orphans is correlated with male-biased gene expression.

Orphans with male-biased gene expression in D. pseudoobscura were grouped into classes according to expression bias strength. The fraction of conserved orphans in each bin shows a significant positive correlation with expression bias (Spearman’s rho = 0.811, p=0.02692). This correlation suggests that orphans with a more pronounced male-biased expression tend to persist longer than less male-biased orphans. No similar trend was seen for female-biased orphans (Spearman’s rho = 0.78262, p=0.1176).

https://doi.org/10.7554/eLife.01311.023
Comparison of strength of sex-biased gene expression for conserved and lost orphans in D. miranda.

A sex-biased gene expression larger than zero indicates a higher gene expression intensity in males than in females (male-biased gene expression). Conserved orphans have significantly higher male-biased expression than lost orphans (Mann–Whitney test, p=0.03158).

https://doi.org/10.7554/eLife.01311.024

Discussion

Our study provides the missing link to understand orphan dynamics. Until now, orphan evolution was primarily studied on long phylogenetic branches. Although this approach is well suited to discover new orphans, it does not allow tracing the evolution of orphans. Previous studies showed a high rate of orphan gain, which is not reflected in an increase in gene number. To resolve this apparent paradox, it has been postulated that orphans must be lost at a high rate as well (Tautz and Domazet-Loso, 2011). In this study, we used the framework of closely related species in the obscura group to study the patterns of orphan gain and losses. We show that orphans not only emerge at high rates, but that they are also rapidly lost (Figure 10). Interestingly, most losses (∼76%) were due to disabling mutations rather than deletions of the orphan gene. Although under equilibrium conditions the number of losses balances the number of orphan gains, here, we observed a surplus of orphan gains (Figure 6). We caution that this discrepancy probably does not imply an increase of gene number, but rather reflects the limited evolutionary time to acquire mutations. Using a rather conservative criterion for disabling mutations, either premature stop codons or frameshift indels, we have probably not identified all orphans that have lost their function. Furthermore, we do not account for the possibility of loss of function due to changes in gene regulation.

Importantly, codon usage bias, dN/dS values and sequence conservation clearly suggest that orphan genes are functionally constrained and these constraints do not differ among orphans that are conserved in the obscura group and those that lost function in at least one species of the group. Hence, it may be possible that orphan loss is stochastic and reflects weak purifying selection. Nevertheless, lost orphans differ in some aspects from conserved ones. Orphans that are lost contain more microsatellite stretches and have a lower, less sex-biased gene expression than retained ones. Furthermore, we also found that the rate of orphan loss decreases with orphan age, a result consistent with orphans serving a functional role only temporarily. Previous work suggested that orphans are important for adaptation to novel environments (Khalturin et al., 2009; Colbourne et al., 2011), but it is also possible that orphans contribute to stabilize new connections in gene networks (Capra et al., 2010; Warnefors and Eyre-Walker, 2011) and become obsolete once such new connections have been optimized. Our data suggest that orphans become quickly functional, which is reflected in their codon usage bias, dN/dS ratio and sequence conservation.

The chromosomal translocation resulting in the neo-X chromosome provides another interesting perspective on the evolution of orphan genes. Despite the fact that the neo-X is now fully dosage compensated (Abraham and Lucchesi, 1974), and has obtained a similar base composition as the XL (Gallach et al., 2007), we noted that the translocation resulted in a preferential loss of orphan genes on the neo-X. Since this pattern is restricted to orphans that most likely originated before the chromosomal fusion, we argue that the change in chromosomal environment has affected the function of orphan genes, most likely via expression differences. We speculate that the selective advantage conferred by these orphans has diminished, which resulted in a higher loss rate. Interestingly, the elevated rate of orphan loss after the neo-X formation seems to be still ongoing. This differential loss of orphan genes point in a similar direction as the observation that the gene composition of the neo-X has been altered by gene duplication (Meisel et al., 2009). Hence, both (orphan) gene loss and duplication contribute to fast gene content remodeling on a newly formed sex chromosome in Drosophila.

Materials and methods

Species data collection

An individual species sample of D. affinis (stock number 140120141.02) was ordered from the Drosophila Species Stock Center (https://stockcenter.ucsd.edu/info/welcome.php) and sequenced on the Illumina GAIIx following the paired-end library preparation protocol (version Illumina 1.7) in two runs (run 1: read length = 101 bp, insert size = 230 bp; run 2: read length = 101 bp, insert size = 550 bp). Short genomic reads for D. lowei (accessions SRX091466 and SRX091467) and D. persimilis (accession SRX091471) were downloaded from the Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra). The genome of D. miranda was downloaded from NCBI (GenBank Assembly ID GCA_000269505.1). The genome of D. pseudoobscura was downloaded from FlyBase (release 2.23).

Assembly of the obscura species

Reads for D. affinis, D. lowei, and D. persimilis were trimmed using the Perl script trim_fastq.pl (parameters –quality-threshold 20 −−min-length 40) from PoPoolation (Kofler et al., 2011). For each species, a de novo assembly (parameters: min-contig-length 200) was performed using CLC Genomics Workbench 4.6 (http://www.clcbio.com/products/clc-genomics-workbench/), followed by scaffolding with nucmer (parameters: –c 30 –g 1000 –b 1000 –l 15) against the D. pseudoobscura genome. Average coverage per assembled genome was calculated by realigning the reads against the contigs of the respective species with Bowtie 2.1.0 (parameters: --very-fast) and selecting only reads with mapping quality >20.

Annotation of the obscura species

The annotation of D. affinis, D. lowei, D. miranda, and D. persimilis is based on orthology to D. pseudoobscura using Exonerate 2.2.0 (parameters: -model protein2genome–bestn 1 -showtargetgff), by aligning the longest isoform of D. pseudoobscura proteins extracted from a recent re-annotation of D. pseudoobscura (Palmieri et al., 2012) to the genomes of D. affinis, D. lowei, D. miranda, and D. persimilis. For each gene, the best unambiguous hit was retained. To remove non-informative hits, we also required a minimum fraction of the gene to be recovered. Since the sequence conservation of orthologs decreases with divergence time, the expected length of the ortholog depends strongly on the phylogenetic distance between query and subject sequence. To apply consistent criteria for all species, we empirically determined the expected fraction of a gene with sequence homology. Based on genes that are conserved between D. pseudoobscura and the 10 Drosophila species outside the obscura clade (old genes) (Figure 6—figure supplement 1), we determined the distribution of the fraction of the genes that could be aligned. As cutoff the value we used the 5th percentile of the distribution of aligned protein length of old genes. This resulted in a threshold of 47% for D. affinis, 52% for D. lowei, 59% for D. miranda and 53% for D. persimilis. Hence, only orphan orthologs that showed a fraction of aligned coding sequence higher than the empirically determined cutoffs were retained. In addition to this ortholog set, we generated an alternative, more conservative ortholog set. For this one, at least one of the flanking genes of D. pseudoobscura was required to be in synteny with the respective orthologs in D. affinis, D. miranda and D. persimilis. D. lowei was not considered in the synteny analysis since most of the genes in this species are flanked by genomic gaps, due to the shorter contig length of the D. lowei assembly (Table 1), which caused many contigs to contain only a single gene (Table 2), thus precluding proper synteny assignments. Assembly and annotation of all the species are available at http://popoolation.at/affinis_genome, http://popoolation.at/lowei_genome, http://popoolation.at/miranda_genome and http://popoolation.at/persimilis_genome. Detailed annotation statistics for each gene are available at http://dx.doi.org/10.5061/dryad.hq564 (Palmieri et al., 2014).

Table 1

De novo assembly statistics

https://doi.org/10.7554/eLife.01311.025
D. affinisD. loweiD. persimilis
Number of contigs28,946106,46517,387
N759,4781,21810,359
N5025,1603,23024,172
N2549,0627,35749,047
Minimum length121162147
Maximum length216,90387,164204,742
Average length5,1831,3887,736
Total bp150,030,247147,756,871134,501,523
Average coverage51 X92 X44 X
  1. The D. miranda genome was available at NCBI, thus no de novo assembly was made for this species.

Table 2

Orthology annotation statistics

https://doi.org/10.7554/eLife.01311.026
D. affinisD. loweiD. mirandaD. persimilis
Total genes14,28714,95215,28214,995
Genes with frameshifts/PTC*1,2331,2661,171898
Mean number of genes per contig3.41.63.4
Median number of genes per contig212
Maximum number of genes per contig352437
  1. *

    PTC = Premature termination codons.

Detection of orphan genes

D. pseudoobscura proteins corresponding to the longest isoform for each gene were aligned using BLASTP (E < 10−4) and TBLASTN (E < 10−4) against the published proteomes and genomes of 10 Drosophila species outside the obscura group (D. melanogaster, D. simulans, D. sechellia, D. erecta, D. yakuba, D. ananassae, D. willistoni, D. mojavensis, D. virilis, D. grimshawi). Genes without BLAST hits and without annotated orthologs in FlyBase (gene orthologs release 09-2011) were classified as orphans.

Polymorphism analysis

Illumina reads for 45 D. pseudoobscura strains were downloaded from NCBI (Sequence Read Archive, accession SRP017196). Reads were trimmed using PoPoolation (Kofler et al., 2011) and a total of 3.5 million reads was randomly extracted for each strain and combined into a single FASTQ file. The combined reads were treated as a Pool-Seq dataset and mapped to the FlyBase D. pseudoobscura genome release 2.23 with BWA (Li and Durbin, 2009) (parameters -o 1 -n 0.01 -l 200 -e 12 -d 12) on a hadoop cluster using DistMap (Pandey and Schlötterer, 2013). From the resulting BAM file, PCR duplicates were removed with Picard (http://picard.sourceforge.net) using the tool MarkDuplicates.jar (parameters REMOVE_DUPLICATES = true, VALIDATION_STRINGENCY = SILENT). Proper-pairs with mapping quality >20 were extracted with samtools (version 0.1.18) (Li et al., 2009). Indels were detected with PoPoolation using the script identify-genomic-indel-regions.pl (parameters--min-count 2 --indel-window 5) and masked from the reference genome prior to SNP calling. Coverage was subsampled to 50X for all the chromosomes. Only SNPs on the 3rd chromosome were considered in all analyses, since a balancer chromosome was used to extract the 3rd chromosome, precluding an unbiased polymorphism analysis for the remaining chromosomes. SNPs were called with the PoPoolation script Variance-sliding.pl (parameters--min-coverage 10 --min-count 2 --max-coverage 500 --min-qual 20 --window-size 500 --step-size 500 --fastq-type sanger--pool-size 45).

Calculation of orphan gains and losses

Orphan gains and losses (pseudogenizations) were inferred by Dollo parsimony. Based on the phylogenetic tree of Figure 6, a gene was assigned as gained at a given node if an intact ortholog was present in both external branches of the subtree corresponding to that node. For example, a gene having an intact ORF in D. lowei but not in D. affinis was classified as gained at node 3 (Figure 6). A gene was considered to be lost at a terminal branch if at least one ORF-disrupting mutation (frameshift/premature stop codon) was present in the gene at that branch and two intact ORFs were detected at both external leaves (Wang et al., 2006). The relatively high coverage of our assemblies (Table 1) makes unlikely that disrupting mutations are sequencing errors. In D. affinis for instance, only 8 genes had an average coverage lower than 20x.

A gene was considered as completely deleted in a species if no ortholog was detected in that species and no BLASTP (E < 10−4) or TBLASTN (E < 10−4) hit was found. Deletions were not considered into analyses of gene turnover, since they cannot be distinguished from missing annotations.

Expression analysis

Four RNA-Seq datasets of D. pseudoobscura males and females (strains ps94 and ps88 from the ArrayExpress database—accession E-MTAB-1424), together with two RNA-Seq samples of D. miranda males and females from the Sequence Read Archive (accessions SRX106024, SRX106025), were used for expression analysis. For each sample, reads were trimmed using PoPoolation (Kofler et al., 2011) and aligned to the genome of the respective species with GSNAP version 2012-07-12 (Wu and Nacu, 2010) (parameters: –N 1). Only proper pairs mapping unambiguously to one position were retained. Expression in FPKM was calculated with Cufflinks version 1.2.1 (parameters: -F 0.10 –j 0.15 –I 300000). For D. pseudoobscura sex-bias was calculated using the package DESeq (Anders and Huber, 2010), treating the strains as two biological replicates for each sex and applying an FDR = 0.1. Differential expression between D. miranda males and females was calculated for both species using the log2 fold change on the normalized expression counts using the normalization protocol implemented in the R package DESeq (Anders and Huber, 2010) version 1.10.1.

Codon usage bias

Codon usage bias was calculated using the R package seqinr (function cai) based on the D. pseudoobscura codon usage table downloaded from http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=7237.

Evolutionary rates

Coding sequences of D. pseudoobscura and D. miranda orthologs without frameshifts/stop codons were aligned using PRANK (Loytynoja and Goldman, 2005) (parameters: –codon). To test for purifying selection on orphans, dN/dS was compared between orphans and a set of randomly selected intergenic regions. This set was generated as follows: (1) we identified the intergenic regions from the D. pseudoobscura annotation from Palmieri et al. (2012), (2) for each CDS belonging to an orphan gene we extracted all the intergenic regions longer than that CDS, (3) we randomly selected one intergenic region and we extracted from that a random subregion with the same length of a given orphan CDS, (4) this procedure was repeated for all orphan CDS, resulting in a set of intergenic regions with the same length distribution as orphan CDSs. These regions were aligned with BLASTN (cutoff 10−5) to the D. affinis genome and for each region the best hit was kept and realigned with PRANK (default parameters) to the D. pseudoobscura query sequence. Each alignment was truncated at the 5’end to get an alignment length, which is a multiple of 3. Internal stop codons were replaced by Ns. The ratio of the rates of nonsynonymous and synonymous substitutions per gene (dN/dS) was measured using Markov models of codon evolution and maximum likelihood methods implemented in PAML (Yang, 2007).

Comparison of genomic features among old-X, neo-X, and autosomes

To shed light on the differences in orphan number between XL and XR, different features were compared among old-X, neo-X, and autosomes in D. pseudoobscura (unassembled contigs were not considered in this analysis): (A) GC content was calculated with the R package seqinr for 100 kb sliding windows along each chromosome (B) microsatellite density was calculated using SciRoKo 3.4 (Kofler et al., 2007) (parameters: -mode mmfp–l 15 –r 3 –s 15 –p 5 –seedl 8 –seedr 3 –mmao 3) for 100 kb sliding windows along each chromosome; (C) transposon density was estimated with RepeatMasker 3.2.9 (parameters: –q–gff -nolow–norna–species drosophila) for 100 kb sliding windows along each chromosome; (D) length of intergenic regions were calculated using BEDTools (-complement) by interval subtraction between genome and gene coordinates; (E) recombination rates for different windows were taken from McGaugh et al. (2012).

Microsatellite detection

Microsatellites were detected on the transcript sequences of the longest isoform for each D. pseudoobscura gene using the tool SciRoKo 3.4 (Kofler et al., 2007) (parameters: -mode mmfp–l 15 –r 3 –s 15 –p 5 –seedl 8 –seedr 3 –mmao 3).

Transposons detection

Genomic annotation of transposons was performed in D. pseudoobscura using RepeatMasker 3.2.9 (parameters: –q–gff -nolow–norna–species drosophila). Only transposons longer than 50 bp and not overlapping with microsatellites (see ‘Microsatellite detection’) were retained. We required for an orphan to contain a full transposon sequence in one of its exons in order to classify it as associated with a transposon.

References

  1. 1
    Dosage compensation of genes on the left and right arms of the X chromosome of Drosophila pseudoobscura and Drosophila willistoni
    1. I Abraham
    2. JC Lucchesi
    (1974)
    Genetics 78:1119–1126.
  2. 2
  3. 3
    Relationships in the Drosophila obscura species group, inferred from mitochondrial cytochrome oxidase II sequences
    1. AT Beckenbach
    2. YW Wei
    3. H Liu
    (1993)
    Molecular Biology and Evolution 10:619–634.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
    Evolution of genes and genomes on the Drosophila phylogeny
    1. AG Clark
    2. MB Eisen
    3. DR Smith
    4. CM Bergman
    5. B Oliver
    6. TA Markow
    7. TC Kaufman
    8. M Kellis
    9. W Gelbart
    10. VN Iyer
    11. DA Pollard
    12. TB Sackton
    13. AM Larracuente
    14. ND Singh
    15. JP Abad
    16. DN Abt
    17. B Adryan
    18. M Aguade
    19. H Akashi
    20. WW Anderson
    21. CF Aquadro
    22. DH Ardell
    23. R Arguello
    24. CG Artieri
    25. DA Barbash
    26. D Barker
    27. P Barsanti
    28. P Batterham
    29. S Batzoglou
    30. D Begun
    31. A Bhutkar
    32. E Blanco
    33. SA Bosak
    34. RK Bradley
    35. AD Brand
    36. MR Brent
    37. AN Brooks
    38. RH Brown
    39. RK Butlin
    40. C Caggese
    41. BR Calvi
    42. A Bernardo de Carvalho
    43. A Caspi
    44. S Castrezana
    45. SE Celniker
    46. JL Chang
    47. C Chapple
    48. S Chatterji
    49. A Chinwalla
    50. A Civetta
    51. SW Clifton
    52. JM Comeron
    53. JC Costello
    54. JA Coyne
    55. J Daub
    56. RG David
    57. AL Delcher
    58. K Delehaunty
    59. CB Do
    60. H Ebling
    61. K Edwards
    62. T Eickbush
    63. JD Evans
    64. A Filipski
    65. S Findeiss
    66. E Freyhult
    67. L Fulton
    68. R Fulton
    69. AC Garcia
    70. A Gardiner
    71. DA Garfield
    72. BE Garvin
    73. G Gibson
    74. D Gilbert
    75. S Gnerre
    76. J Godfrey
    77. R Good
    78. V Gotea
    79. B Gravely
    80. AJ Greenberg
    81. S Griffiths-Jones
    82. S Gross
    83. R Guigo
    84. EA Gustafson
    85. W Haerty
    86. MW Hahn
    87. DL Halligan
    88. AL Halpern
    89. GM Halter
    90. MV Han
    91. A Heger
    92. L Hillier
    93. AS Hinrichs
    94. I Holmes
    95. RA Hoskins
    96. MJ Hubisz
    97. D Hultmark
    98. MA Huntley
    99. DB Jaffe
    100. S Jagadeeshan
    101. WR Jeck
    102. J Johnson
    103. CD Jones
    104. WC Jordan
    105. GH Karpen
    106. E Kataoka
    107. PD Keightley
    108. P Kheradpour
    109. EF Kirkness
    110. LB Koerich
    111. K Kristiansen
    112. D Kudrna
    113. RJ Kulathinal
    114. S Kumar
    115. R Kwok
    116. E Lander
    117. CH Langley
    118. R Lapoint
    119. BP Lazzaro
    120. SJ Lee
    121. L Levesque
    122. R Li
    123. CF Lin
    124. MF Lin
    125. K Lindblad-Toh
    126. A Llopart
    127. M Long
    128. L Low
    129. E Lozovsky
    130. J Lu
    131. M Luo
    132. CA Machado
    133. W Makalowski
    134. M Marzo
    135. M Matsuda
    136. L Matzkin
    137. B McAllister
    138. CS McBride
    139. B McKernan
    140. K McKernan
    141. M Mendez-Lago
    142. P Minx
    143. MU Mollenhauer
    144. K Montooth
    145. SM Mount
    146. X Mu
    147. E Myers
    148. B Negre
    149. S Newfeld
    150. R Nielsen
    151. MA Noor
    152. P O’Grady
    153. L Pachter
    154. M Papaceit
    155. MJ Parisi
    156. M Parisi
    157. L Parts
    158. JS Pedersen
    159. G Pesole
    160. AM Phillippy
    161. CP Ponting
    162. M Pop
    163. D Porcelli
    164. JR Powell
    165. S Prohaska
    166. K Pruitt
    167. M Puig
    168. H Quesneville
    169. KR Ram
    170. D Rand
    171. MD Rasmussen
    172. LK Reed
    173. R Reenan
    174. A Reily
    175. KA Remington
    176. TT Rieger
    177. MG Ritchie
    178. C Robin
    179. YH Rogers
    180. C Rohde
    181. J Rozas
    182. MJ Rubenfield
    183. A Ruiz
    184. S Russo
    185. SL Salzberg
    186. A Sanchez-Gracia
    187. DJ Saranga
    188. H Sato
    189. SW Schaeffer
    190. MC Schatz
    191. T Schlenke
    192. R Schwartz
    193. C Segarra
    194. RS Singh
    195. L Sirot
    196. M Sirota
    197. NB Sisneros
    198. CD Smith
    199. TF Smith
    200. J Spieth
    201. DE Stage
    202. A Stark
    203. W Stephan
    204. RL Strausberg
    205. S Strempel
    206. D Sturgill
    207. G Sutton
    208. GG Sutton
    209. W Tao
    210. S Teichmann
    211. YN Tobari
    212. Y Tomimura
    213. JM Tsolas
    214. VL Valente
    215. E Venter
    216. JC Venter
    217. S Vicario
    218. FG Vieira
    219. AJ Vilella
    220. A Villasante
    221. B Walenz
    222. J Wang
    223. M Wasserman
    224. T Watts
    225. D Wilson
    226. RK Wilson
    227. RA Wing
    228. MF Wolfner
    229. A Wong
    230. GK Wong
    231. CI Wu
    232. G Wu
    233. D Yamamoto
    234. HP Yang
    235. SP Yang
    236. JA Yorke
    237. K Yoshida
    238. E Zdobnov
    239. P Zhang
    240. Y Zhang
    241. AV Zimin
    242. J Baldwin
    243. A Abdouelleil
    244. J Abdulkadir
    245. A Abebe
    246. B Abera
    247. J Abreu
    248. SC Acer
    249. L Aftuck
    250. A Alexander
    251. P An
    252. E Anderson
    253. S Anderson
    254. H Arachi
    255. M Azer
    256. P Bachantsang
    257. A Barry
    258. T Bayul
    259. A Berlin
    260. D Bessette
    261. T Bloom
    262. J Blye
    263. L Boguslavskiy
    264. C Bonnet
    265. B Boukhgalter
    266. I Bourzgui
    267. A Brown
    268. P Cahill
    269. S Channer
    270. Y Cheshatsang
    271. L Chuda
    272. M Citroen
    273. A Collymore
    274. P Cooke
    275. M Costello
    276. K D’Aco
    277. R Daza
    278. G De Haan
    279. S DeGray
    280. C DeMaso
    281. N Dhargay
    282. K Dooley
    283. E Dooley
    284. M Doricent
    285. P Dorje
    286. K Dorjee
    287. A Dupes
    288. R Elong
    289. J Falk
    290. A Farina
    291. S Faro
    292. D Ferguson
    293. S Fisher
    294. CD Foley
    295. A Franke
    296. D Friedrich
    297. L Gadbois
    298. G Gearin
    299. CR Gearin
    300. G Giannoukos
    301. T Goode
    302. J Graham
    303. E Grandbois
    304. S Grewal
    305. K Gyaltsen
    306. N Hafez
    307. B Hagos
    308. J Hall
    309. C Henson
    310. A Hollinger
    311. T Honan
    312. MD Huard
    313. L Hughes
    314. B Hurhula
    315. ME Husby
    316. A Kamat
    317. B Kanga
    318. S Kashin
    319. D Khazanovich
    320. P Kisner
    321. K Lance
    322. M Lara
    323. W Lee
    324. N Lennon
    325. F Letendre
    326. R LeVine
    327. A Lipovsky
    328. X Liu
    329. J Liu
    330. S Liu
    331. T Lokyitsang
    332. Y Lokyitsang
    333. R Lubonja
    334. A Lui
    335. P MacDonald
    336. V Magnisalis
    337. K Maru
    338. C Matthews
    339. W McCusker
    340. S McDonough
    341. T Mehta
    342. J Meldrim
    343. L Meneus
    344. O Mihai
    345. A Mihalev
    346. T Mihova
    347. R Mittelman
    348. V Mlenga
    349. A Montmayeur
    350. L Mulrain
    351. A Navidi
    352. J Naylor
    353. T Negash
    354. T Nguyen
    355. N Nguyen
    356. R Nicol
    357. C Norbu
    358. N Norbu
    359. N Novod
    360. B O’Neill
    361. S Osman
    362. E Markiewicz
    363. OL Oyono
    364. C Patti
    365. P Phunkhang
    366. F Pierre
    367. M Priest
    368. S Raghuraman
    369. F Rege
    370. R Reyes
    371. C Rise
    372. P Rogov
    373. K Ross
    374. E Ryan
    375. S Settipalli
    376. T Shea
    377. N Sherpa
    378. L Shi
    379. D Shih
    380. T Sparrow
    381. J Spaulding
    382. J Stalker
    383. N Stange-Thomann
    384. S Stavropoulos
    385. C Stone
    386. C Strader
    387. S Tesfaye
    388. T Thomson
    389. Y Thoulutsang
    390. D Thoulutsang
    391. K Topham
    392. I Topping
    393. T Tsamla
    394. H Vassiliev
    395. A Vo
    396. T Wangchuk
    397. T Wangdi
    398. M Weiand
    399. J Wilkinson
    400. A Wilson
    401. S Yadav
    402. G Young
    403. Q Yu
    404. L Zembek
    405. D Zhong
    406. A Zimmer
    407. Z Zwirko
    408. DB Jaffe
    409. P Alvarez
    410. W Brockman
    411. J Butler
    412. C Chin
    413. S Gnerre
    414. M Grabherr
    415. M Kleber
    416. E Mauceli
    417. I MacCallum
    (2007)
    Nature 450:203–218.
    https://doi.org/10.1038/nature06341
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
    An algorithm for progressive multiple alignment of sequences with insertions
    1. A Loytynoja
    2. N Goldman
    (2005)
    Proceedings of the National Academy of Sciences of the United States of America 102:10557–10562.
    https://doi.org/10.1073/pnas.0409137102
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
    Data from: the life cycle of Drosophila orphan genes
    1. N Palmieri
    2. C Kosiol
    3. C Schlötterer
    (2014)
    Dryad Digital Repository, 10.5061/dryad.hq564.
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63

Decision letter

  1. Diethard Tautz
    Reviewing Editor; Max Planck Institute for Evolutionary Biology, Germany

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for sending your work entitled “The life cycle of Drosophila orphan genes” for consideration at eLife. Your article has been favorably evaluated by a Senior editor and 3 reviewers, one of whom is a member of our Board of Reviewing Editors.

The Reviewing editor and the other reviewers discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments to help you prepare a revised submission.

The authors discuss an interesting topic on life cycle of Drosophila orphan genes. They investigated the features of orphan genes in the D. pseudoobscura group based on an extensive new genomic dataset and by further analyses on substitution rate, gain/loss, and expression pattern of them they found that the orphan genes experienced purifying selection, and highly expressed orphan genes with strong male-bias are more likely to be retained. They proposed that orphan gene loss reflects lineage specific functional requirements.

The dataset consists of one newly sequenced genome as well as new assemblies of other genomes from read archives. The analysis is mostly solid and convincing. Overall, the paper represents a significant step forward in the field, but it requires further clarification of technical issues before it can be considered for acceptance in eLife.

The following comments are condensed from the comments of the three referees.

Referee 1:

My major comment concerning the methodology concerns the procedure to classify a “loss”. They say in the Methods section: “ Based on a set of orthologs in the obscura group, pseudogenes were identified by the presence of an intact ORF”, but Figure 9–figure supplement 2 shows that a large number of such disruptions occur in the last 10% of the reading frame (this is actually an interesting result that should not be relegated to the supplements). These might not be function-disabling mutations. I suggest to check relative conservation levels along the ORFs to assess this question or treat these genes differently in the overall comparisons.

My major comment concerning data analysis is towards the exclusive use of D. miranda for calculating divergence levels in Figure 1. Why is this not done for all species in the tree? And why D. miranda, rather than the newly sequenced D. affinis that is at the base of the comparisons?

All figures that support relevant data discussed in the text should be in the paper – supplementary figures should only be used for further clarifications.

Referee 2:

1) Loss of orphan genes is a key point of this study. To make the results convincing, calculation on gene gains and loses should be described in detail.

2) How exactly was the “intergenic” dataset constructed?

Referee 3:

1) When the average or median values are compared between gene classes, this comparison hides the potential fact that the distribution might include genes far from the mean. For example, even if Ka/Ks is on average <1 some orphan genes might be fast evolving and actually not easy to recognize in other species. These orphans are not “real” orphans but evolving to quick for the threshold of similarity to catch them. Authors should plot and compare the distributions not only means or medians.

In this sense, I do not think the authors have a conservative definition of orphans because they do not ask for Ka/Ks to be significantly smaller than 1. This would be conservative. Some might evolve as non-coding regions. Their definition is annotated genes of D. pseudoobscura that are expressed but show not similarity in other species at a particular threshold. However, whenever possible functionality should be supported with a Ka/Ks significantly <1 or Ka/Ks significantly >1 for the particular orphans to provide additional evidence of functionality. They should add a significance value in the supplementary table and comment on it in the text.

Is there any polymorphism data available for D. pseudoobscura to perform M-K test to test the functionality of the 225 new orphan genes in that genome?

2) Figure 5 should contain divergence times in the nodes. Figure 5 should be discussed in the light of the quality of every genome considered to make sure that losses or gains depend on the specific rates of orphan gene evolution per lineage and not on the quality of sequencing and assembly.

3) In the Results section I would like the authors to describe whether the loss of XR orphans is pseudogenization or actually a relocation. Can they comment if the gene has moved whenever possible to perform these analyses? They comment about the Meisel et al. work (2009). In that work some male-specific duplicates were found to relocate to autosomes. Is this the case for these orphans.

4) How were GAPs in the assembly treated? Can it be that some losses are actually GAPs in the assembly?

5) Can the high rate of orphans on the X chromosome be explained by the X being sequenced at a lower coverage if males and females were sequenced and this affecting the presence of genes in different assemblies? Please mention what individuals (males, females or 50%/50%) were sequenced for every genome.

6) In the Materials and methods section, an orphan is considered to be disabled if a single disrupting mutations is present. What is the probability that this single disabling mutation is actually a sequencing error? Please comment.

7) Intron turnover should be analyzed in more detail. How are those introns gained or lost? Are they part of the initially annotated coding region? This could reveal disablements. Some annotation software will put an intron in regions with disablements and still annotate a coding region.

8) In the Results section the authors mention that their definition of orphan is more conservative than in previous work. Please explain in detail in what sense the definition is conservative in this work.

9) Are orphans related (similar in sequence) to other orphans?

https://doi.org/10.7554/eLife.01311.027

Author response

Referee 1:

My major comment concerning the methodology concerns the procedure to classify a “loss”. They say in the Methods section: “ Based on a set of orthologs in the obscura group, pseudogenes were identified by the presence of an intact ORF”, but Figure 9–figure supplement 2 shows that a large number of such disruptions occur in the last 10% of the reading frame (this is actually an interesting result that should not be relegated to the supplements). These might not be function-disabling mutations. I suggest to check relative conservation levels along the ORFs to assess this question or treat these genes differently in the overall comparisons.

We agree that disruptions present at the end of the ORFs might not be considered as real pseudogenizations events. To account for this, we looked at the correlation between gene age and loss by considering only the disrupting mutations occurring in the first half of the ORF. We obtained qualitatively similar results, as shown in Figure 13.

My major comment concerning data analysis is towards the exclusive use of D. miranda for calculating divergence levels in Figure 1. Why is this not done for all species in the tree? And why D. miranda, rather than the newly sequenced D. affinis that is at the base of the comparisons?

We calculated now the divergence levels for all the species in the tree and obtained qualitatively similar results (Figure 1–figure supplement 2).

All figures that support relevant data discussed in the text should be in the paper - supplementary figures should only be used for further clarifications.

We moved all the figure supplements from Figure 9 (now Figure 10) to the main figures.

Referee 2:

1) Loss of orphan genes is a key point of this study. To make the results convincing, calculation on gene gains and loses should be described in detail.

We rewrote the methods section “Calculation of orphan gain and losses” providing a more detailed description.

2) How exactly was the “intergenic” dataset constructed?

We added the details in the method section “Evolutionary rates”.

Referee 3:

1) When the average or median values are compared between gene classes, this comparison hides the potential fact that the distribution might include genes far from the mean. For example, even if Ka/Ks is on average <1 some orphan genes might be fast evolving and actually not easy to recognize in other species. These orphans are not “real” orphans but evolving to quick for the threshold of similarity to catch them. Authors should plot and compare the distributions not only means or medians.

We plotted the distribution of dN/dS for orphans in Figure 1–figure supplement 1 and showed that most of the orphans have a dN/dS lower than 1.

In this sense, I do not think the authors have a conservative definition of orphans because they do not ask for Ka/Ks to be significantly smaller than 1. This would be conservative. Some might evolve as non-coding regions. Their definition is annotated genes of D. pseudoobscura that are expressed but show not similarity in other species at a particular threshold. However, whenever possible functionality should be supported with a Ka/Ks significantly <1 or Ka/Ks significantly >1 for the particular orphans to provide additional evidence of functionality. They should add a significance value in the supplementary table and comment on it in the text.

The reason why we called our orphan identification procedure conservative is that we included an additional filtering step that has not been used in previous publications. The procedure suggested by the reviewer is certainly a truly conservative one, but given the short evolutionary time scale of our study, too few mutations have occurred to identify genes with a statistically significant difference of Ka/Ks from one. We tested all orphans and none was found to have a Ka/Ks significantly different from 1. Hence, we did not further pursue the suggested strategy.

Since the reviewer feels that the wording “conservative” is a problem, we have removed it from the manuscript since we do not think that it is of major importance.

Is there any polymorphism data available for D. pseudoobscura to perform M-K test to test the functionality of the 225 new orphan genes in that genome?

We agree that the inclusion of polymorphism data could be interesting and we analyzed a set of polymorphism dataset for D. pseudoobscura from 45 individuals. Since the MK test is designed to identify genes with different patterns of evolution for divergence and polymorphism data, it is not well-suited for the identification of purifying selection. Hence, we used the polymorphism data set to calculate pN/pS for orphans (Figure 2) and compared them to conserved genes and intergenic regions. In support of our other data, also the polymorphism data supports the hypothesis of purifying selection acting on orphan genes.

2) Figure 5 should contain divergence times in the nodes. Figure 5 should be discussed in the light of the quality of every genome considered to make sure that losses or gains depend on the specific rates of orphan gene evolution per lineage and not on the quality of sequencing and assembly.

We added a new figure supplement (Figure 6–figure supplement 1), which includes the other Drosophila species and this figure includes scale with divergence times.

We understand the concern of the reviewer that the rates of gains and losses may also influenced by assembly quality, but we feel that this is difficult to put into this figure. Rather, we would like to refer the reviewer to Figure 10, which shows the percentage of lost non-orphan genes (age class 5). We caution, however, that it is not clear to what extent premature stops or frameshift mutations are segregating in these species or are the outcome of sequencing errors. We think that the high coverage of the genomes (Table 1) will result only in a small fraction of sequencing errors.

3) In the Results section I would like the authors to describe whether the loss of XR orphans is pseudogenization or actually a relocation. Can they comment if the gene has moved whenever possible to perform these analyses? They comment about the Meisel et al. work (2009). In that work some male-specific duplicates were found to relocate to autosomes. Is this the case for these orphans.

We performed this analysis by extracting the orphans pseudogenized on the XR in D. miranda and D. lowei and by aligning them to the respective genomes to find evidence for gene duplication. We found no support for this hypothesis, as described at in the Results section.

4) How were GAPs in the assembly treated? Can it be that some losses are actually GAPs in the assembly?

We understand the issue that GAPs in the assembly might confound the detection of genuine gene losses if the definition of loss is merely based on the presence/absence of the gene sequence in a given species. To avoid this problem we considered only pseudogenization events (i.e., premature stops or frameshift mutations), as specified in the methods. Thus GAPs in the assembly do not affect our definition of gene loss. At the same time, we also recognize that our definition of gene loss is an underestimate of the true number of losses.

5) Can the high rate of orphans on the X chromosome be explained by the X being sequenced at a lower coverage if males and females were sequenced and this affecting the presence of genes in different assemblies? Please mention what individuals (males, females or 50%/50%) were sequenced for every genome.

All the genome assemblies used in our study come from female individuals. This ensures similar coverage between autosomes and X chromosome. Thus, the high rate of orphans on the X cannot be explained by the X being sequenced at lower coverage. The strain details / accession numbers are given in the Methods. Nevertheless, we would point out that the coverage of the genome assemblies was high enough (Table 1), such that even a 50% lower coverage of males would not have resulted in a substantial increase in sequencing errors due to low coverage.

6) In the Materials and methods section, an orphan is considered to be disabled if a single disrupting mutations is present. What is the probability that this single disabling mutation is actually a sequencing error? Please comment.

Overall, our assemblies were not based on low coverage genomes (Table 1), thus sequencing errors are not very likely to contribute significantly to disrupting mutations, since this is expected for a coverage of 3 reads or less. However, we also think that a fraction of the disrupting mutations may still be segregating. We did a more careful analysis of disrupting mutations in D. mauritiana (a species for which we had assembled the genome ourselves and had also high quality polymorphism data available) and found consistent with other results in the literature that several of the interrupting mutations were not fixed in the population.

7) Intron turnover should be analyzed in more detail. How are those introns gained or lost? Are they part of the initially annotated coding region? This could reveal disablements. Some annotation software will put an intron in regions with disablements and still annotate a coding region.

The reviewer brings up a very important aspect. This made us recognize that intron turnover was not explored in sufficient detail in our study. Since we feel that this is not the main focus of our work and it would require in depth analyses, we removed the section about intron turnover.

8) In the Results section the authors mention that their definition of orphan is more conservative than in previous work. Please explain in detail in what sense the definition is conservative in this work.

Our definition of orphan is more conservative, since we used an additional filtering using TBLASTN. We explained how orphans were detected in previous studies at the beginning of the result section and how our filtering makes our definition more conservative. In any case, we do not consider this an essential aspect of the manuscript and have removed the statement that our definition of orphans is conservative.

9) Are orphans related (similar in sequence) to other orphans?

We clustered the protein sequences of orphan genes using the CD-hit software, widely used for gene family analyses. Generally, we find that only about 7% of the orphans show similarity to other orphans. We summarized the results of the orphan gene family analysis in the revised manuscript.

https://doi.org/10.7554/eLife.01311.028

Article and author information

Author details

  1. Nicola Palmieri

    Institut für Populationsgenetik, Vetmeduni Vienna, Vienna, Austria
    Contribution
    NP, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  2. Carolin Kosiol

    Institut für Populationsgenetik, Vetmeduni Vienna, Vienna, Austria
    Contribution
    CK, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  3. Christian Schlötterer

    Institut für Populationsgenetik, Vetmeduni Vienna, Vienna, Austria
    Contribution
    CS, Conception and design, Acquisition of data, Drafting or revising the article
    For correspondence
    christian.schloetterer@vetmeduni.ac.at
    Competing interests
    The authors declare that no competing interests exist.

Funding

Austrian Science Funds (FWF) (P22834)

  • Christian Schlötterer

The funder had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Ram Vinay Pandey for help with programming and Viola Nolte for support with genome annotations and library preparation. We are grateful to the members of the institute, in particular A Betancourt and R Kofler, for valuable discussion and comments on the manuscript.

Reviewing Editor

  1. Diethard Tautz, Reviewing Editor, Max Planck Institute for Evolutionary Biology, Germany

Publication history

  1. Received: July 30, 2013
  2. Accepted: January 19, 2014
  3. Version of Record published: February 19, 2014 (version 1)

Copyright

© 2014, Palmieri et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,870
    Page views
  • 267
    Downloads
  • 34
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    2. Genes and Chromosomes
    Wahid A Mulla et al.
    Research Article Updated
    1. Neuroscience
    Stefanie Hampel et al.
    Research Advance Updated