Phylostratigraphic analysis of tumor and developmental transcriptomes reveals relationship between oncogenesis, phylogenesis and ontogenesis

Multicellular stability is the outcome of two antagonistic processes: one promotes individual cell fitness and the other promotes tissue (collective) fitness. The activities of oncogenes and tumor suppressor genes are central to the equilibrium of these two processes. Overt proliferation and cell-autonomy can result from over-activity of oncogenes that stimulate the cell division cycle, or from the reduced activity of the tumor suppressor genes. Both scenarios can be interpreted as biological functions that support single cellular fitness and can result from activating mutations in the former or inactivating mutations in the latter.

Accordingly, what we now call tumor suppressor genes have been proposed to have played a critical role in the evolution of multicellularity [50–54], as part of the complex machinery of top-down controls necessary for establishing organismal coherence and homeostatic stability. On the other hand, retaining oncogenes offers a less apparent general advantage to multicellular organisms. Given that cell proliferation could be a default state in unicellular organism [55], as long as there is sufficient energy supply and space, one possible reason to retain oncogenes is that extant oncogenes evolved from genes involved in the machinery that implements this default state, or are derived from ancestral 'contingency genes' of unicellular organism that originally evolved to enable populations of cells to respond to stressful environments by promoting proliferation. In multicellular organism then such genes would have evolved to promote tissue expansion during development (providing the somatic cells) and regeneration (repairing the somatic cells)—both would serve the propagation of the germline. In any case, one can expect that evolution of a fraction of oncogene precursor genes predate the appearance of multicellularity or at least, that of tumor suppressor genes.

We first used existing knowledge bases to define three set of genes: all human genes, genes that have been labeled 'tumor suppressor genes' and genes that have been labeled 'oncogenes' [56, 57]. In parallel, we categorized all genes in five evolutionary age groups based on the eggNOG database [58], as shown in figure 1(a), from old to young: LUCA (last unicellular common ancestor), Eukaryota, Metazoa, Vertebrata and Primata (for details about the determination of gene ages see Material and Method). When multicellular life appeared on earth, the cellular traits that are heir to unicellular protozoa were suppressed or tightly controlled by regulatory mechanisms which evolved to promote collective over individual fitness [59]. Since evolution tends to build layers of phenotype on top of existing ones by amending developmental processes and suppressing primitive elements [47, 60], a concept that explains the existence of atavism (such as tails in human) [61], cancer could be depicted as a form of cellular atavism: a reversion to single cell behavior due to reactivation of hidden (normally suppressed) ancient 'genetic programs' (figure 1(b)). This reactivation could be triggered by mutations (or other disturbances) that disrupt the functional mechanisms that protect multicellularity [41, 43, 45]. Recent work [43] showed that ancient mutational responses to stress, co-opted in multicellular life to maintain diversity in the germline and immune system, can cause cells to revert to a primordial form that displays the cancer phenotype, including elevated genetic instability and accelerated somatic adaptation.

Zoom In Zoom Out Reset image size

In a first analysis, we counted the number of genes with the labels 'tumor suppressor' or 'oncogene' and computed their relative frequency among the gene classes of the different evolutionary age and compared it with the frequency of all genes as the background, i.e. we determined the relative enrichment and depletion of suppressor and oncogenes in each of the evolutionary age groups. Figure 2(a) presents a simple comparison of the frequencies of occurrence in each evolutionary age group for the three sets of genes, all human genes, tumor suppressor genes and oncogenes. In figure 2(b) we use an enrichment score based on an adjusted residual analysis (see Material and Methods), comparing suppressor genes and oncogenes against all human genes as the background. Both plots show that the frequency of the tumor suppressor genes are depleted for genes at the Vertebrata level while they are enriched in the age groups of the Eukaryota and Metazoa level. Oncogenes, however, also display a stronger enrichment for genes of the Eukaryota and Metazoa levels. By contrast, the frequency of oncogenes are specifically depleted for genes in Vertebrata and Primata levels. Thus in general both gene sets follow a very similar pattern: both groups of genes mutated in cancer typically evolved before or during the evolution of multicellularity, but oncogenes seem to have a stronger pre-metazoan component earlier, appearing during the evolution of eukaryotes, as suggested by the fact that their frequency is actually depleted in LUCA.

Zoom In Zoom Out Reset image size

Since genes function as interacting gene networks, we also examined the gene sets involved in cancer based on gene ontology (GO) annotations which places genes with the same function into groups based on shared GO labels. In figure 2(c) we separately curated all GO terms and cancer-related GO terms and then calculated the average age of each one using our gene age database. The cancer-related GO terms display a strong enrichment in the age groups between Eukaryota and Metazoa, again suggesting that cancer-related genes evolved before or during the evolution of multicellularity.

Since there is no significant difference in the evolutionary age enrichment between oncogenes and tumor suppressors found in tumor genomes, it has been suggested that considering pathways as the 'units' targeted by mutations may be more appropriate [62]. Therefore, we next curated sets of genes involved in cancer-related pathways from the KEGG database [63] and analyzed the gene age enrichments by comparing the age distribution of such gene sets with the background distribution of non-cancer genes. We analyzed 25 pathways associated with cancer (see table 1 (stacks.iop.org/CSPO/4/025002/mmedia)) as well as 14 housekeeping pathways (see table 2). The latter set can be used to establish a benchmark for the comparison of the patterns observed in cancer.

In each case, cancer pathways were grouped according to their KEGG-BRITE orthologous groups, which spans categories ranging from those responsible for essential cell functions, such as metabolism, genetic information processing and cellular processes, to those of organismal-level processes, such as environmental information processing, organismal systems and a few categories involved in human diseases.

Among the cancer pathways (figure 3(a)), our analysis shows that pathways associated with cellular growth and genetic information processing, which could be considered cell-level processes, are typically enriched in the ancient LUCA and Eukaryota gene groups, while most pathways controlling the organismal-level processes are enriched in the group of genes that appeared with Metazoa and Vertebrata. House-keeping pathways (figure 3(b)) related to metabolism are enriched in the group of LUCA genes while the pathways involved in signal transduction again, are enriched in the groups of Metazoa and Vertebrata genes. Thus, the functional distinction between cell-level and organismal-level functions is consistent with the gene age profile.

Zoom In Zoom Out Reset image size

It is important to mention that enrichment in this analysis only pertains to the number of genes (see details in table 1) in a pathway that falls in each age category, and thus this is a crude analysis: pathways constitute of genes connected by a regulatory relationship and as such, different subsets of genes do not necessarily have the same relevance in the dynamics and functionality of each pathway. Therefore, the enrichments that we observe are just a testament of the overall trend in the fraction of genes that evolved before or after multicellularity. The changes in the activities and topologies of these pathways could of course be equally important to a full understanding of the evolution of their biological functions, but we still obtain fairly consistent results from the gene membership annotation alone. Another important observation is that cancer-related and housekeeping pathways are not mutually exclusive, as they share six common pathways (VEGF, mTOR, PI3K-Akt, WNT, Hedgehog and Adherens). Firstly, we must emphasize that there is no aggregation of statistics across pathways, and thus results for each pathway are independent of pathway family selection. Second, these six overlapping pathways regard essential biological functions accessed or modified by mechanisms relevant to cancer progression. But whether we contrast the features displayed by cancer pathways with each or against housekeeping pathways in general, our results indicate a consistent trend that depends more strongly on the evolutionary makeup of the pathway than on the cancer or normal nature of it. This observation is manifested by the contrast of positive and negative enrichment scores between cellular and environmental information processes.

Beyond looking at the evolutionary age profile of genes and pathways from curated knowledge bases that predetermines relationship to cancer, we next analyzed the characteristic gene expression profiles directly from cancer tissue samples as has recently been done by Trigos et al [45]. We used the currently largest cancer omics database, TCGA, which contains thousands of primary tumor RNA-Seq samples and their corresponding normal tissue samples for various cancer types.

To define a gene as differentially expressed between cancer and normal tissues in view of the vast heterogeneity of the transcriptomes even among the samples within the same tumor type ('inter-patient heterogeneity') we first performed a standard t-test on each cancer sample against the population of corresponding normal tissues. If a gene expression value in a cancer sample is three standard deviations ${{\sigma}^{n}}$ above the population mean $\overline{{{x}^{n}}}$ for the corresponding normal tissue, this gene is deemed 'overexpressed', i.e. if ${{x}^{c}}>\overline{{{x}^{n}}}+3{{\sigma}^{n}}$, then gene ${{x}^{c}}$ is 'overexpressed' in this cancer tissue sample; conversely, ${{x}^{c}}<\overline{{{x}^{n}}}-3{{\sigma}^{n}}$, gene ${{x}^{c}}$ is considered 'suppressed' in cancer tissue sample. Second, the overlap frequency $F$, as the percentage of samples that share the commonly suppressed or overexpressed gene, was used as another criterion for categorization. For instance the set of up-regulated genes with $F>10\%$ consists of all the genes that are overexpressed in more than 10 percent of the studied samples and so forth, with higher % values indicating more robust alteration of expression in cancer. We then proceeded to study the evolution age distribution enrichment of overexpressed and suppressed genes with different overlap frequencies $F$.

We first performed this analysis for breast cancer samples and the corresponding normal tissue samples. As shown in figure 4(a), with regard to the evolutionary age group, the frequency of the genes that exhibited most pronounced differential suppression in breast cancer were significantly enriched for Metazoa and Vertebrata evolutionary ages while slightly depleted in the Eukaryota age group. This result suggests that malignancies tend to suppress the expression of post-metazoan genes, those likely involved in multicellular functions such as intercellular signaling. At the same time, the more primordial genes, likely involved in essential cell functions, are typically not down-regulated in tumors. Similarly in figure 4(b), the frequencies of genes overexpressed in cancer were significantly enriched with genes of the LUCA level age group and depleted of genes of the Vertebrata age group. Thus, the overall trend observed was that genes up-regulated in breast cancer are typically ancient, although not consistently up-regulated across all samples.

Zoom In Zoom Out Reset image size

A similar analysis was performed for nine other cancer types (see details about cancer types in Material and Method Section) for genes with a sample overlap frequency taken at F  >  30% in all cases in order to focus on genes with consistent overexpression or suppression patterns across the cancer data. The frequencies of genes of a given evolutionary age that were suppressed in cancer were typically higher for the post-metazoan genes and reduced for pre-metazoan genes in several of the cancer types studied (see figure 5(a)). Notable exceptions were Glioblastoma (GBM) and Kidney Renal Clear Cell Carcinoma (KIRC), which were depleted of the Vertebrata age group. Colon Adenocarcinoma (COAD), Head-Neck Squamous Cell Carcinoma (HNSC) and Uterine Corpus Endometrial Carcinoma (UCEC) displayed no significant evolutionary age enrichment. Interestingly, genes suppressed in Glioblastoma (GBM) were enriched in the Eukaryota age group. Analysis of genes in underrepresented age groups using the DAVID Gene Functional Annotation Clustering Tool [64] suggests that the inversion in the trend of GBM and KIRC compared to other cancer types is driven itself by an enrichment of genes functionally involved with production of immunoglobulins and membrane functions, which are typically system-level functions that could be associated with tissue organization and interactions with the environment. This result suggests that genes involved in these functions are better regulated for these particular tissues. Conversely in figure 5(b) shows that overexpressed genes typically present a more regular behavior: enrichment in the LUCA age group and depletion in the post-metazoan age groups for most cancer types. Kidney Renal Clear Cell Carcinoma (KIRC) was again an outlier displaying enrichment in the Vertebrata age group. Functional analysis using DAVID shows that this set of genes was in turn enriched with genes involved in cell adhesion, calcium ion binding and cadherin, which are also tissue-level functions suggesting perhaps that the organ function enforces particular susceptibility to the overexpression of the associated genes.

Zoom In Zoom Out Reset image size

In conclusion, the enrichment of cancer-associated genes with respect to their evolutionary age in ten cancer types from the TCGA RNA-Seq transcriptome database suggests that cancer cells tend to suppress the expression of post-metazoan genes and overexpress pre-metazoan genes. These results are in line with the analysis of oncogenes versus suppressor genes and suggest that cancer cells acquire a state that is close to a primordial state defined by expression of genes that existed in pre-metazoan stage of evolution, as predicted by the cancer atavism hypothesis.

Besides knowing the most frequent evolutionary age of genes that are overexpressed and suppressed in cancer tissues, it is also desirable to identify their biological functions of cancer-associated genes for which we had analyzed the evolutionary age. We computed the enrichment of GO terms for the differentially expressed genes with overlap frequency $F>30\%$ for each cancer type. From figure 6(a) it is evident that overexpressed genes in cancer are involved with basic cell processes such as cell division, cell cycle (cell cycle phase, mitotic cell cycle etc) and organelle fission (including spindle, cytoskeleton and nuclear division). On the other hand, the suppressed genes in cancer are related to multicellular development, embryonic organ development, morphogenesis, cell differentiation and inter-cellular signal transmission, as shown in figure 6(b).

Zoom In Zoom Out Reset image size

This functional analysis, which agrees with a vast body of work on GO enrichment among cancer genes in the past [65–71], is in line with our results of the evolution age enrichment analysis: cancer cells overexpress old genes mainly involved with essential cell functions, such as cell cycle and cell division, and suppresses young genes typically associated with multicellular functions, such as cell differentiation and inter-cellular signaling. This agrees with the observation that an early characteristic of cancer cells the loss the differentiated cell identify and reversion to a stem-cell-like state [72]. Specific mutational patterns of oncogenes and tumor suppressors are not necessary conditions for cells to achieve a cancer cell state, explaining why patients with the same disease can have totally different spectra of gene mutations and overlap minimally with each other. Since the gene regulatory network is a complex system, there are multiple ways to reach the cancer cell state from many distinct combinations of gene mutations or disrupted regulation controls. Therefore, genes that are mutated in cancer display a large diversity while the implementations of the cancerous cell states share a typical core functionality, namely, suppressing the specific functions of cell differentiation and inter-cellular signaling while enhancing the functions of cell division cycle.

Genes differentially expressed between embryo stem cells and the differentiated breast cells show the divergent enrichment of gene evolution age

If tumorigenesis is a process in which cells revert to a well-programmed early-evolved stem-cell-like state, it is very tempting to look at the normal cell differentiation process and align the oncogenesis-phylogenesis axis of analysis performed above with that of the long recognized oncogenesis-ontogenesis axis. If a similar mechanism is at work, we should see that the differentiation process runs opposite to that of tumorigenesis: differentiated cells will suppress the evolutionarily 'old' genes involved with essential cell functions such as cell cycle and cell division but overexpress 'young' genes which are associated to the multicellular processes and inter-cellular signaling.

We curated microarray gene expression profiles from 123 samples of embryonic stem cells (ESC) and from 159 breast epithelial cell lines of primary cell cultures from GEO (see details in Material and Methods). After applying quantile normalization to all microarray data, we used the same methodology to select the set of genes exhibiting differential expression between breast epithelium and ESC: if a gene ${{x}_{{\rm B}}}$ in the breast cell transcriptome is expressed at three standard deviations ${{\sigma}_{{\rm ESC}}}$ above the population mean $\overline{{{x}_{{\rm ESC}}}}$ of the ESC, this gene is considered as overexpressed. Thus, if ${{x}_{{\rm B}}}>\overline{{{x}_{{\rm ESC}}}}+3{{\sigma}_{{\rm ESC}}}$, gene ${{x}_{{\rm B}}}$ is overexpressed. Conversely if ${{x}_{{\rm B}}}>\overline{{{x}_{{\rm ESC}}}}+3{{\sigma}_{{\rm ESC}}}$, gene ${{x}_{{\rm B}}}$ is suppressed. The overlap frequency of samples $F$ is the percentage of samples that share the commonly suppressed/overexpressed gene.

As shown in figure 7(a), the frequency of the genes suppressed in breast cells were significantly enriched in the LUCA and Eukaryota age groups and depleted in the Metazoa, Vertebrata and Primata age groups. By contrast, genes overexpressed in these epithelial cells were enriched in the Metazoa and Vertebrata age groups while depleted in pre-Metazoan age groups (figure 7(b)). Thus, comparing epithelial cells (mature) to embryonic cells (immature) as opposed to comparing cancer cells (presumed immature) to their normal counterpart (mature) yielded opposite results for the evolutionary age of the genes suppressed/overexpressed in these two comparison pairs. Since here the transcriptomes were mostly measured from the same or similar cells, the diversity of cells does not play as much a role as was the case for the comparison between cancer and normal tissues. This may explain why the enrichment scores vary little with respect to different overlap frequencies $F$. In conclusion, human breast cell differentiation seems to have the opposite effect on evolutionary age of the genes expressed compared to from tumorigenesis. To reach the fully differentiated state, cells tend to suppress the pre-Metazoan genes and overexpress post-Metazoan genes, underscoring the known alignment of ontogenesis with phylogenesis.

Zoom In Zoom Out Reset image size

The evolution ages of human genes were extracted from the eggNOG 4.0 database [58]. This database has over 7.7 million proteins categorized into 1.7 million orthologous gene families derived by functional annotations and genetic sequence overlaps, comprising more than 3,600 species. Initially, we identified the orthologous groups that contain human genes and downloaded the pre-computed phylogenetic trees. If the orthologous group included a single human gene, then the phylostratigraphic age was the most ancient phyla of species represented. If there are multiple human species, the phylogenetic trees where parsed at the last common ancestor node for each individual human genes. Each sub-tree is then used to assign the phylostratigraphic age.

The gene families provided by eggNOG are covered in eleven main taxonomic levels for the human lineage: Last Universal Common Ancestor (LUCA), Eukaryota, Opisthokonta, Metazoa, Bilateria, Chordata, Vertebrata, Mammalia, Euarchontoglires, Primata, and Hominidae. In order to have a good number of species per level (and have better statistical significance) we reduce these 11 levels into five basic levels (see the mapping table in table 4):

• 1.  
  Last Universal Common Ancestor (LUCA)
• 2.  
  Eukaryota
• 3.  
  Metazoa
• 4.  
  Vertebrata
• 5.  
  Primata

Genes from intermediate levels were included at a correspondingly lower accepted level. As a result of this methodology, we obtain the taxonomic level for each human gene in the database (19,177 genes). We refer to these levels in order to indicate the rank of ancestry of genes, pathways and functions involved in cancer, motivated by the observation that orthologous genes seem more likely to retain ancestral gene functions [75].

We studied 1,895 genes whose mutations are causally implicated in cancer. These genes are annotated as either oncogenes or tumor suppressor genes accordingly if the evidence indicates that the mutation promotes or suppresses cancer cell growth. This list was compiled by merging the well documented Catalogue Of Somatic Mutations In Cancer (COSMIC) [56], the more recent list of cancer driver genes published by Vogelstein's team [73] and the lists of tumor suppressor genes [76] and oncogenes [57] published by M Zhao and collaborators. From these gene lists, only those genes that have associated EGGNog phylostratigraphic ages were chosen to yield 980 oncogenes and 915 tumor suppressor genes.

The gene expression data of RNA-Seq in TPM (Transcripts Per Million), which used MapSplice to do the alignment and RSEM to quantify the gene expression leves [77], were retrieved from The Caner Genome Atlas (TCGA) project data portal (https://gdc.nci.nih.gov/). We focused our analysis on ten cancer types: breast invasive carcinoma (BRCA), Uterine Corpus Endometrial Carcinoma (UCEC), Colon adenocarcinoma (COAD), Head and Neck squamous cell carcinoma (HNSC), Kidney renal clear cell carcinoma (KIRC), Kidney renal papillary cell carcinoma (KIRP), Liver hepatocellular carcinoma (LIHC), Lung adenocarcinoma (LUAD), Lung squamous cell carcinoma (LUSC), Thyroid carcinoma (THCA).

We curated microarray gene expression profiles from 50 samples of embryonic stem cells (ESC) from GEO (www.ncbi.nlm.nih.gov/geo/). Their GSE accession numbers are GSE6561, GSE7234, GSE7896, GSE9086, GSE9440, GSE9709, GSE9832, GSE9865, GSE9940, GSE12390, GSE13828, GSE14711, GSE14897, GSE15148, GSE16654 and GSE18679. We curated microarray gene expression profiles from 73 samples of iPS cells from GEO. Their GSE accession numbers are GSE9832, GSE9709, GSE9865, GSE12390, GSE12583, GSE13828, GSE14711, GSE15148, GSE16654 and GSE18111. We also curated microarray gene expression profiles from 159 samples of breast epithelial cell lines of primary cell cultures from GEO. Their GSE accession numbers are GSE10780, GSE9649, GSE12917 and GSE13671.

The cancer-related pathway data was obtained from the KEGG PATHWAY Database [63] (www.kegg.jp/). We downloaded lists of human genes involved in pathways associated with cancer (all pathways included under the KEGG category ID 05200 Pathways in Cancer) and 14 essential housekeeping pathways. These two sets of pathwasy are not mutually exclusive. The list of pathways is shown in tables 1 and 2. The sets of genes were further subset to aged genes, thus all the human genes involved in each pathway that have a corresponding EGGNog age level are ultimately considered in our analysis.

We define gene age enrichment score for an observed gene list as the adjusted residual [78] calculated in the following way:

• 1.  
  Calculate the distribution of ages of the background gene list, as the frequency ${{F}_{{\rm bg}}}\left(a \right)$ of genes in each age group $a$
• 2.  
  Calculate the expected frequency of observed ages according to the usual null hypothesis in a contingency table:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{F}_{{\rm exp}}}\left(a \right)={{F}_{{\rm bg}}}\left(a \right)\times {{N}_{{\rm obs}}}\nonumber \end{align}$$

  where ${{N}_{{\rm obs}}}$ is the number of genes in the observed list (e.g. tumor suppressor genes, oncogenes, etc)
• 3.  
  For each set 1,000 uniformly random samples of the same size as the observed list are drawn without replacement from the background list. Age frequencies $F_{{\rm trial}}^{\{i\}}\left(a \right)$ are calculated for each trial $i$.
• 4.  
  We use the distribution of frequency values for each age group across sample trials ${{G}_{{\rm a}}}\left(F \right)$ to statistically test the likelihood of observed frequencies.
• 5.  
  In particular, it's true that the expected value of the frequency for each age group ${{F}_{{\rm exp}}}\left(a \right)=\langle {{G}_{{\rm a}}}\left({{F}_{i}} \right) \rangle$, with $\langle X \rangle$ the mean across trial samples.
• 6.  
  The adjusted residual is the deviation of observed from expected frequencies normalized by 2 times the standard error of the bootstrapped trial samples ${{G}_{{\rm a}}}\left({{F}_{i}} \right)$ in each age group $a$: $\Delta \left(a \right)=\frac{1}{2}\left(\frac{{{F}_{{\rm bg}}}\left(a \right)-~{{F}_{{\rm exp}}}(a)}{S(a)} \right),$ with $S\left(a \right)=\sqrt{\frac{1}{N}}~\times ~\sigma \left({{G}_{{\rm a}}}\left({{F}_{i}} \right) \right)$ the standard error across trial samples and $N=1000$ the number of trials.

This standardized measure takes into account the individual sizes of the samples [78]. The value of the adjusted residual can be directly interpreted in a probabilistic way, a value of  ±1 corresponds to the standard 2-sigma test with p-value ~0.05. Hence adjusted residuals larger than 1 indicate that the observed value is larger than expected with statistical significance, i.e. the corresponding age group is over-represented, and equivalently adjusted residuals smaller than  −1 indicate that the corresponding age group is under-represented. Thus, the adjusted residual is used as a statistical score for gene age enrichment.

Additionally for each gene set we performed a direct comparison of the age distribution with the background set (e.g. all human genes) using a standard Kolmogorov-Smirnov (KS) test. When the associated p-values obtained from the KS-statistic is  <0.05, the gene set is not likely to have the same age distribution as the background.