Towards an evolutionary model of transcription networks

Abstract

DNA evolution models made invaluable contributions to comparative genomics, although it seemed formidable to include non-genomic features into these models. In order to build an evolutionary model of transcription networks (TNs), we had to forfeit the substitution model used in DNA evolution and to start from modeling the evolution of the regulatory relationships. We present a quantitative evolutionary model of TNs, subjecting the phylogenetic distance and the evolutionary changes of cis-regulatory sequence, gene expression and network structure to one probabilistic framework. Using the genome sequences and gene expression data from multiple species, this model can predict regulatory relationships between a transcription factor (TF) and its target genes in all species, and thus identify TN re-wiring events. Applying this model to analyze the pre-implantation development of three mammalian species, we identified the conserved and re-wired components of the TNs downstream to a set of TFs including Oct4, Gata3/4/6, cMyc and nMyc. Evolutionary events on the DNA sequence that led to turnover of TF binding sites were identified, including a birth of an Oct4 binding site by a 2nt deletion. In contrast to recent reports of large interspecies differences of TF binding sites and gene expression patterns, the interspecies difference in TF-target relationship is much smaller. The data showed increasing conservation levels from genomic sequences to TF-DNA interaction, gene expression, TN, and finally to morphology, suggesting that evolutionary changes are larger at molecular levels and smaller at functional levels. The data also showed that evolutionarily older TFs are more likely to have conserved target genes, whereas younger TFs tend to have larger re-wiring rates.

Figure 1. A model for TN evolution.

(A) An example of evolving TNs. Transcription factors and target genes are depicted in purple and blue nodes. Conserved and species-specific regulatory relationships are depicted with red and black arrows. (B) The regulatory states are modeled as a continuous time Markov chain, and the regulatory sequence and the gene expression data are emitted from the hidden regulatory states.

Figure 2. Conservation and rewiring of Oct4 downstream TN.

(A) Venn diagram of model-inferred Oct4 target genes in three species. Each node is a gene and each edge represents a transcriptional regulatory link. When a set of genes are linked to a regulatory hub gene, this set of genes is suppressed into one node, annotated with the name of the hub gene followed by the number of other genes transcriptionally linked to the hub. (B) Inferred conservation levels for DNA, TFBS, gene expression, TN, and morphology. (C) Comparison of model-inferred re-wiring rates with (Exp+Seq+Phy) and without (Exp+Phy) DNA sequence data. The re-wiring rate of TF was calculated as the percentage of non-conserved TF-target relationships among all TF-target relationships. A background distribution (red curve) of the re-wiring rate is derived from randomly permuting the rows and columns of the Oct4 position specific score matrix (PSSM) and feeding the permutated PSSM to the model with unperturbed sequence and gene expression data.

Figure 3. Examples of model-inferred species-specific target genes of Oct4 in humans (A) and mice (B).

Oct4 is zygotically expressed (4–8 cells), and it is strongly increased at the late stages of pre-implantation development, including morula (M) and blastocyst (B). The target genes' expression is zygotically activated in a species-specific manner (Left panel). The species-specific binding of Oct4 to target genes is observed in ChIP-seq experiments (TFBS track, right panel). The genomic sequences of the Oct4 binding regions are shown in yellow, with Oct4 binding motifs shown in red boxes.

Figure 4. Relationship between conservation and binding affinities.

(A) The distribution of the average binding affinities between Oct4 and orthologous regulatory sequences in three species. (B) Interspecies difference in binding affinities of the 20 kb upstream sequences of orthologous genes. The orthologous genes are put into four categories of conservation. Error bars show 95% confidence intervals.

Figure 5. Re-wiring rate against divergence time.

(A) Re-wiring rate against phylogenetic distance. Circles: Gata3/Gata4/Gata6 genes. Cross: cMyc(MYC)/nMyc(MYCN) genes. TreeFam branch length of a gene is the distance between the speciation event of the gene and the first duplication event in the paralogous family. (B, C) Phylogenetic trees of the Gata and the Myc families in mammals. Blue dot: gene speciation event. Red dot: gene duplication event. Branch lengths are estimated for the consensus tree of bootstrapped trees using HKY model. Numbers on branching events are the support numbers to the consensus tree in 100 bootstrapped trees.