BACKGROUND
rna-seq allows the measurement of the gene expression levels in a manner far more precise than previous methods.
studies using this approach have already altered our view of the extent and complexity of the eukaryotic transcriptomes  <cit> .
important analyses have been performed based on transcriptomic data obtained from multiple human tissues.
the fantom project, for example, presents the gene expression profiling of  <dig> human healthy tissues associated with the functional annotation of the mammalian genomes  <cit> .
the human protein atlas applies rna-seq to  <dig> human tissues to find the correlation between gene expression and protein presence in the characterization of several parts of the human proteome: membrane proteome, druggable proteome, cancer proteome, secretome and proteome involved in metabolic processes  <cit> .
the genotype-tissue expression   <cit>  provides an extensive resource of human transcriptomic data permitting the inference of different patterns across human tissues and individuals  <cit> , and the construction of tissue-specific gene co-expression networks  <cit> .
other studies are more focused on specific cells or tissue types, like the article entitled “a comprehensive analysis of the human placenta transcriptome” that characterizes the transcriptome of placenta from  <dig> healthy women with uncomplicated pregnancies using rna-seq  <cit> .
the assembly of comprehensive maps of the human transcriptome is essential for a clear identification of the functional elements of our genome and to reveal the molecular constituents of different cells and tissues  <cit> .
despite the accomplishment of many transcriptomic studies, little has been reported until recently about the evolutionary determinants of human cell identity, particularly from a joint perspective of protein evolution and gene expression  <cit> .
the evolutionary history of a gene can be very informative about its function.
gene age  is an important piece of information that can be inferred in different ways and has been used in some genome-scale studies and in some studies on gene families  <cit> .
phylostratigraphy is the usual methodology applied to find the origin and emergence of genes  <cit> .
previous phylogenetic studies showed that the evolutionary history of different coding parts of the genome have relations with diseases  <cit> , codon usage  <cit> , essentiality, interactions  <cit> , stemness and self-renewal  <cit> .
other studies have shown that older genes evolve slower  <cit> , encode longer proteins, present higher expression levels, possess higher intron density and are subject to stronger purifying selection  <cit> .
several of these studies approach the question of gene age in different ways, but most of them are not focused on human genes or do not applied phylostratigraphy using large-scale genomic data.
in this work, we address the key question about the age of the human genes using a combination of complementary genome-wide data and databases that are analysed and integrated to achieve a map of the human genes on the evolutionary time-scale.
this integrative analysis uses a comprehensive human tissue rna-seq dataset that allows deep expression profiling of protein-coding genes  <cit> ; a database of orthologous proteins that allows finding the oldest relatives to each human protein along different species  <cit> ; the taxonomy mapping of genes to lineage clades from the ncbi taxonomy database ; and the time-scale mapping provided by timetree resource   <cit> .
the generation of clusters of orthologous proteins built along multiple species is more accurate for phylogenetic studies than simple sequence homology search  <cit> , because it does not look for singular best homologous but it implies a conservation along the evolutionary tree and can provide a way to date the origin of protein modules involved in specific functions or in specific biological pathways  <cit> .
finally, to complete the view of the protein-coding gene phylostratigraphy, we use the genome-wide expression data to produce a human gene network based on a coexpression analysis of the transcriptomic rna-seq profiles along multiple tissues and identify the genes that can be considered house-keeping  or tissue-enriched .
the allocation of these gene subsets  on the evolutionary time map shows a clear difference in gene age, indicating that house-keeping genes are older.
methods
gene expression data from human normal tissues
the genome-wide expression dataset used in this work corresponds to a series of rna-seq analyses performed with illumina hiseq  <dig> paired end sequencing on cdna libraries prepared from samples of  <dig> human individuals from  <dig> different tissues  <cit>  .
the data provided reads for  <dig>  genes detected per sample, and  <dig>  of these genes showed relevant expression signal corresponding to mean ≥  <dig> in all the selected tissues.
after normalization and comparative analysis of the expression distributions of the samples from this dataset, we selected a total of  <dig> samples with two to five biological replicates for the following  <dig> tissues: adrenal gland , appendices , bone marrow , brain , colon rectum , duodenum , endometrium , esophagus , fallopian tube , adipose tissue , gallbladder , heart , kidney , liver , lung , lymph node , ovary , pancreas , placenta , prostate , rectum , salivary gland , skeletal muscle , skin , small intestine , smooth muscle , spleen , stomach , testis , thyroid , tonsil  and urinary bladder .
expression profiling and coexpression data analyses
rna-seq expression data from all the tissue samples, taken as normalized fpkm  from  <cit> , were log <dig> transformed to obtain the final expression values as: log <dig>  normalized expression distributions of these samples can be seen in additional file 1: figures s <dig> and s <dig>  unsupervised clustering of the samples based in all genes expression was done using agglomerative hierarchical clustering  and calculating the distances based on: .
this clustering was done for all the  <dig> samples  and just for the  <dig> tissues using the average expression of the replicates .
the coexpression dataset was built calculating the pair-wise spearman correlation coefficient  of all the genes along the  <dig> samples and only selecting, as positive gene-pairs, the ones with a correlation coefficient ≥  <dig> .
crossvalidation of these correlation values was applied by a random selection of two sample replicates from each tissue  and recalculating again the spearman correlation for these random subsets of the data.
this sampling was done  <dig> times, annotating for each gene-pair the number of times that its r coefficient was ≥  <dig> .
only the gene-pairs validated  <dig> times in this sampling were selected.
in this way, a final set of highly correlated gene-pairs was produced that included  <dig> genes and  <dig>  coexpression interactions.
this coexpression dataset is provided as additional file  <dig>  indicating the names of all the gene-pairs and their correlation value.
a gene coexpression network derived from the coexpression data was built using cytoscape  and we applied the mcode algorithm to identify clusters inside the network.
this algorithm performs an analysis of the topology of the network to find densely connected regions that define modules.
the coexpression network built with cytoscape including all the subnetworks found , as well as the parameters derived from the graph analysis, is provided in additional file  <dig> fig.
<dig> heatmap presenting the comparison of the transcriptomic profiles of  <dig> human tissues.
unsupervised clustering analysis of the global expression correlation along  <dig> samples of  <dig> human normal tissues.
the main plot presents the color heatmap and dendrogram produced by the comparison of all the expressed genes  along all the samples.
the pair-wise distances between samples are calculated using the values of:  <dig> – spearman correlation coefficient.
the color bar with scales shows dark-red corresponding to minimum distances  and dark-blue to maximum distances
evolutionary analyses
for the evolutionary analysis and determination of lowest common ancestor , we used a database of orthologous proteins: orthologous matrix   <cit> .
oma includes a database and resource with methods for the inference of orthologous among complete genomes  <cit> .
we downloaded the oma database into a local mysql database and created python scripts to search for the ensembl id's from our transcriptomic data into this local database and to calculate the lca into each respective ortholog group.
in this way we obtained a table with the number of protein-coding genes assigned to each clade in the human taxonomy lineage, as defined by the taxonomy resource from ncbi .
furthermore to produce a time-scale, we integrated these data with the timetree of life   <cit> , that includes the tree of living species calibrated to time.
the analysis of the number of genes placed along the evolutionary time-scale allowed visualization of the profile of human genes origin for the whole genome  or for specific subsets of genes.
in this genes/time profile we calculated the number of protein-coding genes that had lca corresponding to each clade  in the taxonomy lineage –that for human includes  <dig> consecutive levels– and we identified certain levels where major changes occur.
these are taken as most significant stages and proposed as key evolutionary hallmarks including specific sets of the human protein-coding genes that are identified.
functional enrichment analysis and identification of gene modules
for the functional enrichment analysis we used david   <cit>  and geneterm-linker   <cit>  bioinformatic tools with the list of genes from each evolutionary stage level of the human lineage.
in all cases the enrichment analyses were done using a hypergeometric test and adjusting the p-values for multiple testing with the benjamini-hochberg procedure  <cit> .
in the same way, we also investigated the functional enrichment of the subnetworks generated by the clusters and modules found in the analysis of the gene coexpression network.
general calculations and statistics
all the data analyses and graphics have been produced in the r statistic environment.
general functions as boxplot, image, qplot  or wilcox.test, have been applied over the different data presented.
some specific methods or algorithms are cited along different sections of this manuscript.
RESULTS
human global transcriptome profile based on expression in multiple tissues
the unsupervised clustering analysis of the global expression correlation along  <dig> samples of  <dig> human normal tissues is presented in fig.
<dig>  this plot shows the color heatmap and dendrogram produced by the comparison of all the expressed genes  along all samples.
the pair-wise distances were calculated as: ; that is a non-parametric approach to calculate similarity and distance.
in this way, the heatmap shows a clear relationship between the samples from the same tissue  and also presents the proximity between the tissues that have strong biological and physiological links, such as: spleen, lymph nodes and tonsils ; or stomach, duodenum, small intestine, colon and rectum .
a color bar with scales for the heatmap is included in fig.
<dig>  to indicate that dark-red corresponds to minimum distance  and dark-blue to maximum distance .
white color corresponds to medium values and the distributions inside the color bars show that the major part of the comparisons have values around r =  <dig> – <dig>  .
the results also reveal that some tissues have transcriptomic profiles very different to the rest, being testis the most different one that produces a clearly separated branch in the dendrogram.
house-keeping and tissue-enriched genes
the transcriptomic data allowed creation of two relevant subsets of genes: house-keeping  and tissue-enriched  genes.
figure 2a shows the number of genes expressed, with fpkm equal or higher than  <dig>  per number of tissues.
this calculation allowed the identification of  <dig> ubiquitous genes present in at least one replicate of all tissues, and  <dig> genes expressed in all the  <dig> samples .
the intersection of these  <dig> genes with a curated dataset of  <dig> house-keeping genes created by eisenberg et al.
<cit>  gave a total of  <dig> hkg , indicating a large overlap of 93 %.
to identify how relevant was this, we calculated the odds ratio  of such overlap and obtained a very significant value of or =  <dig>  .
for the tissue-enriched genes we explored the other side of the data in fig.
2a and considered just the genes that were expressed in only one, two or three tissues .
we did not take only one, but also two or three tissues, because some tissues are physiologically very related and in fact presented high correlation between them, for example: colon and rectum; small intestine and duodenum, etc.
these observations indicated that such tissues share a large number of common genes in their expression profile.
finally, a global comparison of the expression distributions of hkg versus teg indicated that the tissue-enriched genes showed significantly lower expression values than the house-keeping genes .
to demonstrate this difference, we made two statistical tests, t-test and wilcoxon.
in both cases, the p-values were very low  and the difference between the mean expressions of teg and hkg was  <dig>  .
thus, we can conclude that this difference is not by chance but due to a real different expression regulation that should be considered in further analyses.
these boxplots also show that the global variability of teg  was larger than the global variability of hkg.
this observation indicates that the teg present a larger range of expression values along different tissues and the hkg have a tighter regulation.fig.
<dig> number of genes expressed along  <dig> tissues derived from rna-seq transcriptomic data.
a plot showing the number of genes expressed, with fpkm equal or higher than  <dig>  per number of tissues.
the house-keeping genes are labeled as hkg and tissue-enriched genes as teg.
b comparison of the expression distributions of hkg versus teg.
c venn diagrams showing the intersection of  <dig> genes  with the dataset of  <dig> house-keeping genes obtained from eisenberg and levanon   <cit>
human gene hallmarks on the evolutionary time-scale
for the evolutionary analysis, we used a phylostratigraphic method for reconstruction of macro-evolutionary trends based on the principle of “founder gene” formation  <cit> .
typically, these methods first identify the homologues of a given gene and then use the divergence between the two most distant to determine the gene age.
historically, such studies have been using blast  <cit>  for homology searches.
however, this approach was shown to introduce some biases into the analyses  <cit> .
another approach is to use orthologous groups to determine the age of a gene.
orthologues are believed to be functionally more similar than paralogues  <cit>  and by definition they trace back to an ancestral gene that was present in a common ancestor of the compared species  <cit> .
the parameters used for clustering orthologous groups affect the age estimations for a gene; for example, restrictive parameters tend to limit the set of possible progenitors  <cit> .
nevertheless, both approaches used for dating the gene origin depend on the correct identification of homologues and/or orthologues, but in the second case the accurate reconstruction of orthologous families imposes a higher stringency and it implies a conservation along the evolutionary clades.
for our analysis, we identified the group of orthologues that corresponded to each of the  <dig>  genes detected in the transcriptomic study, mapping them to the corresponding human protein-coding genes in the oma database  <cit>  , and then assigning the lowest common ancestor  to each human protein according to its orthologous family.
the use of oma in comparison with blast homology approach gives us a more detailed view of gene origin since, as we indicated above, it uses a more restrictive grouping method.
thereby, the number of genes dated on ancient clades is lower.
once we identified the lca for each human protein-coding gene we assigned such protein/gene to the corresponding taxonomy level in the human lineage as defined in ncbi , that includes  <dig> taxonomic groups as consecutive clades from the first one, named cellular organisms, to the last one homo sapiens.
figure  <dig> presents these  <dig> taxonomic clades placed along the time-scale  from the origin to present.
for each one of these taxonomy levels we represented the cumulative percentage of proteins that are dated at such level.
this is done in the following way: first, plotting all the genes mapped to oma proteins ; second, the same plot is produced but including only the proteins that correspond to house-keeping genes ; third, plot including only the tissue-enriched genes .fig.
<dig> evolutionary hallmarks of human protein-coding genes along time-scale.
plot presenting the number of human protein-coding genes  that are assigned to each of the  <dig> taxonomic clades .
for each one of these taxonomy levels the graph represents the cumulative percentage of protein-coding genes that are dated at such level.
in this way, the  <dig> taxonomic clades are placed as dots along the time-scale  from the origin to present.
line in black includes all the  <dig>  proteins derived from the mapping of expressed genes in oma.
line in blue includes only the hk genes:  <dig>  line in red includes only the te genes: 2157
the analyses of these plots obtained with the phylostratigraphic method revealed the presence of some major differential steps on the emergence of protein-coding genes along the evolutionary time-scale from origin to present.
looking at all the expressed coding genes, we can see the global evolutionary profile of the organism , but along this profile we can identify some more prominent steps in the accumulated relative number of genes  along time.
for example, looking at the plot of the hk a large increase is observed at the start of the curve, from the first taxonomy level  to second  taxonomy level .
by contrast, the te have the major emergence of genes much later  .
the complete timeline includes  <dig> phylogenetic clades , but by analysing these points it was possible to identify eight major steps or stage levels  that appear on the human gene profile along evolutionary time.
moreover, we could assign the number of human protein-coding genes that emerged along each one of these eight stages for the three categories reported: all the expressed protein-coding genes mapped to oma, the hk and the te.
the eight stage levels identified are: st1)
cellular organisms ; st2)
cellular organisms to eukaryota; st3)
eukaryota to metazoa; st4)
metazoa to vertebrata; st5)
vertebrata to euteleostomi; st6)
euteleostomi to mammalia; st7)
mammalia to primates; and st8)
primates to homo sapiens.
all the numbers corresponding to the human protein-coding genes assigned to each of these eight evolutionary hallmarks are included in fig.
<dig>  that indicates how many are in each stage either considering the complete human gene set or just the hk or the te.
all the information about each one of the  <dig>  human protein-coding genes including the assignment to stages is also provided as additional file  <dig> fig.
<dig> evolutionary hallmarks of human protein-coding genes along time-scale.
table showing the numbers corresponding to the human proteins found at each of the eight major evolutionary steps depicted from placing them on the time-scale.
these eight stages provide the evolutionary hallmarks of the human proteome.
the stage levels correspond to: st1)
cellular organisms ; st2)
cellular organisms to eukaryota; st3)
eukaryota to metazoa; st4)
metazoa to vertebrata; st5)
vertebrata to euteleostomi; st6)
euteleostomi to mammalia; st7)
mammalia to primates; and st8)
primates to homo sapiens.
the numbers are presented for all the  <dig>  protein-coding genes, for the  <dig> hk genes and for the  <dig> te genes
the analysis of the hallmarks also reveals that the hk genes are more ancient than te genes.
the plots in fig.
<dig> and the data in fig.
<dig> show clear difference between them.
the hk genes present a major increase or expansion in stage  <dig> , with  <dig> genes and a change of ≈ 30 % with respect to the total.
by contrast, the te genes show a major increase in stage  <dig>  with  <dig> genes and a change of ≈ 37 %.
these observations seem to indicate that house-keeping genes emerged early in evolution and are older in age, knowing that they reflect more essential and constitutive functions.
the idea that gene essentiality is associated with older genes has been reported in several studies, for example on yeast  <cit>  and mammalian genes  <cit> .
by contrast, the observations that human tissue-specific genes had emerged later in evolution may reveal that human specific cellular or physiological roles are implemented at molecular level by the appearance of newer genes.
gene age data comparison
as we indicated above, there are some studies that use the phylostratigraphic method to explore the age of human genes, but most of these studies use sequence similarity search  to look for the oldest homologues to the human  <cit> .
to compare the results on human gene age assignment done in this work with other available age assignments, we took the published data from domazet-lošo   <cit>  and from neme   <cit> , and we represented the information about allocation to lowest common ancestor  of the human genes in phylogenetic clades of the evolutionary tree.
the assignments were done using  <dig> common phylostratum to allow the comparison of the data.
the results of this comparison are included in additional file 1: figure s <dig> and they show a general similarity but some important differences.
the most significant difference corresponds to the fact that both domazet-lošo  <cit>  and neme  <cit>  placed a very large number of genes on the first stage of the evolutionary time-scale that goes from the origin of life to first cellular organisms :  <dig> of  <dig>    <cit>  and  <dig> of  <dig>    <cit> .
this result denotes a bias that, as we indicated above, can be due to the methodology of using the homology search approach.
in any case, the idea suggested that one third of the human proteoma may have emerged in evolutionary time before the origin of eukaryotic cells needs deeper studies and it is not what we observed in our analyses.
another important difference is that the proposed age mapping allocates the largest number of genes, first, to the chordata-vertebrata-euteleostomi phylostratums  and, second, to the mammalia-eutheria  .
from the evolutionary point of view these results make a lot of sense since the time-scale of life  <cit>  reveals two large expansions of the species precisely around the vertebrates time  and around the time of the mammals appearance .
these expansions are well reflected in our time-scale profile .
finally, it is important to indicate that the age mapping presented in our study only considers human protein-coding genes that are included in orthologous families  and in this way it has a lower coverage over human genes than the other reported studies which include more than  <dig>  genes in each case  <cit> .
functional enrichment of the genes at different evolutionary hallmarks
we performed functional enrichment analyses of the sets of protein-coding genes included in each one of eight major stages found in the evolutionary study.
the full results of these analyses are provided as additional file  <dig>  in all the stages the functional enrichment makes clear biological sense and provides a strong support to the allocation of many biological processes in evolutionary time.
in brief, we comment and discuss below some interesting functions enriched in each stage.
stage- <dig>  from the origin of life to first cellular organisms .
this stage comprises the genes occurring over two major domains of life: archaea and bacteria.
determining the lca, our data shows that human has  <dig>  %  of the protein-coding genes assigned to prokaryotic age.
prokaryotes are organisms that lack both membrane-bound organelles and nucleus.
functional enrichment analysis showed that this stage involved many basic metabolic processes like glycolysis , the krebs cycle , and lipid oxidation .
the enrichment also shows the appearance of the oldest cellular organelle, the mitochondria, and the oldest macromolecular machine, the ribosome, that are well reported to be dated to prokaryotic times.
stage- <dig>  cellular organisms to eukaryota.
according to basic literature, the defining feature of eukaryotic cells is that they have membrane-bound organelles, especially the nucleus, which contains the genetic material, and is enclosed by the nuclear envelope.
protists, fungi, animals, and plants all consist of eukaryotic cells.
eukaryotic cells also contain other membrane-bound organelles such as the golgi apparatus.
eukaryotic organisms can be unicellular or multicellular.
the functional enrichment analysis for the  <dig> genes that emerged along this stage showed well the formation of the principal complexes expected in eukaryotes.
in this way it is remarkable the enrichment on nuclear pore proteins, nuclear import proteins, nucleosome and chromatin proteins, as well as many proteins involved dna and rna activity: mrna and rrna processing, mrna splicing, dna unwinding, dna polymerase, dna/rna helicase.
this stage also marks in time the appearance and biogenesis of some major molecular complexes: the proteasome , the spliceosome, and the ribosome .
stage- <dig>  eukaryota to metazoa.
the third stage comprises organisms from opisthokonta and metazoan clades with  <dig> protein-coding genes .
the opisthokonts are a broad group of eukaryotes, including both the animal and fungi kingdoms, sometimes referred to as the “fungi/metazoan group”  <cit> .
this stage comprises metazoan, fungal and protistan taxa, and other multicellular taxa   <cit> .
they also include known fungi and/or parasites of plants like: ascomycota, basidiomycota, chytridiomycetes, glomeromycota, microsporidia, urediniomycetes, ustilaginomycetes and zygomycota  <cit> .
according to our functional enrichment analysis, this stage involves different genes responsible for signal transduction like the gtpases.
some of the enriched terms are: pyrophosphatase activity, nucleoside-triphosphatase activity, gtp binding, transferring phosphorus-containing groups.
all these functions indicate that it may be the time when phosphorus and phospate acquired a key role in protein function and regulation.
other enriched functions, like post-translational protein modification, calcium-binding ef-hand, protein transport and localization also indicate cellular protein regulation.
stage- <dig>  metazoa to vertebrata.
this stage includes organisms from eumetazoa, bilateria, deuterostomia, chordata, craniata and vertebrata with  <dig> genes .
the main novelties of this stage are the appearance of protein kinase activity, and the presence of growth factors and some specific signaling proteins like wnt.
all biochemically characterized members of the wnt superfamily encode enzymes that transfer organic acids, typically fatty acids, onto hydroxyl groups of membrane-embedded targets  <cit> .
other enriched terms in this stage, like sarcomere and contractile fiber part, may indicate the emergence of the muscular structures present in vertebrates  <cit> .
stage- <dig>  vertebrata to euteleostomi.
this stage represents the largest step in the human lineage according to the number of protein-coding genes assigned , that correspond to a 29 % of the total.
the stage comprises organisms from gnathostomata , teleostomi  and euteleostomi   <cit> .
the enrichment analysis shows a large functional expansion including new biological systems, like the neural and the vascular-circulatory systems, represented in enriched terms like: neurogenesis, neuron differentiation, axogenesis, voltage-gated channels, neuromuscular junction development, blood vessel development, vasculature development, mesenchymal cell development and differentiation, etc.
many other genes are assigned to biological regulation and regulation of cellular processes; including cell death and apoptosis.
finally, the appearance of the large family of homeobox proteins seems to be placed at this stage.
stage- <dig>  euteleostomi to mammalia.
in this stage there are organisms from sarcopterygii   <cit> , dipnotetrapodomorpha , tetrapoda , amniota  and up to mammalia clades.
with  <dig> genes at this stage, the human lineage achieves 80 % of its gene composition.
the most relevant enriched terms are related to the hematologic system, marking the appearance of the leukocytes and the lymphocytes.
previous phylogenetic analyses based on gene expression data also placed the date of many proteins from leukocytes around the time of the mammals’ clade  <cit> .
stage- <dig>  mammalia to primates.
this stage comprises clades of theria, eutheria, boreoeutheria, euarchontoglires and primates, representing organisms that give birth to live young without using a shelled egg up to placental mammals  <cit> .
there are  <dig> genes emerged on this stage, adding up to  <dig>  % of the cummulative profile in the human gene lineage.
a large amount of these genes are enriched in the terms: regulation of gene expression and transcription.
other more specific terms are related to the skin  or with the sexual reproductive system .
this stage also includes a family of cytochrome p <dig> proteins  and the mammalian defensins : defa1b, defa <dig>  defa <dig>  defa <dig>  defa <dig>  and defb4a.
defensins are a family of antimicrobial peptides and vital contributors to host immune response.
being constitutive or inducible expressed genes, they have been shown to contribute to innate host defense via direct bactericidal activity, as well as to adaptive immunity through effector and regulatory functions  <cit> .
stage- <dig>  primates to homo sapiens.
the last stage of human development, with  <dig> genes, presents a group of quite specific functions played by specific protein families, such as: somatotropin hormone, cytochrome p <dig>  gtpase activator activity, defense response to fungus and bacterium provided by histatins.
his <dig> and his <dig>  have been found only in saliva of humans, macaques and some other primates but not in any other mammals  <cit> .
they are a family of histidine-rich polypeptides that probably function as part of the non-immune host defense system and appeared very late in evolution  <cit> .
cytochromes p <dig> constitute a superfamily of proteins that existed in virtually all species from prokaryotes to humans.
most of these proteins in the cyp <dig>  cyp <dig>  cyp <dig> and cyp <dig> families encode enzymes involved in the metabolism and elimination of potential toxic compounds like drugs or foreign xenobiotics, and are inducible by various environmental stimuli  <cit> .
this last stage includes a small subset of six cytochromes p <dig> that seem to be very specific of the primates-human clades: cyp2a <dig>  cyp2c <dig>  cyp2d <dig>  cyp2j <dig>  cyp2s <dig> and cyp3a <dig>  the appearance late in evolution of some of these genes may reflect their functional specificity and it is known that they play a key role in human health  <cit> .
network analysis reveals evolutionary age conservation of coexpressed proteins
the global coexpression analysis of the human protein-coding genes allowed the construction of a network including highly correlated protein pairs.
the integration of these data with the data from the evolutionary analysis –that provided the identification of the eight stages along evolutionary time– did allow mapping the age of the genes on the network according to such stages.
these results are presented in fig.
<dig> that shows a complex network –like a galaxy– that includes  <dig> protein nodes associated by  <dig>  interactions.
this network corresponds to a subset of the coexpression data, build as indicated in methods, which included  <dig> proteins and  <dig>  interactions .
the subset is done with only the groups that had at least five linked proteins, since we wanted to provide in fig.
<dig> a visible representation of the network with a clear color mapping of the eight stages that were identified in the evolutionary study.
the colors of the stages are also presented in the illustrated table, fig.
<dig>  to allow better identification of the number and % of proteins at each age hallmark.
the analysis of the network  done with the algorithm mcode revealed the existence of  <dig> major subnetworks, which can be considered as major constellations in the galaxy of relational nodes.
the color code of this large graph indicates that there is enrichment in certain colors in each subnetwork.
to prove this in a more accurate way, we built a graphic representation for each one of the  <dig> subnetworks found to show the proportion of proteins assigned to each of the eight evolutionary stages with their corresponding color code .
this graphic is presented in additional file 1: figure s <dig>  that shows each subnetwork with its specific color pattern, indicating that there are always some predominant colors: subnetwork  <dig> is the oldest with red predominant colors and subnetwork  <dig> is the newest with blue predominant colors.
as a conclusion, these results revealed that in the groups of highly coexpressed proteins there is a tendency to include proteins of the same evolutionary age.fig.
<dig> human coexpression network mapping the evolutionary age on highly correlated nodes.
representation of the coexpression complex network –like a galaxy– that includes  <dig> protein nodes related with  <dig>  interactions.
this network corresponds to a subset of the larger coexpression network produced, which included  <dig> proteins and  <dig>  interactions .
the subset is done to include only the groups that had at least five linked proteins.
the color mapping of the nodes correspond to the eight stages that were identified in the evolutionary study .
the network also includes numbers for  <dig> major subnetworks –clusters of closely related proteins that include more than  <dig> nodes– considered as major constellations in the galaxy of nodes.
three panels on the right show an enlarged view of three subnetworks corresponding to: ribosomal proteins , mitochondrial proteins  and angiogenesis proteins
finally, we did a functional enrichment analysis of the proteins forming the  <dig> subnetworks which again showed a coherent biological enrichment in specific functions:  immune response;  cell cycle;  cytoskeleton;  rna splicing;  ribosome;  extracellular matrix;  muscle and contraction;  gametes and reproductive process;  cell junction and cell adhesion;  mitochondria and atp synthesis;  angiogenesis and vasculogenesis.
more detailed results for this analysis are included in additional file 1: figure s <dig>  thus, we observed that age-related proteins are predisposed to present expression coregulation and to have close functional links.
combining age data, functional data and coexpression data can provide a deeper view about the links and roles of the human protein-coding genes.
in this way, we observed for example that subnetwork  <dig>  presented as expected an overwhelming majority of ancient genes from the prokaryotic or eukaryotic age .
on the other hand, subnetworks  <dig> ,  <dig>  and  <dig>  showed a higher proportion of recent genes dated after vertebrata .
these results agree with studies based on yeast protein physical interaction networks, arguing that proteins preferentially interact with proteins of same age and origin  <cit> .
moreover, it was previously shown that coexpression networks can be conserved over the evolutionary history, and these genes tend to be functionally related and provide selective advantages  <cit> .
it has been also reported that coexpression networks are found associated to functions like cell adhesion, cell cycle, dna replication and dna repair  <cit> , and this is in agreement with functions found enriched in subnetworks of our analyses: subnetwork  <dig> and  <dig> .
CONCLUSIONS
the transcriptomic study of the human gene expression distributions and profiles along  <dig> tissues provided a global mapping of the activity of most human genes and of the links between them, showing the expected association of samples from common physiological regions, i.e.
: the gastrointestinal tract , the hematopoietic and lymphatic system  and the muscle .
the evolutionary study of the human protein-coding genes placed them in the time-scale of the living species and revealed eight distinct hallmarks along such time-scale , showing that the hk genes are more ancient than the te genes.
the hk genes present the major emergence in stage  <dig> of the evolutionary profile, while the te genes have the major emergence in stage  <dig>  the functional enrichment study found coherent groups of terms and annotations assigned to the genes placed at each evolutionary stage.
for example, in stage  <dig> there were many functional terms on essential metabolic processes, like aerobic respiration and mitochondrial activity; and in stage  <dig> there were enriched functions related to the nucleus and genome regulation, like chromatin and nucleosome assembly, dna replication, mrna processing.
finally, the study of the pair-wise correlation of the gene expression profiles along tissues allowed building human gene coexpression networks and find modules with functional and biological meaning.
the mapping of the age of the protein-coding genes on these networks demonstrated the existence of tight links between age-related proteins.
additional files
additional file 1: additional pdf file that includes six supplementary figures.
additional file 2: table  including the values of the gene pair-wise spearman correlation and the cross-validation from the coexpression analysis, as well as the gene ids to allow the reconstruction of a human coexpression network with  <dig> proteins and  <dig>  interactions.
additional file 3: cytoscape file  that includes the complete coexpression network that is produced in this work and presented in fig.
<dig> of the main article.
this is a file to be opened only with cytoscape program.
additional file 4: table  including the data and ids of the  <dig>  human protein-coding genes included in each of the eight evolutionary stages.
additional file 5: table  including the results of the functional enrichment analysis of the  <dig>  human protein-coding genes included in each of the eight evolutionary stages.
the results for each of the eight stages are presented in different spreadsheets.

 rna-seq allows the measurement of the gene expression levels in a manner far more precise than previous methods. studies using this approach have already altered our view of the extent and complexity of the eukaryotic transcriptomes. important analyses have been performed based on transcriptomic data obtained from multiple human tissues. the fantom project, for example, presents the gene expression profiling of human healthy tissues associated with the functional annotation of the mammalian genomes. the human protein atlas applies rna-seq to human tissues to find the correlation between gene expression and protein presence in the characterization of several parts of the human proteome: membrane proteome, druggable proteome, cancer proteome, secretome and proteome involved in metabolic processes. the genotype-tissue expression provides an extensive resource of human transcriptomic data permitting the inference of different patterns across human tissues and individuals, and the construction of tissue-specific gene co-expression networks. other studies are more focused on specific cells or tissue types, like the article entitled “a comprehensive analysis of the human placenta transcriptome” that characterizes the transcriptome of placenta from healthy women with uncomplicated pregnancies using rna-seq. the assembly of comprehensive maps of the human transcriptome is essential for a clear identification of the functional elements of our genome and to reveal the molecular constituents of different cells and tissues. despite the accomplishment of many transcriptomic studies, little has been reported until recently about the evolutionary determinants of human cell identity, particularly from a joint perspective of protein evolution and gene expression. the evolutionary history of a gene can be very informative about its function. gene age is an important piece of information that can be inferred in different ways and has been used in some genome-scale studies and in some studies on gene families. phylostratigraphy is the usual methodology applied to find the origin and emergence of genes. previous phylogenetic studies showed that the evolutionary history of different coding parts of the genome have relations with diseases, codon usage, essentiality, interactions, stemness and self-renewal. other studies have shown that older genes evolve slower, encode longer proteins, present higher expression levels, possess higher intron density and are subject to stronger purifying selection. several of these studies approach the question of gene age in different ways, but most of them are not focused on human genes or do not applied phylostratigraphy using large-scale genomic data. in this work, we address the key question about the age of the human genes using a combination of complementary genome-wide data and databases that are analysed and integrated to achieve a map of the human genes on the evolutionary time-scale. this integrative analysis uses a comprehensive human tissue rna-seq dataset that allows deep expression profiling of protein-coding genes; a database of orthologous proteins that allows finding the oldest relatives to each human protein along different species; the taxonomy mapping of genes to lineage clades from the ncbi taxonomy database; and the time-scale mapping provided by timetree resource. the generation of clusters of orthologous proteins built along multiple species is more accurate for phylogenetic studies than simple sequence homology search, because it does not look for singular best homologous but it implies a conservation along the evolutionary tree and can provide a way to date the origin of protein modules involved in specific functions or in specific biological pathways. finally, to complete the view of the protein-coding gene phylostratigraphy, we use the genome-wide expression data to produce a human gene network based on a coexpression analysis of the transcriptomic rna-seq profiles along multiple tissues and identify the genes that can be considered house-keeping or tissue-enriched. the allocation of these gene subsets on the evolutionary time map shows a clear difference in gene age, indicating that house-keeping genes are older. methods gene expression data from human normal tissues the genome-wide expression dataset used in this work corresponds to a series of rna-seq analyses performed with illumina hiseq paired end sequencing on cdna libraries prepared from samples of human individuals from different tissues. the data provided reads for genes detected per sample, and of these genes showed relevant expression signal corresponding to mean ≥  in all the selected tissues. after normalization and comparative analysis of the expression distributions of the samples from this dataset, we selected a total of samples with two to five biological replicates for the following tissues: adrenal gland, appendices, bone marrow, brain, colon rectum, duodenum, endometrium, esophagus, fallopian tube, adipose tissue, gallbladder, heart, kidney, liver, lung, lymph node, ovary, pancreas, placenta, prostate, rectum, salivary gland, skeletal muscle, skin, small intestine, smooth muscle, spleen, stomach, testis, thyroid, tonsil and urinary bladder. expression profiling and coexpression data analyses rna-seq expression data from all the tissue samples, taken as normalized fpkm from, were log transformed to obtain the final expression values as: log normalized expression distributions of these samples can be seen in additional file 1: figures s and s unsupervised clustering of the samples based in all genes expression was done using agglomerative hierarchical clustering and calculating the distances based on:. this clustering was done for all the samples and just for the tissues using the average expression of the replicates. the coexpression dataset was built calculating the pair-wise spearman correlation coefficient of all the genes along the samples and only selecting, as positive gene-pairs, the ones with a correlation coefficient ≥ . crossvalidation of these correlation values was applied by a random selection of two sample replicates from each tissue and recalculating again the spearman correlation for these random subsets of the data. this sampling was done times, annotating for each gene-pair the number of times that its r coefficient was ≥ . only the gene-pairs validated times in this sampling were selected. in this way, a final set of highly correlated gene-pairs was produced that included genes and coexpression interactions. this coexpression dataset is provided as additional file indicating the names of all the gene-pairs and their correlation value. a gene coexpression network derived from the coexpression data was built using cytoscape and we applied the mcode algorithm to identify clusters inside the network. this algorithm performs an analysis of the topology of the network to find densely connected regions that define modules. the coexpression network built with cytoscape including all the subnetworks found, as well as the parameters derived from the graph analysis, is provided in additional file fig. heatmap presenting the comparison of the transcriptomic profiles of human tissues. unsupervised clustering analysis of the global expression correlation along samples of human normal tissues. the main plot presents the color heatmap and dendrogram produced by the comparison of all the expressed genes along all the samples. the pair-wise distances between samples are calculated using the values of: – spearman correlation coefficient. the color bar with scales shows dark-red corresponding to minimum distances and dark-blue to maximum distances evolutionary analyses for the evolutionary analysis and determination of lowest common ancestor, we used a database of orthologous proteins: orthologous matrix. oma includes a database and resource with methods for the inference of orthologous among complete genomes. we downloaded the oma database into a local mysql database and created python scripts to search for the ensembl id's from our transcriptomic data into this local database and to calculate the lca into each respective ortholog group. in this way we obtained a table with the number of protein-coding genes assigned to each clade in the human taxonomy lineage, as defined by the taxonomy resource from ncbi. furthermore to produce a time-scale, we integrated these data with the timetree of life, that includes the tree of living species calibrated to time. the analysis of the number of genes placed along the evolutionary time-scale allowed visualization of the profile of human genes origin for the whole genome or for specific subsets of genes. in this genes/time profile we calculated the number of protein-coding genes that had lca corresponding to each clade in the taxonomy lineage –that for human includes consecutive levels– and we identified certain levels where major changes occur. these are taken as most significant stages and proposed as key evolutionary hallmarks including specific sets of the human protein-coding genes that are identified. functional enrichment analysis and identification of gene modules for the functional enrichment analysis we used david and geneterm-linker bioinformatic tools with the list of genes from each evolutionary stage level of the human lineage. in all cases the enrichment analyses were done using a hypergeometric test and adjusting the p-values for multiple testing with the benjamini-hochberg procedure. in the same way, we also investigated the functional enrichment of the subnetworks generated by the clusters and modules found in the analysis of the gene coexpression network. general calculations and statistics all the data analyses and graphics have been produced in the r statistic environment. general functions as boxplot, image, qplot or wilcox.test, have been applied over the different data presented. some specific methods or algorithms are cited along different sections of this manuscript.  human global transcriptome profile based on expression in multiple tissues the unsupervised clustering analysis of the global expression correlation along samples of human normal tissues is presented in fig.  this plot shows the color heatmap and dendrogram produced by the comparison of all the expressed genes along all samples. the pair-wise distances were calculated as:; that is a non-parametric approach to calculate similarity and distance. in this way, the heatmap shows a clear relationship between the samples from the same tissue and also presents the proximity between the tissues that have strong biological and physiological links, such as: spleen, lymph nodes and tonsils; or stomach, duodenum, small intestine, colon and rectum. a color bar with scales for the heatmap is included in fig.  to indicate that dark-red corresponds to minimum distance and dark-blue to maximum distance. white color corresponds to medium values and the distributions inside the color bars show that the major part of the comparisons have values around r =  –. the results also reveal that some tissues have transcriptomic profiles very different to the rest, being testis the most different one that produces a clearly separated branch in the dendrogram. house-keeping and tissue-enriched genes the transcriptomic data allowed creation of two relevant subsets of genes: house-keeping and tissue-enriched genes. figure 2a shows the number of genes expressed, with fpkm equal or higher than per number of tissues. this calculation allowed the identification of ubiquitous genes present in at least one replicate of all tissues, and genes expressed in all the samples. the intersection of these genes with a curated dataset of house-keeping genes created by eisenberg et al. gave a total of hkg, indicating a large overlap of 93 %. to identify how relevant was this, we calculated the odds ratio of such overlap and obtained a very significant value of or = . for the tissue-enriched genes we explored the other side of the data in fig. 2a and considered just the genes that were expressed in only one, two or three tissues. we did not take only one, but also two or three tissues, because some tissues are physiologically very related and in fact presented high correlation between them, for example: colon and rectum; small intestine and duodenum, etc. these observations indicated that such tissues share a large number of common genes in their expression profile. finally, a global comparison of the expression distributions of hkg versus teg indicated that the tissue-enriched genes showed significantly lower expression values than the house-keeping genes. to demonstrate this difference, we made two statistical tests, t-test and wilcoxon. in both cases, the p-values were very low and the difference between the mean expressions of teg and hkg was. thus, we can conclude that this difference is not by chance but due to a real different expression regulation that should be considered in further analyses. these boxplots also show that the global variability of teg was larger than the global variability of hkg. this observation indicates that the teg present a larger range of expression values along different tissues and the hkg have a tighter regulation.fig. number of genes expressed along tissues derived from rna-seq transcriptomic data. a plot showing the number of genes expressed, with fpkm equal or higher than per number of tissues. the house-keeping genes are labeled as hkg and tissue-enriched genes as teg. b comparison of the expression distributions of hkg versus teg. c venn diagrams showing the intersection of genes with the dataset of house-keeping genes obtained from eisenberg and levanon human gene hallmarks on the evolutionary time-scale for the evolutionary analysis, we used a phylostratigraphic method for reconstruction of macro-evolutionary trends based on the principle of “founder gene” formation. typically, these methods first identify the homologues of a given gene and then use the divergence between the two most distant to determine the gene age. historically, such studies have been using blast for homology searches. however, this approach was shown to introduce some biases into the analyses. another approach is to use orthologous groups to determine the age of a gene. orthologues are believed to be functionally more similar than paralogues and by definition they trace back to an ancestral gene that was present in a common ancestor of the compared species. the parameters used for clustering orthologous groups affect the age estimations for a gene; for example, restrictive parameters tend to limit the set of possible progenitors. nevertheless, both approaches used for dating the gene origin depend on the correct identification of homologues and/or orthologues, but in the second case the accurate reconstruction of orthologous families imposes a higher stringency and it implies a conservation along the evolutionary clades. for our analysis, we identified the group of orthologues that corresponded to each of the genes detected in the transcriptomic study, mapping them to the corresponding human protein-coding genes in the oma database, and then assigning the lowest common ancestor to each human protein according to its orthologous family. the use of oma in comparison with blast homology approach gives us a more detailed view of gene origin since, as we indicated above, it uses a more restrictive grouping method. thereby, the number of genes dated on ancient clades is lower. once we identified the lca for each human protein-coding gene we assigned such protein/gene to the corresponding taxonomy level in the human lineage as defined in ncbi, that includes taxonomic groups as consecutive clades from the first one, named cellular organisms, to the last one homo sapiens. figure  presents these taxonomic clades placed along the time-scale from the origin to present. for each one of these taxonomy levels we represented the cumulative percentage of proteins that are dated at such level. this is done in the following way: first, plotting all the genes mapped to oma proteins; second, the same plot is produced but including only the proteins that correspond to house-keeping genes; third, plot including only the tissue-enriched genes.fig. evolutionary hallmarks of human protein-coding genes along time-scale. plot presenting the number of human protein-coding genes that are assigned to each of the taxonomic clades. for each one of these taxonomy levels the graph represents the cumulative percentage of protein-coding genes that are dated at such level. in this way, the taxonomic clades are placed as dots along the time-scale from the origin to present. line in black includes all the proteins derived from the mapping of expressed genes in oma. line in blue includes only the hk genes: line in red includes only the te genes: 2157 the analyses of these plots obtained with the phylostratigraphic method revealed the presence of some major differential steps on the emergence of protein-coding genes along the evolutionary time-scale from origin to present. looking at all the expressed coding genes, we can see the global evolutionary profile of the organism, but along this profile we can identify some more prominent steps in the accumulated relative number of genes along time. for example, looking at the plot of the hk a large increase is observed at the start of the curve, from the first taxonomy level to second taxonomy level. by contrast, the te have the major emergence of genes much later. the complete timeline includes phylogenetic clades, but by analysing these points it was possible to identify eight major steps or stage levels that appear on the human gene profile along evolutionary time. moreover, we could assign the number of human protein-coding genes that emerged along each one of these eight stages for the three categories reported: all the expressed protein-coding genes mapped to oma, the hk and the te. the eight stage levels identified are: st1) cellular organisms; st2) cellular organisms to eukaryota; st3) eukaryota to metazoa; st4) metazoa to vertebrata; st5) vertebrata to euteleostomi; st6) euteleostomi to mammalia; st7) mammalia to primates; and st8) primates to homo sapiens. all the numbers corresponding to the human protein-coding genes assigned to each of these eight evolutionary hallmarks are included in fig.  that indicates how many are in each stage either considering the complete human gene set or just the hk or the te. all the information about each one of the human protein-coding genes including the assignment to stages is also provided as additional file fig. evolutionary hallmarks of human protein-coding genes along time-scale. table showing the numbers corresponding to the human proteins found at each of the eight major evolutionary steps depicted from placing them on the time-scale. these eight stages provide the evolutionary hallmarks of the human proteome. the stage levels correspond to: st1) cellular organisms; st2) cellular organisms to eukaryota; st3) eukaryota to metazoa; st4) metazoa to vertebrata; st5) vertebrata to euteleostomi; st6) euteleostomi to mammalia; st7) mammalia to primates; and st8) primates to homo sapiens. the numbers are presented for all the protein-coding genes, for the hk genes and for the te genes the analysis of the hallmarks also reveals that the hk genes are more ancient than te genes. the plots in fig.  and the data in fig.  show clear difference between them. the hk genes present a major increase or expansion in stage, with genes and a change of ≈ 30 % with respect to the total. by contrast, the te genes show a major increase in stage with genes and a change of ≈ 37 %. these observations seem to indicate that house-keeping genes emerged early in evolution and are older in age, knowing that they reflect more essential and constitutive functions. the idea that gene essentiality is associated with older genes has been reported in several studies, for example on yeast and mammalian genes. by contrast, the observations that human tissue-specific genes had emerged later in evolution may reveal that human specific cellular or physiological roles are implemented at molecular level by the appearance of newer genes. gene age data comparison as we indicated above, there are some studies that use the phylostratigraphic method to explore the age of human genes, but most of these studies use sequence similarity search to look for the oldest homologues to the human. to compare the results on human gene age assignment done in this work with other available age assignments, we took the published data from domazet-lošo and from neme, and we represented the information about allocation to lowest common ancestor of the human genes in phylogenetic clades of the evolutionary tree. the assignments were done using common phylostratum to allow the comparison of the data. the results of this comparison are included in additional file 1: figure s and they show a general similarity but some important differences. the most significant difference corresponds to the fact that both domazet-lošo and neme placed a very large number of genes on the first stage of the evolutionary time-scale that goes from the origin of life to first cellular organisms: of and of. this result denotes a bias that, as we indicated above, can be due to the methodology of using the homology search approach. in any case, the idea suggested that one third of the human proteoma may have emerged in evolutionary time before the origin of eukaryotic cells needs deeper studies and it is not what we observed in our analyses. another important difference is that the proposed age mapping allocates the largest number of genes, first, to the chordata-vertebrata-euteleostomi phylostratums and, second, to the mammalia-eutheria. from the evolutionary point of view these results make a lot of sense since the time-scale of life reveals two large expansions of the species precisely around the vertebrates time and around the time of the mammals appearance. these expansions are well reflected in our time-scale profile. finally, it is important to indicate that the age mapping presented in our study only considers human protein-coding genes that are included in orthologous families and in this way it has a lower coverage over human genes than the other reported studies which include more than genes in each case. functional enrichment of the genes at different evolutionary hallmarks we performed functional enrichment analyses of the sets of protein-coding genes included in each one of eight major stages found in the evolutionary study. the full results of these analyses are provided as additional file in all the stages the functional enrichment makes clear biological sense and provides a strong support to the allocation of many biological processes in evolutionary time. in brief, we comment and discuss below some interesting functions enriched in each stage. stage- from the origin of life to first cellular organisms. this stage comprises the genes occurring over two major domains of life: archaea and bacteria. determining the lca, our data shows that human has  % of the protein-coding genes assigned to prokaryotic age. prokaryotes are organisms that lack both membrane-bound organelles and nucleus. functional enrichment analysis showed that this stage involved many basic metabolic processes like glycolysis, the krebs cycle, and lipid oxidation. the enrichment also shows the appearance of the oldest cellular organelle, the mitochondria, and the oldest macromolecular machine, the ribosome, that are well reported to be dated to prokaryotic times. stage- cellular organisms to eukaryota. according to basic literature, the defining feature of eukaryotic cells is that they have membrane-bound organelles, especially the nucleus, which contains the genetic material, and is enclosed by the nuclear envelope. protists, fungi, animals, and plants all consist of eukaryotic cells. eukaryotic cells also contain other membrane-bound organelles such as the golgi apparatus. eukaryotic organisms can be unicellular or multicellular. the functional enrichment analysis for the genes that emerged along this stage showed well the formation of the principal complexes expected in eukaryotes. in this way it is remarkable the enrichment on nuclear pore proteins, nuclear import proteins, nucleosome and chromatin proteins, as well as many proteins involved dna and rna activity: mrna and rrna processing, mrna splicing, dna unwinding, dna polymerase, dna/rna helicase. this stage also marks in time the appearance and biogenesis of some major molecular complexes: the proteasome, the spliceosome, and the ribosome. stage- eukaryota to metazoa. the third stage comprises organisms from opisthokonta and metazoan clades with protein-coding genes. the opisthokonts are a broad group of eukaryotes, including both the animal and fungi kingdoms, sometimes referred to as the “fungi/metazoan group”. this stage comprises metazoan, fungal and protistan taxa, and other multicellular taxa. they also include known fungi and/or parasites of plants like: ascomycota, basidiomycota, chytridiomycetes, glomeromycota, microsporidia, urediniomycetes, ustilaginomycetes and zygomycota. according to our functional enrichment analysis, this stage involves different genes responsible for signal transduction like the gtpases. some of the enriched terms are: pyrophosphatase activity, nucleoside-triphosphatase activity, gtp binding, transferring phosphorus-containing groups. all these functions indicate that it may be the time when phosphorus and phospate acquired a key role in protein function and regulation. other enriched functions, like post-translational protein modification, calcium-binding ef-hand, protein transport and localization also indicate cellular protein regulation. stage- metazoa to vertebrata. this stage includes organisms from eumetazoa, bilateria, deuterostomia, chordata, craniata and vertebrata with genes. the main novelties of this stage are the appearance of protein kinase activity, and the presence of growth factors and some specific signaling proteins like wnt. all biochemically characterized members of the wnt superfamily encode enzymes that transfer organic acids, typically fatty acids, onto hydroxyl groups of membrane-embedded targets. other enriched terms in this stage, like sarcomere and contractile fiber part, may indicate the emergence of the muscular structures present in vertebrates. stage- vertebrata to euteleostomi. this stage represents the largest step in the human lineage according to the number of protein-coding genes assigned, that correspond to a 29 % of the total. the stage comprises organisms from gnathostomata, teleostomi and euteleostomi. the enrichment analysis shows a large functional expansion including new biological systems, like the neural and the vascular-circulatory systems, represented in enriched terms like: neurogenesis, neuron differentiation, axogenesis, voltage-gated channels, neuromuscular junction development, blood vessel development, vasculature development, mesenchymal cell development and differentiation, etc. many other genes are assigned to biological regulation and regulation of cellular processes; including cell death and apoptosis. finally, the appearance of the large family of homeobox proteins seems to be placed at this stage. stage- euteleostomi to mammalia. in this stage there are organisms from sarcopterygii, dipnotetrapodomorpha, tetrapoda, amniota and up to mammalia clades. with genes at this stage, the human lineage achieves 80 % of its gene composition. the most relevant enriched terms are related to the hematologic system, marking the appearance of the leukocytes and the lymphocytes. previous phylogenetic analyses based on gene expression data also placed the date of many proteins from leukocytes around the time of the mammals’ clade. stage- mammalia to primates. this stage comprises clades of theria, eutheria, boreoeutheria, euarchontoglires and primates, representing organisms that give birth to live young without using a shelled egg up to placental mammals. there are genes emerged on this stage, adding up to  % of the cummulative profile in the human gene lineage. a large amount of these genes are enriched in the terms: regulation of gene expression and transcription. other more specific terms are related to the skin or with the sexual reproductive system. this stage also includes a family of cytochrome p proteins and the mammalian defensins: defa1b, defa defa defa defa and defb4a. defensins are a family of antimicrobial peptides and vital contributors to host immune response. being constitutive or inducible expressed genes, they have been shown to contribute to innate host defense via direct bactericidal activity, as well as to adaptive immunity through effector and regulatory functions. stage- primates to homo sapiens. the last stage of human development, with genes, presents a group of quite specific functions played by specific protein families, such as: somatotropin hormone, cytochrome p gtpase activator activity, defense response to fungus and bacterium provided by histatins. his and his have been found only in saliva of humans, macaques and some other primates but not in any other mammals. they are a family of histidine-rich polypeptides that probably function as part of the non-immune host defense system and appeared very late in evolution. cytochromes p constitute a superfamily of proteins that existed in virtually all species from prokaryotes to humans. most of these proteins in the cyp cyp cyp and cyp families encode enzymes involved in the metabolism and elimination of potential toxic compounds like drugs or foreign xenobiotics, and are inducible by various environmental stimuli. this last stage includes a small subset of six cytochromes p that seem to be very specific of the primates-human clades: cyp2a cyp2c cyp2d cyp2j cyp2s and cyp3a the appearance late in evolution of some of these genes may reflect their functional specificity and it is known that they play a key role in human health. network analysis reveals evolutionary age conservation of coexpressed proteins the global coexpression analysis of the human protein-coding genes allowed the construction of a network including highly correlated protein pairs. the integration of these data with the data from the evolutionary analysis –that provided the identification of the eight stages along evolutionary time– did allow mapping the age of the genes on the network according to such stages. these results are presented in fig.  that shows a complex network –like a galaxy– that includes protein nodes associated by interactions. this network corresponds to a subset of the coexpression data, build as indicated in methods, which included proteins and interactions. the subset is done with only the groups that had at least five linked proteins, since we wanted to provide in fig.  a visible representation of the network with a clear color mapping of the eight stages that were identified in the evolutionary study. the colors of the stages are also presented in the illustrated table, fig.  to allow better identification of the number and % of proteins at each age hallmark. the analysis of the network done with the algorithm mcode revealed the existence of major subnetworks, which can be considered as major constellations in the galaxy of relational nodes. the color code of this large graph indicates that there is enrichment in certain colors in each subnetwork. to prove this in a more accurate way, we built a graphic representation for each one of the subnetworks found to show the proportion of proteins assigned to each of the eight evolutionary stages with their corresponding color code. this graphic is presented in additional file 1: figure s that shows each subnetwork with its specific color pattern, indicating that there are always some predominant colors: subnetwork is the oldest with red predominant colors and subnetwork is the newest with blue predominant colors. as a conclusion, these results revealed that in the groups of highly coexpressed proteins there is a tendency to include proteins of the same evolutionary age.fig. human coexpression network mapping the evolutionary age on highly correlated nodes. representation of the coexpression complex network –like a galaxy– that includes protein nodes related with interactions. this network corresponds to a subset of the larger coexpression network produced, which included proteins and interactions. the subset is done to include only the groups that had at least five linked proteins. the color mapping of the nodes correspond to the eight stages that were identified in the evolutionary study. the network also includes numbers for major subnetworks –clusters of closely related proteins that include more than nodes– considered as major constellations in the galaxy of nodes. three panels on the right show an enlarged view of three subnetworks corresponding to: ribosomal proteins, mitochondrial proteins and angiogenesis proteins finally, we did a functional enrichment analysis of the proteins forming the subnetworks which again showed a coherent biological enrichment in specific functions: immune response; cell cycle; cytoskeleton; rna splicing; ribosome; extracellular matrix; muscle and contraction; gametes and reproductive process; cell junction and cell adhesion; mitochondria and atp synthesis; angiogenesis and vasculogenesis. more detailed results for this analysis are included in additional file 1: figure s thus, we observed that age-related proteins are predisposed to present expression coregulation and to have close functional links. combining age data, functional data and coexpression data can provide a deeper view about the links and roles of the human protein-coding genes. in this way, we observed for example that subnetwork presented as expected an overwhelming majority of ancient genes from the prokaryotic or eukaryotic age. on the other hand, subnetworks, and showed a higher proportion of recent genes dated after vertebrata. these results agree with studies based on yeast protein physical interaction networks, arguing that proteins preferentially interact with proteins of same age and origin. moreover, it was previously shown that coexpression networks can be conserved over the evolutionary history, and these genes tend to be functionally related and provide selective advantages. it has been also reported that coexpression networks are found associated to functions like cell adhesion, cell cycle, dna replication and dna repair, and this is in agreement with functions found enriched in subnetworks of our analyses: subnetwork and.  the transcriptomic study of the human gene expression distributions and profiles along tissues provided a global mapping of the activity of most human genes and of the links between them, showing the expected association of samples from common physiological regions, i.e.: the gastrointestinal tract, the hematopoietic and lymphatic system and the muscle. the evolutionary study of the human protein-coding genes placed them in the time-scale of the living species and revealed eight distinct hallmarks along such time-scale, showing that the hk genes are more ancient than the te genes. the hk genes present the major emergence in stage of the evolutionary profile, while the te genes have the major emergence in stage the functional enrichment study found coherent groups of terms and annotations assigned to the genes placed at each evolutionary stage. for example, in stage there were many functional terms on essential metabolic processes, like aerobic respiration and mitochondrial activity; and in stage there were enriched functions related to the nucleus and genome regulation, like chromatin and nucleosome assembly, dna replication, mrna processing. finally, the study of the pair-wise correlation of the gene expression profiles along tissues allowed building human gene coexpression networks and find modules with functional and biological meaning. the mapping of the age of the protein-coding genes on these networks demonstrated the existence of tight links between age-related proteins. additional files additional file 1: additional pdf file that includes six supplementary figures. additional file 2: table including the values of the gene pair-wise spearman correlation and the cross-validation from the coexpression analysis, as well as the gene ids to allow the reconstruction of a human coexpression network with proteins and interactions. additional file 3: cytoscape file that includes the complete coexpression network that is produced in this work and presented in fig.  of the main article. this is a file to be opened only with cytoscape program. additional file 4: table including the data and ids of the human protein-coding genes included in each of the eight evolutionary stages. additional file 5: table including the results of the functional enrichment analysis of the human protein-coding genes included in each of the eight evolutionary stages. the results for each of the eight stages are presented in different spreadsheets. 
understanding the relational landscape of the human protein-coding genes is essential for interpreting the functional elements and modules of our active genome. moreover, decoding the evolutionary history of the human genes can provide very valuable information to reveal or uncover their origin and function. however, until recently, little has been reported about the evolutionary marks present in human protein-coding genes, particularly from the combined perspective of gene expression and protein evolution. human protein evolutionhuman gene evolutiontranscriptomicsrna-seqtissue transcriptomicsprotein familiesgene coexpressiongene house-keepinggene tissue-enriched6th solbio international conference soibio 2016riviera maya, mexico 22- april http://www.soibio.org/en/icmexico2016issue-copyright-statement© the author 2016 the evolutionary analysis of the human proteins is performed combining the information from: a database of orthologous proteins, the taxonomy mapping of genes to lineage clades and the evolution time-scale mapping provided by timetree. the development of large-scale technologies for quantitative transcriptomics has enabled comprehensive analysis of the gene expression profiles in complete genomes. rna-seq allows the measurement of gene expression levels in a manner far more precise and global than previous methods. in this respect, multiple efforts have been done to determine and analyse the gene expression patterns of human cell types in different conditions, either in normal or pathological states. we present a combined analysis of human protein-coding gene expression profiling and time-scale ancestry mapping, that places the genes in taxonomy clades and reveals eight evolutionary major steps, that include clusters of functionally coherent proteins.
 the development of large-scale technologies for quantitative transcriptomics has enabled comprehensive analysis of the gene expression profiles in complete genomes. rna-seq allows the measurement of gene expression levels in a manner far more precise and global than previous methods. studies using this technology are altering our view about the extent and complexity of the eukaryotic transcriptomes. in this respect, multiple efforts have been done to determine and analyse the gene expression patterns of human cell types in different conditions, either in normal or pathological states. however, until recently, little has been reported about the evolutionary marks present in human protein-coding genes, particularly from the combined perspective of gene expression and protein evolution.  we present a combined analysis of human protein-coding gene expression profiling and time-scale ancestry mapping, that places the genes in taxonomy clades and reveals eight evolutionary major steps, that include clusters of functionally coherent proteins. the human expressed genes are analysed using a rna-seq dataset of samples from tissues. the evolutionary analysis of the human proteins is performed combining the information from: a database of orthologous proteins, the taxonomy mapping of genes to lineage clades and the evolution time-scale mapping provided by timetree. the human protein-coding genes are also placed in a relational context based in the construction of a robust gene coexpression network, that reveals tighter links between age-related protein-coding genes and finds functionally coherent gene modules.  understanding the relational landscape of the human protein-coding genes is essential for interpreting the functional elements and modules of our active genome. moreover, decoding the evolutionary history of the human genes can provide very valuable information to reveal or uncover their origin and function. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords human protein evolutionhuman gene evolutiontranscriptomicsrna-seqtissue transcriptomicsprotein familiesgene coexpressiongene house-keepinggene tissue-enriched6th solbio international conference soibio 2016riviera maya, mexico 22- april http://www.soibio.org/en/icmexico2016issue-copyright-statement© the author 2016 