Genes Encoding Intracellular Signaling Proteins in Animals Originated Along with Metazoa and Chordata: Chance or Necessity?

Abstract

In this work, we investigate whether the construction of signaling pathways during evolution follows a deterministic law through a study of the eventual link between age of appearance in the tree of life and position in the signaling pathway of genes involved in these pathways. We use the 47 human signaling pathways described in the Kyoto Encyclopedia of Genes and Genomes and investigate the orthologs of these genes in 315 animal species plus a yeast taxon, representing 15 large clades. Many genes appear on two key branches: those between the last common ancestor of Opisthokonta and Metazoa and between Deuterostomia and Chordata. We look for a link between the age of appearance of an upstream A gene and that of its downstream B partner. We observe that for all the interactions of two partners, only 20.6% of the corresponding genes arose simultaneously in the tree of life, 40.7% being called “backward” (i.e. B appearing before A) and 38.7% “forward” (A appearing before B). For 16 of the 47 pathways, there is a positive correlation between the age rank difference between interacting partner genes and the position of the corresponding proteins in the pathway: the more upstream a protein is involved in the pathway, the greater the rank difference is (the correlation, positive or negative, is not significant for 30 pathways). For the sole insulin signaling pathway, this correlation is negative. Moreover, by permutation test, we find that 14 of the 47 observed pathway contained larger modules (subset respecting a homogeneous appearance pattern) than expected by chance alone. Finally, for 20 of the 47 pathways, the construction scenario appears to be random, as these pathways do not validate any of our statistical tests (permutation tests on interaction direction and module sizes as well as correlation test on pathway position and age rank). Given that only 14.9% of the tests are significant and that significant effects are different among pathways, we conclude that there is no deterministic rule in the establishment of the pathways herein studied or that the patterns have been obscured by subsequent transformations.

evolution, genome, signaling pathways

Significance

This work shows that the majority of genes coding for intracellular signaling proteins appeared early and are associated with the origin of Metazoa and Chordata. It also shows that for 20 of the 47 signaling pathways, there is no significant correlation between the age rank difference between interacting partner genes and the position of the corresponding protein in the signaling pathway. Surprisingly, the age of appearance of the genes that code proteins in these pathways appears to have occurred randomly, rather than in an order reflecting position in the proteins in the pathway.

Introduction

The co-evolution of genes encoding interacting molecules has interested many scientists (Fraser et al. 2004; Rand et al. 2004; Andreani and Guerois 2014; Lynch and Hagner 2015) because of the intriguing question of the modes of mutation and selection that act on two molecules simultaneously. In particular, co-evolution of the binding motif has been well investigated (Lewis et al. 2010), focusing on the fitness (Williams et al. 2001; Bloom et al. 2004), conservation of interaction (Wuchty et al. 2003; Lovell and Robertson 2010; Kachroo et al. 2015), or evolution of the residues at the interface of the molecules (Mintseris and Weng 2005; Jack et al. 2016; Echave and Wilke 2017). These studies on the co-evolution of binding partners often require the integration of different disciplines (chemistry, evolutionary and molecular biology, etc.). Paradoxically, establishment of the interaction in a phylogenetic context, which involves fewer fields, is less studied. Little is known, for example, about the age of appearance and evolution of the different partners prior to their first interaction. Does the emergence of one partner favor the emergence of the second partner, and if so, in which order?

In the case of interacting molecules, the appearance of genes coding for molecules included in a complex is intricate (Kauffman 1993). For two molecules that will eventually interact, the appearance of one may depend on the appearance and conservation of the other. This may be the case, for example, when the presence of the first molecule is not advantageous if its partner has not yet appeared.

The existence of interacting proteins without partners (“orphan protein”) has been described frequently (Howard et al. 2001), even though it is sometimes difficult to assess whether an interacting protein is a true orphan or whether its ligand is simply unknown (Benoit et al. 2006). The relative order of appearance of genes encoding protein partners is thus an open question. Furthermore, several types of interactions can be observed in living organisms, with different numbers of interacting partners (Koshland 1958; Albelda and Buck 1990; Maslov and Sneppen 2002; Sullivan and Holyoak 2008), varying affinities (Kent et al. 1980), or different durations for the interactions (Nooren and Thornton 2003) through evolution, making the problem even more complex.

Concerning the pairs of ligands and their receptor(s), Thornton (2001) has shown that the first steroid receptor, present in lamprey and believed to be present in the last common ancestor of vertebrates, was an estrogen receptor and that several duplications led to other steroid receptors specialized in other functions with other ligands. However, more recent investigations suggested that the ancestral ligand for the ancestral steroid receptor was a molecule with a structure distinct from modern estrogen (an aromatized steroid with a side chain) called paraestrol A (Markov et al. 2017). In a previous work, we found that 41% of the membrane receptors and their respective first ligands appeared on the same branch of the evolutionary tree, which is 2.5-fold more often than expected by chance, thus suggesting an evolutionary dynamic of interdependence and conservation between these partners (Grandchamp and Monget 2018). In contrast, 21% of the receptors appeared after their ligand (i.e. 3-fold less often than expected by chance) and 38% appeared before their first ligand (i.e. as much as expected by chance). These results suggest that selective pressure is exerted on ligands and receptors once they appear, which would eliminate molecules whose partner does not appear quickly.

Using both descriptive and hypothesis testing approaches, we investigate the age rank of appearance of genes encoding proteins involved in signaling pathways downstream of membrane receptors in the animal tree of life (Fig. 1). We also study the relationship between the position of the protein in the signaling pathway relative to its partners, and their relative age rank of appearance in the tree of life. In other words, does the order of appearance of genes in the tree of life reflect that of their coded protein in the signaling pathway (upstream or downstream of the pathway), with a direct, reverse or random order (Fig. 2)? We also test whether the size of modules, defined as subsets respecting a homogeneous appearance pattern, departs from randomness.

Fig. 1.

Simplified animal tree of life and sampled clades. The rectangles represent the branches on which proteins appeared, with estimated times (ages of the nodes to the right of the rectangles) and can include one to several lower-ranking taxa included in this study. The circles represent the terminal taxa (A–Y) included in this study; these are sorted from the most distantly related to the most closely related to humans. The nodes available on Ensembl Metazoa (or Genomicus Metazoa) appear in italics, and the nodes available on Ensembl Vertebrata (or Genomicus Vertebrata) appear in bold.

Fig. 2.

Three hypotheses about the possible scenarios (A–C) linking the age rank of appearance of a gene to the pathway position of the corresponding protein. The colors represent the age of genes (dated by their appearances on branches below internal nodes 1–15). We arbitrarily chose 10 proteins/positions per pathway for the illustration. In scenario A (“backward”), the pattern of protein activation in the pathway represents the reverse of the pattern of appearance of genes in the tree. Scenario B (“forward”) is characterized by the opposite relationship (the pattern between protein position in the pathway matches the appearance of the genes coding these proteins). In scenario C (“random”), there is no link between the position of the protein in the pathway and the pattern of appearance of the corresponding gene in evolution.

Three mutually exclusive hypotheses can be formulated about the possible relationships between gene position in a pathway and its relative age of origin (which we call below “scenarios”), based on principles of molecular biology, evolution, and paleontology. The first is that the genes expressed downstream in the sequence appeared first because this is where the interactions with the genome occur, in the nucleus (Fig. 2A). The second hypothesis, which is diametrically opposed, is that the genes expressed upstream appeared first because the first life forms most likely had an outer cell membrane, but they almost certainly lacked a cell nucleus. Indeed, recent paleontological research suggests that crown eukaryotes (hence, potentially, the nucleus) may have appeared later than usually thought (Porter, 2020), and this may add support for this hypothesis because a greater proportion of the evolution of life may have occurred prior to the appearance of the nucleus (Figure 2B). Also, the outer cell membrane is where interactions with the environment occur. The third hypothesis, to an extent complementary with the first two, is that some genes are recruited in a pathway (an example of exaptation) randomly, and their age of appearance bears no logical relationship with the age of other genes in the cascade (Figure 2C). With time, this phenomenon, a kind of evolutionary tinkering (Jacob, 1977), could erase the order originally created by successive gene appearance following a backward or forward rule. The signal loss generated by this process is vaguely reminiscent of the well-known saturation in phylogenetic signal in gene sequences that makes fast-evolving genes unsuitable to infer ancient divergence events, which was evidenced long ago (Hillis and Huelsenbeck, 1992).

Results

Age of Appearance of Genes Encoding Proteins of Signaling Pathways

Our results, based on 2,298 unique genes involved in 47 human intracellular signaling pathways available in KEGG database V104.0 (see Materials and Methods), reveal interesting patterns on the distribution of ages of appearance for the genes encoding proteins involved in animal signaling pathways. Two parts of the tree of life are overrepresented in gene origins (Fig. 3): at the base of the tree, on the branch leading to Metazoa (before 765 mya), and at the base of Chordata, which includes Cephalochordata, Tunicata, and Vertebrata, around 635 to 558 mya (Kumar et al. 2022). Overall, in cumulative data, 29% of genes present in humans from our dataset had appeared before the origin of Metazoa, 45% before diversification of Deuterostomia, 72% before Chordata, and 92% before the origin of Vertebrata.

Fig. 3.

Cumulative distribution of genes encoding proteins involved in the 47 signaling pathways. The colored bars represent the cumulative percentage of genes involved in a signaling pathway progressively appearing during evolution. The light gray bars represent the raw data.

The temporal relationship of gene appearances among interacting partners in a pathway also reveals interesting patterns. Of the 2,666 interactions studied, 20.6% of interactions are “simultaneous” (i.e. both partners in the interaction appeared on the same branch), ranging from 9% for the FoxO and chemokine pathway to 35% for the VEGF pathway. Therefore, 79.4% are asynchronous, of which 40.7% are “backward”, with partner B (downstream in the pathway) appearing first in the tree of life (ranging from 18% for toll-like receptor pathways to 79% for the ovarian steroidogenesis pathway) and 38.7% are “forward,” with partner A (upstream in the pathway) appearing first (ranging from 11% for ovarian steroidogenesis to 65% for IL-17; Fig. 4). Some pathways have a greater proportion of forward relationships, such as the toll-like receptor, IL-17, and T-cell receptor pathways; it is the opposite for other pathways, such as the PPAR, cAMP, and ovarian steroidogenesis pathways (Fig. 4).

Fig. 4.

Forward, backward, or simultaneous interaction for each signaling pathway. Each horizontal column represents a signaling pathway (described on the left). On each column, for an A × B interaction, the red color corresponds to a backward interaction (A was born after B), the orange color corresponds to a simultaneous interaction (A was born at the same age as B), and the green color corresponds to a forward interaction (A was born before B).

The age rank difference or delta (i.e. the age rank of appearance of the first interactant and that of the second interactant) is highly variable among the 47 pathways (supplementary Data S1, Supplementary Material online). For eight pathways (at the top of supplementary Data S1, Supplementary Material online), the delta median is > 0 (between 1 and 4; more backward than forward); for nine pathways, the median is < 0 (between −2.5 and −0.5; more forward than backward); and for 30 pathways, the median is 0. Moreover, some pathways have specific profiles. More precisely, some pathways exhibit a particularly eccentric median and biased distribution compared with 0, such as the prolactin pathway (min = −13, Q1 = −7, median = −2.5, Q3 = 2, and max = 8) and the ovarian steroidogenesis pathway (min = −2, Q1 = 0, median = 4, Q3 = 7, and max = 12), where Q1 and Q3 are the first and third quantiles of a dataset (supplementary Data S1, Supplementary Material online, columns Q–U in supplementary Data S4, Supplementary Material online).

For each pathway, the delta distribution in relation to the age rank of appearance of the first interactant displays surprising patterns (supplementary Data S2 and S3, Supplementary Material online). Within the Hippo pathway, for example, several proteins appearing in either Opisthokonta or Metazoa (in blue) have “waited” a long time for their partner (delta = −13), whereas proteins appearing with Vertebrata (in orange) appeared after their partner. For this Hippo pathway, more than half of partners (that appeared on a different branch) present a delta of −6 to 13 (21.1% of delta = 0; supplementary Data S2, Supplementary Material online). For the ovarian steroidogenesis and cAMP pathways, most of the deltas are > 0 (median = 4 and 1, respectively) and the first partner of the interaction appeared after their direct downstream partner, the interaction being A (position n in the pathway) → B (position n+ 1) (supplementary Data S2, Supplementary Material online). All the pathways are represented in supplementary Data S2, Supplementary Material online and supplementary Data S3, Supplementary Material online. These examples suggest that these pathways settled with genes that have been replaced since.

Relationship Between the Delta of Age Rank of Appearance of the Gene Encoding a Partner Involved in Each Pathway and the Upstream or Downstream Position in the Pathway

For 16 of the 47 pathways studied here (from the AMPK to mTOR pathway; Table 1), there is a positive correlation (P ranges from 2.63 × 10⁻¹¹ to 0.04) between the delta of a gene pair and the position of the corresponding protein in the pathway, which means that the more upstream a protein is involved in the pathway, the older is the upstream gene compared with downstream gene of the pair (hypothesis B of Fig. 2, Table 1, and Fig. 5). These include the NOD-like pathway, among others. For a single pathway (the insulin signaling pathway), there is a negative correlation (r = −0.18, P = 0.03), which means that the more upstream a protein is involved in the pathway, the younger is the activator compared with the downstream gene of the pair (hypothesis A of Fig. 2 and Table 1). For 30 pathways, such as the Hippo and Notch pathways, we do not find any significant correlation between the age of appearance and the position in the pathway (Table 1).

Fig. 5.

Relative age (rank) of appearance of genes encoding proteins involved in each pathway according to the pathway position of the protein. Each row corresponds to a pathway. The upstream proteins of the pathway (close to the cell membrane) are on the left, and the downstream proteins of the pathway (close to the nucleus) are on the right. Depending on the signaling pathway, there are 3–16 positions on the abscissa. For a given pathway position, each rectangle represents the distribution of age rank of appearance for each of the proteins occupying the pathway position. For example, for the GnRH pathway at position 1, all the proteins appeared at the base of Osteichthyes. The pathways are sorted by correlation between pathway position and delta age rank (Table 1).

Table 1

Correlation between rank age differences of appearance and pathway position

path	categories	nb interacts	nb genes	nb subpath	max subpath	r	pvalue
AMPK	Signal transduction	109	120	184	10	0.419	9.49E−06*
ErbB	Signal transduction	114	90	309	10	0.381	2.88E−05*
JAK-STAT	Signal transduction	381	183	12,958	12	0.333	2.63E−11*
Adipocytokine	Endocrine system	72	68	54	8	0.304	5.54E−03*
cAMP	Signal transduction	203	206	1966	11	0.255	1.47E−04*
Wnt	Signal transduction	129	126	402	11	0.246	4.91E−03*
HIF-1	Signal transduction	72	77	54	7	0.245	4.78E−02*
Oxytocin	Endocrine system	79	90	51	9	0.243	2.59E−02*
AGE-RAGE	Endocrine and metabolic disease	220	146	358	6	0.202	2.11E−03*
NOD-like receptor	Immune system	210	208	420	9	0.200	3.71E−03*
IL-17	Immune system	176	168	709	9	0.158	3.63E−02*
MAPK	Signal transduction	410	354	2083	11	0.154	1.82E−03*
T cell receptor	Immune system	216	160	376	9	0.153	2.46E−02*
Toll-like receptor	Immune system	200	174	371	10	0.151	3.33E−02*
Neurotrophin	Nervous system	231	193	952	16	0.137	3.81E−02*
mTOR	Signal transduction	293	282	375	14	0.136	2.12E−02*
Insulin	Endocrine system	133	130	255	11	−0.184	2.95E−02*
Relaxin	Endocrine system	115	131	140	13	0.186	5.02E−02
Prolactin	Endocrine system	123	131	63	14	−0.165	6.63E−02
FoxO	Signal transduction	74	78	70	9	0.216	6.64E−02
Ovarian steroidogenesis	Endocrine system	19	29	15	6	−0.323	2.27E−01
cGMP-PKG	Signal transduction	60	79	201	12	−0.138	2.46E−01
Hedgehog	Signal transduction	51	64	41	5	−0.161	2.65E−01
Notch	Signal transduction	35	34	64	4	0.174	3.17E−01
Calcium	Signal transduction	40	59	77	8	−0.156	3.40E−01
Phospholipase D	Signal transduction	55	64	235	11	0.119	3.49E−01
Apelin	Signal transduction	79	86	26	11	0.104	3.88E−01
PPAR	Endocrine system	58	50	61	2	0.110	4.30E−01
NF-kappa B	Signal transduction	122	124	110	6	−0.063	4.91E−01
GnRH	Endocrine system	56	57	16	16	0.088	5.19E−01
Chemokine	Immune system	91	84	183	11	−0.058	5.76E−01
B cell receptor	Immune system	109	104	88	13	−0.054	5.90E−01
PI3K–Akt	Signal transduction	264	267	668	13	−0.033	5.97E−01
Ras	Signal transduction	115	112	232	9	−0.045	6.43E−01
Thyroid hormone	Endocrine system	113	121	102	7	−0.044	6.53E−01
RIG-I-like receptor	Immune system	98	88	374	8	0.039	7.03E−01
TNF	Signal transduction	64	70	27	8	0.048	7.04E−01
Fc epsilon RI	Immune system	49	49	42	10	−0.037	7.81E−01
Rap1	Signal transduction	89	90	108	5	−0.025	8.12E−01
Glucagon	Endocrine system	65	63	40	9	0.029	8.27E−01
Hippo	Signal transduction	90	106	73	6	0.022	8.39E−01
VEGF	Signal transduction	40	41	19	10	0.033	8.42E−01
p53	Cell growth and death	72	71	254	7	−0.022	8.52E−01
C-type lectin receptor	Immune system	159	174	107	8	−0.014	8.59E−01
TGF-beta	Signal transduction	85	93	61	7	−0.011	9.21E−01
Estrogen	Endocrine system	58	69	29	16	−0.010	9.47E−01
Sphingolipid	Signal transduction	56	62	51	6	0.002	9.86E-01

The correlation index (r) quantifies the strength of the relationship: strong between [−1, −0.5] and [0.5, 1], low between [−0.5, 0] and [0, 0.5], and no correlation at r = 0. Significant P-values (*P < 0.05) are presented in bold italics.

Number of Appearance Pattern Occurrences

Among the 47 pathways, we observe that six pathways have significantly more backward pairs than expected by chance alone, while only three are more forward and nine have more simultaneous relationships (Table 2).

Table 2

Summary of the statistical tests about appearance patterns

path	forward	backward	simult	sizeeffct_f	pvalue_f	sizeeffct_b	pvalue_b	sizeeffct_s	pvalue_s
Ovarian steroidogenesis	2	15	2	−5.681	6.00E−03*	7.426	1.00E−03*	−1.745	1.27E−01
JAK-STAT	218	84	79	73.352	6.00E−03*	−59.607	2.50E−02*	−13.745	4.02E−01
IL-17	114	38	24	37.512	1.90E−02*	−36.842	3.60E−02*	−0.670	8.11E−01
Toll-like receptor	124	37	39	34.972	2.80E−02*	−39.441	2.70E−02*	4.469	5.34E−01
VEGF	10	16	14	−6.800	4.70E−02*	−0.636	8.10E−01	7.436	6.00E−03*
RIG-I-like receptor	48	21	29	5.292	4.20E−01	−19.746	3.00E−03*	14.454	3.00E−03*
cAMP	48	133	22	−43.153	7.40E−02	56.890	4.00E−03*	−11.964	3.39E−01
ErbB	33	63	18	−11.158	8.80E−02	18.942	5.00E−03*	−7.784	1.42E−01
NF-kappa B	42	65	15	−9.688	2.01E−01	13.781	4.30E−02*	−3.558	4.00E−01
NOD-like receptor	109	63	38	12.062	3.16E−01	−22.191	4.70E−02*	10.129	1.02E−01
Wnt	37	41	51	−11.751	8.90E−02	−7.755	3.03E−01	19.506	2.00E−03*
TNF	22	21	21	−6.466	8.70E−02	−2.293	5.20E−01	8.759	1.10E−02*
MAPK	151	129	130	−13.537	6.53E−01	−35.085	2.35E−01	48.622	1.20E−02*
B cell receptor	38	33	38	−4.137	6.03E−01	−8.395	3.24E−01	20.834	1.60E−02*
Sphingolipid	16	21	19	−6.023	8.20E−02	−1.285	6.84E−01	7.308	2.00E−02*
C-type lectin receptor	53	68	38	−13.284	7.60E−02	1.824	7.49E−01	11.460	2.10E−02*
p53	36	16	20	5.444	5.33E−01	−14.261	2.12E−01	8.817	3.60E−02*
Hedgehog	18	17	16	−2.569	4.00E−01	−3.738	2.67E−01	6.307	3.60E−02*
AMPK	55	30	24	9.078	4.08E−01	−18.355	1.14E−01	9.277	8.30E−02
Ras	50	38	27	2.420	6.90E−01	−10.048	1.33E−01	7.628	8.70E−02
TGF-beta	35	31	19	5.864	1.73E−01	1.818	6.23E−01	−7.682	1.07E−01
FoxO	30	37	7	−1.370	8.90E−01	5.721	5.37E−01	−4.351	1.84E−01
Phospholipase D	21	24	10	0.160	9.22E−01	5.144	1.93E−01	−4.358	2.03E−01
Oxytocin	41	31	7	5.939	4.21E−01	−0.968	8.87E−01	−4.971	2.09E−01
Estrogen	23	24	11	−1.895	5.61E−01	−0.867	7.35E−01	2.762	2.19E−01
HIF-1	17	40	15	−13.277	1.52E−01	9.380	4.01E−01	3.897	2.35E−01
cGMP-PKG	22	26	12	−5.159	3.83E−01	2.518	7.05E−01	2.641	3.17E−01
Apelin	38	27	14	6.393	2.18E−01	−2.562	4.99E−01	−3.831	3.22E−01
Relaxin	46	55	14	−6.445	3.12E−01	10.315	1.34E−01	−3.870	3.62E−01
PI3K–Akt	120	114	30	11.032	5.88E−01	4.054	8.56E−01	−10.684	4.19E−01
Thyroid hormone	46	45	22	−1.259	8.82E−01	−1.856	8.13E−01	3.115	4.29E−01
Insulin	48	63	22	−4.945	5.56E−01	10.104	2.95E−01	−5.159	4.41E−01
Calcium	19	14	7	1.651	5.38E−01	−2.485	5.43E−01	1.323	5.57E−01
Fc epsilon RI	17	22	10	0.817	7.12E−01	1.720	5.37E−01	0.169	5.93E−01
GnRH	26	23	7	2.154	5.71E−01	−0.969	7.28E−01	−1.185	6.59E−01
AGE-RAGE	101	87	32	7.630	7.80E−01	−3.293	9.32E−01	−4.337	7.47E−01
PPAR	8	38	12	−15.834	1.46E−01	14.454	1.29E−01	1.380	7.56E−01
Chemokine	39	35	17	−2.644	4.68E−01	1.418	6.20E−01	1.226	7.60E−01
Adipocytokine	30	26	16	0.737	8.02E−01	−1.598	6.81E−01	0.861	7.70E−01
Prolactin	65	38	20	14.013	2.73E−01	−12.916	3.83E−01	−1.097	8.42E−01
Notch	16	10	9	3.160	3.33E−01	−3.082	3.19E−01	−0.078	8.63E−01
T cell receptor	126	58	32	32.827	9.80E−02	−30.713	1.29E−01	−2.114	8.75E−01
Rap1	26	46	17	−10.120	8.90E−02	9.978	9.20E−02	0.142	8.88E−01
Hippo	35	36	19	−0.312	9.40E−01	0.920	7.81E−01	0.119	8.88E−01
Glucagon	18	37	10	−9.534	1.24E−01	10.062	1.65E−01	−0.528	9.02E−01
Neurotrophin	88	98	45	−5.820	6.09E−01	5.238	6.66E−01	0.582	9.05E−01
mTOR	90	139	64	−25.693	4.57E−01	25.742	4.71E−01	−0.049	9.52E−01

The columns forward/backward and simultaneous represent the number of interaction of this direction; Size effect_f/b/s: average number of interaction that meets the “forward”/“backward”/“simultaneous” condition; pvalue_f/b/s: P-value of “forward”/“backward”/“simultaneous” condition. Significant P-values (*P < 0.05) are presented in bold italics.

We find that 18 of the observed pathways (ovarian steroidogenesis, JAK-STAT, IL-17, toll-like receptor, VEGF, RIG-I-like receptor, cAMP, ErbB, NF-kappa B, NOD-like receptor, Wnt, TNF, MAPK, B-cell receptor, Sphingolipid, C-type lectin receptor, p53, and Hedgehog) display a distribution of appearance pattern that has significantly more or less backward/simultaneous/forward pairs than expected by chance alone based on 1,000 of random permutations. For 29 other pathways, the statistical test suggests no deviation from a random gene appearance pattern (Table 2 and supplementary Data S4, Supplementary Material online). Given the large number of tests, it cannot be ruled out that the significant results are random. This further supports the third scenario (Fig. 2c).

Size of Modules

Fourteen of the observed pathways (adipocytokine, C-type lectin receptor, cAMP, ErbB, Hedgehog, IL-17, JAK-STAT, mTOR, ovarian steroidogenesis, RIG-I-like receptor, sphingolipid, T-cell receptor, toll-like receptor, and VEGF) display modules significantly larger than expected by chance alone (Table 3 and supplementary Data S4, Supplementary Material online). Thus, 33 other pathways do not have any deviation from distributions under randomness. In a similar way to the number of appearance patterns, the size of modules is mainly random in the majority of the pathways.

Table 3

Summary of module size statistical tests

	backward or simult		backward		simultaneous		forward		forward or simult
path	sizeeffct	pvalue	sizeeffct	pvalue	sizeeffct	pvalue	sizeeffct	pvalue	sizeeffct	pvalue
IL-17	−1.203	1.10E−02*	−1.120	9.00E−03*	0.126	6.29E−01	0.421	3.34E−01	−0.006	9.76E−01
Ovarian steroidogenesis	0.721	1.40E−02*	0.993	1.10E−02*	−0.468	1.63E−01	−1.099	5.00E−03*	−1.170	1.00E−03*
mTOR	1.267	3.50E−02*	0.162	5.20E−01	0.368	1.70E−01	−0.626	9.60E−02	−0.728	2.37E−01
Adipocytokine	3.020	4.30E−02*	0.738	1.96E−01	0.082	6.91E−01	−0.434	5.13E−01	−1.368	4.06E−01
cAMP	5.520	4.30E−02*	6.378	3.00E−03*	−1.307	1.16E−01	−1.920	1.65E−01	−3.757	1.51E−01
RIG-I-like receptor	−0.234	7.88E−01	−1.088	3.20E−02*	1.163	3.00E−03*	0.309	5.19E−01	1.170	5.00E−02
Sphingolipid	0.610	9.20E−02	−0.314	2.84E−01	1.135	1.00E−03*	−0.325	2.80E−01	0.192	5.21E−01
JAK-STAT	0.557	1.70E−01	−0.563	2.41E−01	0.758	5.00E−03*	−0.564	1.53E−01	−0.596	2.05E−01
Hedgehog	0.271	2.74E−01	−0.202	4.43E−01	0.602	8.00E−03*	−0.251	3.02E−01	0.182	4.89E−01
C-type lectin receptor	0.297	7.90E−02	−0.027	8.59E−01	0.309	2.40E−02*	−0.096	5.12E−01	0.108	5.10E−01
VEGF	2.838	9.10E−02	−0.899	5.58E−01	1.169	4.40E−02*	−1.343	3.04E−01	0.107	6.66E−01
Toll-like receptor	−1.894	5.70E−02	−1.022	8.90E−02	−0.156	5.73E−01	1.553	1.20E−02*	1.560	6.80E−02
T-cell receptor	−2.459	3.57E−01	−1.727	1.66E−01	−0.091	9.57E−01	2.552	4.60E−02*	2.974	1.07E−01
ErbB	−0.396	7.45E−01	0.702	2.33E−01	−0.376	5.26E−01	−1.516	5.10E−02	−2.444	2.70E−02*
TGF-beta	0.637	1.85E−01	0.497	1.01E−01	−0.413	1.76E−01	−0.110	7.61E−01	−0.918	7.40E−02
Apelin	−1.097	2.40E−01	0.005	8.10E−01	−0.427	3.03E−01	0.842	5.70E−02	1.266	8.20E−02
NF-kappa B	0.397	1.38E−01	0.342	1.41E−01	−0.184	3.63E−01	−0.341	1.69E−01	−0.478	8.90E−02
AMPK	−0.160	9.44E−01	−0.609	2.47E−01	0.548	5.40E−02	−0.031	9.48E−01	0.919	8.90E−02
Oxytocin	−0.016	8.26E−01	−0.159	9.13E−01	−0.376	2.29E−01	−0.396	5.18E−01	−1.631	1.30E−01
Chemokine	0.305	6.20E−01	−0.195	8.61E−01	−0.552	1.97E−01	0.168	6.33E−01	−1.908	1.30E−01
p53	−0.325	9.46E−01	−2.089	3.46E−01	0.937	2.30E−01	1.093	3.33E−01	3.059	1.34E−01
Relaxin	0.028	8.69E−01	0.488	9.30E−02	−0.335	1.87E−01	−0.011	9.81E−01	−0.821	1.47E−01
PPAR	2.934	1.99E−01	3.487	1.24E−01	0.699	4.04E−01	−3.458	1.10E−01	−2.825	1.50E−01
Insulin	−0.758	6.02E−01	0.044	8.11E−01	−0.204	6.42E−01	−0.371	6.18E−01	−1.443	2.69E−01
GnRH	−0.909	3.19E−01	−0.198	8.76E−01	−0.443	2.20E−01	−0.333	7.29E−01	−0.870	3.32E−01
NOD-like receptor	−0.005	9.81E−01	−0.424	1.09E−01	0.255	2.11E−01	−0.192	4.82E−01	0.292	3.54E−01
AGE-RAGE	2.226	1.16E−01	0.260	6.58E−01	−0.104	8.98E−01	−0.734	5.68E−01	−1.725	3.99E−01
MAPK	0.249	7.35E−01	0.566	4.55E−01	0.249	5.25E−01	−1.552	1.67E−01	−2.128	4.16E−01
Rap1	0.981	1.48E−01	0.907	5.80E−02	−0.101	7.04E−01	−0.429	4.08E−01	−0.829	4.17E−01
B cell receptor	0.619	3.85E−01	−0.660	3.17E−01	0.665	6.70E−02	−0.137	8.92E−01	0.520	4.48E−01
Neurotrophin	0.494	6.56E−01	0.533	4.47E−01	0.324	3.83E−01	−1.527	1.20E−01	−1.628	5.26E−01
HIF-1	0.957	4.44E−01	0.042	9.15E−01	0.481	3.02E−01	−1.473	2.25E−01	−1.027	5.49E−01
Wnt	4.331	7.60E−02	1.853	8.10E−02	0.814	1.25E−01	−0.817	6.28E−01	0.657	5.52E−01
Fc epsilon RI	−0.296	6.70E−01	−0.126	9.05E−01	0.332	2.19E−01	−0.372	5.52E−01	−0.449	5.55E−01
Prolactin	0.096	8.39E−01	−0.670	2.57E−01	0.304	4.08E−01	0.104	7.12E−01	0.415	5.67E−01
Hippo	0.105	6.66E−01	−0.057	8.05E−01	0.066	7.13E−01	−0.082	6.81E−01	0.083	7.13E−01
PI3K–Akt	−0.476	5.37E−01	−0.513	2.68E−01	0.205	3.96E−01	0.340	3.00E−01	0.086	7.75E−01
FoxO	−1.066	5.56E−01	−0.605	6.66E−01	−0.593	3.39E−01	−0.175	9.82E−01	−0.604	7.93E−01
Phospholipase D	−1.227	2.58E−01	−0.532	3.25E−01	−0.713	1.82E−01	0.431	2.14E−01	0.003	8.26E−01
Thyroid hormone	−0.242	4.84E−01	−0.160	6.13E−01	−0.074	7.76E−01	0.176	5.64E−01	0.060	8.44E−01
Estrogen	0.298	4.10E−01	0.172	5.37E−01	0.116	7.80E−01	−0.002	9.65E−01	−0.004	9.06E−01
Notch	−0.339	6.71E−01	−0.349	6.80E−01	−0.191	6.90E−01	0.493	3.43E−01	−0.057	9.39E−01
Ras	−0.266	9.62E−01	−0.614	5.44E−01	−0.007	8.69E−01	0.066	7.96E−01	−0.345	9.42E−01
Glucagon	0.784	1.07E−01	0.585	1.46E−01	−0.101	7.31E−01	−0.553	2.15E−01	−0.168	9.44E−01
TNF	0.461	2.30E−01	−0.093	7.60E−01	0.287	3.47E−01	−0.164	6.12E−01	−0.078	9.59E−01
Calcium	−0.553	1.71E−01	−0.131	7.01E−01	−0.201	4.36E−01	0.125	6.84E−01	−0.020	9.62E−01
cGMP-PKG	−0.052	9.86E−01	−0.036	9.32E−01	0.158	6.05E−01	−0.181	9.01E−01	−0.105	9.98E−01

Size effect: average size of the module that meets each condition, calculated by averaging the size of the module across 1,000 simulated pathways; pvalue: P-value of each condition. Significant P-values (*P < 0.05) are presented in bold italics.

Summary of Hypotheses Testing

Among the 47 pathways, 30 have no correlation between the age rank differences (delta) with respect to the position in the pathway (Table 1), 29 do not have more backward/simultaneous/forward pairs than randomness (Table 2), and 33 do not have modules larger than obtained by chance alone (Table 3). Interestingly, 20 pathways constitute the intersection of these tests: they can be interpreted as randomly organized for all the variables considered in our study. These 20 pathways are: apelin, B-cell receptor, calcium, cGMP-PKG, chemokine, estrogen, Fc epsilon RI, FoxO, glucagon, GnRH, Hippo, NF-kappa B, Notch, p53, phospholipase D, PI3K–Akt, PPAR, prolactin, Rap1, Ras, relaxin, TGF-beta, thyroid hormone, and TNF. A total of 581 genes are involved only in the pathways for which a correlation and/or permutation statistical test yielded significant results, and 334 genes are exclusively associated with the pathways that provide support for one of these architectures. Finally, 14.9% of the tests are significant and their effects are different among pathways. This low ratio can rule out that scenarios A and B but is compatible with scenario C.

Discussion

Ages of Appearance and Their Patterns in Evolution

For most of the pathways with significant correlations between position of proteins in the pathway and the delta of appearance of their corresponding gene, the correlation is surprisingly positive: the more upstream a protein is involved in the pathway, the older is the upstream gene (Fig. 2b). The sole exception is the insulin pathway that has a negative correlation (Fig. 2a).

Finally, for 24 of the 47 pathways studied here, there is no significant correlation between the position of the proteins in the pathway and the delta of appearance of the corresponding gene in the tree of life, which reflects evolutionary tinkering (Fig. 2c); furthermore, the architecture of proteins and their interactions is not significantly different from 1,000 simulation taken at random (a simulation is a possible configuration randomly generated from a signaling pathway, where the interactions and their directions are maintained, but the ages of the involved proteins are randomly swapped). This result suggests that for these 24 pathways, the organization of proteins and interactions during evolution was not guided by a molecular determinism. One may formulate two hypotheses: that the first interactions may have arisen under determinism/pressure (e.g. 100% backward; Fig. 2a) and may subsequently have been modified by new interactions during evolution (a form of evolutionary tinkering that may improve fitness) that erased the signal in the initial pathway (Fig. 2c) and that several modules arose partly by chance and then further new interactions made the entire pathway necessary. This “chance and necessity” hypothesis was proposed by Monod (1970) in his famous book “Le hasard et la nécessité”. This evolutionary concept has been widely discussed and remains accepted by many (Merlin 2015; Davankov 2021; Xie et al. 2021; Alexander 2022).

Moreover, even if the correlation between the position of proteins in the pathway and the delta of appearance of their corresponding genes is not significant, we observe that almost all signaling pathways have target factors (e.g. factors entering the nucleus) that appeared very early during evolution (dark blue factors to the right in the colored KEGG diagrams; supplementary Data S3, Supplementary Material online). For various pathways, some factors are close to the targets and appeared very early (Opisthokonta or Metazoa), whereas those further upstream in the pathway appeared later (PPAR and cAMP). For other pathways (IL-17 and T-cell receptor), some factors close to the targets appeared later than the factors acting upstream of the pathway. This also supports that pathways or parts thereof were not formed/refined according to a unique deterministic law like scenarios A or B but rather conforms to scenario C (Fig. 2c) and suggests that some evolutionary tinkering took place.

In a previous work (Grandchamp and Monget 2018), we have shown that the genes encoding the pairs of ligand/membrane receptors mainly appeared at the root of Metazoa, Vertebrata, and Osteichthyes. In the present work, we show that the genes encoding proteins of signaling pathways mainly appeared between the last common ancestor of Opisthokonta and Metazoa and between Deuterostomia and Chordata. Moreover, as expected, among all the pathways, at least 119 proteins (from 1 to 23 by pathway) appeared in Opisthokonta, showing, as expected, that numerous genes encoding proteins of signaling pathways appeared long before the membrane receptors and their ligands. Our previous study on membrane receptors and their ligands was done on a tree with only 10 nodes, but we refine it by using 15 nodes in the present paper.

Some signaling pathways are relatively old (i.e. contain more members that appeared in Opisthokonta or Metazoa than in more recent clades), such as the MAPK pathway, and others are more recent, such as the pathways involved in immunity (Cooper and Alder 2006).

Some Interesting Cases

In the signaling pathways studied, numerous proteins interact with their partner(s) in humans, but their appearance was desynchronized in the tree of life (Fig. 2c). This raises the question of the functioning of these pathways without the full set of their members. For example, p53, which appeared in Metazoa, is inhibited by several factors whose genes appeared later, such as MDM2 and MDM4 (Olfactora). Without these inhibitors, mammalian embryos would die in utero, as shown in the mouse (Jones et al. 1995; Montes de Oca Luna et al. 1995). This suggests that another p53 inhibitor existed before Olfactora appeared. In addition, human p53 targets also emerged later in evolution, such as IGFBP3 (Vertebrata) and pro-apoptotic factor BAX (Olfactora). It is these and other factors that explain the tumor suppressor role of p53 in humans (Lehmann-Che et al. 2007) and in Drosophila (Zhou 2019). We can therefore hypothesize a progressive refinement of the mechanisms of action and inhibition of p53 during evolution; for example, a possible first inhibitor of p53 in an ancestor may have been replaced by another during evolution.

Another example is the interaction between PIN1 and IRF3 in the RIG-I-like receptor pathway, characterized by a delta of 8, the PIN1 gene having appeared at the root of Opisthokonta and IRF3 in stem-gnathostomes. PIN1 also activates p53, which is an apomorphy of Metazoa (Berger et al. 2005). Another example of desynchronization concerns the sphingolipid pathway, in which CTSD interacts only with BID, whereas BID interacts with CTSD and BAX in humans. The CTSD gene appeared in Metazoa, BAX in Olfactora, and BID in Mammalia. Moreover, although the yeast genome does not contain genes encoding BCL2 proteins, the heterologous expression of mammalian BAX in yeast induces a suppressible lethal phenotype that is associated with characteristics of metazoan apoptosis, strongly suggesting that its targets are present in yeast and thus hark back to Opisthokonta (Zha et al. 1996; Khoury and Greenwood 2008).

Limits and Future Improvements

This study could have been completed by an analysis based on known pathways in yeast to verify certain interactions. However, only one pathway from our list, MAPK, is described in KEGG for Saccharomyces cerevisiae. Unsurprisingly, this pathway is widely documented (Zou et al. 2008; Saito 2010; Chavel et al. 2014), including in plants (Meng and Zhang 2013). In the literature, several other pathways (not studied here) are documented in yeast, such as cAMP (Tamaki 2007; Portela and Rossi 2020), mTOR (Powers et al. 2004), Ras (Tisi et al. 2014), or sphingolipid (Montefusco et al. 2014). Thus, using the KEGG database, it is not possible to establish an exhaustive evolutionary bridge of signaling pathways between yeast and humans. Moreover, certain taxa lack specific pathways; for example, the yeast S. cerevisiae does not have the NF-kappa B pathway (Ho et al. 2017; Saleski et al. 2017), and yet, some elements of the pathway are present in S. cerevisiae, such as casein kinase 2 (CSNK2A1/2/3 and CSNK2B). In yeast, these caseins are known to be essential for mitophagy (Kanki et al. 2013).

Another problem concerns the limitations of our method of dating the appearance of a gene. In this, we are limited by the various versions of the Genomicus and Genomicus Metazoa databases. In particular, the species common to these two databases are Drosophila melanogaster and Caenorhabditis elegans. Indeed, if for a given gene the oldest ortholog found is not from these two species, or if this gene is lost in both species but present in taxa more distantly related to humans, our method does not allow us to adequately use the Genomicus Metazoa trees. Genomicus trees are modified trees of the Ensembl database, and therefore, the limit comes from its origin, Ensembl. Furthermore, depending on the versions of Ensembl, there may be mega-trees with all the paralogs of a gene in the same tree or in sub-trees with one paralog per sub-tree. This change from mega-trees into several sub-trees in Ensembl V94 and further (Emily 2018) makes the attribution of appearance times more complex, to an extent that depends on the size of the paralog family.

Not all of the 15 clades that we have sampled have been studied to the same extent. As shown in supplementary Data S5, Supplementary Material online, we can see that yeast are represented by a single species (S. cerevisiae) out of about 64,000 known species, whereas for Aves, we were able to collect data for 13 species, out of a little more than 11,000 currently recognized. However, our data show that appearance of the genes is concentrated on two nodes (Opisthokonta and Vertebrata), and the genes that appeared at the base of these taxa represent 75.7% of all the genes involved in the KEGG pathways.

It must also be considered that the signaling pathways noted on KEGG (or databases more broadly) are human representations and that simplifications have been made to facilitate reading and understanding. However, as given in supplementary Data S3, Supplementary Material online, some proteins, such as Shc, Grb2, SOS, Ras, Raf1, MEK, and ERK, are involved in multiple pathways. In fact, interactions within these proteins are involved (partially or entirely) in 22 pathways, including ErbB, estrogen, GnRH, insulin, prolactin, relaxin, B-cell receptor, chemokine, F epsilon RI, T-cell receptor, neurotrophin, cAMP, FoxO, JAK-STAT, MAPK, mTOR, phospholipase D, PI3K–Akt, Rap1, Ras, sphingolipid, and VEGF. The RTK/RAS/ERK component is known to be common to Drosophila, Nematoda, and humans (Ashton-Beaucage and Therrien 2010).

In a previous study, we described cases where the genes coding for membrane receptors appeared before their ligand (Grandchamp and Monget 2020). We had studied more precisely the 30 cases for which 3D structures were known in the Protein Data Bank, to formulate hypotheses on the plausible scenarios of evolution of the amino acids involved in the binding pocket of the receptor until the ligand appeared. In future, a similar study of the partners that appeared before the proteins with which they interact would be worthwhile. Such a study could be performed for Ras/SOS complexes, for example, on Grb2/SOS, Akt/Gsk3, or Akt/mTor, for which the 3D structure is available.

Materials and Methods

Implementation of the Database

Signaling Pathways by KEGG

We follow the method described in our previous paper (Picolo et al. 2023). We use the KEGG pathways because they are annotated in humans and therefore its genes are also annotated in humans, and we use parsimony and the trees to locate the branch on which each gene originated. We retrieved a list of 2,298 genes encoding proteins involved in the 47 human intracellular signaling pathways available on the KEGG V104.0 database (https://www.genome.jp/kegg) (Bader et al. 2006) using the keywords “signaling pathway” and “human” (columns C and D in supplementary Data S4, Supplementary Material online). KEGG is one of the most referenced and used databases listing signaling pathways.

Each signaling pathway was retrieved from the KEGG pathway analysis tool, and the information was compiled into an Extensible Markup Language (XML) file. These files were read using the R XML (Lang and Kalibera 2023) library, and pathways were represented graphically using the R library igraph (Csárdi and Nepusz 2006; Csárdi et al. 2023).

The gene products of KEGG pathways are inscribed in rectangular blocks that we call “labels”; for the XML file, this is the label “name,” and to keep it simple for our readers, these labels can cover several paralogs (supplementary Data S3, Supplementary Material online). Some paralogous genes appear under different labels, as in the PPAR pathway for which the proteins PPARα, PPARγ, and PPARβ/δ have similar interactions (with FABP3, thiolase B, AP2, and UCP-1) and more specific interactions (e.g. PPARα → HMGCS2 and PPARγ → GyK). Each pathway is constructed in such a way that there can be one or more inputs (e.g. a ligand) and one or more outputs (e.g. transcription factor). Each subpathway is analyzed, which means that the same gene can be involved in several subpathways (e.g. A → B → C and A → D → E).

Tree of Life From Genomicus

To determine the time of origin of genes, two phylogenetic trees from Genomicus (Nguyen et al. 2018) were used: the vertebrate tree and the metazoan tree. The vertebrate tree is the V109 tree (https://www.genomicus.bio.ens.psl.eu/genomicus-109.01), comprising 199 animal species and covering a range of vertebrate species, as well as D. melanogaster and C. elegans, and to this is added the yeast species S. cerevisiae. The metazoan tree is that of V51 (https://www.genomicus.bio.ens.psl.eu/genomicus-metazoa-51.01), comprising 116 animal species and covering all main metazoan clades except for Vertebrata. Drosophila melanogaster and C. elegans are the two species present in both of these databases (supplementary Data S5, Supplementary Material online).

From these Genomicus trees, we determined 15 clades on a similar and “enriched” model of a previous work (Grandchamp and Monget 2018): (i) Opisthokonta (∼1,300 my; taxon age here and below refers crown groups and to the divergence between the various taxa (S. cerevisiae, in this case) and Homo sapiens, unless stated otherwise); (ii) Metazoa (∼765 my); (iii) Eumetazoa (∼743 my); (iv) Bilateria (∼708 my); (v) Deuterostomia (∼635 my); (vi) Chordata (∼588 my); (vii) Olfactora (∼570 my); (viii) Vertebrata (∼563 my); (ix) Gnathostomata (∼462 my); (x) Osteichthyes (462 my); (xi) Sarcopterygii (∼415 my); (xii) Tetrapoda (∼350 my); (xiii) Amniota (∼330 my); (xiv) Mammalia (∼166 my); and (xv) Theria (∼98 my). To provide consistent ages and given that several of the oldest nodes are poorly constrained by the fossil record, we have used molecular time estimates throughout (de Vienne 2016; Kumar et al. 2022; Nguyen et al. 2022). These clades are our references for the present study. The gene is considered to appear on the branch that leads to the clade for which parsimony unambiguously indicates its presence ancestrally.

Dating the Appearance of Genes

In order to identify the origin of a gene, we trace its orthologs along the evolutionary tree, starting with humans (our reference species because our knowledge of these pathways is human-centered). We identify the deepest evolutionary node that includes H. sapiens where orthologs are already present as the age of appearance. For example, if we identify a teleost orthologue of a human gene, but we do not find one for more distantly related taxon, we infer that the node of emergence of this gene is probably Osteichthyes. (https://github.com/florianepicolo/origin-gene).

Some cases require a more complex analysis. If, under a label in the KEGG pathways, several paralogs are listed, then the latest possible date of origin of this gene is the oldest node at which paralogs occur. If several paralogs are present in a pathway but under different labels, we consider them independently from other paralogs. In the case of a protein complex, we consider the complex to be a group, so we consider the node of origin of the complex to be the deepest node for which all the genes of the complex are present, because the complex cannot be functional unless all the proteins are present.

All elements in a pathway are assigned a position in the pathway, and multiple positions are possible for the elements that are present in several pathways. We have not considered the “components” of the pathways present on KEGG in the allocation of the positions if they are not proteins and therefore do not involve genes, such as Ca²⁺ or lactate (e.g. A₁ → B₂ → Ca²⁺ → C₃, with pathway position in index).

Furthermore, in an interaction A → B (A is closer to the receptor than B), the interaction is said to be “backward” if gene A was born after B, “forward” if gene A was born before B, and “simultaneous” if A and B were born on the same branch (this condition may reflect extinction of the relevant taxa or insufficient taxonomic resolution of our study, which prevents us from untangling the pattern of appearance).

Analyses and Statistics

To explore the relative age of genes in the signaling pathways, we used several variables with both descriptive and hypothesis testing approaches (supplementary Data S6, Supplementary Material online). Indeed, for each pathway, we considered the pathway position of genes (supplementary Data S6-b, Supplementary Material online), rank age of appearance (supplementary Data S6-c, Supplementary Material online), appearance pattern in an interaction (supplementary Data S6-d, Supplementary Material online), rank age difference between genes in an interaction (supplementary Data S6-e, Supplementary Material online), and size of module sharing the same appearance pattern. A module is a subset of a pathway that achieves a specific condition, in our case a pattern of appearance for a given interaction (e.g. A -> B): forward (A before B), backward (B before A), or simultaneous (A at the same time as B) (supplementary Data S6-f, Supplementary Material online). As the distribution of variable values in networks is usually not normal and implies complex structures of dependences, we used a permutation algorithm to test the hypotheses (supplementary Data S6-g, Supplementary Material online). The algorithm works as follows for each pathway: (i) compute the statistics of interest; (ii) randomly permute the rank ages of origin of all genes from the pathway; (iii) compute again the statistics on the pathway with random rank ages; (iv) redo steps 2 and 3 a large number of times (here, 1,000 times); (v) if the statistics have more than one value per pathway, compute the medians; and (vi) use the values of the statistics (or the medians) as the distribution, under the null hypothesis that there is no relationship between the rank ages of origin and the statistics under scrutiny. Comparing the values from the raw data (step 1) to the distribution of values under the null hypothesis (step 6) allows P-value (here, two-sided) to be computed. In other words, we evaluate how probable the values of the statistics on the raw data are if we accept the null hypothesis. For clarity, the nomenclature of the two types of pathway subsets (subpaths and modules) is distinguished in supplementary Data S6-f and S6-i, Supplementary Material online.

Finally, we tested the correlation between the differences of rank ages (delta age; supplementary Data S6-e, Supplementary Material online) and the position of the upstream gene in the pathway (supplementary Data S6-b, Supplementary Material online). Because pathways are networks, cycles imply several pathway position attributions for some genes (supplementary Data S6-h, Supplementary Material online). We were not able to choose one pathway position hypothesis over another in networks with cycles. Therefore, all hypotheses were included, and each was down-weighted proportionally to the number of alternative hypotheses.

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Acknowledgments

The authors thank Éric Reiter and Pascale Crépieux for helpful and encouraging discussion.

Data Availability

All data are available free of charge online or on request.

Author notes

Floriane Picolo and Jérémie Bardin equally contributed to the present work.