1. Computational and Systems Biology
  2. Evolutionary Biology

Evidence for a common evolutionary rate in metazoan transcriptional networks

  1. Anne-Ruxandra Carvunis
  2. Tina Wang
  3. Dylan Skola
  4. Alice Yu
  5. Jonathan Chen
  6. Jason F Kreisberg
  7. Trey Ideker  Is a corresponding author
  1. University of California, San Diego, United States
https://doi.org/10.7554/eLife.11615
Share this article
Cite this article
  1. Anne-Ruxandra Carvunis
  2. Tina Wang
  3. Dylan Skola
  4. Alice Yu
  5. Jonathan Chen
  6. Jason F Kreisberg
  7. Trey Ideker
(2015)
Evidence for a common evolutionary rate in metazoan transcriptional networks
eLife 4:e11615.
https://doi.org/10.7554/eLife.11615
  2. Download BibTeX
  3. Download .RIS

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
Download asset Open asset
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
Download asset Open asset
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
Download asset Open asset
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
Figure 3—source data 1

Accession numbers used in RNA-seq analyses.

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

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
Download asset Open asset
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
Figure 4—source data 1

Accession numbers used in ChIP-seq analyses.

https://doi.org/10.7554/eLife.11615.013
Figure 4—source data 2

648 segment-based ChIP analyses.

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

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
Download asset Open asset
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. (c, d) 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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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=Mappedwithmotif−ShuffledMappedwithmotifMapped*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).

Data availability

The following previously published data sets were used
    1. Malone JH
    2. Artieri CG
    3. Sturgill D
    4. Zhang Y
    5. Oliver B
    (2011) mRNA-Seq of whole flies from Drosophila
    Publicly available at the NCBI Gene Expression Omnibus (Accession no: GSE28078).
    http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28078
    1. Wang C
    2. Auerbach SS
    3. Shi L
    4. Tong W
    (2014) SEQC Toxicogenomics Study: RNA-Seq data set
    Publicly available at the NCBI Gene Expression Omnibus (Accession no: GSE55347).
    http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55347
    1. Malone JH
    2. Artieri CG
    3. Sturgill D
    4. Zhang Y
    5. Oliver B
    (2011) mRNA-Seq of whole flies from Drosophila
    Publicly available at the NCBI Gene Expression Omnibus (Accession no: GSE28078).
    http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28078
    1. Li N
    2. Huang Y
    (2013) Genome and transcriptome of the duck
    Publicly available at the NCBI Gene Expression Omnibus (Accession no: GSE22967).
    http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22967
    1. He Q
    2. Bardet AF
    3. Patton B
    4. Purvis J
    5. Johnston J
    6. Paulson A
    7. Gogol M
    8. Stark A
    9. Zeitlinger J
    (2011) cross_species
    Publicly available at the Sequence Read Archive + (Accession no: ERP000357).
    http://sra.dnanexus.com/studies/ERP000357

References

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

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 "This ORCID iD identifies the author of this article:" 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.

Version 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,866
    views
  • 607
    downloads
  • 21
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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)

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

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

  1. Anne-Ruxandra Carvunis
  2. Tina Wang
  3. Dylan Skola
  4. Alice Yu
  5. Jonathan Chen
  6. Jason F Kreisberg
  7. Trey Ideker
(2015)
Evidence for a common evolutionary rate in metazoan transcriptional networks
eLife 4:e11615.
https://doi.org/10.7554/eLife.11615

Share this article

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

Categories and tags

Research organisms

Further reading

    1. Computational and Systems Biology
    2. Evolutionary Biology

    Comparative Genomics: One for all

    Duncan T Odom
    Insight

    Using a common analysis pipeline to compare data from three major lineages of complex eukaryotes reveals that transcription seems to evolve at a common rate.

    1. Computational and Systems Biology
    2. Genetics and Genomics

    Haplotype function score improves biological interpretation and cross-ancestry polygenic prediction of human complex traits

    Weichen Song, Yongyong Shi, Guan Ning Lin
    Tools and Resources

    We propose a new framework for human genetic association studies: at each locus, a deep learning model (in this study, Sei) is used to calculate the functional genomic activity score for two haplotypes per individual. This score, defined as the Haplotype Function Score (HFS), replaces the original genotype in association studies. Applying the HFS framework to 14 complex traits in the UK Biobank, we identified 3619 independent HFS–trait associations with a significance of p < 5 × 10−8. Fine-mapping revealed 2699 causal associations, corresponding to a median increase of 63 causal findings per trait compared with single-nucleotide polymorphism (SNP)-based analysis. HFS-based enrichment analysis uncovered 727 pathway–trait associations and 153 tissue–trait associations with strong biological interpretability, including ‘circadian pathway-chronotype’ and ‘arachidonic acid-intelligence’. Lastly, we applied least absolute shrinkage and selection operator (LASSO) regression to integrate HFS prediction score with SNP-based polygenic risk scores, which showed an improvement of 16.1–39.8% in cross-ancestry polygenic prediction. We concluded that HFS is a promising strategy for understanding the genetic basis of human complex traits.

    1. Computational and Systems Biology

    Genome-scale annotation of protein binding sites via language model and geometric deep learning

    Qianmu Yuan, Chong Tian, Yuedong Yang
    Tools and Resources

    Revealing protein binding sites with other molecules, such as nucleic acids, peptides, or small ligands, sheds light on disease mechanism elucidation and novel drug design. With the explosive growth of proteins in sequence databases, how to accurately and efficiently identify these binding sites from sequences becomes essential. However, current methods mostly rely on expensive multiple sequence alignments or experimental protein structures, limiting their genome-scale applications. Besides, these methods haven’t fully explored the geometry of the protein structures. Here, we propose GPSite, a multi-task network for simultaneously predicting binding residues of DNA, RNA, peptide, protein, ATP, HEM, and metal ions on proteins. GPSite was trained on informative sequence embeddings and predicted structures from protein language models, while comprehensively extracting residual and relational geometric contexts in an end-to-end manner. Experiments demonstrate that GPSite substantially surpasses state-of-the-art sequence-based and structure-based approaches on various benchmark datasets, even when the structures are not well-predicted. The low computational cost of GPSite enables rapid genome-scale binding residue annotations for over 568,000 sequences, providing opportunities to unveil unexplored associations of binding sites with molecular functions, biological processes, and genetic variants. The GPSite webserver and annotation database can be freely accessed at https://bio-web1.nscc-gz.cn/app/GPSite.