Evolutionary Rate and Duplicability in the Arabidopsis thaliana Protein–Protein Interaction Network

Protein–Protein Interaction Data

The PIN was obtained by merging the A. thaliana networks available from the BioGRID database v3.1.81 (Stark et al. 2011) and from Arabidopsis Interactome Mapping Consortium (2011). Only physical interactions among pairs of A. thaliana proteins were considered. All interactions used in the current analysis had been determined experimentally in A. thaliana (i.e., the data set does not contain interactions inferred computationally or derived from other species).

Impact of Natural Selection

For each gene in the network, we attempted to identify its 1:1 ortholog in the A. lyrata genome (Hu et al. 2011) using a best reciprocal Basic Local Alignment Search Tool (BLAST) approach (using BLASTP and an E-value cut-off of 10⁻¹⁰). Each pair of A. thaliana–A. lyrata orthologous protein sequences was aligned using ProbCons 1.12 (Do et al. 2005), and the resulting alignments were used to guide the alignment of the corresponding coding sequences. Estimates of d_N, d_S, and ω were obtained using the one-ratio model M0 from the PAML 4.4 package (Yang 2007).

Paralogs Identification

Each A. thaliana protein was used as query in a BLASTP (Altschul et al. 1997) search against the A. thaliana proteome. Genes were classified as singleton if no hit was obtained with an E value < 0.1 or as duplicated if at least a homolog was found that met the following criteria: 1) E value ≤ 10⁻¹⁰; 2) the aligned region length (L) was ≥80% of the length of the query sequence; and 3) amino acid identity was ≥30% if L > 150 amino acids or 0.06 + 4.8L^{−0.32[1 + exp(−L/1,000)]} otherwise (Rost 1999). Duplicated genes were further classified as whole-genome duplication (WGD) genes if classified as such in Blanc et al. (2003); otherwise, they were classified as resulting from small-scale genome duplication (SSD).

Gene Expression Level and Breadth

Expression data from 79 A. thaliana tissues were obtained from Schmid et al. (2005). For each gene and tissue, values were averaged across the three replicates, and the gene was considered to be expressed if it was annotated as “present” in at least two of the replicates. For each gene, expression level was computed as the median across the 79 tissues, and expression breadth was computed as the number of tissues in which the gene is expressed. For genes matching multiple probe sets, the set yielding a higher expression level was used. Probe sets matching multiple genes were discarded.

Functional Information

Each A. thaliana gene was assigned to one (or sometimes a few) eukaryotic orthologous group (KOG) categories using the eggNOG v2 database (Muller et al. 2010).

Ribosomal Proteins

Proteins were considered to be ribosomal if identified as such in Barakat et al. (2001), if their description contained the text “ribosomal protein” or if they were assigned to the Gene Ontology (Ashburner et al. 2000) terms “large ribosomal subunit” or “small ribosomal subunit.”

Age of Genes

Each of the 5,789 A. thaliana network proteins was used as query in a BLASTP (Altschul et al. 1997) search against the nr database (obtained from the National Center for Biotechnology Information on January 2012). An E-value cut-off of 10⁻⁶ was used. Only hits whose aligned region length was ≥80% the length of the query sequence were considered. Genes presenting hits in prokaryotes or in nonplant eukaryotes (non-Embryophyta) were deemed as “ancient.” The remaining genes were considered to be “land plant-specific.”

Genes Acting at the Center of the A. thaliana PIN Are More Selectively Constrained Than Those Acting at the Network Periphery

We assembled a PIN for A. thaliana by merging all interactions available in the BioGRID database (Stark et al. 2011) and in Arabidopsis Interactome Mapping Consortium (2011). The resulting network consisted of 5,789 unique proteins connected by 14,368 physical binary interactions. For each protein in the network, we computed three centrality measures: 1) the number of interactors (degree); 2) the number of shortest paths between all pairs of proteins of which the protein is part (betweenness; Freeman 1977); and 3) the inverse of the average shortest distance to all the other proteins in the network (closeness). Using a best reciprocal BLAST approach, we found 1:1 A. lyrata orthologs for 3,868 of the 5,789 A. thaliana network genes. For each pair of orthologs, the impact of natural selection was inferred from the nonsynonymous to synonymous divergence ratio (ω = d_N/d_S). Values of ω = 1 are expected for genes evolving neutrally, whereas ω < 1 is indicative of the action of purifying selection preserving the sequence of the encoded proteins, and ω > 1 in a number of codons is indicative of the action of positive selection (adaptive evolution) driving the fixation of nonsynonymous substitutions. The estimated ω values exhibit a median of 0.128, indicating that these genes evolved under relatively high levels of selective constraint.

We found that the 3,868 genes encoding proteins that are represented in the network (i.e., those with described interactors) exhibit lower rates of evolution than those not represented in the network (median ω values of 0.128 and 0.183, respectively; Mann–Whitney test, P < 10⁻¹⁵). In addition, among these 3,868 network genes, ω values negatively correlate with betweenness (Spearman’s rank correlation coefficient, ρ = −0.053, P = 0.001; fig. 1) and closeness (ρ = −0.034, P = 0.035): the larger the values of betweenness and closeness for a protein, the stronger are the selective constraints acting on that protein. Although the correlation is only marginally significant for degree (ρ = −0.030, P = 0.063), the 1,467 genes with degree >1 show significantly lower ω values than the 2,401 genes with degree = 1 (median ω values of 0.125 and 0.135, respectively; Mann–Whitney test, P = 0.005). Taken together, these results indicate that the most central genes in the A. thaliana network are subject to stronger levels of purifying selection than those acting at the network periphery. Similar results were obtained for d_N but not for d_S (table 1), indicating that the distribution across the network of selection at the amino acid level is the main responsible for the observed trend.

The Relationship between Centrality and Evolutionary Rate Is Independent of Gene Expression Level and Breadth, Gene Function, and the Length of the Encoded Products

Having established an association between proteins’ centralities and their rates of evolution, we sought to determine whether this association reflected a direct effect of network position or, on the contrary, it was the result of the distribution across the network of a number of factors that correlate both with rates of evolution and network centralities.

Levels of selective constraint acting on a gene are known to depend on a number of factors, among which the level and breadth of gene expression seem to be the most important (Duret and Mouchiroud 2000; Pál et al. 2001; Krylov et al. 2003; Subramanian and Kumar 2004; Wright et al. 2004; Drummond et al. 2005, 2006; Ingvarsson 2007; Slotte et al. 2011; Yang and Gaut 2011). In agreement with previous reports, we observed that the ω and d_N values exhibit a strong negative correlation with both expression level and breadth in our data set (table 1). Additionally, we observed that, similar to the PINs of other organisms (Bloom and Adami 2003; Lemos et al. 2005), the most central genes in the Arabidopsis PIN are more highly and broadly expressed than those acting at the periphery: expression levels positively correlate with betweenness and closeness and expression breadths with betweenness (table 1). Combined, these observations raise the possibility that the correlation between centrality measures and evolutionary rates could be driven by the higher levels and/or breadths of expression of genes acting at the center of the PIN. However, partial correlation analysis (supplementary table S1, Supplementary Material online) shows that the association between ω and betweenness is independent of expression level and breadth (ρ = −0.038, P = 0.026).

Another potential confounding factor is protein length, as it correlates positively with ω (i.e., genes encoding shorter proteins tend to be more selectively constrained; Subramanian and Kumar 2004; Ingvarsson 2007; table 1) and negatively with closeness (i.e., genes acting at the center of the network tend to encode shorter proteins, in agreement with previous observations in Drosophila; Lemos et al. 2005; table 1). However, the correlation between ω and both betweenness (ρ = −0.052, P = 0.001) and closeness (ρ = −0.032, P = 0.048) remains significant when protein lengths are controlled for (supplementary table S1, Supplementary Material online), indicating that the tendency for central genes to evolve under strong levels of selective constraint is also independent of protein length. Furthermore, the correlation between ω and betweenness remains significant when protein length and expression level and breadth are simultaneously controlled for (ρ = −0.037, P = 0.028; supplementary table S1, Supplementary Material online).

Genes performing different functions are subject to different levels of selective constraint (e.g., Castillo-Davis et al. 2004; Alvarez-Ponce and McInerney 2011) and encode proteins with different network centralities (Kunin et al. 2004; Cotton and McInerney 2010; Alvarez-Ponce and McInerney 2011). Consistently, we found that genes in different categories exhibit different ω and centrality values (Kruskal–Wallis test, ω: P < 10⁻¹⁵; degree: P < 10⁻¹⁵; betweenness: P = 1.17 × 10⁻¹¹; closeness: P = 7.32 × 10⁻¹⁰). Therefore, the observed correlation between centrality and evolutionary rate (table 1) could conceivably be the result of these differences. To discard this possibility, the correlation between evolutionary rate and degree was evaluated separately for genes involved in each of the KOG functional categories (Tatusov et al. 2003) (similar results were obtained for betweenness and closeness; data not shown). The analysis was restricted to the 20 KOG categories comprising more than 25 A. thaliana network genes. Remarkably, the correlation coefficient was negative for 18 of the 20 categories (table 2). The correlation was significant for only five of these categories, probably as a result of the reduced statistical power resulting from partitioning the data set. These observations indicate that, in spite of the fact that genes with different functions exhibit different ω and centrality values, the negative correlation between centrality and ω is largely independent of gene function.

Remarkably, correlation coefficients are highly variable across the different KOG categories (table 2). First, pairwise comparison shows that 10 pairs of categories exhibit significantly different correlation coefficients (supplementary table S2, Supplementary Material online). Second, for each functional category, we compared the correlation coefficient for genes within that category with the correlation coefficient for all other genes (i.e., all network genes not belonging to that category); this analysis shows that categories “signal transduction mechanisms” (P = 0.009) and “chromatin structure and dynamics” (P = 0.025) exhibit significantly higher correlation coefficients than the rest of the network (table 2). Taken together, these results indicate that the extent to which sequence evolution is affected by protein centrality is dependent on gene function. In fact, some categories exhibit strong negative ω-degree correlations that are comparable in strength to (or even stronger than) the ω-expression level and the ω-expression breadth correlations (tables 1 and 2).

Finally, because genes encoding ribosomal proteins are highly conserved, and presumably highly connected, the association between centrality and evolutionary rate observed here could potentially be due to these genes. However, when these genes are eliminated, the correlation between ω and betweenness remains significant (ρ = −0.062, P = 1.34 × 10⁻⁴) and the correlation between ω and degree—which was not significant in the complete data set (table 1)—reaches significance (ρ = −0.041, P = 0.011). Furthermore, these correlations remain significant when expression level, expression breadth, and protein length are simultaneously controlled for (degree: ρ = −0.038, P = 0.029; betweenness: ρ = −0.042, P = 0.014).

Taken together, these observations point to a direct effect of the position of proteins in the A. thaliana PIN on the rates of evolution of their encoding genes.

Genes Encoding Physically Interacting Proteins Are Subject to Similar Levels of Selective Constraint

We considered whether interacting proteins in the A. thaliana network are subject to relatively similar levels of selective constraint compared with random protein pairs. For that purpose, we computed the normalized absolute difference between the evolutionary rates of all pairs of interacting genes in the network (X):

Here, m is the number of interactions in the network, and ω_i1 and ω_i2 are the ω values of the two genes involved in interaction i. Self-interactions (i.e., interactions among proteins encoded by the same gene, a total of 383 in the data set) were not considered in this analysis, as PINs are enriched in such interactions, which can inflate the average similarity between interacting genes (Ispolatov et al. 2005; Alvarez-Ponce and McInerney 2011). The observed value (X = 0.892) was compared with a null distribution obtained from a collection of random networks with the same proteins, number of interactions, and degree for each node, which we generated from the original network by repeatedly switching pairs of edges (as in Luisi et al. 2012). Out of 10,000 random networks, only 426 showed an X value lower or equal to the observed one, indicating that interacting proteins exhibit rates of evolution that are more similar than expected from a random network (P = 0.0426; fig. 3). Random networks exhibit an average X value of 0.900, indicating that the ω values for interacting genes are 0.8% more similar than random pairs of proteins.

Ribosomal proteins might also represent a source of confounding bias, as they are subject to strong levels of purifying selection and are highly connected to each other. However, significant results were obtained when ribosomal proteins were also removed from the network (X = 0.889, P = 0.039).

Centrality and Duplicability in the A. thaliana PIN

Finally, we studied the effect of proteins’ network position on the duplicability of their encoding genes. Among the 5,789 genes encoding the A. thaliana PIN, 883 were found to be singleton, and 3,532 were deemed as duplicated based on similarity searches (the rest of the genes remained unclassified; see Materials and Methods). Among duplicated genes, 1,540 are the result of one of the WGD events that took place in the Arabidopsis lineage (Blanc et al. 2003; De Bodt et al. 2005), and 1,992 were deemed as the result of SSD events. Unexpectedly, proteins encoded by duplicated genes have more interactors in the network (average degree of 5.18) than those encoded by singleton genes (average degree of 4.10) (Mann–Whitney test, P = 0.008; fig. 2). These differences remain significant when self-interactions and interactions among proteins encoded by paralogs are removed from the analysis (P = 0.049), or when genes are classified as duplicated or singleton based on whether they present annotated paralogs in the Ensembl plants database release 14 (Kersey et al. 2012) (P = 1.62 × 10⁻⁴). Furthermore, degree positively correlates with the number of paralogs annotated in Ensembl plants (ρ = 0.089, P = 1.02 × 10⁻¹¹), even when only duplicated genes are considered (ρ = 0.086, P = 5.85 × 10⁻⁹). This observation is in agreement with the patterns observed in the human interactome (Liang and Li 2007; D'Antonio and Ciccarelli 2011; Doherty et al. 2012) but in contrast with those observed in E. coli, yeast, worm, and fly, whose duplicated genes tend to encode lowly connected proteins compared with singleton genes (Hughes and Friedman 2005; Prachumwat and Li 2006; Makino et al. 2009; D'Antonio and Ciccarelli 2011). Among duplicated genes, those resulting from WGDs are more highly connected in the Arabidopsis network than those resulting from SSD events (average degree of 5.70 and 4.78, respectively; Mann–Whitney test, P = 5.63 × 10⁻⁶; fig. 2).

D’Antonio and Ciccarelli (2011) found that the relationship between centrality and duplicability was different among ancient and metazoan-specific human genes: among genes of premetazoan origin, duplicated genes are less central than singleton genes, whereas among genes that are specific to metazoans, duplicated genes tend to occupy more central positions than singleton genes. We considered whether the relationship between centrality and duplicability observed in Arabidopsis also depended on the age of the genes. For that purpose, genes were classified as “ancient” (if they had homologs in prokaryotes or in nonplant eukaryotes) or as “land plant-specific” (otherwise). A total of 3,114 network genes were classified as ancient, and the remaining ones were deemed as plant specific. We found that, among plant-specific proteins, those encoded by duplicated genes are more highly connected than those encoded by singleton genes (average degrees of 5.39 and 3.85, respectively; Mann–Whitney test, P = 4.73 × 10⁻⁶). This difference, however, is not significant among ancient genes (average degrees for proteins encoded by duplicated and singleton genes: 5.06 and 4.41 interactions, respectively; P = 0.484). Therefore, the relationship between centrality and duplicability observed in Arabidopsis seems to be specific to plant-specific genes.

Because duplicated Arabidopsis genes are more selectively constrained than singleton genes (Yang and Gaut 2011), the lower evolutionary rates of genes central to the A. thaliana PIN could potentially be a by-product of their enrichment in duplicated genes. However, the correlation between ω and both betweenness (ρ = −0.050, P = 0.015) and closeness (ρ = −0.049, P = 0.016) remains significant when only duplicated genes are considered. Although these correlations are not significant for singleton genes (ρ = −0.044, P = 0.288 for betweenness; ρ = 0.043, P = 0.293 for closeness), this might result from the small number of genes in this category (n = 883).

Analysis of a High-Quality Subnetwork Provides Consistent Results

Currently available interactomes are the result, to a high extent, of the application of high-throughput techniques of protein–protein interaction discovery. As a result, interactomic data sets are subject to high rates of false positives and negatives (Bader et al. 2004; Deeds et al. 2006). Given the potential that this could be affecting our observations, we repeated our analyses in a high-quality subset of the A. thaliana interactome, containing only reliable interactions. Similar to D'Antonio and Ciccarelli (2011), we filtered our data set, retaining only interactions that had been determined using low-throughput techniques (i.e., more accurate analysis of protein interactions on a one-by-one basis), and those that had been identified by two or more high-throughput analyses independently. The filtered network consisted of 4,808 interactions connecting 2,798 proteins, out of which 1,842 are encoded by genes with 1:1 orthologs in A. lyrata. This represents a dramatic decrease in the amount of data when compared with the full data set (only 48.3% of the proteins and 33.5% of the interactions were present in the high-quality subnetwork), which potentially involves a decrease in statistical power but also an increase in the signal-to-noise ratio.

Although the correlation between degree and ω is not significant in this reduced data set (ρ = 0.003, P = 0.907), ω values significantly correlate with both betweenness and closeness, with correlation coefficients that are higher in magnitude than those observed in the entire data set (betweenness: ρ = −0.061, P = 0.008; closeness: ρ = −0.108, P = 3.22 × 10⁻⁶). The correlation between ω and closeness remains significant when the effects of expression level and breadth, and protein length, are simultaneously controlled for (ρ = −0.061, P = 0.013). Furthermore, genes encoding proteins that are represented in the network are more selectively constrained than those that are not represented in the network (median ω values of 0.120 and 0.176, respectively; Mann–Whitney test, P < 10⁻¹⁵).

Consistent with the observations in the full interactome, duplicated genes are more highly connected than singleton genes (average degrees of 3.58 and 2.83, respectively) in the high-quality subnetwork. Although the degrees of both groups are not significantly different (Mann–Whitney test, P = 0.132), the duplicated/singleton degree ratio is equivalent to that observed in the full data set (1.27), suggesting that the lack of significance in the analysis of the high-quality subnetwork may result from the reduced statistical power resulting from trimming the data set.