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
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
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.
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).
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).
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).
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).
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.
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).
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).
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).
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.
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.
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.
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).
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.
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).
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).
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:
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).
MEME SUITE: tools for motif discovery and searchingNucleic Acids Research 37:W202–W208.https://doi.org/10.1093/nar/gkp335
A computational pipeline for comparative ChIP-seq analysesNature Protocols 7:45–61.https://doi.org/10.1038/nprot.2011.420
Trimmomatic: a flexible trimmer for illumina sequence dataBioinformatics 30:2114–2120.https://doi.org/10.1093/bioinformatics/btu170
Conservation of core gene expression in vertebrate tissuesJournal of Biology 8:33.https://doi.org/10.1186/jbiol130
Mapping complex disease traits with global gene expressionNature Reviews Genetics 10:184–194.https://doi.org/10.1038/nrg2537
Tempo and mode of regulatory evolution in drosophilaGenome Research 24:797–808.https://doi.org/10.1101/gr.163014.113
R: A Language and Environment for Statistical ComputingVienna, Austria: R Foundation for Statistical Computing.
On the immortality of television sets: "function" in the human genome according to the evolution-free gospel of ENCODEGenome Biology and Evolution 5:578–590.https://doi.org/10.1093/gbe/evt028
The Timetree of LifeNew York: Oxford University Press.
Variation in transcription factor binding among humansNew Biotechnology 27:S81–S235.https://doi.org/10.1016/j.nbt.2010.01.227
Defining functional DNA elements in the human genomeProceedings of the National Academy of Sciences of the United States of America 111:6131–6138.https://doi.org/10.1073/pnas.1318948111
Design and analysis of ChIP-seq experiments for DNA-binding proteinsNature Biotechnology 26:1351–1359.https://doi.org/10.1038/nbt.1508
ChIP-seq guidelines and practices of the ENCODE and modENCODE consortiaGenome Research 22:1813–1831.https://doi.org/10.1101/gr.136184.111
Comparison of the transcriptional landscapes between human and mouse tissuesProceedings of the National Academy of Sciences of the United States of America 111:17224–17229.https://doi.org/10.1073/pnas.1413624111
Can ENCODE tell us how much junk DNA we carry in our genome?Biochemical and Biophysical Research Communications 430:1340–1343.https://doi.org/10.1016/j.bbrc.2012.12.074
The nearly neutral theory of molecular evolutionAnnual Review of Ecology and Systematics 23:263–286.https://doi.org/10.1146/annurev.es.23.110192.001403
EMBOSS: the european molecular biology open software suiteTrends in Genetics 16:276–277.https://doi.org/10.1016/S0168-9525(00)02024-2
Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomesGenome Research 15:1034–1050.https://doi.org/10.1101/gr.3715005
Genome-wide analysis of chromatin states reveals distinct mechanisms of sex-dependent gene regulation in male and female mouse liverMolecular and Cellular Biology 33:3594–3610.https://doi.org/10.1128/MCB.00280-13
Temporal patterns of fruit fly (drosophila) evolution revealed by mutation clocksMolecular Biology and Evolution 21:36–44.https://doi.org/10.1093/molbev/msg236
Transcriptome sequencing reveals the character of incomplete dosage compensation across multiple tissues in flycatchersGenome Biology and Evolution 5:1555–1566.https://doi.org/10.1093/gbe/evt114
Evolution of transcription factor binding in metazoans — mechanisms and functional implicationsNature Reviews Genetics 15:221–233.https://doi.org/10.1038/nrg3481
Duncan T OdomReviewing 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.
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.
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).
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.
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).
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.
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
- Trey Ideker
- Trey Ideker
- Tina Wang
- Anne-Ruxandra Carvunis
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
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.
- Duncan T Odom, Reviewing Editor, University of Cambridge / Cancer Research UK, United Kingdom
© 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.