A comparative evolutionary study of transcription networks.

We use a simple duplication-divergence model for the growth of the transcription network tesikirill ; EI-06 ; CLJ+07 to guide the data analysis. The model is formulated as follows. At each step all, or a fraction of the nodes are duplicated. The duplicate nodes inherit all the in- and outgoing edges of the original (ancestral) node. In other words one supposes that before divergence, the binding sites on the proteins and the regulatory regions on DNA remain identical. Subsequently, the new edges are removed with a certain probability, that might depend on the status of the node (regulator or target, subject or regulating a new or an old gene, etc.).

A qualitative study of the model leads to the following schematic resultstesikirill . First, by definition there is no possibility of de novo addition of edges by rewiring. Hence no edges can be created from an homology class that initially does not regulate another one. In the simplest case, if feedback is initially absent, it cannot arise spontaneously. Moreover, no hierarchical layer can be added to the network, and duplicates selected for fixation will lie in the same layer. On the contrary, in the presence of feedback, and ARs in particular, the model behaves differently. Duplication of ARs can give rise to higher order feedback and to new ARs. However, this feedback will be strictly confined among members of the same homology class. Analogously, the degree distribution is piloted by the degree distribution of duplicates. The in- and out- degree distributions can be decoupled by choosing different removal probabilities for old/new and new/old edges. One can obtain power-law out-degree and compact in-degree distributions that are in close relation with the phylogenetic conservation of the network proteins EI-06 ; EI07 : conserved nodes must exhibit at least one between the in- and out-degree sequence with a scale-free distribution, while proteins with only exponential degrees cannot be phylogenetically conserved under a general duplication-divergence model. According to this model, bacterial transcription networks with scale-free out-degrees and exponential in-degrees are thus consistent with a phylogenetically conserved set of transcription factors and a non-phylogenetically conserved set of target genes (possibly due to abundant horizontal tranfers of target genes).

In a realistic situation, the above observations may not apply because of the possibility of edge rewiring, that may be able to shuffle the hierarchy, create new feedback, affect the degree distribution, and mix edges among homology classes. Note, however, that some regulatory edges between apparently different homology classes may have also resulted from very old edges within ancestral homology classes that are no longer classified as single homology classes due to their gradual divergence. Besides, while actual edge rewiring may have occurred after homology class separation, it is clear that edge deletion must typically happen after gene duplication (or else extant regulatory networks would be densely cross-edged graphs.) Hence, this model, which focuses on the necessary deletion of duplicated edges, can assess the specific role of gene duplication and edge deletion in the growth of the network and also underline their possible shortcomings to delineate the role of other evolutionary processes such as edge rewiring.

Having in mind the above qualitative results we examine the topological roles in empirical networks of nodes coming from the same common ancestor. We define proteins that are likely to share a common ancestor through structural domain architectures TB04 ; MBT03 ; BD92 . These domains allow for the definition of larger classes than sequence comparison alone TB04 . The database enables to associate an ordered sequence of domains, or “domain architecture” to each protein. We define protein homologs as proteins whose domain architectures are identical neglecting domain repeats. This corresponds to a conservative view of homology where no new domains are acquired or lost after duplication. The results are tested using different definitions of homologs sellerio002 . Homology allows us to construct sets, or classes, of proteins which have then supposedly evolved from a common ancestral gene. We have analyzed the distribution of regulatory edges between and within classes of the likely duplicate genes. The statistical significance of the analysis in terms of homology classes is established BLA+04 by comparison with random shuffling of genes (TFs and TGs separately) between classes. This analysis is limited by the number of nodes for which homology classes can be constructed, and thus less sensitive to the particular data set.

In all cases we find that the motifs of duplicated TGs, regulated by a single or duplicate TFs, and duplicate TFs regulating a common target (Fig. 2), are significantly overrepresented TB04 , which can be seen as a validation of the model. The signal for this is smaller in the smaller data sets (such as B. subtilis), because of poor statistics.

A first question to address is whether gene duplications are able to pilot the observed degree sequences of the network, as predicted by the model. For the case of yeast, one can perform this test on the above defined homology classes, and also, on the doubly conserved genes in the yeast whole-genome duplication proved by Kellis et al. KBL04 . Precisely, one can measure the degree distributions restricted to targets or regulator falling in the same homology class, and verify whether it scales in a similar way as the unrestricted distribution.

Our results, considering the Balaji data set for yeast and RegulonDB 5.5 for the E. coli, are shown in Fig. 2. The homologs follow a qualitatively similar distribution to the entire network, exponential for the in-degree, and power-law-like for the out-degree. For the in-degree distributions we find a better agreement between the decay of the entire network and the restriction to homology classes. This confirms the general idea that the degree distributions are strongly dependent upon duplications. However, in the case of yeast, the out-degree of the entire network seems to follow a broader distribution than that of the duplicates (top-left panel of Fig. 2). This can be seen as an indication that other processes, such as rewiring, concur in defining the targets of a TF. The behavior of the duplicates from the whole-genome duplication is conditioned by the small size of the sample, and its degree distribution drops with a faster decay for both the in- and the out-degree.

Our analysis shows that the role of autoregulatory interactions is rather different in the bacterial data sets compared to yeast. The fraction of self-regulated TFs is above $50\%$ in all three bacterial data sets, while being around $10\%$ only in the two S. cerevisiae data sets. If we measure directly the distribution of ARs in classes, in the E. coli datasets we find a consistent signal that the population of ARs in homology classes is more dense and more variable than in randomized instances CLJ+07 . This is a direct evidence of duplication and conservation of ARs, and agrees with the conservation of a hierarchical structure during evolutionary growth. In B. subtilis and yeast, we find that the population of ARs in classes falls closer to (although systematically above) the average of the randomized samples. This means that, in these cases, the pure analysis of AR population in homology classes is not sufficient to assess that ARs in the same class are the result of duplication and inheritance of self-edges. In the smaller data-sets (such as DBTBS MNO+04 ) this absence of signal could be likely due to small size of the sample, so it is possible that this is another trend differentiating bacteria from eukaryotes.

On the other hand, in all the data sets we find evidence for the relevance of AR duplications if we consider the distribution of the subgraphs that can stem from AR duplication (Fig. 3) within homology classes (Table 1). These subgraphs are compatible with AR duplication conserving crosstalks. In order to quantitatively evaluate the relevance of the subgraphs, we define an observable that counts the occurrence of the following events

1. 1.

   An AR regulates a homologous, non-self-regulatory transcription factor.

2. 2.

   An AR is regulated by a homologous, non-self-regulatory transcription factor.

3. 3.

   Two homolog ARs possess reciprocal cross regulation.

We indicate with $N_{c}$ the total number of homology classes, with $M_{i}$ the total number of the subgraphs of the three kinds found in homology class $i$, and by $L_{i}$ the number of edges in the class, and define

$$a=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\frac{M_{i}}{L_{i}}\quad.$$

This observable represents the ratio of subgraphs to the total number of (self- and non-self-) edges present in the homology class.

We find that, although the fixation of cross-talks within homology classes is typically less frequent that the retention of self-edges in bacteria CLJ+07 , the AR duplication subgraphs are systematically overrepresented, as shown by the positive signal. Specifically (see Supplementary Table 1), in E. coli and B. subtilis, this signal is dominated by duplicated self-interactions followed by loss of crosstalks, while for the S. cerevisiae network, subgraphs including crosstalks are more abundant, compatibly with the observed more complex hierarchy CLJ+07 . In addition, in the significantly overrepresented AR duplication subgraphs (type ii in figure 3), only one of the two possible self-regulations of duplicate TFs survives, indicating that ARs do not proliferate by duplication. This observation indicates that duplications have been important to shape the non-self-regulatory edges within homology classes, thereby redefining the TF hierarchy, besides leading to the fixation of new ARs in the bacterial networks. This is especially the case in yeast, indicating that circuits with crosstalks or feedback stemming from AR duplication are not negligible, and might have contributed to the build-up of the existing feedback.

To gain further insight into this point, we evaluated the distribution of interactions within and across homology classes. Considering the transcription factors homology classes, we constructed a “collapsed” weighted network as follows. The nodes are homology classes, and weighted oriented edges represent the number of members of a class regulating members of another class. In absence of rewiring, this would be a likely “ancestral” TF-TF network, inherited by duplication of “primitive” self- and heterologous edges. Our model dictates that this ancestral graph can only be hierarchical or maintain the same amount of feedback of the initial configuration. In particular, if an ancestral edge is sent out from a class to a second one and the reverse is not true, there will be no edges coming from the second class to the first in the evolved network. This corresponds to an asymmetric regulatory control of the first homology class over the second homology class. Conversely, the appearance of symmetric regulations may be a signature of rewiring, if one assumes that retaining both ancestral crosstalks produced by the duplication of an ancestral AR was then as unlikely as it appears to be between recent AR duplicates. Hence, symmetric regulation between homology classes in the extant network actually suggests that the likely initial asymmetry conserved by duplication was progressively smoothed down by edge rewiring or other evolutionary mechanisms with an equivalent effect.

We define the observable $R$, which represents a measure of the strength of asymmetric regulation across the classes. We indicate as $L_{ij}$ the number of regulatory edges between homology class $i$ and class $j$. Let then $C_{i}$ represent the number of elements in homology class $i$. Then we define

$$R=\sum_{{i\neq j}}\frac{|L_{ij}-L_{ji}|}{\sqrt{C_{i}^{2}+C_{j}^{2}}}\quad,$$

where the sum is performed over all the pairs of homology classes $i$ and $j$, without considering self-interactions of classes. The normalization factor in the denominator is chosen to be a linear function of the number of edges. The linear relation is suggested by the experimental observation that the number of edges in these TF-TF networks is comparable to the number of nodes. We introduced the square root normalization factor in order to compensate for the bias generated by the different size of classes $i$ and $j$ ¹¹1The important fact here is that the normalization factor should compare to the number of nodes, as does the total number of edges in the network, alternative choices such as $C_{i}+C_{j}$ or $sqrt(C_{i}C_{j})$ lead to the same results. .

Note that the definition the collapsed network absorbs the presence of crosstalks or feedback likely inherited by AR duplication, and this is reflected by $R$. To understand this, let us assume that the probability of a rewiring event is smaller than the probability of keeping inherited interactions. Then, homology classes, generated by duplication of ancestral genes, would tend to maintain their inherited regulations, thus measuring an high value of $R$. In other words, $R$ measures the extent to which feedback is retained within homology classes. This likely comes from AR duplication, or duplication of existing higher-order feedbacks between homologs.

To sum up, the larger is $R$ (compared to randomizations), the less likely rewiring has shaped the network, and the more the hierarchy between the homology classes is maintained. Evaluation of empirical data (Table 2) shows that in the bacterial data sets, the estimated rewiring is consistent with, or slightly lower than, typical randomized values. Perhaps unexpectedly, we find an even stronger signal for the yeast TF-TF network, where, as we have shown, the feedback is more common. Most of this feedback involves homolog TFs, and the network is likely to have been shaped by a small amount of rewiring, or by rearrangements that keep into account the homology relations of transcription factors. Note that, empirically, one alternative possibility is that a duplication of ARs with inherited crosstalk is followed by a subsequent split the homology class, so that the two homologs cease to be classified as such. This phenomenon would lead to a false positive for rewiring in our interpretation. However, this problem is less relevant in presence of a negative signal, such as the one we find. Another way to put the same observation, is that most of the feedback observed in yeast was shaped by duplication with retained crosstalks, rather than by rewiring.

An issue where we found notable evolutionary differences between the evolutionary transcriptional architecture of E. coli and S. cerevisiae is the duplication of homo- and heterodimeric TFs²²2This analysis was not possible for B. subtilis due to the lack of large-scale interaction data..

It is known from sparse observations that yeast tends to use more heterodimeric TF pairs than prokaryotes IYM+05 ; LPC+06 ; PLK+07 . A systematic analysis using large protein interaction datasets confirms this trend. In E. coli (out of 150 TFs) we find 21 homodimers, and only 5 heterodimer pairs, all formed with the histone-like protein HU, which has a rather special status. Though no systematic data is available, the common opinion is that it likely that these figures strongly underestimate the number of homodimers. A computational evaluation of the TFs available in the 3DComplex database LPC+06 is not able to enrich the sample. On the other hand, it shows that 37/42 TF entries in PDB are scored as homodimers, corroborating the common belief. Conversely, in yeast, where more systematic data is available, we find (out of 157 TFs) 45 homodimers, and 91 heterodimer pairs.

Inspection of the homology classes (Table 3 and Supplementary Figure) shows that in S. cerevisiae, heterodimeric TFs tend to cluster in homology classes, and also to co-occur with homodimers. We interpret this as a strong indication that heterodimers stem from duplications of other dimeric TFs. This is measured by evaluating the overrepresentation of number of homodimers and heterodimers in homology classes, and the product of homodimers and heterodimers in the same class. The same process seems to be forbidden in E. coli, possibly because of the same selective pressure erasing crosstalks from duplicate ARs. Due to insufficient data, we could not show that homodimeric TFs in his bacterium are likely to form classes of homodimers.

Finally, we proceed to evaluate the long and short path hierarchy conservation in evolution. The model predicts that they will be conserved, in absence of rewiring and multiple AR duplications. In particular, we expect to find many members of the same homology class falling in the same computational layer, defined by longest paths to a root or shortest paths. Our observations are reported in Table 4. In E. Coli both data sets indicate a strong tendency to have duplicates in the same layer, with both definitions of hierarchy. In B. subtilis we find a similar trend only for the long-path hierarchy (the signal is much weaker, possibly due to the small size of the sample). Finally, in yeast we also find evidence for hierarchy conservation all data sets. In the larger dataset, the long-path hierarchy was evaluated for the hierarchical component only.

A simple model of network evolution through duplication-divergence was considered. At each time step all the nodes of the graph are duplicated, while the number of edges rises fourfold. This happens for the following reason: for each edge connecting two original (old) nodes (old-old edge), duplication of interaction gives rise to edges between the two old nodes and the two duplicate nodes. The original old-old edge therefore generates the four old-old, old-new, new-old, new-new edges. Duplication of the graph is followed by erasing of edges with prescribed probabilities tesikirill ; CLJ+07 . One can formulate the model with partial or global duplications, and including or not the duplication and removal of self-edges (in this case, we considered the probability of retaining a self-edge equal to that of any other edge). The behavior of this model was compared, through a set of observables, with the observed trends of the experimental data.

We considered the following data sets for the transcription networks. For E. coli, the Shen-orr data-set and the larger and more recent RegulonDB5.5 SMM+02 ; SGP+06 ; SSG+06 . For B. subtilis, DBTBS MNO+04 . For S. cerevisiae, the Guelzim GBB+02 data-set and the more recent Balaji BBI+06 data set. The Balaji data set was modified to include auto regulating interactions taken from the literature. In the case of yeast, probably due to the much larger size of the more recent data, there is no compatibility between the two sets of data. Domain architecture data are taken from the SUPERFAMILY database GKH+01 ; WMV+07 , versions 1.61 and 1.69, as in the data sets in BLA+04 . Homodimers and Heterodimers forS. cerevisiae and E. coli respectively, were obtained from the SGD database CWB+04 and from ref. BPL+05 .

We used the leaf-removal algorithm CLJ+07 on the data-sets (including ARs) and their randomized counterparts. This algorithm prunes the input and output tree-like components of a directed network leaving with a “feedback core” of nodes, where each node is involved in at least one feedback loop. Each iteration of this pruning algorithm defines a hierarchical layer. The main observables we considered were the size of the core and the number of iterations to reach the core. For the evaluation of shortest paths, we used the Dijkstra algorithm Dij59 , considering the distribution of shortest-path lengths and the longest paths. The results for the empirical networks were compared to randomizations obtained using a standard Markov Chain Monte Carlo (MCMC) algorithm RJB96 or an Importance Sampling Monte Carlo algorithm FBJ+07 that preserve the degree sequence (marginals of the adjacency matrix). We considered random counterparts of the networks where the only constraints come from the in- and out- degree sequences. In particular, we chose not to conserve the number of self-regulators. Our previous work CLJ+07 ; FBJ+07 shows that in general this null model yields qualitatively different results. In particular, in the case of the Shen-Orr data set the number of layers falls in the average of the random ensemble that conserves self-regulations. For the scopes of this work, we did not consider this kind of randomization, and we hypothesized that, as in the model for the graph growth by duplication, auto-regulatory edges have the same status of other edges.

We used domain architectures from the SUPERFAMILY database to build protein architecture databases, one for each specie. We then constructed classes of homologous genes using similarity criteria between these architectures, as was done in TB04 . Two genes are considered homologs if they share the same domains in the same order, neglecting domain repeats. For this analysis, proteins coded by the same operon were considered as separate entities. Since the definition of homology is rather arbitrary it is rather natural to test different definition and observe the stability of results. All these results were filtered for consistency using stricter or looser homology criteria sellerio002 . In the case of the S. cerevisiae, and only with respect to the network degree distributions (as shown in Fig. 2), we also considered the notion of homology descending from the gene pairs defined by blocks of doubly conserved syntheny described in KBL04 .

We assessed the evolution of the Auto Regulated transcription factors by studying the distribution of the subgraphs shown in Fig. 3 within homology classes. Quantitative measure was given by the evaluation of observable $a$, as described in the main text, which was computed on both the experimental data set and on the randomized instances of the null model described later.

We studied the distribution of regulating interaction within and across homology classes. A number of observables were implemented to describe the relationship standing between the architecture data set and the transcription network data set. To quantify the effect of rewiring we introduced the observable $R$ as described in main text. Again, this observable was evaluated upon the experimental data set and upon the null model.

We studied the evolutionary properties of homodimers and heterodimers with respect to our duplication and divergence model. We assessed the co-existence of homo- and hetero- dimers within the homology classes by evaluating the sum (over all the homology classes) of the number of all dimeric pairs found in each homology class. The same analysis were performed upon the randomized instances of the null model, and the result is shown in Fig. 2.

Most of the measures performed upon the experimental data sets were compared with a null model of homology. Keeping fixed the homology classes, we randomly shuffled the architecture associations to the gene names. This randomizes all the interaction existing between homology properties (which are responsible of class generation) and all the other data sets (adjacency list of transcription network, homodimers, heterodimer pairs). The advantage of this null model lies in the fact that no experimental database is actually randomized, only their interactions. This is important as this randomization does not destroy the homology information nor the global (large scale) properties of the network. Finally, gene name shuffling was done separately for TFs and TGs, as was done in TB04 , due to their inherently different DNA-binding properties (which depend on their domains). The data shown in this paper correspond to $10^{5}$ randomized instances of the null model, allowing us to estimate P-values larger than $1\times 10^{-5}$.