1. Computational and Systems Biology
  2. Genomics and Evolutionary Biology
Download icon

Evidence for a common evolutionary rate in metazoan transcriptional networks

Short Report
Cited
4
Views
2,235
Comments
0
Cite as: eLife 2015;4:e11615 doi: 10.7554/eLife.11615

Abstract

Genome sequences diverge more rapidly in mammals than in other animal lineages, such as birds or insects. However, the effect of this rapid divergence on transcriptional evolution remains unclear. Recent reports have indicated a faster divergence of transcription factor binding in mammals than in insects, but others found the reverse for mRNA expression. Here, we show that these conflicting interpretations resulted from differing methodologies. We performed an integrated analysis of transcriptional network evolution by examining mRNA expression, transcription factor binding and cis-regulatory motifs across >25 animal species, including mammals, birds and insects. Strikingly, we found that transcriptional networks evolve at a common rate across the three animal lineages. Furthermore, differences in rates of genome divergence were greatly reduced when restricting comparisons to chromatin-accessible sequences. The evolution of transcription is thus decoupled from the global rate of genome sequence evolution, suggesting that a small fraction of the genome regulates transcription.

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

eLife digest

The genetic information that makes each individual unique is encoded in DNA molecules. Cells read this molecular instruction manual by a process called transcription, in which proteins called transcription factors bind to DNA in specific places and regulate which sections of the DNA will be expressed. These 'transcripts' are active molecules that determine the cell’s – and ultimately the individual’s – characteristics. However, it is not well understood how alterations in the DNA of different individuals or species can lead to changes in where the transcription factors bind, and in which transcripts are expressed.

Carvunis, Wang, Skola et al. set out to determine if there is a relationship between how often DNA changes and how often transcription changes during the evolution of animals. The experiments examined the abundance of transcripts in the cells of a variety of animal species with close or distant evolutionary relationships. For example, the house mouse was compared to a close relative called the Algerian mouse, to another species of rodent (rat) and to humans.

The experiments show that the changes in transcript abundances are happening at similar rates in mammals, birds and insects, even though DNA changes at very different rates in these groups of animals. This similarity was also observed for other aspects of transcription, such as in changes to where transcription factors bind to DNA. The next challenges are to find out what makes transcription evolve at such similar rates in these groups of animals, and whether these findings extend to other species and to other processes in cells.

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

Introduction

A long-standing question in biology is what fraction of the genome regulates transcription (ENCODE Project Consortium, 2012; Graur et al., 2013; Niu and Jiang, 2013; Kellis et al., 2014). Recent studies of chromatin structure have implicated half of the human genome in regulatory interactions (ENCODE Project Consortium, 2012). Comparative genomic studies, however, have shown that less than 10% of the human genome is evolutionarily conserved (Siepel et al., 2005), suggesting that many of the experimentally-detected interactions are not functional (Graur et al., 2013). Recent studies have measured the association between sequence changes and changes in transcript levels, epigenetic modifications or binding of transcription factors regulating specific gene sets (gene-specific transcription factors, GSTF) (Cookson et al., 2009; McVicker et al., 2013; Kasowski et al., 2010; 2013; Heinz et al., 2013; Villar et al., 2014; Wong et al., 2015; Brem et al., 2002; Chan et al., 2009; Shibata et al., 2012). These experiments demonstrated that genomic sequences can influence transcription even in the absence of evolutionary conservation. For instance, some repetitive elements previously thought to be 'junk' DNA have been shown to effectively regulate gene expression (Rebollo et al., 2012). The rapid evolution of repetitive and other rapidly-evolving sequences could cause pervasive rewiring of transcriptional networks through creation and destruction of regulatory motifs (Villar et al., 2014). Such rapid transcriptional evolution would set mammals apart from other metazoans like birds or insects, whose genomes contain far fewer repetitive elements (Taft et al., 2007) and tend to be more constrained (Siepel et al., 2005; Zhang et al., 2014).

A few studies have attempted to assess whether transcriptional networks evolve more rapidly in mammals than in insects from the fruit fly genus Drosophila. These studies have reached conflicting conclusions. When examining the evolution of GSTF binding, chromatin immuno-precipitation (ChIP) studies in mammalian livers have generally described faster divergence rates than similar studies in fly embryos (Villar et al., 2014; Stefflova et al., 2013). However, divergence rates were estimated with different analytical methods in the different ChIP studies (Supplementary file 1) (Villar et al., 2014; Bardet et al., 2012). Another study found that gene expression levels may diverge at a slower rate in mammals than in flies, by comparing genome-wide correlations of mRNA abundances estimated by RNA sequencing (RNA-seq) for mammals but by a mixture of technologies for flies including microarrays (Coolon et al., 2014). Although the inconsistencies between these conclusions may indicate that the evolution of transcriptional networks is fundamentally different in mammals and insects, they may also reflect a sensitivity of evolutionary rate estimations to technical methodology.

Here, we jointly examined the evolution of gene expression levels and the underlying genome-wide changes in GSTF binding and cis-regulatory sequences using consistent methodologies both within and across various animal lineages.

Results

We assembled a comparative genomics platform encompassing >40 publicly available datasets spanning >25 organisms representative of the Mammalia (mammals), Aves (birds) and Insecta (insects) phylogenetic classes (Figure 1—figure supplement 1). We designed a statistical framework to objectively compare the rates of divergence of these various datasets across lineages. In brief, an exponential model describing evolutionary divergence under a common, lineage-naïve rate was evaluated against a lineage-aware model, accounting for both statistical significance and effect size (Figure 1; Materials and methods). We assessed the power of this statistical framework using simulations and found that it could detect differences in divergence rates with high sensitivity (Materials and methods; Figure 1—figure supplement 2).

As a baseline, we first performed a comparative analysis of the evolution of genome sequences. We randomly sampled genomic segments from designated reference genomes: Mus musculus domesticus (C57BL/6) for mammals, Gallus gallus for birds and Drosophila melanogaster for insects. The rates at which genomic segments that retained homologs with the other species within each lineage accumulate nucleotide substitutions were then estimated and compared using our statistical framework. Segments retaining homologs displayed high sequence conservation across all three lineages, although our framework detected a slightly but significantly faster divergence in insects than in mammals or birds (P<0.05; Figure 2—figure supplement 1). Next, we compared the rates at which randomly sampled genomic segments lost homology with the other species within each lineage. We observed a much larger difference in evolutionary rates across lineages using this measure (P<0.05; Figure 2; Figure 2—figure supplement 2). For instance, after 100 million years (Myrs) of evolution, only ~30% of mammalian segments retained homology, whereas >60% of bird and insect segments did. These findings recapitulated previous observations according to which genome sequences are less constrained in mammals than in insects (Siepel et al., 2005) or birds (Zhang et al., 2014).

Figure 1 with 2 supplements see all
Statistical framework to evaluate differences in evolutionary rates of change.

Throughout this study we frequently evaluated whether the rate of evolutionary divergence of a given layer of transcriptional regulation differs between lineages. Our approach is equivalent to asking: if the lineage labels were hidden, would one be able to tell that the data points correspond to several lineages or would they seem equally likely to belong to a common distribution? (a, b) Depict an example of statistically indistinguishable evolutionary rates. Without lineage labels (a), the similarity data are modeled by an exponential decay as well as with lineage labels (b). Adding lineage labels does not significantly improve the fit. (c, d) Depict an example of statistically different evolutionary rates. Adding lineage labels (d) significantly improves the fit of an exponential decay model over unlabeled data (c).

https://doi.org/10.7554/eLife.11615.003
Figure 2 with 2 supplements see all
Genomic sequences evolve more rapidly in mammals than in birds and insects.

The evolutionary retention of 5000 randomly sampled 75 bp segments was averaged over 20 trials. Organisms compared to reference species are as follows: M. musculus domesticus (AJ), M. musculus castaneus, M. spretus, rat, guinea pig, rabbit, human, chimpanzee and dog for Mammalia; turkey, zebrafinch and flycatcher for Aves; D. simulans, D. erecta, D. yakuba, D. ananassae, D. pseudoobscura, D. virilis, D. willistoni and D. grimshawi for Insecta. Colored dashed lines: lineage-specific exponential fits, here and in all following displays. The trends were robust to variations in segment length and sequence similarity filters (Figure 2—figure supplement 2).

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

We then studied the evolution of gene expression levels, using exclusively RNA-seq datasets. In mammals and birds, these datasets were generated from adult livers; in insects, they were from whole bodies of adult female fruit flies (Materials and methods; Figure 3—source data 1). After determining expression levels for each gene in each species using a common data processing pipeline, we correlated the expression levels of genes in the reference species with the expression levels of their one-to-one orthologs in all other species within the same lineage (Materials and methods). We found that correlations of gene expression levels decreased over time at similar rates that were statistically indistinguishable: a lineage-naïve model describing the evolution of gene expression levels under a common rate fitted the data as well as a lineage-aware model (Figure 3). This result was robust to changes in correlation metrics or inclusion/exclusion of poorly expressed genes (Figure 3—figure supplement 1).

Figure 3 with 1 supplement see all
Gene expression levels diverge at a common rate in mammals, birds and insects.

Gene expression levels were derived independently from two RNA-seq experiments for each reference species and then correlated against each other and against gene expression levels derived from individual experiments in other species within the same lineage. Black dashed line: lineage-naïve exponential fit of all the data, without differentiating the lineages, here and in all following displays. Organisms compared to reference species are as follows: M. musculus castaneus, M. spretus, rat, human and gorilla for Mammalia; turkey, duck and flycatcher for Aves; D. simulans, D. yakuba, D. ananassae and D. pseudoobscura for Insecta.

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

Several lines of evidence suggest that gene expression levels can remain relatively stable even as the genomic locations bound by GSTFs change rapidly over time (Wong et al., 2015; Chan et al., 2009; Paris et al., 2013). Therefore, we next examined the evolution of GSTF binding patterns. We considered all GSTFs that were profiled using ChIP followed by massively parallel sequencing (ChIP-seq) in at least three related species, where separate ChIPs were performed per species. GSTFs meeting these requirements were Twist and Giant in fruit fly embryos, and CEBPA, FOXA1 and HNF4A in mammalian livers (Materials and methods; Figure 4—source data 1; Supplementary file 1). We aimed to measure cross-species similarity in GSTF occupancy with a unified analytical method across all of these datasets. Despite the widespread use of ChIP-seq, there is no consensus on the appropriate analytical method (Wilbanks and Facciotti, 2010). ChIP-seq analysis pipelines typically discretize continuous occupancy profiles into a set of occupied segments ('peaks'), but this step requires choosing a signal processing algorithm (a peak caller) and associated parameters (Figure 4a). Further comparison of occupied segments across species requires additional analytical choices (Figure 4a), some of which can strongly influence downstream findings (Bardet et al., 2012).

Figure 4 with 1 supplement see all
GSTF occupancy diverges at a common rate in mammals and insects.

(a) Estimating shared GSTF occupancy across species requires multiple parameter choices. This diagram summarizes the main steps involved in comparing GSTF-occupied segments across species, showing a representative sample of choices at each step (steps represented by purple shapes, specific choices by the first letter bolded). The detailed methods and specific choices illustrated here and implemented in panels b–d are described in Materials and methods. (b, c) An example of different analytical choices leading to different results despite starting from the same underlying data. Organisms compared to reference species are as follows: M. musculus domesticus (AJ), M. musculus castaneus, M. spretus, rat, human and dog for Mammalia; D. simulans, D. erecta, D. yakuba, D. ananassae and D. pseudoobscura for Insecta. (d) Most combinations of choices yield indistinguishable evolutionary rates of GSTF binding patterns across lineages. The comparison of Twist and CEBPA is enlarged to show the color labels corresponding to the statistical interpretation regarding relative evolutionary rates. (e) A genome-wide comparison of GSTF occupancy profiles at single-nucleotide resolution shows indistinguishable evolutionary rates for CEBPA, HNF4A and FOXA1 in mammals, and for Twist and Giant in insects. PCC: Pearson correlation coefficient. (f) CTCF occupancy is highly conserved in mammals. Transparent points and lines are identical as in panel e. Hexagons correspond to cross-species correlations of CTCF occupancy at single-nucleotide resolution.

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

To explore the impact of these choices, we processed all ChIP-seq data using systematic combinations of parameters representative of, and expanding from, previous studies (Supplementary file 1) (Landt et al., 2012). In total, we executed 108 analytical pipelines to compare divergence rates across 6 pairs of GSTFs (2 in insects each compared with 3 in mammals), the occupancy profiles of which were examined in 3–7 species per lineage (Materials and methods). The values of the estimated rates varied greatly from one combination of parameters to the next (Figure 4b,c). However, in the majority of cases (56–78% over the 6 comparisons), GSTF binding patterns diverged at statistically indistinguishable rates in mammals and insects (Figure 4d; Figure 4—source data 2). Although the computed divergence rates were sensitive to technical methodology (Figure 4—figure supplement 1), for a given method the results were generally similar across lineages for all of the five GSTFs investigated.

To substantiate these findings, we devised a method to compare genome-wide occupancy profiles at single-nucleotide resolution without discretization. We correlated occupancy profiles between pairs of species across all nucleotides where genomes aligned, after accounting for the differences in sequencing depth, read length and fragment size across datasets (Materials and methods). Again, we found indistinguishable divergence rates, regardless of which GSTF or lineage was examined (Figure 4e). After 100 Myrs of evolution, the correlation of GSTF occupancy profiles was 0.10 in mammals and 0.13 in insects. As a control, we also applied this method to CTCF, a pleiotropic DNA-binding protein that acts as chromatin insulator and looping factor (Ohlsson et al., 2010). In mammals, patterns of DNA occupancy have been shown to be more conserved for CTCF than for GSTFs using unified analytical methods (Schmidt et al., 2012). In contrast, CTCF DNA occupancy was shown to diverge rapidly in insects, perhaps due to the existence of other insulator proteins (Villar et al., 2014; Ni et al., 2012). Our analysis successfully recapitulated this difference (Figure 4f), demonstrating that the common evolutionary rate observed among GSTFs (Figure 4e) was not an artifact of our method for profile correlation.

The similarity of divergence rates observed across lineages for gene expression levels (Figure 3) and GSTF binding patterns (Figure 4) was unexpected given the rapid evolution of genomic sequences in mammals relative to insects (Siepel et al., 2005) or birds (Zhang et al., 2014) (Figure 2). We therefore further examined these trends at the level of cis-regulatory sequences. First, we considered the DNA sequence motifs thought to be specifically recognized by the mammalian and insect GSTFs included in the previous ChIP-seq analysis (Figure 4). We identified locations with significant matches to these motifs throughout the genomes of the reference species and estimated how frequently these loci retained the same motifs relative to background expectations (Materials and methods). We found similar, indistinguishable retention rates in mammals and insects (Figure 5a). Next, we studied the evolution of a broader set of motifs corresponding to GSTFs shared between M. musculus and D. melanogaster. We found that these motifs were retained at similar rates across lineages relative to background expectations in 8 out of 12 cases (one example shown in Figure 5b; all other cases in Figure 5—figure supplement 1).

Figure 5 with 2 supplements see all
Regulatory sequences diverge at similar rates across lineages.

(a) The motifs for CEBPA, HNF4A and FOXA1 in mammals and for Twist and Giant in insects are retained at a common rate. Organisms compared to reference species are the same as Figure 4. (b) The motifs for GSTFs shared in mammals and insects are retained at common rates. One example is shown here for the motifs corresponding to PHO (FBgn0002521) in D. melanogaster and YY1 (ENSMUSG00000021264) in M. musculus, which are orthologous GSTFs. Eleven other cases of motif evolution for shared GSTFs conserved in mammals and insects are shown in Figure 5—figure supplement 1. Organisms compared to reference species are as in Figure 4. (cd) Chromatin-accessible sequences are retained at similar rates in mammals, birds and insects. Analyses were performed as in Figure 2, limiting sampling to the inaccessible (c) and accessible (d) portions of the intergenic regions. Organisms compared to reference species are the same as Figure 2. The trends were robust to variations in segment length and sequence similarity filters (Figure 5—figure supplement 2).

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

Most active cis-regulatory sequences are located in genomic regions with accessible chromatin (Hesselberth et al., 2009). A recent study showed that chromatin-accessible sequences were significantly more conserved between human and mouse than expected by chance (Yue et al., 2014). We expanded this analysis to a wide range of species by using chromatin-accessible sequences identified by DNAse I hypersensitivity in M. musculus livers, D. melanogaster embryos and G. gallus MSB-1 cells (Materials and methods). We performed the segment sampling procedure described previously (Figure 2), after excluding genes and promoter regions since they typically are highly conserved (Materials and methods). Whereas inaccessible segments lost homology much faster in mammals than in insects and birds (P<0.05; Figure 5c), accessible segments retained homologs at more similar rates in the three lineages (Figure 5d; Figure 5—figure supplement 2). We still detected statistically significant differences across lineages (P<0.05), but the effect sizes were considerably smaller than for inaccessible segments. For instance, ~60% of segments retained homology after 100 Myrs in birds and insects, independently of accessibility, whereas ~50% of chromatin-accessible segments and only ~20% of inaccessible segments did so in mammals.

Discussion

To our knowledge, the analyses presented here represent the most comprehensive study conducted to date on the evolution of transcriptional networks across animal lineages. By applying unified analytical methods to data from different lineages, we were able to glean novel insights into the evolution of transcription in animals. We observed that gene expression levels, GSTF binding patterns, regulatory motifs and chromatin-accessible sequences each diverged at rates that were similar across mammals, birds and insects. These unexpected results reconcile previously conflicting findings (Villar et al., 2014; Coolon et al., 2014), highlighting the importance of unified study methodologies and providing evidence for a common evolutionary rate in metazoan transcriptional networks.

Most functional genomics studies have focused on humans and model organisms such as D. melanogaster or M. musculus, which are distantly related to each other. However, data on closely related species, like those which we collected in this study, are needed to investigate the dynamics of molecular network evolution. Unfortunately, such data remain scarce, leading to important limitations of our work. We only investigated three lineages and six to twelve organisms per lineage with non-uniform coverage over evolutionary time. In addition, we only examined a small number of tissues for each lineage and a total of five GSTFs (none in birds). The generalizability of our observations thus remains to be further evaluated as more data becomes available. Despite these limitations, our finding that transcriptional networks evolve at a common rate per year across animal lineages was strikingly robust across data layers.

The underlying mechanisms responsible for this concordance of evolutionary rates are unclear. Mammals, birds and insects exhibit wide differences in the features that are traditionally associated with evolutionary rates, such as generation times and breeding sizes. Populations with small breeding sizes, such as mammals, are thought to be more prone to genetic drift (Ohta, 1992). This theory accounts for the abundance of repetitive elements and the rapid evolution of genomic sequences in mammals relative to insects, which have much larger breeding sizes. If the same theoretical principles also governed the evolution of transcriptional networks, we would have expected that transcription would evolve more rapidly in mammals than in insects. Instead, our results show that the evolution of transcriptional networks, whether slow (e.g. transcript levels) or fast (e.g. GSTF binding), is decoupled from the lineage-specific features that govern genome sequence evolution.

One potential model could be that repetitive and rapidly-evolving sequences, which make up the majority of the mammalian genome (Siepel et al., 2005; Taft et al., 2007), play a negligible role in the global regulation of gene expression. Rather, chromatin-accessible regions may represent the only portion of the mammalian genome that effectively regulates transcription. We observed that chromatin-accessible regions diverge much more slowly than other non-coding sequences in mammals, consistent with previous findings (Yue et al., 2014). These differences in divergence rates, however, were not found in birds and insects. As a result, chromatin-accessible regions in mammals are conserved at levels similar to those in birds and insects, in contrast to the genome as a whole. According to this model, the similar rates of evolution of chromatin-accessible sequences would constrain the dynamics of transcriptional evolution to be similar across lineages. The regulatory potential of repetitive and other rapidly-evolving elements could be rendered functionally inconsequential by silencing, or could be concentrated on controlling the expression of genetic elements that we did not investigate, such as non-coding RNAs or species-specific genes (Sundaram et al., 2014).

An alternative model could be that the sequences that control transcriptional regulation in birds and insects evolve particularly rapidly within otherwise stable genomes. In these organisms, transcriptional networks would diverge under the action of natural selection, through specific single nucleotide substitutions resulting in rapid compensatory turnover (He et al., 2011a). In mammals, transcriptional networks would diverge in a largely neutral fashion driven for instance by transposable elements (Sundaram et al., 2014). In this case, similar rates of transcriptional divergence across lineages would arise through very different evolutionary processes.

Importantly, none of the aforementioned models account for the differences in generation times between lineages. Evolutionary changes occurring based on chronological time and not generation time has also been observed for many protein-coding sequences. Observations such as these led to the molecular clock theory (Kumar, 2005). The mechanisms through which environmental forces entrain these chronological evolutionary clocks remain to be elucidated (Kumar, 2005).

Materials and methods

Genome and annotation sources

We downloaded genome sequences for organisms belonging to three metazoan lineages: mammals, birds and insects. The mammalian and insect genome sequences were downloaded from the UCSC Genome Bioinformatics website (Rosenbloom et al., 2015): mm9 for Mus musculus domesticus, rn5 for Rattus norvegicus and hg19 for Homo sapiens; dm3 for Drosophila melanogaster, droSim1 for Drosophila simulans, droEre2 for Drosophila erecta, droYak2 for Drosophila yakuba, droAna3 for Drosophila ananassae and dp4 for Drosophila pseudoobscura. Genomes for mice strains and species not available from the UCSC Genome Bioinformatics site (M. musculus domesticus [AJ], M. musculus castaneus and M. spretus) were downloaded from (Stefflova et al., 2013). We downloaded bird genome sequences from Ensembl version 80 BioMart (Cunningham et al., 2015): galGal4 for Gallus gallus, Turkey_2.01 for Meleagris gallopavo, taeGut3.2.4 for Taeniopygia guttata and FicAlb_1.4 for Ficedula albicollis. Protein-coding gene names and symbols along with associated transcripts sequences were obtained from FlyBase (dos Santos et al., 2015) for insect species (dmel-r5.46, dsim-r1.4, dere-r1.3, dyak-r1.3, dana-r1.3 and dpse-r2.30), from Ensembl version 80 BioMart for bird species and from Ensembl version 59 BioMart for mammalian species (Cunningham et al., 2015). For M. spretus and M. musculus castaneus, we used the same transcript annotations as for M. musculus. Within the genomes of our designated reference organisms (M. musculus domesticus, G. gallus and D. melanogaster), we defined promoters as the region 0-2 kb upstream of transcription start site and delineated intergenic regions as regions that did not overlap annotated genes or promoters. Chromatin accessibility tracks used in Figure 5c,d and Figure 5—figure supplement 2 were downloaded from the UCSC bioinformatics website (Rosenbloom et al., 2015) for M. musculus domesticus and D. melanogaster, and obtained from (He et al., 2014) for G. gallus. We restricted our analyses to the sequences or annotations in, or homologous to, the well-defined chromosome scaffolds of the reference organism. Specific reference chromosomes analyzed are as follows: G. gallus (1–28, Z, W), D. melanogaster (2L, 2R, 3L, 3R, 4, X) and M. musculus (1–19, X, Y).

Homology and evolutionary relationships

We obtained orthology relationships between protein-coding genes using Ensembl COMPARA (Vilella et al., 2009), matching the Ensembl versions used for protein coding genes for each species described above. These relationships were used in Figure 3, Figure 3—figure supplement 1, Figure 5b and Figure 5—figure supplement 1. Homology between genomic segments was assigned using the LiftOver tool (Rosenbloom et al., 2015), for all analyses presented in Figures 2, 4 and 5 and associated figure supplements, with the exception of the nucleotide-resolution analysis of GSTF occupancy profiles presented in Figure 4e,f. We used pre-computed chain files from UCSC matching the genome versions listed above when chains were readily available (Rosenbloom et al., 2015). When chain files were not available, we built chain files to map the UCSC M. musculus C57BL/6 mm9 to the genomes of M. musculus domesticus AJ, Mus musculus castaneus and Mus spretus, as well as to map the Ensembl 80 galGal4 to the genomes of M. gallopavo, F. albicollis and T. guttata (Figure 1—figure supplement 1). These chains were constructed by following the steps recommended by UCSC (Supplementary file 2) (http://genomewiki.ucsc.edu/index.php/Whole_genome_alignment_howto).

For the nucleotide-resolution analysis of GSTF occupancy profiles, we assigned homology relationships using the chain files, or, in the case of mice strains, using genome mapping tables from (Stefflova et al., 2013). We filtered the chain files to obtain one-to-one unambiguous mappings by retaining only highest scoring alignment for each position. These filtered mappings were then used to transfer data to from any organism onto the corresponding reference genome. Regions in the reference species genome lacking one-to-one unambiguous mappings were excluded from analysis.

To define evolutionary distances separating species in Myrs, we chose published estimates generated as homogenously as possible within each lineage using a combination of sequence alignments and fossil records. All distances between insect species were taken from (Tamura et al., 2004); all distances between bird species were taken from (Lu et al., 2015); distances between mammalian species were taken from (Stefflova et al., 2013) and TimeTree (Hedges, 2009).

Data sources

For RNA-seq analyses (Figure 3; Figure 3—figure supplement 1), sequencing data for the reference species corresponding to two experiments performed independently by different research groups, and, when possible, representing different genotypes, were downloaded from public repositories. For M. musculus domesticus, we used data from (Goncalves et al., 2012; Sugathan and Waxman, 2013), for G. gallus we used data from (Brawand et al., 2011) and (Coble et al., 2014), for D. melanogaster we used data from (ENCODE Project Consortium, 2012; Chen et al., 2014). Other species included were M. musculus castaneus (Goncalves et al., 2012), M. spretus (Wong et al., 2015), R. norvegicus (Gong et al., 2014), H. sapiens (ENCODE Project Consortium, 2012; Lin et al., 2014), G. gorilla (Brawand et al., 2011), D. simulans (Chen et al., 2014), D. yakuba (Chen et al., 2014), D. ananassae (Chen et al., 2014), D. pseudoobscura (Chen et al., 2014), M. gallopavo (Monson et al., 2014), A. platyrhynchos (Huang et al., 2013) and F. albicollis (Uebbing et al., 2013). Specific accession numbers are listed in Figure 3—source data 1.

For ChIP-seq analyses (Figure 4), we downloaded data for FOXA1 in M. musculus domesticus (C57BL/6) (Stefflova et al., 2013), M. musculus domesticus (AJ) (Stefflova et al., 2013), M. musculus castaneus (Stefflova et al., 2013), M. spretus (Stefflova et al., 2013) and R. norvegicus (Stefflova et al., 2013); HNF4A and CEBPA in M. musculus domesticus (C57BL/6) (Stefflova et al., 2013), M. musculus domesticus (AJ) (Stefflova et al., 2013), M. musculus castaneus (Stefflova et al., 2013), M. spretus (Stefflova et al., 2013), R. norvegicus (Stefflova et al., 2013), H. sapiens (Schmidt et al., 2010) and C. familiaris (Schmidt et al., 2010); Twist in D. melanogaster (He et al., 2011b), D. simulans (He et al., 2011b), D. erecta (He et al., 2011b), D. yakuba (He et al., 2011b), D. ananassae (He et al., 2011b) and D. pseudoobscura (He et al., 2011b); Giant in D. melanogaster (Paris et al., 2013; Bradley et al., 2010), D. yakuba (Bradley et al., 2010) and D. pseudoobscura (Paris et al., 2013). We also gathered data for CTCF in M. musculus domesticus (C57BL/6) (Schmidt et al., 2012), R. norvegicus (Schmidt et al., 2012), H. sapiens (Schmidt et al., 2012), C. familiaris (Schmidt et al., 2012), D. melanogaster (Ni et al., 2012), D. simulans (Ni et al., 2012), D. yakuba (Ni et al., 2012) and D. pseudoobscura (Ni et al., 2012). Accession numbers corresponding to the specific experimental replicates and control samples are listed in Figure 4—source data 1.

For motif analyses (Figure 5a,b; Figure 5—figure supplement 1), we gathered known position-weight matrixes from the JASPAR database (Mathelier et al., 2014) and the Fly Factor survey (Zhu et al., 2011). We focused on the motifs corresponding to Twist and Giant in D. melanogaster, to CEBPA, HNF4A and FOXA1 in M. musculus domesticus, and on a set of 12 other motifs corresponding to GSTFs conserved across mammals and insects. This set was constructed by downloading all Core A vertebrata motifs from JASPAR (Mathelier et al., 2014), identifying those corresponding to conserved GSTFs with one-to-one orthologs between M. musculus domesticus and D. melanogaster using COMPARA (Vilella et al., 2009), and filtering the list down to those 12 instances where a position-weight matrix was also described in Fly Factor (Zhu et al., 2011) and were not already analyzed.

Comparing evolutionary rates

We developed a statistical framework to compare evolutionary rates between lineages, and implemented it in R (Development Core Team, 2011). This framework takes as inputs: measures of pairwise cross-species similarity (e.g. correlation of gene expression or sequence conservation), pairwise cross-species evolutionary distances and lineage labels. Conceptually, the framework estimates both a statistical significance and an effect size to determine whether rates of evolutionary divergence are indistinguishable or different between lineages (Figure 1).

In practice, we model evolutionary divergence by an exponential decay in log-linear space. First, the nls function in R is applied to the log-transformed cross-species similarity data as a function of evolutionary distances to derive the following linear models:

  • a lineage-naïve model that estimates a shared intercept and slope for all the data without specifying the lineage labels

  • a lineage-aware model that estimates a shared intercept for all the data and lineage-specific slopes based on lineage labels

  • lineage-specific models that estimate intercept and slope individually for each lineage

Second, an R function written in-house to handle nls model structures estimates the significance level of an ANOVA with a likelihood ratio test comparing the lineage-naïve and the lineage-aware model. Third, we define the effect size as the predicted absolute difference in similarity between lineage pairs after 100 Myrs of divergence as estimated from the lineage-specific models. We consider that the framework detected a difference between evolutionary divergence rates when the significance level is <0.05 and the effect size is >5%.

We chose to use an exponential decay function because it is the simplest evolutionary model that fit all our input measures of cross-species similarity reasonably well. We chose to model the exponential decay in log-linear space because we noted that a simple exponential decay in linear space failed to capture the conservation observed between distant species (mouse versus human at 91 Myrs and dog at 97.4 Myrs) when analyzing the evolutionary dynamics of GSTF binding (Figure 4) and motif retention (Figure 5). We hypothesize that these data layers likely follow a more complex decay model, but we did not want to explore this with our current data set to avoid over-fitting.

The power of our statistical framework was assessed by simulating data for two lineages with measure of cross-species similarity decaying exponentially at different rates over time (Figure 1—figure supplement 2). We fixed one lineage to decay at set rates: −0.007, −0.005 and −0.003. We fixed the second lineage to be faster by a range of given differences. Over 1000 simulations, we sampled two values from a normal distribution centered on the expected values from the set exponential decay rates corresponding to the evolutionary distances shown in Figure 4b, with standard deviations set at 0.5% or 5%. Our framework detected an absolute rate difference of 0.001 in 39.3% of simulations and an absolute rate difference of 0.003 in 88.9% of simulations when the standard deviation was high (5%). When the standard deviation was low (0.5%), our framework detected an absolute rate difference of 0.001 in 25.7% of simulations and an absolute rate difference of 0.003 in 100% of simulations.

Gene expression evolutionary rates (related to Figure 3)

Analysis of gene expression evolutionary rates was performed in four steps. First, we preprocessed the raw RNA sequencing data downloaded for public data sources. Second, we quantified the abundance of all annotated transcripts corresponding to protein-coding genes. Third, we estimated cross-species similarity by correlating transcript abundances at the genome-scale. Finally, we used these cross-species similarity estimates as input to our statistical framework to evaluate a common model against a lineage-aware model.

RNA sequencing data were first preprocessed using FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Trimmomatic (Bolger et al., 2014). In order to quantify transcript abundances, we then used the program Sailfish (Patro et al., 2014) (1) to build transcriptome indices for each species using the transcriptome sequences described above, using the parameters '-p 8 -k 20'; and (2) to quantify transcript abundance using the transcriptome indices with the parameters '-p 8 -l "T=PE:O=><:S=U"' for samples with paired-end reads and '-p 8 –l "T=SE:S=U”' for samples with single-end reads. The bias-corrected transcripts per million (TPM) abundances estimated by Sailfish were then summed over the transcripts corresponding to the same gene locus.

To estimate cross-species similarities in gene expression levels, for each lineage, we used R (Development Core Team, 2011) to build a matrix containing the gene expression values for all the protein-coding genes of the reference organism and their one-to-one orthologs across other organisms within each lineage. We discarded instances where the abundance of a particular gene locus was less than or equal to 5 TPM. We then calculated the Spearman’s rank correlation for the expression of all genes between the reference and all other organisms within each lineage and plotted these correlations as against the evolutionary distance separating each organism pair (Figure 3). We also repeated the calculations using Kendall’s rank correlation coefficient and Pearson’s product-moment correlation on log2-transformed expression values (Figure 3—figure supplement 1a,b). Finally, we calculated Spearman’s correlations among all genes including those with less than 5 TPM (Figure 3—figure supplement 1c). All these scenarios were evaluated using our statistical framework. None indicated that a lineage-aware model described the data better than a common model.

GSTF occupancy – segment-resolution (related to Figure 4a–d)

The first step of all our occupancy analyses was to align the ChIP-seq reads to the corresponding genomes in order to obtain occupancy profiles (Figure 4a). For each accession (Figure 4—source data 1), the sequencing reads were aligned to reference genomes using Bowtie2 version 2.2.4 (Langmead and Salzberg, 2012) with the parameters '-very-sensitive -N 1.' Reads containing the 'XS:' field (multi-mappers) were removed. Reads having the same start site were presumed to be PCR duplicates and removed using the 'rmdup' command of SAMtools version 1.1 (Li et al., 2009). The filtered reads were then converted to tagAlign format. The tagAlign files corresponding to CEBPA, HNF4A, FOXA1, Twist and Giant were then processed using 108 different segment-resolution methods and one nucleotide-resolution method; the tagAlign files corresponding to CTCF were only processed using the nucleotide-resolution method. The nucleotide-resolution method is described below and relates to Figure 4e,f.

The aim of our segment-resolution analyses was to examine how robust the evolution of GSTF binding patterns was across 108 different analysis pipelines (Figure 4a–d). We implemented all these pipelines, which follow the same general framework and differ only in the choice of 5 parameters, described and underlined below.

First, the occupancy profiles in the tagAlign files were discretized into candidate occupied segments using a peak caller algorithm that aims at identifying segments where the ChIP sample is enriched in reads relative to the control sample. We implemented two peak callers: MACS version 2 (M) (Zhang et al., 2008) and SPP (S) (Kharchenko et al., 2008).

The occupied segments were then selected from the candidate set using a quality filter: stringent (S), lenient (L) or asymmetric (A). When using MACS2 (Zhang et al., 2008) as a peak caller, lenient segments were called using a p-value cutoff of 10−5 (default) and merged across replicates when available using the merge function in BEDTools (Quinlan and Hall, 2010). Stringent segments were called using a p-value cutoff of 10−22 and intersected across replicates when replicates were available. The intersection procedure, inspired from (Stefflova et al., 2013), used BEDTools (Quinlan and Hall, 2010) to implement the following two steps: (1) merge the two replicates and (2) select the merged segments corresponding to at least one segment in each original replicate. When using SPP (Kharchenko et al., 2008) as a peak caller, lenient segments were called using a q-value of 10−2 (default), and merged across replicates when available (Quinlan and Hall, 2010). Stringent segments were called by selecting all candidate segments assigned to the lowest possible q-value in the sample, then intersected across replicates when available using the same intersection procedure. The asymmetric quality filter, inspired by (Bardet et al., 2012; He et al., 2011b), indicates that segments were called stringently in the reference species and leniently in the other organism.

The coordinates of the occupied segments called in the reference organism were projected onto the other organism’s genome using the LiftOver tool from the UCSC genome browser (Rosenbloom et al., 2015) and specifying a sequence similarity filter through the minMatch parameter. We used 3 different minMatch thresholds: stringent (S: 0.95 default), lenient (L: 0.5), and none (N: 0.001).

After cross-species coordinate projection, a reference subset was chosen to define the set of reference-occupied segments that would be further analyzed. Three choices were implemented: all reference-occupied segments independently of whether they map to any other species (A); for each pair of species, only reference-occupied segments with a homolog in the second species (P); only reference-occupied segments that had homologs across all the other species considered within the lineage (S).

The projected coordinates of the reference subset were then overlapped with the coordinates of the occupied segments in the other species using the intersect function in BEDTools (Quinlan and Hall, 2010). The overlap requirement was either lenient (L; default parameter of 1 bp) or stringent (S; required a reciprocal overlap of half of the segments length: '-f 0.5 -r').

We systematically executed all combinations of the aforementioned 2 peak callers, 3 quality filters, 3 sequence similarity filters, 3 reference subsets, and 2 overlap requirements, yielding a total of 108 pipelines. The output of each pipeline was the fraction of reference subset segments that overlapped segments occupied in the others species (i.e. segments retaining occupancy between the two species). This output was used as a cross-species similarity measure for GSTF binding patterns. We analyzed these similarity measures for 6 pairs of GSTFs (Twist and Giant were each compared to FOXA1, CEBPA and HNF4A) using our statistical framework. Two GSTFs were considered to diverge differently from each other over time when 1) the significance of the test was less than 0.05 and 2) the effect size was greater than 5%. In summary we found that the choice of parameters greatly influenced what the evolutionary dynamics of a given GSTF looked like (Figure 4b,c) but that in general the rate of divergence of mammal and insects GSTFs were statistically indistinguishable (Figure 4d). The results of these tests for all GSTF pairs considered across 108 pipelines are reported in Figure 4—source data 2 and summarized as pie-charts in Figure 4. Observations about general trends of parameters and evolutionary divergence are further elaborated in Figure 4—figure supplement 1.

As a control we also conducted an analysis between FOXA1 and CEBPA since FOXA1 lacks data past 20 Myrs of evolutionary divergence, whereas for all others GSTFs we have broader coverage across the 100 Myrs range. We applied the same statistical framework to the within-lineage comparison between FOXA1 and CEBPA and detected that FOXA1 evolves faster than CEBPA in 74/108 instances. We believe that most of these detected differences are artifacts because the conservation of binding patterns for FOXA1 and CEBPA is in fact highly correlated throughout all combinations of parameters when restricting analyses to data points up to 20 Myrs (Pearson’s r = 0.96). We suspect that this type of artifact also affects the results of comparing FOXA1 with Twist or Giant (Figure 4d).

GSTF occupancy – nucleotide-resolution (related to Figure 4e,f)

In order to compare occupancy profiles directly without discretizing them into occupied segments and unoccupied segments, we correlated sets of imputed fragment density vectors across species. The inputs to this method were the tagAlign files described above. To generate these vectors we first estimated the mean fragment size using a method adapted from (Kharchenko et al., 2008), whereby the mean fragment size is computed as the number of base pairs of offset between the positive and negative strands that maximizes the Pearson’s correlation coefficient of their mapped read density. We used a modified approach that considered only the density of 5' read start sites on each strand, rather than the density of the entire read. The first peak of the cross-correlation values was identified by approximating the first derivative by the finite difference method, smoothing the derivative values with a Gaussian kernel of bandwidth 10, and identifying the first downward zero-crossing of the curve. This position was used as the estimated mean fragment size L. We created imputed fragments by extending each read start site by L base pairs in the 3' direction. We then calculated a fragment density vector for each chromosome as the number of such imputed fragments that overlap each genomic position. When multiple replicates were available, replicates were merged by adding the fragment density vectors.

In order to minimize bias introduced by the presence of unmappable regions, we implemented a masking scheme that adaptively normalizes each dataset depending on the read length and estimated fragment size of each sequencing run. First, all possible error-free reads of a given length were generated synthetically and aligned back to the genome using Bowtie2 2.2.4 with the following parameters: '-r -N 0 -D 0 -R 0 --dpad 0 --score-min “C,0,-1”'. Any multi-mapping reads with the ‘XS:’ flag were removed and the 5’ and 3’-most positions of the remaining read alignments recorded. The imputed fragment densities computed from the ChIP data were then normalized by dividing the density at each position by the fraction of positions within L base pairs upstream that were covered by the start site (5’ for positive-strand density and 3’ for negative-strand density) of a uniquely-mapped genomic read. Positions with 0 uniquely-mappable read start sites within L base pairs upstream were excluded from further analysis.

In order to compare between species, we transferred data from query organisms to the reference genome using the one-to-one filtered chain files described previously, and calculated the Pearson’s correlation between the concatenated chromosome vectors of reference and reference-mapped query data. The evolution of the correlation was modeled and compared using the statistical framework described above.

Genome sequence evolutionary rates (related to Figure 2 and Figure 5c,d)

We calculated the percentage of randomly sampled segments retaining homology. Within the genomes of the reference species, we delineated the boundaries of the regions from which to sample: whole genome (Figure 2; Figure 2—figure supplement 1), intergenic regions in accessible chromatin and intergenic regions in inaccessible chromatin (Figure 5; Figure 5—figure supplement 2). We used the BEDTools shuffle command (Quinlan and Hall, 2010) to randomize the locations of 5000 segments of 75 bp length within the delineated boundaries using the option '-noOverlapping'. The resulting 5000 shuffled segments were then mapped across species using the LiftOver tool with minMatch parameter 0.001 (Rosenbloom et al., 2015). We then calculated the percentage of segments that were successfully mapped (i.e. retained homology), excluding segments that mapped to a region longer than 1000 bp. The entire simulation was repeated 20 times, starting each time with different sets of 5000 segments. The percentages of segments retaining homology were recorded for each of the 20 simulations, and averaged for each pair of species. These averages were plotted and used as inputs for our statistical framework. Varying the minMatch parameter of the LiftOver tool to 0.5 and segment length to 150 bp allowed us to verify that the observed trends were robust to sequence similarity thresholds and length sampled (Figure 2—figure supplement 2; Figure 5—figure supplement 2).

Nucleotide substitution rate within retained genomic segments (related to Figure 2—figure supplement 1)

The nucleotide sequences of the genomic segments from Figure 2 that retained enough homology to undergo a pairwise alignment were extracted using the getfasta function of BEDTools (Quinlan and Hall, 2010). These sequences were then pairwise aligned using EMBOSS suite’s implementation of Smith-Waterman local alignment (Rice et al., 2000). Default values for gap open penalty (10), gap extend penalty (0.5) and scoring matrix (EDNAFULL) were used to dynamically choose the best local alignment between reference and query sequences. For each cross-species comparison, we calculated the average percent identity of the ungapped alignments of all the segments across 20 randomizations. This procedure yielded values similar to those described previously for the mouse / human (Waterston et al., 2002) and D. melanogaster / D. pseudoobscura comparisons (Richards et al., 2005). The average percent identity of ungapped alignments were used as inputs for our statistical framework, revealing that a model that incorporates lineage labels significantly improved fit to the data relative to a common model (P<0.05; Figure 2—figure supplement 1).

Motif evolutionary rates (related to Figure 5a,b)

Using the FIMO tool (Grant et al., 2011) in the MEME suite (Bailey et al., 2009), the genomes of D. melanogaster and M. musculus domesticus were scanned for matches to experimentally-determined position-weight matrixes corresponding to the GSTFs of interest. Motif matches were called significant according to the default threshold of FIMO, P<10−4. The genomic coordinates of significant motif matches were mapped to the other species within the same lineage using LiftOver (minMatch 0.001). The corresponding coordinates (Mapped) were then extended by 50 bp, and the resulting segments were scanned for motif occurrence (Mappedwithmotif). In order to estimate background expectation, we randomly shuffled the locations of the Mapped segments and scanned these shuffled segments for motifs (ShuffledMappedwithmotif). The percentage of motifs retained relative to background was calculated as:

F=MappedwithmotifShuffledMappedwithmotifMapped*100

The percentages F were then used as measures of cross-species similarity to estimate whether a lineage-aware model would describe the evolution of DNA binding motifs better than a common model (Figure 5—figure supplement 1).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
    R: A Language and Environment for Statistical Computing
    1. R Development Core Team
    (2011)
     Vienna, Austria: R Foundation for Statistical Computing.
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
    The Timetree of Life
    1. SB Hedges
    (2009)
    New York: Oxford University Press.
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
    Initial sequencing and comparative analysis of the mouse genome
    1. RH Waterston
    2. K Lindblad-Toh
    3. E Birney
    4. J Rogers
    5. JF Abril
    6. P Agarwal
    7. R Agarwala
    8. R Ainscough
    9. M Alexandersson
    10. P An
    11. SE Antonarakis
    12. J Attwood
    13. R Baertsch
    14. J Bailey
    15. K Barlow
    16. S Beck
    17. E Berry
    18. B Birren
    19. T Bloom
    20. P Bork
    21. M Botcherby
    22. N Bray
    23. MR Brent
    24. DG Brown
    25. SD Brown
    26. C Bult
    27. J Burton
    28. J Butler
    29. RD Campbell
    30. P Carninci
    31. S Cawley
    32. F Chiaromonte
    33. AT Chinwalla
    34. DM Church
    35. M Clamp
    36. C Clee
    37. FS Collins
    38. LL Cook
    39. RR Copley
    40. A Coulson
    41. O Couronne
    42. J Cuff
    43. V Curwen
    44. T Cutts
    45. M Daly
    46. R David
    47. J Davies
    48. KD Delehaunty
    49. J Deri
    50. ET Dermitzakis
    51. C Dewey
    52. NJ Dickens
    53. M Diekhans
    54. S Dodge
    55. I Dubchak
    56. DM Dunn
    57. SR Eddy
    58. L Elnitski
    59. RD Emes
    60. P Eswara
    61. E Eyras
    62. A Felsenfeld
    63. GA Fewell
    64. P Flicek
    65. K Foley
    66. WN Frankel
    67. LA Fulton
    68. RS Fulton
    69. TS Furey
    70. D Gage
    71. RA Gibbs
    72. G Glusman
    73. S Gnerre
    74. N Goldman
    75. L Goodstadt
    76. D Grafham
    77. TA Graves
    78. ED Green
    79. S Gregory
    80. R Guigó
    81. M Guyer
    82. RC Hardison
    83. D Haussler
    84. Y Hayashizaki
    85. LW Hillier
    86. A Hinrichs
    87. W Hlavina
    88. T Holzer
    89. F Hsu
    90. A Hua
    91. T Hubbard
    92. A Hunt
    93. I Jackson
    94. DB Jaffe
    95. LS Johnson
    96. M Jones
    97. TA Jones
    98. A Joy
    99. M Kamal
    100. EK Karlsson
    101. D Karolchik
    102. A Kasprzyk
    103. J Kawai
    104. E Keibler
    105. C Kells
    106. WJ Kent
    107. A Kirby
    108. DL Kolbe
    109. I Korf
    110. RS Kucherlapati
    111. EJ Kulbokas
    112. D Kulp
    113. T Landers
    114. JP Leger
    115. S Leonard
    116. I Letunic
    117. R Levine
    118. J Li
    119. M Li
    120. C Lloyd
    121. S Lucas
    122. B Ma
    123. DR Maglott
    124. ER Mardis
    125. L Matthews
    126. E Mauceli
    127. JH Mayer
    128. M McCarthy
    129. WR McCombie
    130. S McLaren
    131. K McLay
    132. JD McPherson
    133. J Meldrim
    134. B Meredith
    135. JP Mesirov
    136. W Miller
    137. TL Miner
    138. E Mongin
    139. KT Montgomery
    140. M Morgan
    141. R Mott
    142. JC Mullikin
    143. DM Muzny
    144. WE Nash
    145. JO Nelson
    146. MN Nhan
    147. R Nicol
    148. Z Ning
    149. C Nusbaum
    150. MJ O'Connor
    151. Y Okazaki
    152. K Oliver
    153. E Overton-Larty
    154. L Pachter
    155. G Parra
    156. KH Pepin
    157. J Peterson
    158. P Pevzner
    159. R Plumb
    160. CS Pohl
    161. A Poliakov
    162. TC Ponce
    163. CP Ponting
    164. S Potter
    165. M Quail
    166. A Reymond
    167. BA Roe
    168. KM Roskin
    169. EM Rubin
    170. AG Rust
    171. R Santos
    172. V Sapojnikov
    173. B Schultz
    174. J Schultz
    175. MS Schwartz
    176. S Schwartz
    177. C Scott
    178. S Seaman
    179. S Searle
    180. T Sharpe
    181. A Sheridan
    182. R Shownkeen
    183. S Sims
    184. JB Singer
    185. G Slater
    186. A Smit
    187. DR Smith
    188. B Spencer
    189. A Stabenau
    190. N Stange-Thomann
    191. C Sugnet
    192. M Suyama
    193. G Tesler
    194. J Thompson
    195. D Torrents
    196. E Trevaskis
    197. J Tromp
    198. C Ucla
    199. A Ureta-Vidal
    200. JP Vinson
    201. AC Von Niederhausern
    202. CM Wade
    203. M Wall
    204. RJ Weber
    205. RB Weiss
    206. MC Wendl
    207. AP West
    208. K Wetterstrand
    209. R Wheeler
    210. S Whelan
    211. J Wierzbowski
    212. D Willey
    213. S Williams
    214. RK Wilson
    215. E Winter
    216. KC Worley
    217. D Wyman
    218. S Yang
    219. SP Yang
    220. EM Zdobnov
    221. MC Zody
    222. ES Lander
    223. Mouse Genome Sequencing Consortium
    224. Mouse Genome Sequencing Consortium
    (2002)
    Nature 420:520–562.
    https://doi.org/10.1038/nature01262
  64. 64
  65. 65
  66. 66
    A comparative encyclopedia of DNA elements in the mouse genome
    1. F Yue
    2. Y Cheng
    3. A Breschi
    4. J Vierstra
    5. W Wu
    6. T Ryba
    7. R Sandstrom
    8. Z Ma
    9. C Davis
    10. BD Pope
    11. Y Shen
    12. DD Pervouchine
    13. S Djebali
    14. RE Thurman
    15. R Kaul
    16. E Rynes
    17. A Kirilusha
    18. GK Marinov
    19. BA Williams
    20. D Trout
    21. H Amrhein
    22. K Fisher-Aylor
    23. I Antoshechkin
    24. G DeSalvo
    25. L-H See
    26. M Fastuca
    27. J Drenkow
    28. C Zaleski
    29. A Dobin
    30. P Prieto
    31. J Lagarde
    32. G Bussotti
    33. A Tanzer
    34. O Denas
    35. K Li
    36. MA Bender
    37. M Zhang
    38. R Byron
    39. MT Groudine
    40. D McCleary
    41. L Pham
    42. Z Ye
    43. S Kuan
    44. L Edsall
    45. Y-C Wu
    46. MD Rasmussen
    47. MS Bansal
    48. M Kellis
    49. CA Keller
    50. CS Morrissey
    51. T Mishra
    52. D Jain
    53. N Dogan
    54. RS Harris
    55. P Cayting
    56. T Kawli
    57. AP Boyle
    58. G Euskirchen
    59. A Kundaje
    60. S Lin
    61. Y Lin
    62. C Jansen
    63. VS Malladi
    64. MS Cline
    65. DT Erickson
    66. VM Kirkup
    67. K Learned
    68. CA Sloan
    69. KR Rosenbloom
    70. B Lacerda de Sousa
    71. K Beal
    72. M Pignatelli
    73. P Flicek
    74. J Lian
    75. T Kahveci
    76. D Lee
    77. W James Kent
    78. M Ramalho Santos
    79. J Herrero
    80. C Notredame
    81. A Johnson
    82. S Vong
    83. K Lee
    84. D Bates
    85. F Neri
    86. M Diegel
    87. T Canfield
    88. PJ Sabo
    89. MS Wilken
    90. TA Reh
    91. E Giste
    92. A Shafer
    93. T Kutyavin
    94. E Haugen
    95. D Dunn
    96. AP Reynolds
    97. S Neph
    98. R Humbert
    99. R Scott Hansen
    100. M De Bruijn
    101. L Selleri
    102. A Rudensky
    103. S Josefowicz
    104. R Samstein
    105. EE Eichler
    106. SH Orkin
    107. D Levasseur
    108. T Papayannopoulou
    109. K-H Chang
    110. A Skoultchi
    111. S Gosh
    112. C Disteche
    113. P Treuting
    114. Y Wang
    115. MJ Weiss
    116. GA Blobel
    117. X Cao
    118. S Zhong
    119. T Wang
    120. PJ Good
    121. RF Lowdon
    122. LB Adams
    123. X-Q Zhou
    124. MJ Pazin
    125. EA Feingold
    126. B Wold
    127. J Taylor
    128. A Mortazavi
    129. SM Weissman
    130. JA Stamatoyannopoulos
    131. MP Snyder
    132. R Guigo
    133. TR Gingeras
    134. DM Gilbert
    135. RC Hardison
    136. MA Beer
    137. B Ren
    138. Mouse ENCODE Consortium
    (2014)
    Nature 515:355–364.
    https://doi.org/10.1038/nature13992
  67. 67
    Comparative genomics reveals insights into avian genome evolution and adaptation
    1. G Zhang
    2. C Li
    3. Q Li
    4. B Li
    5. DM Larkin
    6. C Lee
    7. JF Storz
    8. A Antunes
    9. MJ Greenwold
    10. RW Meredith
    11. A Odeen
    12. J Cui
    13. Q Zhou
    14. L Xu
    15. H Pan
    16. Z Wang
    17. L Jin
    18. P Zhang
    19. H Hu
    20. W Yang
    21. J Hu
    22. J Xiao
    23. Z Yang
    24. Y Liu
    25. Q Xie
    26. H Yu
    27. J Lian
    28. P Wen
    29. F Zhang
    30. H Li
    31. Y Zeng
    32. Z Xiong
    33. S Liu
    34. L Zhou
    35. Z Huang
    36. N An
    37. J Wang
    38. Q Zheng
    39. Y Xiong
    40. G Wang
    41. B Wang
    42. J Wang
    43. Y Fan
    44. RR da Fonseca
    45. A Alfaro-Nunez
    46. M Schubert
    47. L Orlando
    48. T Mourier
    49. JT Howard
    50. G Ganapathy
    51. A Pfenning
    52. O Whitney
    53. MV Rivas
    54. E Hara
    55. J Smith
    56. M Farre
    57. J Narayan
    58. G Slavov
    59. MN Romanov
    60. R Borges
    61. JP Machado
    62. I Khan
    63. MS Springer
    64. J Gatesy
    65. FG Hoffmann
    66. JC Opazo
    67. O Hastad
    68. RH Sawyer
    69. H Kim
    70. K-W Kim
    71. HJ Kim
    72. S Cho
    73. N Li
    74. Y Huang
    75. MW Bruford
    76. X Zhan
    77. A Dixon
    78. MF Bertelsen
    79. E Derryberry
    80. W Warren
    81. RK Wilson
    82. S Li
    83. DA Ray
    84. RE Green
    85. SJ O'Brien
    86. D Griffin
    87. WE Johnson
    88. D Haussler
    89. OA Ryder
    90. E Willerslev
    91. GR Graves
    92. P Alstrom
    93. J Fjeldsa
    94. DP Mindell
    95. SV Edwards
    96. EL Braun
    97. C Rahbek
    98. DW Burt
    99. P Houde
    100. Y Zhang
    101. H Yang
    102. J Wang
    103. ED Jarvis
    104. MTP Gilbert
    105. J Wang
    106. C Ye
    107. S Liang
    108. Z Yan
    109. ML Zepeda
    110. PF Campos
    111. AMV Velazquez
    112. JA Samaniego
    113. M Avila-Arcos
    114. MD Martin
    115. R Barnett
    116. AM Ribeiro
    117. CV Mello
    118. PV Lovell
    119. D Almeida
    120. E Maldonado
    121. J Pereira
    122. K Sunagar
    123. S Philip
    124. MG Dominguez-Bello
    125. M Bunce
    126. D Lambert
    127. RT Brumfield
    128. FH Sheldon
    129. EC Holmes
    130. PP Gardner
    131. TE Steeves
    132. PF Stadler
    133. SW Burge
    134. E Lyons
    135. J Smith
    136. F McCarthy
    137. F Pitel
    138. D Rhoads
    139. DP Froman
    140. Z Wang
    141. L Jin
    142. P Zhang
    143. H Hu
    144. W Yang
    145. J Hu
    146. J Xiao
    147. Z Yang
    148. Y Liu
    149. Q Xie
    150. H Yu
    151. J Lian
    152. P Wen
    153. F Zhang
    154. H Li
    155. Y Zeng
    156. Z Xiong
    157. S Liu
    158. L Zhou
    159. Z Huang
    160. N An
    161. J Wang
    162. Q Zheng
    163. Y Xiong
    164. G Wang
    165. B Wang
    166. J Wang
    167. Y Fan
    168. RR da Fonseca
    169. A Alfaro-Núñez
    170. M Schubert
    171. L Orlando
    172. T Mourier
    173. JT Howard
    174. G Ganapathy
    175. A Pfenning
    176. O Whitney
    177. MV Rivas
    178. E Hara
    179. J Smith
    180. M Farré
    181. J Narayan
    182. G Slavov
    183. MN Romanov
    184. R Borges
    185. JP Machado
    186. I Khan
    187. MS Springer
    188. J Gatesy
    189. FG Hoffmann
    190. JC Opazo
    191. O Håstad
    192. RH Sawyer
    193. H Kim
    194. KW Kim
    195. HJ Kim
    196. S Cho
    197. N Li
    198. Y Huang
    199. MW Bruford
    200. X Zhan
    201. A Dixon
    202. MF Bertelsen
    203. E Derryberry
    204. W Warren
    205. RK Wilson
    206. S Li
    207. DA Ray
    208. RE Green
    209. SJ O'Brien
    210. D Griffin
    211. WE Johnson
    212. D Haussler
    213. OA Ryder
    214. E Willerslev
    215. GR Graves
    216. P Alström
    217. J Fjeldså
    218. DP Mindell
    219. SV Edwards
    220. EL Braun
    221. C Rahbek
    222. DW Burt
    223. P Houde
    224. Y Zhang
    225. H Yang
    226. J Wang
    227. ED Jarvis
    228. MT Gilbert
    229. J Wang
    230. Avain Genome Consortium
    (2014)
    Science 346:1311–1320.
    https://doi.org/10.1126/science.1251385
  68. 68
  69. 69

Decision letter

  1. Duncan T Odom
    Reviewing Editor; University of Cambridge / Cancer Research UK, United Kingdom

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your work entitled "Evidence for a common evolutionary rate in metazoan transcriptional networks" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors (Duncan Odom). The evaluation has been overseen by Naama Barkai as the Senior Editor.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission.

Summary:

This was unanimously evaluated as an important analytical contribution to regulatory evolution, and will help rectify a number of discordant results currently in the literature. Carvunis et al. confirm that the rate of gene expression evolution is comparatively stable in both fruit flies and mammalian clades, and newly show that using a uniform analysis approach strongly suggests that underlying TF binding evolution also appears to occur at a comparable rate. They provide some evidence towards the hypothesis that this rate is mediated by the different rates of change found in euchromatic versus heterochromatic components of the mammalian genome.

Essential revisions:

1) The methods should be more carefully documented and the main text and supporting materials sections better cross-referenced (Reviewer #1).

2) The Discussion should be somewhat increased, and the major points brought up by all three reviewers addressed.

3) A few of the results should be revisited with new analysis (see Reviewer #3).

Reviewer #1:

This is an extremely important manuscript in regulatory evolution that, in essence, harmonizes apparently contradictory results from a number of previously published studies by using a carefully designed and uniform analysis methodology. There has been considerable debate around how rapidly transcription and transcriptional regulation evolve between species. In particular, fruit flies and mammals have been reported to have surprisingly different tempos of gene-specific TF binding divergence. Carvunis et al. extract the raw data from a number of studies, and then apply rigorously a combination of analysis protocols with many combinations of widely used tools. They reveal that when analyzed together, gene expression and TF binding evolve similarly in drosophila and mammals. Finally and importantly, Carvunis et al. reveal that a likely culprit for the disconnect between rapid sequence evolution of mammals and the slower SPTF binding and gene expression evolution is the rapid divergence of chromatin inaccessible DNA in mammals. This last point is the key intellectual (as opposed to technical/analytical) insight within this paper.

Aside from more clearly outlined computational handling and a few minor revisions I suggest for consideration below, this straightforward and timely paper will be an excellent contribution to the ongoing discussion of evolutionary genetics.

Major considerations:

1) Materials and methods. Starting from paragraph two of the Results, the authors should direct the reader to the Methods section, and write a concise, but carefully laid out, description of how each step was performed. For instance from paragraph two of the Results, the inter-species normalization of gene expression is a notoriously tricky bit of analysis, per the Brawand et al. Nature 2011 study. I think this really must be carefully walked through. Similarly for all the other sections later.

2) Points not mentioned in this version that could be additionally dissected in the Discussion include:

The impact of effective breeding sizes on how chromatin accessible and inaccessible DNA is handled.

How the breeding rate and absolute time (MY) could interplay in driving evolution (in other words, 40 MY for flies is a LOT more generations than for mice).

Reviewer #2:

This study draws on published datasets to study the evolution of transcriptional regulation in three clades (mammals, birds, and diptera). The basic finding is that rates of evolutionary divergence in transcript abundance, transcription factor binding site occupancy, open chromatin sequence evolution, and transcription factor motif sequences are similar in the three clades. This result is unexpected, as previous studies have reported rather different rates of evolution for some of these features. The authors use this observation to argue that only a small fraction of the genome is involved in transcriptional regulation.

A notable strength of this study is that the authors applied uniform methods of analysis to similar kinds of data, demonstrating that the differences in evolutionary rates reported in the earlier publications are at least in part an artifact of the different methods of analysis that they used. (The impact of different technology platforms, which is also likely a contributing factor, was not addressed.) These results provide an important cautionary lesson, namely that it is essential to work with closely comparable data and to apply uniform methods of analysis before drawing conclusions about biological differences or similarities based on functional genomic datasets. The study is valuable for this reason alone.

Where it is less successful is in providing new insights:

1) The plots showing similar overall rates of evolution (Figure 3, Figure 4, Figure 5) are striking and I'm convinced that the rates are similar among the three clades. But I have no idea what this tells us about the evolution of transcription. I couldn't find any mention in the Results or Discussion sections about how to interpret the observation of rate constancy.

2) The only conclusion presented in the Discussion concerns the fraction of the genome that regulates transcription. Unfortunately, the bulk of the evolutionary comparisons have no bearing on this conclusion. The only one that does is the rate of sequence evolution within chromatin-accessible regions of the genome after removing core promoters. This comparison is relevant, but not because of the rate constancy among clades. It is relevant because the rate in this fraction of the genome is notably slower than the genome as a whole, and specifically within mammals, which have the greatest proportion of noncoding sequences among the three clades. Thus, most of the rate comparisons among the three clades contribute little or nothing to the conclusion that a small fraction of the genome regulates gene expression.

3) The observation that regulatory sequences evolve more slowly than other noncoding sequences is by no means a novel observation. This has been noted previously not just in Sundaram et al. (2014) (as the authors point out), but in several others, including Yue et al. (2014) and papers from Odom's group and Crawford's group. So the basic conclusion of the study was already apparent. Indeed, this is the logic that Graur et al. (2013) used to argue that the ENCODE team overestimated the fraction of noncoding sequences that are functional.

4) The argument that sequence conservation implies that a small fraction of sites regulate transcription is not as straightforward as implied in the Discussion. The ENCODE team wrote a thoughtful rebuttal of the Graur et al. paper (Kellis et al. 2014 PNAS 111:6131) that should be cited in this regard. Beyond the narrow debate about the ENCODE claims, it has been clear for quite some time that some enhancers and transcription factor binding sites turn over quite rapidly in evolution even though they are known to play a role in transcriptional regulation. See Yue et al. (2014) for a recent example, but there are many others. It is clear that negative selection is not the only evolutionary mechanism that operates on regulatory sequences. Using conservation as a criterion for function will thus overlook some functional sites.

Reviewer #3:

Carvunis et al. present an interesting analysis comparing evolutionary rate of genomic sequences and transcriptional regulation, by examining gene expression, transcription factor binding, and open chromatin marks across multiple species of mammals, birds, and drosophila. The authors' analyses recapitulate previous results that sequence gain and loss are much more prominent in mammals than in birds and insects, while rates of gene expression level changes are similar. The authors show that the rates of regulatory changes, such as gain and loss of orthologous TF binding and open chromatin events, are indistinguishable between lineages, providing an answer to conflicting reports on the different rates of evolution for genomic sequences and transcription regulation.

This is a well-written study with important insights. The authors draw a significant conclusion. However, the study has several limitations in their results, reducing their supports to a rather strong conclusion. The authors should either be very clear with the assumptions they made, or presenting stronger and more comprehensive evidence. But I am very excited about this work.

First, the authors use fraction of ortholog sequence retained as a metric to measure rate of genomic sequence evolution. However, the study should consider sequence evolution in the context of rate of substitution, in addition to what the authors already provide. Retained sequences can evolve at different rate, which is directly related to the arguments the authors make.

Second, the authors choose to use the exponential decay model to fit this data. While this is a useful model for evolutionary analysis, this is not the only one. Would their conclusion be solid if they fit the data with a different model?

Third, the species chosen for the study does not provide enough coverage across the 100myrs span. For mammals, there was one data point at 20myrs, and many more data points at either >5myrs, or more than 90myrs. This is a serious concern, because the skewed distribution could potentially bias any model fitting. This is reflected by the TF data for the 20myr comparison. The datapoints at 20myrs were consistently below the fitted curve, often appears quite significantly deviated from the curve.

Lastly, the number of TFs included in the study is too small to support a quite general and significant conclusion. I don't expect many TF data will become available for the study, so perhaps tuning down some of the claims would be a good response.

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

Author response

Essential revisions: 1) The methods should be more carefully documented and the main text and supporting materials sections better cross-referenced (Reviewer #1).

We agree with Reviewer #1 that our methods should be more clearly laid out than they were in the original submission. We thus have more carefully documented our methods and improved cross-referencing throughout the manuscript. For example, the Materials and methods section related to our RNA-seq analyses is now introduced by a succinct paragraph summarizing the various steps of the analysis. Furthermore, each Materials and methods section is now explicitly linked to a figure panel or a set of figure panels.

2) The Discussion should be somewhat increased, and the major points brought up by all three reviewers addressed.

The revised manuscript ends with a considerably longer Discussion. Notably, we now examine our results in light of the differences in breeding sizes and breeding rates (generation times) among the lineages examined. We also list a series of limitations that will need to be addressed in future work to determine how universal our results are, and present a more detailed explanation of why our results suggest that only a small fraction of the genome contributes to transcription.

3) A few of the results should be revisited with new analysis (see Reviewer #3).

Following Reviewer #3’s suggestions, we have made the following improvements to our analyses:

We have investigated evolution of genome sequences using rate of substitution as a metric in addition to the analyses already provided in the original submission.

To increase the number of data points and the quality of fits in our analyses of ChIP-seq data, we have incorporated an additional species (dog) for which data were available in the literature.

We have improved the sensitivity of our statistical framework.

Detailed explanations of each of these improvements are given below in a point-by-point response to Reviewer #3.

Reviewer #1:

1) Methods. Starting from paragraph two of the Results, the authors should direct the reader to the Methods section, and write a concise, but carefully laid out, description of how each step was performed. For instance from paragraph two of the Results, the inter-species normalization of gene expression is a notoriously tricky bit of analysis, per the Brawand et al. Nature 2011 study. I think this really must be carefully walked through. Similarly for all the other sections later.

Following the reviewer’s suggestion, we have added cross-references to the Materials and methods section and to the corresponding supplementary figures and files throughout the manuscript. The description of each step performed for the gene expression analysis is presented in the subsection “Gene expression evolutionary rates (related to Figure 3)” of the revised manuscript and relies on the orthology relationships described in the subsection “Genome and Annotation Sources”. The revised manuscript now includes a paragraph clarifying the steps involved. Similar clarifications were also introduced in other Materials and methods sections. Regarding normalization of RNA-seq data, the algorithm we used for transcript quantification, Sailfish, performs several bias correction steps that are described in the corresponding publication (Patro Nature Biotech, 2014). We note that Brawand and colleagues (Nature, 2011) implemented additional normalization steps in order to search for differentially expressed genes, which we did not attempt in this work.

2) Points not mentioned in this version that could be additionally dissected in the Discussion include:

The impact of effective breeding sizes on how chromatin accessible and inaccessible DNA is handled.

How the breeding rate and absolute time (MY) could interplay in driving evolution (in other words, 40 MY for flies is a LOT more generations than for mice).

To address this comment, we have significantly lengthened our Discussion. Among other improvements, the revised manuscript discusses our results in light of the differences in breeding sizes and breeding rates (generation times) among the lineages examined.

Reviewer #2:

1) The plots showing similar overall rates of evolution (Figure 3, Figure 4, Figure 5) are striking and I'm convinced that the rates are similar among the three clades. But I have no idea what this tells us about the evolution of transcription. I couldn't find any mention in the Results or Discussion sections about how to interpret the observation of rate constancy.

We fully agree with the reviewer that our original submission lacked a thorough Discussion section. Encouraged by this and the other reviewers’ suggestions, we have now included in the revised manuscript a longer Discussion section, including speculative interpretations.

2) The only conclusion presented in the Discussion concerns the fraction of the genome that regulates transcription. Unfortunately, the bulk of the evolutionary comparisons have no bearing on this conclusion. The only one that does is the rate of sequence evolution within chromatin-accessible regions of the genome after removing core promoters. This comparison is relevant, but not because of the rate constancy among clades. It is relevant because the rate in this fraction of the genome is notably slower than the genome as a whole, and specifically within mammals, which have the greatest proportion of noncoding sequences among the three clades. Thus, most of the rate comparisons among the three clades contribute little or nothing to the conclusion that a small fraction of the genome regulates gene expression.

To address this comment, we have revised the Discussion section such that the hypothesis that a small fraction of the mammalian genome regulates transcription is presented as one of several possible interpretations.

In the revised Discussion, we claim that the evolutionary rate of transcription is decoupled from the lineage-specific parameters that govern the evolution of genome sequence. We present the hypothesis that a small fraction of the mammalian genome contributes to transcriptional regulation as a plausible explanation given our observations within chromatin-accessible regions. We also present alternative hypotheses where rapidly evolving portions of the bird and insect genomes may make significant contributions to transcriptional regulation.

3) The observation that regulatory sequences evolve more slowly than other noncoding sequences is by no means a novel observation. This has been noted previously not just in Sundaram et al. (2014) (as the authors point out), but in several others, including Yue et al. (2014) and papers from Odom's group and Crawford's group. So the basic conclusion of the study was already apparent. Indeed, this is the logic that Graur et al. (2013) used to argue that the ENCODE team overestimated the fraction of noncoding sequences that are functional.

We fully agree with the reviewer that the fact that many regulatory sequences are conserved across species was known before our submission. We specifically quoted Yue et al. (2014) in this context as it pertains specifically to comparing chromatin-accessible sequences to random non-coding regions, similar to what we do in the corresponding paragraph. While the slow divergence of chromatin-accessible regions relative to other non-coding regions in mammals may have been easily predicted from pre-existing literature, there was no previously published indication that such difference would be practically nonexistent in birds and insects, leading to a shared divergence rate across the three clades. We have now added text raising and explaining these points in the manuscript Discussion.

4) The argument that sequence conservation implies that a small fraction of sites regulate transcription is not as straightforward as implied in the Discussion. The ENCODE team wrote a thoughtful rebuttal of the Graur et al. paper (Kellis et al. 2014 PNAS 111:6131) that should be cited in this regard. Beyond the narrow debate about the ENCODE claims, it has been clear for quite some time that some enhancers and transcription factor binding sites turn over quite rapidly in evolution even though they are known to play a role in transcriptional regulation. See Yue et al. (2014) for a recent example, but there are many others. It is clear that negative selection is not the only evolutionary mechanism that operates on regulatory sequences. Using conservation as a criterion for function will thus overlook some functional sites.

We wholeheartedly agree with the reviewer that conservation is by no means the only indicator of regulatory function. We have included the Kellis (2014) reference in the Introduction of our revised manuscript. We never meant to imply that only conserved sites can be functional. Our study is showing that, at different levels of the transcription process, evolutionary rates are shared across mammals, birds and insects. Our point is that these rates could not have been predicted from genome-wide patterns of DNA evolution, as they can be slow (transcript levels, chromatin-accessible sequences) or fast (TF binding, motifs), irrespective of how slow or fast the genome evolves as a whole. We have clarified our language in the revised manuscript.

Reviewer #3:

First, the authors use fraction of ortholog sequence retained as a metric to measure rate of genomic sequence evolution. However, the study should consider sequence evolution in the context of rate of substitution, in addition to what the authors already provide. Retained sequences can evolve at different rate, which is directly related to the arguments the authors make.

To address this comment, we have measured how much retained sequences changed in nucleotide identity over time. We found that retained sequences were highly conserved at the nucleotide level in mammals, birds and insects. However, insect sequences diverged slightly but significantly faster than bird and mammal sequences. The results of this analysis are shown in Figure 2—figure supplement 1 and discussed in the revised manuscript.

Second, the authors choose to use the exponential decay model to fit this data. While this is a useful model for evolutionary analysis, this is not the only one. Would their conclusion be solid if they fit the data with a different model?

To address the reviewer’s comment, we conducted an additional analysis aimed at comparing an exponential decay with a linear decay and with a power law decay. Models with more than two parameters were not investigated due to the higher risk of over-fitting.

To compare the quality of these various decay models, we encoded them in R in a unified regression framework relying on non-linear least squares (nls) – as opposed to the linear least square regression (lm) method used in the original submission. We quantitatively estimated fit qualities across all the data layers and lineages used throughout our manuscript using the AICc criteria, a small-sample-size corrected version of Akaike information criterion. We found that the data are almost always better represented by an exponential decay than by either a linear or a power law decay (as determined by the lowest AICc). The only data type better modelled by a power law decay than by an exponential decay was the nucleotide identity of retained genomic segments, a data type introduced to the manuscript as per the above suggestion of the same reviewer. We felt that introducing additional models that do not adequately fit the majority of the data would not make our results more robust. Therefore, we did not specifically investigate if linear or power law decays would give different results than an exponential decay when comparing lineages with our statistical framework, and we kept the exponential model as the basis of our work.

While performing this analysis, we discovered that model fitting using nls was much more sensitive to detect changes in decay rates than fitting with lm. In spite of this increase in sensitivity, the vast majority of our results and all of our major conclusions remained unchanged. We thus opted to use nls throughout our revised manuscript.

Third, the species chosen for the study does not provide enough coverage across the 100myrs span. For mammals, there was one data point at 20myrs, and many more data points at either >5myrs, or more than 90myrs. This is a serious concern, because the skewed distribution could potentially bias any model fitting. This is reflected by the TF data for the 20myr comparison. The datapoints at 20myrs were consistently below the fitted curve, often appears quite significantly deviated from the curve.

We agree that the distribution of mammalian species with fully sequenced genomes across the 100 Myrs range is less than ideal for model fitting. We were able to include a large number of species in some of our analyses (e.g. Figure 2 includes rabbit and guinea pig, separated from mice by 75 and 86 Myrs, respectively), but the lack of ChIP data for species separated from mouse by more than 20 but less than 90 Myrs is indeed a concern (reflected specifically in Figure 4 and Figure 5a).

In these specific figure panels, the data points at 20 Myrs (corresponding to mouse-rat comparison) are often below the fitted curve while the data points at 91 Myrs (mouse-human) are reasonably well captured. Given the data presented in the submitted manuscript, it is unclear whether the mouse-rat or the mouse-human comparison is the outlier. To remedy this gap to the best of our ability given the available ChIP data, we have added dog samples to our analyses of HNF4A, CEBPA and CTCF. Since dog is very distant from mouse (97 Myrs), adding these samples improved fit confidence at large evolutionary distances. In addition to this improvement, we also added to the revised manuscript a discussion of the limitations of our study posed by the non-uniform coverage of datasets across evolutionary time.

Lastly, the number of TFs included in the study is too small to support a quite general and significant conclusion. I don't expect many TF data will become available for the study, so perhaps tuning down some of the claims would be a good response.

We agree that the number of TFs with occupancy data investigated is too small to make a generalization on TF occupancy evolution, and we have therefore toned down our language in the corresponding paragraphs accordingly.

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

Article and author information

Author details

  1. Anne-Ruxandra Carvunis

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    A-RC, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Tina Wang and Dylan Skola
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-6474-6413
  2. Tina Wang

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    TW, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Anne-Ruxandra Carvunis and Dylan Skola
    Competing interests
    The authors declare that no competing interests exist.
  3. Dylan Skola

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    DS, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Anne-Ruxandra Carvunis and Tina Wang
    Competing interests
    The authors declare that no competing interests exist.
  4. Alice Yu

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    AY, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  5. Jonathan Chen

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    JC, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  6. Jason F Kreisberg

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    JFK, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  7. Trey Ideker

    Department of Medicine, University of California, San Diego, La Jolla, United States
    Contribution
    TI, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    tideker@ucsd.edu
    Competing interests
    The authors declare that no competing interests exist.

Funding

National Institute of General Medical Sciences (R01 GM084279)

  • Trey Ideker

National Institute of General Medical Sciences (P50 GM085764)

  • Trey Ideker

National Institute of General Medical Sciences (T32 GM008666)

  • Tina Wang

National Institute of General Medical Sciences (K99 GM108865)

  • Anne-Ruxandra Carvunis

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

Acknowledgements

This work was supported by the National Institutes of Health: R01 GM084279 and P50 GM085764 awarded to TI, T32 GM008666 awarded to TW and the Pathway to Independence K99 GM108865 awarded to A-RC. The authors are grateful to Drs. Pollard KS, Wittkopp PJ, Burke M, Charloteaux B and Coolon JD for review of the manuscript prior to submission and to the editors and reviewers for their help in improving the manuscript after submission.

Reviewing Editor

  1. Duncan T Odom, Reviewing Editor, University of Cambridge / Cancer Research UK, United Kingdom

Publication history

  1. Received: September 15, 2015
  2. Accepted: December 17, 2015
  3. Accepted Manuscript published: December 18, 2015 (version 1)
  4. Version of Record published: February 11, 2016 (version 2)

Copyright

© 2015, Carvunis et al.

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

Metrics

  • 2,235
    Page views
  • 508
    Downloads
  • 4
    Citations

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

Comments

Download links

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

Downloads (link to download the article as PDF)

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

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

Further reading

    1. Computational and Systems Biology
    2. Genomics and Evolutionary Biology
    Matthew Hur et al.
    Tools and Resources Updated
    1. Neuroscience
    Chia-Wei Chang et al.
    Research Article Updated