GBE Adaptation and Conservation throughout the Drosophila melanogaster Life-Cycle Marta Coronado-Zamora1,†, Irepan Salvador-Martınez2,3,†, David Castellano4, Antonio Barbadilla1,*, and Isaac Salazar-Ciudad1,2,5,* 1 noma de Barcelona, Cerdanyola del Vallès, Genomics, Bioinformatics and Evolution, Departament de Genètica i de Microbiologia, Universitat Auto Spain 2 Evo-Devo Helsinki Community, Centre of Excellence in Experimental and Computational Developmental Biology, Institute of Biotechnology, University of Helsinki, Finland 3 Department of Genetics, Evolution and Environment, University College London, United Kingdom 4 Bioinformatics Research Center, Aarhus University, Denmark 5 Centre de Recerca Matematica, Cerdanyola del Vallès, Spain † These authors contributed equally to this work. *Corresponding authors: E-mails: antonio.barbadilla@uab.cat; isaac.salazar@helsinki.fi. Accepted: April 16, 2019 Abstract Previous studies of the evolution of genes expressed at different life-cycle stages of Drosophila melanogaster have not been able to disentangle adaptive from nonadaptive substitutions when using nonsynonymous sites. Here, we overcome this limitation by combining whole-genome polymorphism data from D. melanogaster and divergence data between D. melanogaster and Drosophila yakuba. For the set of genes expressed at different life-cycle stages of D. melanogaster, as reported in modENCODE, we estimate the ratio of substitutions relative to polymorphism between nonsynonymous and synonymous sites (a) and then a is discomposed into the ratio of adaptive (xa) and nonadaptive (xna) substitutions to synonymous substitutions. We find that the genes expressed in mid- and late-embryonic development are the most conserved, whereas those expressed in early development and postembryonic stages are the least conserved. Importantly, we found that low conservation in early devel- opment is due to high rates of nonadaptive substitutions (high xna), whereas in postembryonic stages it is due, instead, to high rates of adaptive substitutions (high xa). By using estimates of different genomic features (codon bias, average intron length, exon number, recombination rate, among others), we also find that genes expressed in mid- and late-embryonic development show the most complex architecture: they are larger, have more exons, more transcripts, and longer introns. In addition, these genes are broadly expressed among all stages. We suggest that all these genomic features are related to the conservation of mid- and late-embryonic development. Globally, our study supports the hourglass pattern of conservation and adaptation over the life-cycle. Key words: adaptation, conservation, natural selection, evo-devo, DFE-alpha, hourglass hypothesis. Introduction there should be a relationship between divergence in species in The relationship between phylogeny and ontogeny has been a a phylogeny and different stages along development. Several hotly debated topic in evolutionary biology since its origins theories and hypotheses about this relationship have been pro- (Darwin 1872; Gould 1977; Raff 1996). Any change in an posed over the years (Gould 1977; Irie and Kuratani 2014). organism’s morphology is first a change in the developmental According to von Baer’s hypotheses or laws (1828), early de- process leading to that morphology (Alberch 1980) and, thus, velopment stages are the most similar between species within a morphological divergence between species in evolution implies phylogenetic group (e.g., vertebrates), whereas late develop- divergence in their development. It seems, thus, intuitive, that ment stages are the most divergent. ß The Author(s) 2019. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Genome Biol. Evol. 11(5):1463–1482. doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1463 Coronado-Zamora et al. GBE Multiple theoretical justifications for this “law” have been explore whether some stages are indeed more conserved proposed (Irie and Kuratani 2011). The most immediate is that than others. The dN/dS ratio, also known as Ka/Ks or x (x changes in early development can have consequences in later from now on), is widely used as a proxy for conservation at development, because late developmental processes are caus- the sequence level (Yang and Nielsen 1998; Hurst 2002). x ally dependent on the correct functioning of earlier develop- quantifies the level of constraint by comparing the rate of mental processes, whereas late developmental processes do nonsynonymous substitutions per nonsynonymous site (dN), not retroactively affect early developmental processes (Arthur which are presumably under selection, to the rate of synon- 1977). As a result, early development should be more ymous substitutions per synonymous site (dS), which are pre- constrained. sumably neutral. Departures of neutrality are detected when The advent of developmental genetics in mouse and x is different of 1. Davis et al. (2005) used this statistic on Drosophila and some observations in comparative embryol- 4,028 genes that were orthologous between Drosophila mel- ogy (Sander 1983) led to a different hypothesis: the hourglass anogaster and Drosophila pseudoobscura. By combining hypothesis of embryonic development evolution (Medawar these analyses with D. melanogaster microarray expression 1954; Slack et al. 1993; Duboule 1994). According to this data through development, they found that genes with the hypothesis, early and late development would be more diver- highest rates of nonsynonymous substitutions were expressed gent between species than intermediate developmental at low levels in late embryonic development and at high levels stages (middevelopment). There is no consensus about which in the larva, pupa, and adult. The genes with the lowest rates would be, exactly, these middevelopment conserved stages, of nonsynonymous substitution (the most conserved genes) but many researchers suggest that it should be some short were expressed at high levels in late embryonic development time after gastrulation (Wilt and Hake 2004). The most con- and at low levels before and after late embryonic develop- served stage within a phylum is called the phylotypic stage ment. This suggests, according to these authors, an hourglass (Sander 1983). There is no consensus either on what is actu- pattern where embryonic stages spanning from 12 to 22 h are ally conserved: embryonic morphology (von Baer 1828), de- highly conserved between D. melanogaster and Drosophila velopmental mechanisms (Raff 1996), or expression patterns pseudoobscura. of specific genes (Slack et al. 1993; Duboule 1994). In another study by Kalinka et al. (2010), the number of There is an extensive literature about the processes that nonsynonymous substitutions, dN, was calculated for 3,019 may lead to an hourglass pattern. Some propose that many genes in six different Drosophila species and their expression whole-body scale interactions take place during middevelop- was measured by microarrays in each species in eight 2- ment, whereas during early and late development interactions h intervals during development. Consistently with the hour- are at a much more restricted spatial scale (Raff 1996), both at glass model, the study found that middevelopment, around the level of mechanical interactions and molecular signaling the 10-h stage, is the period in which gene expression is the between tissues. Accordingly, changes in middevelopment most conserved among the six species. The set of genes would be much more likely to affect the whole embryo and expressed in those stages has also the lowest average dN changes at other stages would have much more spatially re- values across the six species. A similar study by Mensch stricted effects. Other authors argue that the hourglass pat- et al. (2013) estimated x for more than 2,000 genes across tern arises from different selection pressures acting in early six different Drosophila species for three categories of genes: and late development (Slack et al. 1993; Kalinka and maternal genes, genes expressed in early development, and Tomancak 2012). genes expressed in late development. Maternal genes (whose There is also a long ongoing discussion about whether the mRNA is left by the mother in the egg) and late embryonic hourglass hypothesis or the von Baer’s law (or something else) genes show higher x than early expressed genes. In contrast, fit the observed patterns of divergence among developmental another D. melanogaster study by Artieri et al. (2009) found stages in phylogenies (Richardson et al. 1997; Poe and Wake that genes expressed in the adult have higher x than genes 2004; Salazar-Ciudad 2010; Kalinka and Tomancak 2012; Hu expressed in the pupa and those of the pupa have higher x et al. 2017; Yang et al. 2018). In principle, the von Baer and that those expressed in the embryo (for 7,180 analyzed hourglass hypotheses could equally apply at the morpholog- genes), thus favoring von Baer’s law. In another study ical, gene expression, or genomic level and there are indeed (Drost et al. 2015), the x of the genes expressed in each stage studies at these three levels. In here, we briefly explain their are weighted by their expression to provide a stage- main conclusions to clarify why the questions, approaches, conservation measure. This measure supports the hourglass and results in this article significantly add to the understanding hypothesis in the zebrafish, D. melanogaster, and in the plant of the relationship between phylogeny and development. We Arabidopsis thaliana. Similar studies exist for other species of focus, however, on studies combining the transcriptomic and animals (Roux and Robinson-Rechavi 2008; Piasecka et al. genomic levels. 2013; Liu and Robinson-Rechavi 2018b). The latter of such Many studies measure conservation at the DNA sequence studies found that the genes expressed in middevelopment level for genes expressed at different developmental stages to tend to be expressed in more developmental stages and that, 1464 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE thus, they can be seen as more pleiotropic. This higher degree Other methods correcting for the potential biases of the of pleiotropy would then explain the higher degree of con- MKT estimates are based on the estimation of the distribution servation in middevelopment, at least in Caenorhabditis ele- of fitness effects (DFE of mutations from DNA sequence poly- gans, D. melanogaster, Danio rerio, the mouse Mus musculus morphism data) at functional sites (Bustamante et al. 2002, (Liu and Robinson-Rechavi 2018b), and in other eight species 2005; Eyre-Walker 2006; Eyre-Walker and Keightley 2007, of chordates (Hu et al. 2017). 2009; Keightley and Eyre-Walker 2007; Boyko et al. 2008). A limitation of the aforementioned studies is that x is not One of the most popular DFE-based methods is the DFE-alpha able to disentangle adaptive from nonadaptive substitutions. (Eyre-Walker and Keightley 2009). The DFE-alpha method In fact, nonsynonymous substitutions can turn out to be adap- corrects for the segregation of slightly deleterious substitu- tive, neutral, or slightly deleterious (strongly deleterious muta- tions by first estimating the DFE at selected sites by a gamma tions are purged from the population and so they do not distribution and then calculating how many nonadaptive sub- contribute to differences between species). In other words, stitutions are expected to become fixed given the inferred DFE a high x can be the result of relaxed selection (low conserva- from polymorphism data. Additionally, this method attempts tion), a high adaptive substitution rate, or a combination of to correct for possible effects of demography. Although this both. Because adaptive substitutions tend to be fixed quickly method is relatively recent, it has been used to estimate the (McDonald and Kreitman 1991), they will rarely be detected rate of adaptive nonsynonymous substitutions in a number of as polymorphic variants but only as divergent ones between studies (Strasburg et al. 2011; Carneiro et al. 2012; species. Thus, adaptive substitutions can be inferred when Kousathanas et al. 2014; Avila  et al. 2015; Santpere et al. there is an excess of nonsynonymous substitutions between 2015; Cagan et al. 2016; Galtier 2016; James et al. 2016; species relative to nonsynonymous substitutions within a pop- Murray et al. 2017; Steige et al. 2017). ulation, this is an excess of divergence in respect to polymor- There is no reason to expect that the patterns of adaptation phism. Thus, in the McDonald and Kreitman test (MKT) and conservation in the genes expressed over developmental (McDonald and Kreitman 1991), the divergence ratio (Dn/ stages should be directly related to each other. It could be, for Ds) is divided by the polymorphism ratio (Pn/Ps) and adaptive example, that middevelopment would be the most conserved substitutions are inferred if the ratio is larger than 1 [(Dn/Ds)/ stage but that it would still accumulate more adaptive non- (Pn/Ps) > 1]. From that MKT, the proportion of adaptive to synonymous substitutions than other stages of development nonadaptive substitutions that have been fixed, a, can be because nonsynonymous substitutions in the latter would be estimated (Smith and Eyre-Walker 2002). The divergence mostly neutral. We expect, however, that most adaptive sub- and polymorphism in synonymous sites, which are assumed stitutions occur in the adult, the larva, or the late stages of to be neutral, are used to estimate the underlying mutation embryonic or pupal development. The reason for that would rate and the expected polymorphism and divergence under a be that the larva and the adult have simply more functional neutral scenario. aspects of the phenotype to select from (more cell types, To try to circumvent known biases in the MKT, several more tissues, more organs), than developmental stages where modifications have been proposed over the years most of the phenotype is still in the making. In addition, be- (Bustamante et al. 2002, 2005; Smith and Eyre-Walker cause Drosophila is a holometabolous genus, the larva and 2002; Sawyer et al. 2003; Bierne and Eyre-Walker 2004; the adult effectively live in different habitats and, thus, one Mackay et al. 2012; Messer and Petrov 2013). Under station- may exhibit larger rates of adaptive nonsynonymous substitu- ary population size, slightly deleterious mutations lead to an tions than the other. underestimation of a because they tend to contribute more to In this article, we study 1) whether D. melanogaster devel- polymorphism than to divergence (Eyre-Walker 2002). opment follows the hourglass model or the von Baer’s law, 2) Because slightly deleterious substitutions tend to segregate whether there are differences not just in conservation but also at low frequency, its effect can be partially controlled by re- in the rates of adaptive substitutions between stages, and 3) moving low-frequency polymorphisms from the analysis (Fay how these rates are related to specific genomic features such et al. 2001). However, Charlesworth and Eyre-Walker (2008) as gene size, codon usage bias, average intron length, inter- showed that even removing low-frequency variants, a esti- genic distance, Guanine–Cytosine (GC) content, number of mates are always downwardly biased. Mackay et al. (2012) protein–protein interactions (PPI), number of exons and tran- proposed an extension to the MKT, which we will call here scripts, expression bias, or recombination rates. extended MKT (eMKT), in which the count of segregating Question (2) has been recently approached through a sites in nonsynonymous sites is partitioned into the completely different method (Liu and Robinson-Rechavi number of neutral variants and the number of weakly dele- 2018a): the branch-site likelihood test (Zhang et al. 2005). terious variants. This increases the power for detecting This test is an extension of the x ratio (Yang and Nielsen selection and allows the independent estimation of both 1998; Hurst 2002), for the detection of positive selection for adaptive and weakly deleterious selection (see Materials and at least some codons in a phylogenetic branch of interest. Methods). Using this method, Liu and Robinson-Rechavi study three Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1465 Coronado-Zamora et al. GBE phylogenetic branches, the Clupeocephala branch, the 2011). The data set contains the expression data for 30 stages Murinae branch, and the Melanogaster group branch for fo- of the whole life-cycle of D. melanogaster, including 12 em- cal species Danio rerio, Mus musculus, and D. melanogaster, bryonic samples collected at 2-h intervals for 24 h, six larval, six respectively. The test contrasts two hypotheses: no positive pupal, and three sexed adult stages at 1, 5, and 30 days after selection occurred (H0) in the phylogenetic branch of interest eclosion (Graveley et al. 2011). These data were downloaded versus at least some codons experienced positive selection from FlyBase (release 6.06, July 2015, file: “gene (H1). One major advantage of the test is that it allows positive _rpkm_report_fb_2015_03,” last accessed December selection to vary both among codon sites and among phylo- 2015). The downloaded file reports gene expression values genetic branches. In our study, we use tests (MKT, DFE-alpha) as reads per kilobase per million reads (RPKM). RPKM values that use both polymorphic and divergence data. Polymorphic are provided only for exonic regions of the gene (excluding data allow taking into account purifying selection in the di- segments that overlap with other genes), except for genes vergent x ratio, strongly increasing the power of detecting derived from dicistronic/polycistronic transcripts, where all positive selection. These tests covers a shorter time scale, from exon regions were used for the estimation of RPKM expres- the divergence of the outgroup species to the present, low- sion. See Gelbart and Emmert (2013) for details. All RPKM ering the contingencies associated with longer evolutionary values were log2-transformed. processes in the branch-site likelihood test. These include, for We tried different criteria to consider whether a gene is example, the poor sequence alignments due to increased di- expressed in a stage or not. The following five criteria were vergence, which in turn can result in false positives. applied: A “low stringent criterion,” in which a gene is con- In our study, we used the expression data in the sidered expressed in a stage if its RPKM is larger than 0 in that modENCODE project (Graveley et al. 2011) from FlyBase. stage. This criterion is used in many other publications (e.g., This is the most complete gene expression database through Hebenstreit et al. 2011; Wagner et al. 2013; Guillen et al. D. melanogaster life-cycle (it includes 17,788 genes over most 2019). A “low stringent criterion with 2-fold differential developmental and life-cycle stages). We used divergence and expression,” in which a gene is considered expressed in a polymorphism data for the genes expressed in each develop- stage if its RPKM is larger than 0 in that stage and if, in ad- mental stage to estimate four selection statistics a, x, xa, and dition, its maximal gene expression (over all the stages) is at xna. xa is x  a, or the rate of adaptive substitutions per least twice that of its minimal gene expression (also over all nonsynonymous site divided by the synonymous substitution the stages). A “low stringent criterion with a 4-fold differen- rate (dS) (Gossmann et al. 2012), whereas xna is x(1  a), tial expression” that is as the 2-fold criterion but with a 4-fold the rate of nonadaptive substitutions divided by dS (Galtier differential expression criterion. A “medium stringent 2016). In other words, xa is a measure of the rate of adaptive criterion,” in which a gene is considered expressed if its nonsynonymous substitutions, xna is the same measure for RPKM is equal or higher than 2 and a “high stringent nonadaptive substitutions. Thus, xa informs about the overall criterion” in which a gene is considered expressed if its rate of adaptive substitutions occurred in respect to the neu- RPKM is equal or higher than 10. This criterion is also used tral mutation rate, the overall level of adaptation, whereas a as a stringent criterion in RNA-seq analysis by other authors simply measures how much common are adaptive substitu- (e.g., Dezso et al. 2008). Supplementary figure S1, tions in respect to nonadaptive ones (but does not tell if there Supplementary Material online, shows the number of genes has been many of them or not, as xa does). Polymorphism for each criterion and stage. As a result of applying each of data come from the Drosophila Genetic Reference Panel these five criteria on the same modENCODE RNA-seq data, (DGRP) project (Mackay et al. 2012) and divergence comes we obtain five different lists of genes for each developmental from comparing the genomes of D. melanogaster and and life-cycle stages. Drosophila yakuba. Choosing D. yakuba as outgroup has been reported to provide a more reliable estimation of the selection statistics than the closest Drosophila simulans Gene Expression Profile Clustering (Keightley and Eyre-Walker 2012) (see Materials and To identify shared temporal expression patterns among the Methods). genes of the modENCODE RNA-seq experiments, we applied a soft clustering method to the log2-transformed expression values. We used a fuzzy c-means algorithm with the mfuzz() Materials and Methods function of the R package Mfuzz (Futschik 2015). The Mfuzz Gene Expression soft clustering algorithm uses Euclidean distance as distance statistic and requires two main parameters (c ¼ number of Gene Expression through the Developmental and Life- clusters and m ¼ fuzzification parameter). For the clustering, Cycle Stages we z-standardized the log2-transformed expression values, so Gene expression data of 17,875 genes comes from RNA-seq that the average expression value for each gene is 0 and the experiments in the modENCODE project (Graveley et al. standard deviation is 1. The fuzzy soft clustering method is 1466 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE different from hard clustering (like hierarchical clustering), in Additionally, we also removed genes that were exclusively that genes are not uniquely assigned to one cluster. Instead of expressed in one sex but no in the other. See supplementary this, a gene i has gradual degrees of membership lij to a table S4, Supplementary Material online, for the genes cluster j. High membership values indicates a high correlation analyzed. between gene i with the cluster centroid cj (Futschik 2015). To control for the effect of the sex chromosome, we dis- The mfuzz() function uses the fuzzy c-means algorithm, based carded X-linked genes from the analyses. See supplementary on minimization of a weighted square error function with table S5, Supplementary Material online, for the genes which the clusters centroids cj result from the weighted sum analyzed. of all clusters members. The membership value (lij) indicates how well the gene i is represented by cluster j. Soft clustering Maternal, Maternal–Zygotic, and Zygotic Genes is especially useful when clusters are not well defined, as in gene expression time-course data (Futschik and Carlisle A list of maternal, maternal–zygotic, and zygotic genes was 2005). We discarded genes that were constitutively expressed obtained from Thomsen et al. (2010) data using egg and early in all stages. Genes having a cluster membership lower than development microarray analyses. Maternal genes are those 0.8 were excluded. This resulted in 3,819 genes in eight clus- genes whose mRNA is laid in the egg and are never tran- ters (one of the clusters was discarded as it was mostly noise) scribed by the embryo, maternal–zygotic those genes whose for the embryonic development out of the 5,514 genes that mRNA is laid in the egg by the mother but that are also tran- are expressed during the embryo development under the low scribed by the embryo and zygotic genes are genes whose stringent criteria, and 8,167 genes in nine clusters for the life- mRNA is not laid in the egg by the mother. The maternal gene cycle of the 9,241 genes expressed in the whole life-cycle list was obtained joining the original Thomsen’s categories for (discarding females). The values of c and m were optimized not transcribed genes: “maternal decay,” “mixed decay,” using the procedure described in Futschik and Carlisle (2005) “stable,” and “zygotic decay” categories (4,942 genes). and Futschik (2015), resulting in c ¼ 9 for both data sets and The maternal–zygotic gene list was created by joining the m ¼ 1.23 and 1.08 for the embryonic development and life- categories of genes that are both present in the egg and cycle, respectively. Supplementary tables S1 and S2, transcribed later (1,332 genes analyzed): “maternal decay þ Supplementary Material online, show the number of genes transcription” and “stable transcription” categories. Finally, expressed in each cluster for the two analyzes (for the embryo the zygotic gene list equals the “purely zygotic” category development and the whole life-cycle, respectively). (850 genes). In the end, we get three lists of genes, one for the maternal genes, one for the maternal–zygotic, and one for the zygotic genes. See supplementary table S6, Testis and Immune Genes Supplementary Material online, for the genes analyzed. For determining testis- and immune-related genes, we imple- mented the developed pipeline in Castellano et al. (2016). We Population Genomics downloaded the Gene Ontology (GO) terms for our gene set through the R package biomaRt (Durinck et al. 2005) using Selection Statistics the D. melanogaster ENSEMBL database. When a gene was We used three different tests to estimate the selection statis- associated to any term related to the testis or the immune tics on the lists of genes expressed in each life-cycle stage: system (see supplementary table S3, Supplementary Material MKT standard (McDonald and Kreitman 1991), eMKT online, for the related GO terms) we removed it from the data (Mackay et al. 2012), and DFE-alpha (Eyre-Walker and set (a total of 171 out of 2,869 genes were removed). Those Keightley 2009). For the gene expression profiles and for ma- 171 exhibit higher rates of adaptation as it would be expected ternal, maternal–zygotic, and zygotic gene categories, we (permutation test, xa, P ¼ 0.028; x, P < 0.001) when com- used only the DFE-alpha method. pared against our whole data set (6,690 genes). These three methods rely on polymorphism and divergence data to estimate adaptation. For obtaining this data, coding Sex-Biased Genes and Sex-Linked Genes exon annotations from D. melanogaster were retrieved from FlyBase (release 5.50, www.flybase.org , last accessed March To control for the effect of genes with a sex-biased expres- 2013). Genes 1:1 orthologs across D. yakuba–D. mela- sion, we estimated the gene sex-bias ratio with the nogaster were obtained from FlyBase (www.flybase.org, last expression:   accessed March 2013). We used D. yakuba as outgroup spe- RPKM$ cies because, aside from its high coverage (9.1, Clark et al. sex bias ratio ¼ log2 : RPKM# 2007), the time elapsed because its divergence from D. mel- anogaster (7.4 Ma, Tamura et al. 2004) is more suitable for And we discarded genes with a sex-biased above a >j2j estimating adaptation (Keightley and Eyre-Walker 2012). In threshold in any of the three sampled adult points. closely related species (as is the case of D. melanogaster and Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1467 Coronado-Zamora et al. GBE D. simulans, which diverged 4.3 Ma [Cutter 2008]), the esti- of the ancestral state (Hernandez et al. 2007) and because it mated rate of adaptive substitution can be biased due to 1) an performs well for inferring deleterious DFE (Eyre-Walker and erroneous attribution of polymorphism to divergence, 2) an- Keightley 2007; Boyko et al. 2008; Tataru et al. 2017). We cestral polymorphism contributing to divergence, and 3) dif- used a custom-made Perl script to estimate the number of ferences in the rate of fixation of neutral and adaptive short-intron substitutions and to compute the folded SFS. mutations (Keightley and Eyre-Walker 2012). The authors Jukes and Cantor correction for multiple hits was applied (Keightley and Eyre-Walker 2012) find that the adaptive (Jukes and Cantor 1969). For computational reasons, one of rate estimated from closely related species (as in the case of the programs used for estimating adaptation (DFE-alpha, D. melanogaster and D. simulans) may be underestimated by Eyre-Walker and Keightley 2009, see below) needs that all 10% or more. Choosing D. yakuba as outgroup gives a more sites are sampled in the same number of individuals (Eyre- reliable estimation of the adaptive rate than the closest D. Walker and Keightley 2009). Hence, the original data of simulans (Keightley and Eyre-Walker 2012). 168 lines set have been reduced to 128 isogenic lines by ran- The population genomic data comes from 168 inbred lines domly sampling the polymorphisms at each site without re- of D. melanogaster sequenced in the Freeze 1.0 of the DGRP placement. Finally, residual heterozygous sites and sites with project (Mackay et al. 2012). The DGRP population was cre- no quality value were excluded from the analysis. ated collecting gravid females from a single population of Raleigh, NC, followed by the full-sibling inbreeding approach MKT, Standard McDonald and Kreitman Test during 20 generations to obtain full homozygous individuals. The MKT (McDonald and Kreitman, 1991) was developed to After this, the residual heterozygosis in the samples is be applied to protein-coding sequences, combining both di- expected to be 1.4% (inbreeding coefficient F ¼ 0.986). The vergent (D) and polymorphic (P) sites, and categorizing muta- expected 1.4% of residual heterozygosis was true for 90% of tions as synonymous (PS, DS) and nonsynonymous (PN, DN). If the sequenced chromosome lines. DGRP lines showing high all mutations are either strongly deleterious or neutral, then values of residual heterozygosity (>9%) were observed to be DN/DS is expected to roughly equal PN/PS. In contrast, adaptive associated with large polymorphic inversions and, they were mutations rapidly reach fixation and thus contribute relatively not included in our analyses (Huang et al. 2014). The inbreed- more to divergence than to polymorphism when compared ing approach can alter the frequency spectrum of the strongly with neutral mutations, and then DN/DS > PN/PS. Assuming deleterious recessive mutations. However, alternative resour- that adaptive mutations contribute little to polymorphism but ces such as the Drosophila Population Genomics Project substantially to divergence, the proportion of adaptive non- (DPGP) (Langley et al. 2012) would encounter the same prob- synonymous substitutions than have been fixed by selection lem. In fact, previous works comparing adaptation and DFE can be inferred as a ¼ 1  (PN/PS)/(DN/DS) (Smith and Eyre- estimates between DGRP and DPGP data sets have shown no Walker 2002). The significance of effect can be easily quan- differences between those data sets (Castellano et al. 2016, tified using a simple 2  2 contingency table (see table 1). 2018). And given that the DPGP sample size is smaller than the DGRP sample size, it is likely that these mutations contrib- ute very marginally to the estimations of polymorphisms, DFE, eMKT or Extended MKT and adaptation in both databases. The null hypothesis of neutrality is rejected in a MKT not only A multiple genome alignment between the DGRP isogenic when DN/DS > PN/PS inferring adaptation but also when PN/PS lines (Mackay et al. 2012) and the D. yakuba genome using > DN/DS. In this later case, there is an excess of polymorphism the BDGP 5 coordinates (Berkeley Drosophila Genome Project relative to divergence for the nonsynonymous class N, due to 5) was obtained from the publicly available database at http:// 1) slightly deleterious variants segregating at low frequency in popdrowser.uab.cat (Ramia et al. 2012) (last accessed May the population, which contribute to polymorphism but not to 2010). For each gene, we took all nonoverlapping coding divergence or 2) relaxation of selection where sites previously exons, independently of their inclusion levels. When two under strong or weak purifying selection have become neu- exons overlapped, the largest was chosen for subsequent tral, causing an increased level of polymorphism relative to analyses. Only exons without frameshifts, gaps, or early stop divergence. codons were retained. In this way, we tried to avoid potential Adaptive mutations and weakly deleterious selection act in alignment errors that would inflate our mutation rate esti- opposite directions on the MKT, so a will be underestimated mates and create an artifactual positive correlation between when the two selection regimes occur. Because slightly dele- them. Exonic sequences were trimmed in order to contain terious mutations tend to segregate at lower frequencies than only full codons. We calculated the number of substitutions do neutral mutations, they can be partially controlled for by and the folded site frequency spectrum (SFS, Ronen et al. removing low-frequency polymorphisms from the analysis, 2013) of the minor allele frequency (MAF) for short introns, generally the 5% (Fay et al. 2001). However, this method is 0-fold and 4-fold degenerate sites, using an ad hoc Perl script. still expected to lead to biased estimates (Charlesworth and The SFS was folded to avoid difficulties with misidentification Eyre-Walker 2008). To take adaptive and slightly deleterious 1468 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE Table 1 Standard MKT and eMKT Tables Standard MKT Table Number of Segregating Sites by MAF Category Site Class Polymorphism Divergence Site Class MAF < 5% MAF  5% S (neutral class) PS DS S (neutral class) PS MAF < 5% PS MAF  5% N (selected class) PN DN N (selected class) PN MAF < 5% PN MAF  5% mutation into account, PN, the count of segregating sites in of fixed substitutions that are adaptive (a) and the rate of the nonsynonymous class, should be separated into the num- adaptive substitutions relative to the neutral rate (xa, esti- ber of neutral variants and the number of weakly deleterious mated as x  a, Gossmann et al. 2012). Furthermore, in variants, PN ¼ PN neutral þ PN weakly del. If both numbers are our analysis, we include another statistic, xna (estimated as estimated, adaptive and weakly deleterious selection can be x  xa), which represents the proportion of nonadaptive evaluated independently. Consider the following pair of 2  2 substitutions (slightly deleterious and neutral) relative to the contingency tables (table 1). neutral rate (Galtier 2016). Thus, we are able to decompose The table on the left represents the standard MKT table the classical x ratio, into these two statistics: xa and xna, and with the theoretical counts of segregating sites and divergent differentiate whether high rates of x are due to adaptive or sites for each cell. The table on the right contains the count of nonadaptive nonsynonymous substitutions. We estimated PN and PS for two-frequency categories. Mackay et al. (2012) natural selection on coding regions under a two-epoch de- showed using DGRP data that the ideal threshold to separate mographic model. neutral and weakly deleterious mutations is a 5%, and that Finally, to estimate these selection statistics with DFE-alpha, similar results were obtained using a 10% threshold, suggest- it is necessary to concatenate data from several genes because ing that most of the weakly deleterious variants are below a estimates from a single gene cannot be obtained due to the 5%, whereas the mutations above can be considered as lack of segregating (or divergent) sites for some site classes. nearly neutral variants. The estimate of the fraction of sites Thus, we calculated the selection statistics for the genes segregating neutrally within the MAF < 5% is f neutral MAF<5% expressed in each stage and calculated the confidence inter- ¼ PS MAF<5%/PS. The expected number of segregating sites in vals by randomly resampling the genes expressed in each the nonsynonymous class which are neutral within the stage, cluster or gene category (bootstrapping) 100 times MAF < 5% is PN neutral MAF<5% ¼ PN  fneutral MAF<5%. The with replacement (see Bootstrapping section). expected number of neutral segregating sites in the nonsy- nonymous class is PN neutral ¼ PN neutral MAF<5% þ PN MAF5%. Genomic Features To estimate a from the standard MKT table correcting by the segregation of weakly deleterious variants, we have to For the genes in each list and gene expression profile, we substitute the PN by the expected number of neutral segre- measured a number of genomic features. Coding exons gating sites, PN neutral. The correct estimate of a is then a ¼ 1  and short-intron annotations from D. melanogaster were (PN neutral/Ps)/(DN/DS). obtained from FlyBase (release 5.50, www.flybase.org, last accessed March 2013). Average intron length. It is the average distance, in base DFE-alpha pairs, between the successive exons of a gene. The effect of DFE-alpha (Eyre-Walker and Keightley 2009), an extension of average intron length and total intron length into x has been the MKT, also corrects for the segregation of slightly delete- reported to be very similar (Marais et al. 2005). rious alleles, providing a more accurate estimation than the Intergenic distance. It is the average distance, in base pairs, MKT and other methods that do not take polymorphism data between the two closest genes to a given gene. into account. Briefly, this software uses a maximum-likelihood Gene size. Length of the coding DNA of a gene sequence method based on polymorphism data to infer the DFE of new (CDS). mutations. It assumes two classes of sites in the genome: Number of exons and transcripts. Number of different neutral sites (synonymous) and selected sites (nonsynony- exons and transcripts of a gene, respectively. For the latter, mous) and contrasts the SFS at these two classes. As a neutral we count all the transcripts reported for a gene in Fly Base. reference, we used positions 8–30 of short introns (65 bp) This is not the same as the number of different transcripts following Halligan and Keightley (2006) (and 4-fold degener- expressed in a stage. ated sites for some cross-validating analysis). As selected sites, Codon bias. Measured as the frequency of optimal codons, we used 0-fold degenerate sites. Provided the SFS at both Fop. We used the software CodonW for the estimation of Fop neutral and selected sites together with divergence data, (Peden 1999, www.codonw.sourceforge.net, last accessed the DFE-alpha method allows calculating of the proportion June 2012). The index is estimated as the ratio of optimal Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1469 Coronado-Zamora et al. GBE codons to synonymous codons. Values range between 0 (no permutation test was applied. For obtaining a null distribution optimal codons are used) and 1 (only optimal codons are for the differences between gene groups, we shuffled with- used). out replacement 1,000 times the complete list of genes by Expression bias. Proportion of development stages in which means of ad hoc bash and Perl scripts. We estimated the rates a gene is expressed according to modENCODE. Following of adaptive nonsynonymous substitutions and nonadaptive Yanai et al. (2005) and Larracuente et al. (2008), we esti- synonymous substitutions in each randomized list, obtaining mated the expression bias, s, as (eq. 1): an expected null distribution. The two-tailed P value was Pn logðSj Þ obtained by counting the number of replicates below and j¼1 1  logðSmax Þ above the observed difference divided by the total number s¼ ; (1) n1 of replicates (i.e., 1,000), thus obtaining a two-tail P value. Multiple comparisons for each analysis were corrected by the where log(S) is the logarithm of the RPKM and n is the num- false recovery rate approach. ber of developmental stages. s ranges from 0 to 1. Values See in the following supplementary table S7, close to 0 indicate broadly expressed genes and values close Supplementary Material online, a summary of the permuta- to 1 indicate genes with highly biased expression. s ¼ 1 means tion analysis performed in this study and supplementary figure that expression is only detectable in a single sample. S2, Supplementary Material online, for a graphical summary Recombination rates. Recombination rates estimates at of the permutation procedure. 100-kb nonoverlapping windows and microscopically ob- served crossing-over (CD) events were obtained from Spearman’s Rank Correlations Comeron et al. (2012). CG content. We used the software CodonW (Peden 1999, Correlations between temporal profiles were carried www.codonw.sourceforge.net, last accessed June 2012) to out bySpearman’s rank correlations, calculated by using the calculate the GC content of each protein-coding gene. cor.test() function of the R program (R Core Team 2015). PPI. We downloaded PPIs for 10,631 protein-coding genes available at the DroiD database, version DroID_v2018_08 Spearman’s Partial Correlations (Murali et al. 2011, http://www.droidb.org/Downloads.jsp; last accessed January 2019). We integrated the information In order to test for the joint effect of ten genomic features on of six data sets available in DroiD: PPI Curagen yeast two- the selective statistics, we performed three independent mul- hybrid, DPIM co-AP/MS, Finley Lab yeast two-hybrid, FlyBase tiple regression linear models, one for each selection statistic curated, from other databases (BioGRID, IntAct, MINT, and (x, xa, and xna). BIND) and Hybrigenics yeast two-hybrid and counted the total The selection statistics were considered response variables unique protein interactions of each protein-coding gene. The in each model and the ten genomic features as main effects number of genes included in FlyBase curated database was (no interaction terms were included). Some of the ten geno- too small to base our analysis only on it. mic features showed significant pairwise correlations (supple- mentary table S16, Supplementary Material online) as determined using Spearman’s partial correlations (Kim Statistical Analysis 2015). The analysis using multiple regression models were Bootstrapping intended to disentangle their effects on the selection statistics. To assess the relative importance of each analyzed genomic As explained before, to estimate the rate of adaptive evolu- feature in each linear regression model, we used the pmvd tion with the DFE-alpha method (Eyre-Walker and Keightley metric introduced by Feldman (2005) included in the R pack- 2009), it is necessary to pool several genes together. The sam- age relaimpo (Grömping 2006), that averages the R2 of each pling distribution of parameters estimated by the DFE-alpha regressor over all possible orderings using weighted averages. method at each stage, cluster or gene category was estimated by randomly sampling 100 times with replacement the genes at each temporal stage (or cluster or category). This concate- Results nation was not necessary for the MKT and eMKT methods Overall Temporal Pattern of Adaptation and Conservation and, in fact, the same results are obtained with and without concatenation. We calculated four selection statistics (three based on poly- morphism and divergence data xa, xna, and a and one based on divergence data alone, x) for the set of genes expressed in Permutation Test each developmental and life-cycle stage. To assess whether stages, gene clusters, or gene categories We estimated these statistics using three different meth- undergo differential selection compared with the genes not ods, DFE-alpha (Eyre-Walker and Keightley 2009) (fig. 1), the expressed in such stage, gene cluster or gene category, a standard MKT (McDonald and Kreitman 1991) 1470 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE (supplementary fig. S3, Supplementary Material online), and the null distribution, we shuffled without replacement the eMKT (Mackay et al. 2012) (supplementary fig. S3, 1,000 times the complete list of genes expressed during Supplementary Material online), and using different criteria the whole life-cycle (2,869 genes) and randomly assigned to consider whether a gene is expressed in a stage or not the genes into one of the two groups. We estimated the (supplementary figs. S4–S6, Supplementary Material online). selection statistics in each randomized list, obtaining an In these calculations, the mutation rate was estimated using expected null distribution than we then compared with short introns (fig. 1) and 4-fold degenerated sites (supplemen- the observed difference. Multiple comparisons were cor- tary fig. S7, Supplementary Material online). Finally, we also rected by the false recovery rate approach. This analysis explored whether the number of genes in a stage has an shows that mid- and late-embryonic development, the effect on the values of the estimated metrics. For that, we beginning of the larva and genes expressed in female repeated the analyses by sampling 350 genes per stage (with adults show significantly low rates of nonsynonymous replacement) 100 times and calculating the mean values for substitutions (supplementary fig. S9, Supplementary the selection metrics (supplementary fig. S8, Supplementary Material online). The relatively high rates of substitutions Material online). For each combination of methods and ex- in early development and in the larva, pupa, and males are pression criteria, we obtain very similar patterns of change of not significant in this analysis (P values can be consulted xa, xna, x, and a over stages (what we call the temporal on supplementary table S10, Supplementary Material pattern of these selection statistics). online). Figure 1 shows the temporal pattern found when using the However, if the same permutation test is done using all DFE-alpha method, short-intron sites as a proxy for the mu- genes (expressed in all stages or not) the test shows that, in tation rate, and considering all the genes with nonzero ex- addition, early development, late larva, pupa, and male adult pression (excluding the 6,655 genes that were constitutively exhibit significantly high xa, x, and a values (supplementary expressed throughout all stages, see Materials and Methods fig. S10, Supplementary Material online). See the P values of and supplementary table S8, Supplementary Material online, all comparisons in supplementary table S11, Supplementary for the genes analyzed). A total of 2,869 genes are considered Material online. with this criterion. The rate of adaptive (xa) and, especially, Immune system and testis-related genes have been nonadaptive (xna) nonsynonymous substitutions on the set of reported to be under faster rates of adaptation than other genes expressed in each developmental stage (fig. 1) was the genes (Pröschel et al. 2006; Obbard et al. 2009). Immune highest in the first embryonic stage. Both rates gradually de- system genes are expressed mostly in adults, and testis genes crease until the 10-h embryo stage. The next developmental are expressed only in adult males. It could then be that the stages (across mid- and late-embryonic development) show, low rates of conservation in the genes expressed in the male on the contrary, the lowest rates of fixation of substitutions adult stages would be explained by the low rates of conser- (either adaptive or not). At the third larval stage (L3), the rates vation of immune- and testis-related genes. To explore this of adaptive substitution (xa), and to a lesser extent of non- possibility, we repeated our analysis excluding immune sys- adaptive substitution (x), increase and remain high through tem and testis genes. The exclusion of these genes does not all the pupal stages. Finally, in the male adult stage, xa and x modify our results much (supplementary fig. S11, values are very similar to those of the pupa, whereas female Supplementary Material online, see supplementary table adults exhibit lower values. S3, Supplementary Material online, for the GO terms ex- Overall, all the selection statistics show a similar tem- cluded), so the high rates of substitution we observed in poral trend except for nonadaptive substitutions being adult males are not exclusively due to these immune- and seemingly more abundant in the genes expressed in the testis-related genes. We additionally controlled for the effect very early stages as compared with later stages (a differ- of genes with a sex-biased expression (see Materials and ence we further analyze in a coming section). To analyze Methods) by removing genes with a sex-biased ratio above whether these differences between stages were statisti- a ratio of 2 (supplementary fig. S12, Supplementary Material cally significant, we merged stages into eight develop- online). As expected, sex-biased genes are responsible of the mental periods: embryo 0–2 h, embryo 2–6 h, embryo observed differences between males and females adaptation 6–24 h, larva 1–3, larva 4–6, pupa, female, and male rates. When removing sex-biased genes the female stages adults (see supplementary table S9, Supplementary exhibit selection statistics values similar to those of males. Material online, for the genes analyzed). We calculated This indicates that the genes that are overexpressed in by a permutation test the chances that the genes females tend to be more constrained. expressed in a period undergo differential selection com- X-linked genes are known to be under faster rates of evo- pared with the genes not expressed in that period. First, lution compared with autosomal genes (faster-X we calculated the difference in selection statistics be- effect) (Meisel and Connallon 2013). When we excluded X- tween the two groups, that is, the genes of a period linked genes from the analyses (supplementary fig. S13, and the genes not expressed in that period. To obtain Supplementary Material online) the selection pattern Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1471 Coronado-Zamora et al. GBE FIG. 1.—Temporal pattern of the four selection statistics estimated with DFE-alpha (xa, a, xna, and x). (A) xa, the rate of adaptive nonsynonymous substitutions relative to the mutation rate. (B) a, the proportion of base substitutions fixed by natural selection. (C) xna, the rate of nonadaptive non- synonymous substitutions relative to the mutation rate. (D) x, the rate of nonsynonymous substitutions relative to the mutation rate. Each boxplot (A–E, 100 bootstrap replicates per stage) in a plot is calculated for a randomly drawn sample of the set of genes expressed in a stage with replacement. The solid line going through the boxplot is inferred by LOESS. For the male and female stages, the line is simply a linear regression. The dashed line shows the mean value of each statistic for the genes that are expressed in all stages (again with 100 bootstrap replicates). The embryonic stages are named by the hour’s intervals (from 0 to 24 h), the larval stages are the first instar (L1), second instar (L2), and third instar (L3). The L3 stages are subdivided into the first 12 h (L3-12h) and several puff stages (L3-PS1 to L3-PS7). WPP is the white prepupae stage. The pupal stages with RNA-seq are phanerocephalic pupa, 15 h (P5), 25.6 h pupa (P6), yellow pharate, 50.4 h (P8), amber eye-pharate, 74.6 h (P9–10), and green meconium pharate, 96 h (P15). Adult stages are 1, 5, and 30 days after eclosion (1, 5, and 30 days). Number of genes analyzed is in supplementary table S8, Supplementary Material online. The earliest stages show more variation in the selection statistics because they have less genes. remained similar, so the rates of substitutions observed are the high xa in a time period and the expression of a specific not exclusively due to X-linked genes. subset of genes in it. To explore these possibilities, we categorized all the ana- lyzed genes into classes based on their temporal profiles of Gene Expression Profiles Clustering expression. To do that we use an unsupervised soft clustering There are at least three different scenarios that could explain algorithm (Futschik and Carlisle 2005) as explained in the the observed temporal pattern of change in the selection sta- Materials and Methods. Genes within each temporal expres- tistics. In the case of xa, for example, it could be that a subset sion class show relatively similar changes in gene expression of genes with high xa is expressed just with the observed levels over time. We consider eight such classes for embryonic temporal pattern (everywhere but in mid- and late-develop- development (fig. 2A, number of genes analyzed in each clus- ment and in females). Alternatively, it could be that each of ter in supplementary table S1, Supplementary Material online) the time periods with high xa would express a distinct group and nine classes for the whole life-cycle (see supplementary of genes that have high levels of adaptive substitutions. It fig. S14, Supplementary Material online, for the temporal could also be that no simple correspondence exists between profiles, number of genes analyzed in each cluster in 1472 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE supplementary table S2, Supplementary Material online). For content of each gene, and the number of PPIs (Murali embryonic development, clusters 1 and 2 are the ones show- et al. 2011). ing significantly the highest xa and x compared with First of all, we assessed the relationship between the ge- the other clusters using a permutation test (cluster 1, x: nomic features and selection statistics, clumping all expressed P < 0.001; xa: P ¼ 0.008; cluster 2, x: P < 0.001; xa: genes together irrespectively of the stage in which they are P ¼ 0.059) (fig. 2B–E). These clusters correspond to the genes expressed. For that, we categorized the genes in our data set that are expressed at high levels in the earliest development in five different categories based on each genomic feature and that rapidly decrease their expression to very low levels (see supplementary table S14, Supplementary Material online, (fig. 2A). xa, x, and a values in clusters 1 and 2 are larger than for the number of genes considered in each category) and those in the first three developmental stages and, thus, it is resampled with replacement 100 times the genes in each likely that the genes in these clusters, but not the other genes category and estimated the selection statistics in each such that are expressed in these three stages, are responsible for category. We found clear negative correlations between xa the high xa, x, xna, and a values in the earliest development. and gene size, number of exons, number of transcripts, aver- The decline in the values of these selection statistics over early age intron length, codon bias, GC content, and PPIs (similar development would then just be a simple reflection of the correlations were found for xna [supplementary fig. S16, decrease in expression of the genes in those clusters over Supplementary Material online], x, and a [data not shown]) time. Cluster 8 also shows larger x than that observed in and positive correlations between xa and expression bias, re- the other clusters (x: P < 0.001). This cluster is composed of combination rate, and intergenic distance (supplementary fig. genes whose expression increases only in the last stages of S17, Supplementary Material online). Thus, in general, less embryonic development. This high xa cannot be detected complex genes (shorter, with fewer exons and transcripts when directly analyzing the genes in each stage because and short introns), expressed in only a small number of stages, the other genes expressed in these late stages have lower x with lower GC content and/or less PPIs are more likely to values, as it can be seen for cluster 5, which expresses genes accumulate adaptive substitutions (high xa). from the 10-h onward. Thus, the temporal pattern of change We performed linear multiple regression models to test for in the selection statistics seems to come from the temporal the dependence of the selective statistics (x, xa, and xna) dynamics of expression of three different sets of genes (those estimated on each protein-coding gene using the eMKT on of clusters 1, 2, and 8). P values can be consulted on supple- the ten genomic features (supplementary table S15 and fig. mentary table S12, Supplementary Material online. Similar S18, Supplementary Material online). Supplementary table results were found when the clustering was done over the S16, Supplementary Material online, show the matrix of cor- whole life-cycle (supplementary fig. S15 and table S13, relation coefficients from the pairwise correlations between Supplementary Material online). the ten genomic features included in the linear models. The coefficient of determination (multiple R2) was the highest for x, explaining a 21.21% of the total variance Genomic Features (P ¼ 2.2  1016). For xa, it was 3.34% (P ¼ 2.2  1016) To search for sequence properties affecting the tempo- and for xna it was 9.36% (P ¼ 2.2  1016). Average intron ral pattern of the selection statistics, we measured a length, size, and PPI are the only genomic features that do not number of genomic features for each of the genes have a significant effect on any of the selective regimes. expressed in development (see Materials and Expression bias is the genomic feature with the highest effect Methods). These are gene size (the length of the coding on x, as observed in Guillen et al. (2019). sequence of a gene), number of exons (total number of When considered over developmental stages we found exons of a gene; not to confuse with the number of that these genomic features exhibit a temporal pattern that exons a gene expresses in a specific developmental is either very similar to that of xa or its opposite (see fig. 3). To stage), codon usage bias (measured as frequency of op- analyze that we calculated the correlation between each timum codons, Fop), number of transcripts per gene stage’s average xa and the average of each genomic feature (this is not the number of alternative transcripts of a in each stage (see table 2). Gene size, number of exons, codon gene in a stage but the total number of alternative tran- usage bias, and number of transcripts per gene follow a tem- scripts of a gene according to FlyBase annotations), av- poral pattern that is the reverse of that of xa. Average intron erage intron length (measured as the average distance length also shows a temporal pattern contrary to that of xa in base pairs between the exons of a given gene), inter- except that no clear differences between stages are found genic distance (average distance in base pairs between after embryonic development. The intergenic distance shows two adjacent genes), expression bias (a measure of how a temporal pattern similar to that of xa, except that this dis- evenly distributed the expression of a gene is over time), tance is low in the earliest stages in which xa is high. The recombination rate (based in observed cross-overs in average expression bias shows a temporal pattern similar to 100-kb intervals, from Comeron et al. [2012]), the GC that of xa except for an overall increase over time. The same Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1473 Coronado-Zamora et al. GBE FIG. 2.—Clusters of temporal profiles of expression of embryonic development genes and their xa, a, xna, and x estimated using DFE-alpha. (A) Temporal expression profile for all the genes belonging to each cluster. (B) x for each cluster. (C) xa sampling for each cluster. (D) xna for each cluster. (E) a for each cluster. Each point in the plots in (B)–(E) is calculated for a randomly drawn sample of the set of genes in each cluster with replacement (100 bootstrap replicates per cluster). Asterisks represent the significance by a permutation test, the color indicates whether the value was higher (red) or lower (blue) than expected (•, 0.1–0.05, *<0.05, **<0.01, and ***<0.001). Number of genes analyzed in supplementary table S1, Supplementary Material online. Permutation P values are shown in supplementary table S12, Supplementary Material online. 1474 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE FIG. 3.—Nine genomic features over developmental stages. Lines and stages as in figure 1. (A) Size is the coding sequences length of a gene in base pairs. (B) Number of exons is the number of exons for the genes expressed in a stage. (C) Number of transcripts is the number of different isoforms of each gene expressed in a stage. (D) Fop is a measure of codon usage bias: the ratio of optimal codons to synonymous codons. (E) Average intron length is the average distance, in bases, between the exons of a gene. (F) The expression bias is a measure of how much the expression of a gene is restricted to one or few stages estimated as equation (1) (see Materials and Methods). (G) Recombination rate is estimated in 100-kb nonoverlapping windows. (H) CG content of each gene. (I) PPIs are estimated as the number of PPIs of each gene. Mean sampling distribution was obtained by resampling 100 times with replacement the genes from each stage. See supplementary table S8, Supplementary Material online, for the genes analyzed. The same patterns are found when using 4- fold data, see supplementary figure S19, Supplementary Material online. Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1475 Coronado-Zamora et al. GBE correlations are found when we used 4-fold degenerated sites Discussion as a proxy for the mutation rate (see supplementary table S17 and fig. S19, Supplementary Material online). A similar pat- Three main conclusions can be derived from our analyses. tern is found when we analyze the genomic features in the First, the rate of adaptive substitution (xa) measured along gene expression clusters for the embryo development and the life-cycle of D. melanogaster reveals two peak periods: life-cycle (see supplementary fig. S20, Supplementary one encompassing the four initial hours of the embryonic Material online, and P values in supplementary table S18 development and one encompassing from the L3 larval stage and fig. S21, Supplementary Material online, and P values in onward. Drosophila melanogaster, as all holometabolous supplementary table S19, Supplementary Material online, insects, has an indirect development with two active free- respectively). roaming phases, the larva and the adult, and two inactive sessile developmental phases, the embryo and the pupa. Analysis of Maternal, Maternal–Zygotic, and Zygotic Genes The larval and adult phenotypes, especially their morphology, To further explore why xa and xna were high in the ear- arise primarily through the genetic, cellular and tissue inter- liest developmental stages, we analyzed separately mater- actions of embryonic and pupal (metamorphosis) develop- nal, maternal–zygotic, and zygotic genes. For that ment, respectively. Therefore, adaptation in the larva and purpose, we used a microarray study (Thomsen et al. the adult should be reflected, not only in the substitution rates 2010) that categorized developmental genes as maternal, of the genes expressed in the larva and adult but also in those zygotic, and maternal–zygotic by determining which gene expressed during embryonic development (for the larva) and transcripts are already present in the egg and which ones pupal development (for the adult). The observation that genes are not. Maternal genes are defined as genes whose expressed in mid-and late-embryonic development show rel- mRNA is laid within the egg by the mother and are never atively lower rates of nonsynonymous substitution than the transcribed by the embryo. The embryo contains, thus, genes expressed in the larval and pupal stages suggests that mRNAs coming from two different genomes, that of the adaptation has occurred preferentially in the adult rather than mother and that of the embryo. Maternal–zygotic genes in the larva. In contrast to a previous report using a smaller are genes whose mRNA is laid in the egg by the mother gene set (Artieri et al. 2009), we do not find that adults have but that are also transcribed by the embryo. Zygotic genes higher rates of nonadaptive substitutions, neither in x nor in are genes whose mRNA is not laid in the egg by the xna, than the pupa. In a previous study, it was found that the mother. We analyzed the genes reported in Thomsen 150 genes with the highest number of nonsynonymous sub- et al. (2010) for each category with a permutation test, stitutions (this is the highest x) are more intensively expressed testing if the genes in each category undergo differential in the larva and pupa than in the embryo and that their high- selection compared with the genes in the other two cat- est level of expression is in the male adults (Davis et al. 2005). egories (see Materials and Methods for details). We found An important novelty of our analysis is the incorporation of that xa was not significantly different between categories polymorphism data, which allows us to more precisely distin- when assessed with a permutation test (fig. 4). This guish between the rates of nonsynonymous substitutions implies that the large xa of the earliest stages is not due (with conservation indicated by low x), adaptive nonsynon- to any specific gene category. Consistent with the hypoth- ymous substitution (measured by xa), nonadaptive nonsynon- esis of lower efficiency of natural selection on maternal ymous substitution (measured by xna), and the proportion of genes both x (P ¼ 0.024) and xna (P ¼ 0.003) were higher fixed adaptive substitutions (a). In a previous study (Salvador- for maternal genes than for zygotic genes (and interme- Martınez et al. 2018), we have estimated the same selection diate for the maternal–zygotic genes). Zygotic genes show statistics over the embryonic anatomy of D. melanogaster, in lower values of xa and xna than expected from the per- this study, we do it over developmental time. This has allowed mutation test (P ¼ 0.035 and P ¼ 0.036, respectively). us to infer that the relatively lower level of conservation in Supplementary table S20, Supplementary Material online, early pupal and male stages is due to adaptive nonsynony- contains permutation P values. Finally, this analysis was mous substitutions and not to nonadaptive nonsynonymous repeated but using the gene that are in common with substitutions. In the earliest developmental stages, instead, the genes expressed in the first hours according to we can infer that the lower sequence conservation is due, modENCODE to check whether maternal genes account mostly, to nonadaptive substitutions. Our results are, thus, for the high xna in the first stages of the development compatible with those of Liu and Robinson-Rechavi (2018a) (fig. 1). Very similar results are obtained, thus, indicating but the methods we use and the inclusion of polymorphism that the high xna values in the earliest stages are due to data, however, allow us to infer the nonadaptive nature of the the maternal genes in these stages (supplementary fig. nonsynonymous substitutions in the earliest stages of S22, Supplementary Material online). Even so, the xa in development. these earliest stages is still larger than in mid-and late- Our analysis shows that the latter is due, mostly, to mater- development. nal genes. Previous studies have already pointed out that 1476 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE Table 2 Spearman’s Correlations between xa and Genetic Features Genomic Feature Relation with xa Correlation (r2) for Females (P Value) Correlation (r2) for Males (P Value) 6 Average intron length Negative 0.802 (1.12  10 ) 0.808 (1.08  106) Size Negative 0.731 (1.56  106) 0.764 (1.39  106) Number of exons Negative 0.862 (6.53  107) 0.886 (4.82  107) Number of transcripts Negative 0.874 (5.70  107) 0.870 (5.94  107) Fop Negative 0.759 (1.42  106) 0.688 (1.71  106) Expression bias Positive 0.508 (4.89  105) 0.552 (1.58  105) Recombination Positive 0.330 (2.07  103) 0.334 (1.91  103) Intergenic distance N.S. 0.043 (0.299) 0.082 (0.148) CG content Negative 0.905 (3.647  107) 0.772 (1.34  106) PPI Negative 0.152 (0.045) 0.075 (0.165) NOTE.—Spearman’s correlations performed between each stage’s average xa and the average of each genomic feature in each stage. Females and males are separated because their gene expression is measured separately in the last three stages in the modENCODE. Thus, the data considered for females and males are the same for all the stages but the last three. Methodological details about how each genomic feature was calculated are explained in the Materials and Methods. Fop: frequency of optimum codons. selection in maternal genes is less efficient (although these probability because, for any given protein, the codon changes studies do not relate that to the lower conservation of early that improve a protein function would often be different from development [Cruickshank and Wade 2008; Wade et al. the codon changes associated with more efficient codon us- 2009]). This is because, in females, the alleles in the loci of age (Hershberg and Petrov 2008; Presnyak et al. 2015), espe- maternal genes can affect the fitness of the offspring but that cially for highly expressed genes (Pal et al. 2001). The is not the case for males (because males lay no eggs and thus relationship between codon bias and GC content we observe, no maternal mRNAs in their offspring). As a result, for mater- that is, identical temporal trends over development (see nal genes, selection cannot act as effectively on males as on fig. 3), is explained by many optimal codons ending in G or females. Natural selection is then less effective and many non- C in D. melanogaster (Bierne and Eyre-Walker 2006). adaptive variants cannot be eliminated from the population Similarly, it has previously been reported that genes with a (leading to the higher xna). large expression bias tend to be less conserved (low x, The high rate of nonsynonymous nonadaptive substitu- Larracuente et al. 2008) and we found, that in addition, tions in the earliest developmental stages parallels recent find- they also show higher rates of adaptive substitutions. On av- ings in developmental genetics: the genes expressed in the erage, if a gene is expressed in a very specific time window, it earliest developmental stages were found to be different is likely involved in regulating a smaller number of develop- within Diptera, whereas the zygotic genes expressed right mental interactions than if it is expressed in through many after, the gap genes, were found to be the same in all the stages. Changes in such gene are then less likely to interfere Diptera species analyzed so far (Wotton et al. 2015). with many different developmental processes and, thus, these A second conclusion of our study is that the temporal pat- changes are less likely to be deleterious. Similarly, it has been tern of the selection statistics mirrors that of the genomic shown that genes with a high connectivity, as measured with features analyzed. Thus, mid- and late-embryonic develop- the total number of PPIs, tend to be under selective constraint ment express genes that have, on average, more exons, because they are involved in complex functions (Valdar and more different transcripts, a more optimal codon usage, Thornton 2001; Caffrey et al. 2004). We do indeed find that larger introns, and larger gene size. Some of these genomic our selection statistics correlate with PPI, although the tem- features were previously found to correlate with x, and in poral pattern of PPI is not similar to the temporal pattern of some few cases with xa. However, their distribution over de- any of the selection statistics. velopmental stages and their relationship with the von Baer’s Gene size, measured as the protein’s primary structure law or the hourglass model have not been analyzed before, length, has also been shown to negatively correlate with x except for intron length (Liu and Robinson-Rechavi 2018b) (Comeron et al. 1999; Duret and Mouchiroud 1999; and gene length (Liu and Robinson-Rechavi 2018a; Yang Larracuente et al. 2008) and even with a measure of the et al. 2018). This latter study also reports how other measures rate of adaptive substitutions (Liu and Robinson-Rechavi of gene complexity (such as the number of protein domains) 2018a). A possible explanation for this observation could be are higher among developmentally expressed genes. the fact that for a given adaptive mutation, longer genes The negative correlation between codon bias and adaptive would have more adaptive segregating sites competing nonsynonymous substitution at the protein level was already against each other in different haplotypes. This would pro- known (Hershberg and Petrov 2008; Presnyak et al. 2015). duce a kind of intergenic or interexonic Hill–Robertson inter- Such negative correlation should be expected by mere ference (Hill and Robertson 1966). Consistent with this Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1477 Coronado-Zamora et al. GBE FIG. 4.—Selection statistics (xa, x, xna, and a) for maternal, maternal–zygotic and zygotic genes. Maternal genes are those genes whose mRNA are laid by the mother in the egg and are never zygotically transcribed, maternal–zygotic are those genes whose mRNA is present in the egg but that are also transcribed by the zygote. Zygotic genes are genes whose mRNA is not laid in the egg by the mother. (A) xa is not statistically different between these gene categories. (B) x is significantly higher in maternal genes than the other two gene categories (P ¼ 0.024). (C) xna is significantly higher in maternal genes than in the other two gene categories (P value ¼ 0.003). (D) a is marginally lower in maternal genes compared with the other two categories. Each point in a plot (100 bootstrap replicates per group) is calculated for a randomly drawn sample of the set of genes in each gene category. Asterisks represent the significance by a permutation test, the color indicates whether the value was higher (red) or lower (blue) than expected (•, 0.1–0.05, *<0.05, **<0.01, and ***<0.001). The number of genes analyzed is shown in supplementary table S6, Supplementary Material online. P values in supplementary table S20, Supplementary Material online. previous work, we found that gene size correlates negatively features over the life-cycle, or if, on the contrary, these geno- with both x and xa over developmental time. mic features are a consequence of differential adaptation over To our knowledge, the number of exons and transcripts the life-cycle. To establish what leads to what, one would themselves have never been directly associated with neither x have to devise experiments, or at least some form of model- nor xa, as we found in our study. As one would expect, ing, that is beyond the state of the art for fly’s development however, the number of exons correlates with gene size and its underlying genetic bases. There are, however, some (data not shown). In a similar way, the larger the number of intrinsic characteristics of development that can be used to exons, the larger the number of different transcripts that can suggest that the first possibility is more likely. be produced by reshuffling exons by alternative splicing. A possible explanation for that would be that, as suggested Then, the relationship between the number of exons and from a more qualitative evo-devo perspective (Gellon and the number of transcripts with the selection statistics could McGinnis 1998; Kennison 1993), embryonically expressed be either direct or indirect through the correlation of these genes have a more complex regulation than postembryoni- genomic features with each other and gene size. cally expressed genes. Although the former are expressed in It is not possible to establish whether the temporal pattern wider areas of the embryo, their expression changes more in of adaptation is a consequence of differences in genomic time and space than that of postembryonically expressed 1478 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE genes (Salvador-Martınez and Salazar-Ciudad 2015). This rather weak, because there are no major differences in x more complex regulation may require a more complex ge- between embryonic stages after the 8-h, except for genes netic structure, such as manifested by those genes having in cluster 8. During the first 2 h, x is significantly high (hours more exons, more transcripts, larger genes and larger introns 0–2 h: P ¼ 0.032), but from 6–8 h to 22–24 h x is lower than (Kennison 1993; Gellon and McGinnis 1998). The larger av- expected based on the permutation test (P < 0.001). This is erage intron length of developmental genes may also be a also the case for xa (supplementary table S10, Supplementary reflection of complex regulation, but at the level of cis-regu- Material online). In contrast with some previous studies latory elements, because cis-regulatory elements can be lo- (Kalinka et al. 2010; Levin et al. 2016), we do not find that cated within introns too. Larger expression areas, less the latest stages of embryonic development are less conserved temporally restricted expression and more complex gene (this previous study measures only dN and dS). However, structure, may also reflect that mid- and late-development genes, whose expression is high only in late embryonic devel- genes interact with more other genes than genes expressed opment (cluster 8, fig. 2), show a significant high x and mar- later. This would make them, in rough terms, more pleiotropic ginally significant high xa. These genes are only a small and, thus, less likely to adapt. This explanation would be the proportion of the genes expressed in the last embryonic de- reflection at the genomic level of the idea that there are wider velopmental stages and, thus, have a minor effect on our and higher levels of interdependence among body parts in calculations of xa, x, and a of these stages (thus likely explain- middevelopment than in late and early development ing the differences between our study and Kalinka et al. (Duboule 1994; Raff 1996). From this perspective, genomic [2010]). The difference in dN between mid- and late- features would not be a consequence of differential adapta- development stages in Kalinka et al. (2010) is however rather tion over the life-cycle. Instead, the pattern of differential ad- subtle too. Overall, thus, our results are compatible with aptation over the life-cycle would be a consequence of how Kalinka et al. (2010). The hourglass model was proposed on genomic features have to be over development. the basis of what was understood of D. melanogaster and The above arguments would not apply to the earliest em- vertebrate’s development (Slack et al. 1993; Duboule 1994; bryonic stages because of the low efficiency of natural selec- Raff 1996). The life-cycles of the fly and the mouse are quite tion on maternal genes. In addition, it has been suggested different. Mice, as all amniotes, are direct developers, mean- that the small intron and gene size of the genes expressed in ing that development gives rise to a juvenile and later, grad- the first hours of D. melanogaster development (Anderson ually, to an adult. Flies are indirect developers in which 1973; Heyn et al. 2014) would be imposed by the very fast embryonic development gives rise to a free-roaming larva cell divisions occurring in those hours. Because cells divide very and that, by a rather abrupt process of metamorphosis, to rapidly in early D. melanogaster development, there is no time an adult. If the hourglass model is understood for the whole to transcribe, splice, and translate long genes and genes with of the life-cycle then our results are roughly consistent with it long introns. This would preclude the earliest expressed genes at the genetic level: genes expressed during embryonic devel- from having the complex regulatory gene structure (no long opment are highly conserved, except for the genes expressed introns and no long genes) that in the later-expressed genes in the earliest stages, whereas the genes expressed later, from correlates with lower rates of nonsynonymous substitutions. the larval stage L3 onward, show less conservation and more This hypothesis is consistent with the reasoning of the previ- adaptation. On the other hand, this temporal hourglass pat- ous paragraph because, in fact, there are not many genetic tern can also be understood as development generally obey- interactions, signaling or cell movement during these early ing von Baer’s law, but departing from it in the earliest stages. fast cell-division states. In other words, no complex genetic We hypothesize that this departure would arise, in one hand, regulation occurs until the early fast division stages have fin- from the lower efficiency of selection on maternal genes and, ished. This applies not only to the fly but also to many other in the other hand, as a consequence of the reduced gene animal groups, including vertebrates (O’Farrell et al. 2004; structure complexity required for fast nuclei divisions in early Heyn et al. 2014; Siefert et al. 2015). The hourglass model development. has also been suggested for vertebrates and it is then tempt- ing to suggest that, in them, the relative lack of conservation Supplementary Material of early development would also be related to the simpler gene structure required for short cell cycles and, possibly, to Supplementary data are available at Genome Biology and maternal genes. This simpler gene structure would be consis- Evolution online. tent with our results of shorter genes with fewer exons, shorter introns, and less transcripts in the genes expressed in the earliest stages of development. Acknowledgments A third main conclusion of our study is that our results are We thank Roland Zimm, Miquel Marin-Riera, Pascal Hagolani, consistent, roughly, with the hourglass model but not with Rishi Das Roy, and Miguel Brun-Usan for comments. We the von Baer’s law. However, the fit to the hourglass model is thank Carla Giner for helping in developing R scripts. This Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1479 Coronado-Zamora et al. GBE work was supported by the Finnish Academy to IS-C and the Darwin C. 1872. The origins of species by means of natural selection. Spanish Ministerio de Economıa y Empresa [Project MINECO Philadelphia (PA): John Wanamaker. Davis JC, Brandman O, Petrov DA. 2005. Protein evolution in the context CGL2017-89160-P awarded to AB], the Generalitat de of Drosophila development. J Mol Evol. 60(6):774–785. Catalunya [grant 2017SGR 1379], and the Secretaria Dezso Z, et al. 2008. A comprehensive functional analysis of tissue specif- d’Universitats i Recerca del Departament d’Empresa i icity of human gene expression. BMC Biol. 6:49. Coneixement de la Generalitat de Catalunya [FI-DGR 2015 Drost H-G, Gabel A, Grosse I, Quint M. 2015. Evidence for active mainte- to MC-Z]. nance of phylotranscriptomic hourglass patterns in animal and plant embryogenesis. Mol Biol Evol. 32(5):1221–1231. Duboule D. 1994. Temporal colinearity and the phylotypic progression: a Literature Cited basis for the stability of a vertebrate Bauplan and the evolution of Alberch P. 1980. Ontogenesis and morphological diversification. Am Zool. morphologies through heterochrony. Development 1994:135–142. 20(4):653–667. Duret L, Mouchiroud D. 1999. Expression pattern and, surprisingly, gene Anderson DT. 1973. Embryology and phylogeny in annelids and arthro- length shape codon usage in Caenorhabditis, Drosophila, and pods. Oxford, New York: Pergamon Press. Arabidopsis. Proc Natl Acad Sci U S A. 96(8):4482–4487. Arthur W. 1977. The origin of animal body plans: a study in evolutionary Durinck S, et al. 2005. BioMart and Bioconductor: a powerful link between developmental biology. Cambridge: Cambridge University Press. biological databases and microarray data analysis. Bioinformatics Artieri CG, Haerty W, Singh RS. 2009. Ontogeny and phylogeny: molec- 21(16):3439–3440. ular signatures of selection, constraint, and temporal pleiotropy in the Eyre-Walker A. 2002. Changing effective population size and the development of Drosophila. BMC Biol. 7(1):42. McDonald–Kreitman test. Genetics 162(4):2017–2024.  Avila V, Campos JL, Charlesworth B. 2015. The effects of sex-biased gene Eyre-Walker A. 2006. The genomic rate of adaptive evolution. Trends Ecol expression and X-linkage on rates of adaptive protein sequence evo- Evol (Amst). 21(10):569–575. lution in Drosophila. Biol Lett. 11:20150117. Eyre-Walker A, Keightley PD. 2007. The distribution of fitness effects of Bierne N, Eyre-Walker A. 2004. The genomic rate of adaptive amino acid new mutations. Nat Rev Genet. 8(8):610–618. substitution in Drosophila. Mol Biol Evol. 21(7):1350–1360. Eyre-Walker A, Keightley PD. 2009. Estimating the rate of adaptive mo- Bierne N, Eyre-Walker A. 2006. Variation in synonymous codon use and lecular evolution in the presence of slightly deleterious mutations and DNA polymorphism within the Drosophila genome. J Evol Biol. population size change. Mol Biol Evol. 26(9):2097–2108. 19(1):1–11. Fay JC, Wyckoff GJ, Wu C-I. 2001. Positive and negative selection on the Boyko AR, et al. 2008. Assessing the evolutionary impact of amino acid human genome. Genetics. 158(3):1227–1234. mutations in the human genome. PLoS Genet. 4(5):e1000083. Feldman B. 2005. Relative importance and value. Unpublished manuscript Bustamante CD, et al. 2002. The cost of inbreeding in Arabidopsis. Nature (Version 1.1, March 19, 2005). Available from: http//www.gtcenter. 416(6880):531–534. org/Archive/Conf05/Downloads/Conf/Feldman60.pdf (last accessed Bustamante CD, et al. 2005. Natural selection on protein-coding genes in January 2019). the human genome. Nature 437(7062):1153–1157. Futschik M. 2015. Mfuzz: soft clustering of time series gene expression Caffrey DR, Somaroo S, Hughes JD, Mintseris J, Huang ES. 2004. Are data. R package version 2.36.0. Available from: http://www.sysbiolab. protein–protein interfaces more conserved in sequence than the rest eu/software/R/Mfuzz/index.html (last accessed December 2015). of the protein surface? Protein Sci. 13(1):190–202. Futschik ME, Carlisle B. 2005. Noise-robust soft clustering of gene expres- Cagan A, et al. 2016. Natural selection in the great apes. Mol Biol Evol. sion time-course data. J Bioinform Comput Biol. 3(4):965–988. 33(12):3268–3283. Galtier N. 2016. Adaptive protein evolution in animals and the effective Carneiro M, et al. 2012. Evidence for widespread positive and purifying population size hypothesis. PLoS Genet. 12(1):e1005774. selection across the European rabbit (Oryctolagus cuniculus) genome. Gelbart WM, Emmert DB. 2013. FlyBase high throughput expression pat- Mol Biol Evol. 29(7):1837–1849. tern data. FlyBase Analysis. Available from: http://flybase.org/reports/ Castellano D, Coronado-Zamora M, Campos JL, Barbadilla A, Eyre-Walker FBrf0221009.html (last accessed May 2014). A. 2016. Adaptive evolution is substantially impeded by Hill–Robertson Gellon G, McGinnis W. 1998. Shaping animal body plans in development interference in Drosophila. Mol Biol Evol. 33(2):442–455. and evolution by modulation of Hox expression patterns. BioEssays Castellano D, James J, Eyre-Walker A. 2018. Nearly neutral evolution 20(2):116–125. across the Drosophila melanogaster genome. Mol Biol Evol. Gossmann TI, Keightley PD, Eyre-Walker A. 2012. The effect of variation in 35:2685–2694. the effective population size on the rate of adaptive molecular evolu- Charlesworth J, Eyre-Walker A. 2008. The McDonald–Kreitman test and tion in eukaryotes. Genome Biol Evol. 4(5):658–667. slightly deleterious mutations. Mol Biol Evol. 25(6):1007–1015. Gould SJ. 1977. Ontogeny and phylogeny. Cambridge, Massachusetts:The Clark AG, et al. 2007. Evolution of genes and genomes on the Drosophila Belknap Press of Harvard University Press. phylogeny. Nature 450(7167):203–218. Graveley BR, et al. 2011. The developmental transcriptome of Drosophila Comeron JM, Kreitman M, Aguad e M. 1999. Natural selection on synon- melanogaster. Nature 471(7339):473–479. ymous sites is correlated with gene length and recombination in Grömping U. 2006. Relative importance for linear regression in R: the Drosophila. Genetics 151(1):239–249. package relaimpo. J Stat Softw. 17:1–27. Comeron JM, Ratnappan R, Bailin S. 2012. The many landscapes Guillen Y, Casillas S, Ruiz A. 2019. Genome-wide patterns of sequence of recombination in Drosophila melanogaster. PLoS Genet. divergence of protein-coding genes between Drosophila buzzatii and 8(10):e1002905. D. mojavensis. J Hered. 110(1):92–101. Cruickshank T, Wade MJ. 2008. Microevolutionary support for a develop- Halligan DL, Keightley PD. 2006. Ubiquitous selective constraints in the mental hourglass: gene expression patterns shape sequence variation Drosophila genome revealed by a genome-wide interspecies compar- and divergence in Drosophila. Evol Dev. 10(5):583–590. ison. Genome Res. 16(7):875–884. Cutter AD. 2008. Divergence times in Caenorhabditis and Drosophila in- Hebenstreit D, et al. 2011. RNA sequencing reveals two major ferred from direct estimates of the neutral mutation rate. Mol Biol classes of gene expression levels in metazoan cells. Mol Syst Biol. Evol. 25(4):778–786. 7(1):497. 1480 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 Adaptation and Conservation throughout D. melanogaster GBE Hernandez RD, Williamson SH, Bustamante CD. 2007. Context depen- Mensch J, Serra F, Lavagnino NJ, Dopazo H, Hasson E. 2013. Positive dence, ancestral misidentification, and spurious signatures of natural selection in nucleoporins challenges constraints on early expressed selection. Mol Biol Evol. 24(8):1792–1800. genes in Drosophila development. Genome Biol Evol. Hershberg R, Petrov DA. 2008. Selection on codon bias. Annu Rev Genet. 5(11):2231–2241. 42(1):287–299. Messer PW, Petrov DA. 2013. Frequent adaptation and the McDonald– Heyn P, et al. 2014. The earliest transcribed zygotic genes are short, newly Kreitman test. Proc Natl Acad Sci U S A. 110(21):8615–8620. evolved, and different across species. Cell Rep. 6(2):285–292. Murali T, et al. 2011. DroID 2011: a comprehensive, integrated resource Hill WG, Robertson A. 1966. The effect of linkage on limits to artificial for protein, transcription factor, RNA and gene interactions for selection. Genet Res. 8(3):269–294. Drosophila. Nucleic Acids Res. 39(Suppl_1):D736–D743. Hu H, et al. 2017. Constrained vertebrate evolution by pleiotropic genes. Murray GGR, et al. 2017. Natural selection shaped the rise and fall of Nat Ecol Evol. 1(11):1722–1730. passenger pigeon genomic diversity. Science 358(6365):951–954. Huang W, et al. 2014. Natural variation in genome architecture among Obbard DJ, Welch JJ, Kim K-W, Jiggins FM. 2009. Quantifying adaptive 205 Drosophila melanogaster Genetic Reference Panel lines. Genome evolution in the Drosophila immune system. PLoS Genet. Res. 24(7):1193–1208. 5(10):e1000698. Hurst LD. 2002. The Ka/Ks ratio: diagnosing the form of sequence evolu- O’Farrell PH, Stumpff J, Su TT. 2004. Embryonic cleavage cycles: how is a tion. Trends Genet. 18(9):486–487. mouse like a fly? Curr Biol. 14:R35–R45. Irie N, Kuratani S. 2011. Comparative transcriptome analysis reveals ver- Pal C, et al. 2001. Highly expressed genes in yeast evolve slowly. Genetics tebrate phylotypic period during organogenesis. Nat Commun. 2:248. 158:927–931. Irie N, Kuratani S. 2014. The developmental hourglass model: a predictor Peden JF. 1999. Analysis of codon usage. Nottingham (United Kingdom): of the basic body plan? Development 141(24):4649–4655. University of Nottingham. James JE, Piganeau G, Eyre-Walker A. 2016. The rate of adaptive evolution Piasecka B, Lichocki P, Moretti S, Bergmann S, Robinson-Rechavi M. 2013. in animal mitochondria. Mol Ecol. 25(1):67–78. The hourglass and the early conservation models—co-existing patterns Jukes TH, Cantor CR. 1969. Evolution of protein molecules. In: Munro, HN, of developmental constraints in vertebrates. PLoS Genet. editor. Mammalian protein metabolism. New York: Academic Press. 9(4):e1003476. p. 21–32. Poe S, Wake MH. 2004. Quantitative tests of general models for the evo- Kalinka AT, Tomancak P. 2012. The evolution of early animal embryos: lution of development. Am Nat. 164(3):415–422. conservation or divergence? Trends Ecol Evol. 27(7):385–393. Presnyak V, et al. 2015. Codon optimality is a major determinant of mRNA Kalinka AT, et al. 2010. Gene expression divergence recapitulates the de- stability. Cell 160(6):1111–1124. velopmental hourglass model. Nature 468(7325):811–814. Pröschel M, Zhang Z, Parsch J. 2006. Widespread adaptive evolution Keightley PD, Eyre-Walker A. 2007. Joint inference of the distribution of of Drosophila genes with sex-biased expression. Genetics fitness effects of deleterious mutations and population demography 174(2):893–900. based on nucleotide polymorphism frequencies. Genetics R Core Team. 2015. R: A language and environment for statistical com- 177(4):2251–2261. puting. R Foundation for Statistical Computing, Vienna, Austria. Keightley PD, Eyre-Walker A. 2012. Estimating the rate of adaptive mo- Available from: http://www.R-project.org/ (last accessed December lecular evolution when the evolutionary divergence between species is 2015). small. J Mol Evol. 74(1-2):61–68. Raff RA. 1996. The shape of life: genes, development, and the evolution of Kennison JA. 1993. Transcriptional activation of Drosophila homeotic animal form. Chicago (IL): University of Chicago Press. genes form distant regulatory elements. Trends Genet. 9(3):75–79. R amia M, Librado P, Casillas S, Rozas J, Barbadilla A. 2012. PopDrowser: Kim S. 2015. ppcor: an R package for a fast calculation to semi-partial the population Drosophila browser. Bioinformatics 28:595–596. correlation coefficients. Commun Stat Appl Methods. 22:665–674. Richardson MK, et al. 1997. There is no highly conserved embryonic stage Kousathanas A, Halligan DL, Keightley PD. 2014. Faster-X adaptive protein in the vertebrates: implications for current theories of evolution and evolution in house mice. Genetics 196(4):1131–1143. development. Anat Embryol. 196(2):91–106. Langley CH, et al. 2012. Genomic variation in natural populations of Ronen R, Udpa N, Halperin E, Bafna V. 2013. Learning natural selection Drosophila melanogaster. Genetics 192(2):533–598. from the site frequency spectrum. Genetics 195(1):181–193. Larracuente AM, et al. 2008. Evolution of protein-coding genes in Roux J, Robinson-Rechavi M. 2008. Developmental constraints on verte- Drosophila. Trends Genet. 24(3):114–123. brate genome evolution. PLoS Genet. 4(12):e1000311. Levin M, et al. 2016. The mid-developmental transition and the evolution Salazar-Ciudad I. 2010. Morphological evolution and embryonic develop- of animal body plans. Nature 531(7596):637–641. mental diversity in metazoa. Development 137(4):531–539. Liu J, Robinson-Rechavi M. 2018a. Adaptive evolution of animal proteins Salvador-Martınez I, Coronado-Zamora M, Castellano D, Barbadilla A, over development: support for the Darwin selection opportunity hy- Salazar-Ciudad I. 2018. Mapping selection within Drosophila mela- pothesis of evo-devo. Mol Biol Evol. 35:2862–2872. nogaster embryo’s anatomy. Mol Biol Evol. 35(1):66–79. Liu J, Robinson-Rechavi M. 2018b. Developmental constraints on genome Salvador-Martınez I, Salazar-Ciudad I. 2015. How complexity increases in evolution in four bilaterian model species. Genome Biol Evol. development: an analysis of the spatial-temporal dynamics of 1218 10:2266–2277. genes in Drosophila melanogaster. Dev Biol. 405(2):328–339. Mackay TFC, et al. 2012. The Drosophila melanogaster Genetic Reference Sander K. 1983. The evolution of patterning mechanisms: gleanings from Panel. Nature 482(7384):173–178. insect embryogenesis and spermatogenesis. In: Goodwin BB, Holder Marais G, Nouvellet P, Keightley PD, Charlesworth B. 2005. Intron size and N, Wylie CC, editors. Development and evolution. New York: exon evolution in Drosophila. Genetics 170(1):481–485. Cambridge University Press. p. 137–159. McDonald JH, Kreitman M. 1991. Adaptive protein evolution at the Adh Santpere G, et al. 2015. Analysis of five gene sets in chimpanzees suggests locus in Drosophila. Nature 351(6328):652–654. decoupling between the action of selection on protein-coding and on Medawar PB. 1954. The significance of inductive relationships in the de- noncoding elements. Genome Biol Evol. 7(6):1490–1505. velopment of vertebrates. Br J Philos Sci. 2:172–174. Sawyer SA, Kulathinal RJ, Bustamante CD, Hartl DL. 2003. Bayesian anal- Meisel RP, Connallon T. 2013. The faster-X effect: integrating theory and ysis suggests that most amino acid replacements in Drosophila are data. Trends Genet. 29(9):537–544. driven by positive selection. J Mol Evol. 57(Suppl 1):S154–S164. Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019 1481 Coronado-Zamora et al. GBE Siefert JC, Clowdus EA, Sansam CL. 2015. Cell cycle control in the early € von Baer KE. 1828. Uber Entwicklungsgeschichte der Thiere. Beobachtung embryonic development of aquatic animal species. Comp Biochem und Reflexion. Königsberg (Russia): Bei den Gebrüdern Borntr€ager. Physiol C Toxicol Pharmacol. 178:8–15. Wade MJ, Priest NK, Cruickshank TE. 2009. A theoretical overview of Slack JMW, Holland PWH, Graham CF. 1993. The zootype and the phy- genetic maternal effects. In: Maternal effects in mammals. Chicago lotypic stage. Nature 361(6412):490–492. (IL): University of Chicago Press. p. 38–63. Smith NGC, Eyre-Walker A. 2002. Adaptive protein evolution in Wagner GP, Kin K, Lynch VJ. 2013. A model based criterion for gene Drosophila. Nature 415(6875):1022–1024. expression calls using RNA-seq data. Theory Biosci. 132(3):159–164. Steige KA, Laenen B, Reimegård J, Scofield DG, Slotte T. 2017. Genomic Wilt FH, Hake S. 2004. Principles of developmental biology. New York: analysis reveals major determinants of cis-regulatory variation in W.W. Norton. Capsella grandiflora. Proc Natl Acad Sci U S A. 114(5):1087–1092. Wotton KR, et al. 2015. Quantitative system drift compensates for altered Strasburg JL, et al. 2011. Effective population size is positively correlated maternal inputs to the gap gene network of the scuttle fly Megaselia with levels of adaptive divergence among annual sunflowers. Mol Biol abdita. eLife 4:e04785. Evol. 28(5):1569–1580. Yanai I, et al. 2005. Genome-wide midrange transcription profiles reveal Tamura K, Subramanian S, Kumar S. 2004. Temporal patterns of fruit fly expression level relationships in human tissue specification. (Drosophila) evolution revealed by mutation clocks. Mol Biol Evol. Bioinformatics 21(5):650–659. 21(1):36–44. Yang D, et al. 2018. A two-level model for the role of complex and young Tataru P, Mollion M, Gl emin S, Bataillon T. 2017. Inference of distribution genes in the formation of organism complexity and new insights into of fitness effects and proportion of adaptive substitutions from poly- the relationship between evolution and development. Evodevo 9:22. morphism data. Genetics 207(3):1103–1119. Yang Z, Nielsen R. 1998. Synonymous and nonsynonymous rate variation Thomsen S, Anders S, Janga SC, Huber W, Alonso CR. 2010. Genome- in nuclear genes of mammals. J Mol Evol. 46(4):409–418. wide analysis of mRNA decay patterns during early Drosophila devel- Zhang J, Nielsen R, Yang Z. 2005. Evaluation of an improved branch-site opment. Genome Biol. 11(9):R93. likelihood method for detecting positive selection at the molecular Valdar WS, Thornton JM. 2001. Protein–protein interfaces: analysis of level. Mol Biol Evol. 22(12):2472–2479. amino acid conservation in homodimers. Proteins 42(1):108–124. Associate editor: Maria Costantini 1482 Genome Biol. Evol. 11(5):1463–1482 doi:10.1093/gbe/evz086 Advance Access publication April 26, 2019