Introduction
============

Comparative genomics is a powerful tool that allows researchers to infer the genetic changes underlying evolutionary transitions: Correctly inferring the phylogenetic point of origin of the genes involved in those transitions is crucial for that purpose. We illustrate this idea here by considering two important and related biological processes in mammals: placental formation during pregnancy and the exploitation of host genes in cancer development ([@evaa026-B45]), showing how a consideration of the phylogenetic history of the genes involved in them is enlightening. In the first example, the evolutionary origin of mammalian pregnancy is an intriguing process involving reprograming of ancient gene regulatory networks ([@evaa026-B52]), as well as innovations such as the emergence of new cell types ([@evaa026-B81]; [@evaa026-B8]) and the modulation of the immune system to allow for the internal maintenance of an organism genetically distinct from the mother ([@evaa026-B56]). Identifying the group of novel genes that originated at the radiation of the placental mammals would shed light on how reproductive strategies were changed in those mammals by genetic pathway remodeling ([@evaa026-B21]).

Viviparity (internal embryo development and live birth) has evolved multiple times in vertebrates, including in bony and cartilaginous fish, squamates, and mammals ([@evaa026-B77]). However, the degree of maternal--fetal tissue integration varies significantly across these groups ([@evaa026-B68]). The extant mammals contain two subclasses: prototheria, which contains the order monotremata (e.g., platypus), and theria, which contains the metatheria (marsupial) and eutheria (placental) clades ([@evaa026-B70]). Eutherian mammals diverged from marsupials about 160 Ma ([@evaa026-B54]; [@evaa026-B19]; [@evaa026-B64]; [@evaa026-B74]; [@evaa026-B89]). Among therian mammals, marsupials and eutherians are both viviparous, and the extraembryonic membrane in marsupials is homologous to the placenta in eutherian mammals ([@evaa026-B86]). However, the two groups have evolved differing reproductive strategies: marsupials change their milk composition dramatically during the extended lactation window in compensation for their more primitive, underdeveloped, placentae (and hence offspring at birth), whereas eutherians acquired more complicated placentae, allowing the offspring to be more developed at birth ([@evaa026-B48]; [@evaa026-B28]).

Mammalian placentae are formed by close apposition of maternal and fetal tissues ([@evaa026-B57]) and have the vital functions of exchanging nutrients, dissolved gases, and waste during fetal development. Placental morphology and the invasiveness of fetal tissue (how much it erodes the maternal uterine lining) varies greatly across eutherian mammals ([@evaa026-B70]). The least invasive are the epitheliochorial placentae (e.g., in pigs, cows, and horses), whereas the endotheliochorial placentae (e.g., in cats, bats, and elephants) invade by erosion of the trophectoderm in uterine lining. The most invasive form are the hemochorial placentae (e.g., in mice and humans), with deeper fusion of the fetal tissue into the maternal tissue ([@evaa026-B58]; [@evaa026-B87]). Notably, the invasive hemochorial placenta has been reconstructed to be the ancestral state of the eutherian placenta, based on its phylogenetic distribution ([@evaa026-B79]; [@evaa026-B87]).

The evolution of mammalian pregnancy at the genomic level was complex. Rather than a burst of gene evolution at the origin of the placental mammals, a large set of genes contributing to placentation already existed, but their functions shifted or expanded during the origin of eutherian mammals ([@evaa026-B28]), probably through a large-scale rewiring of gene regulatory networks that was at least partly mediated by the insertion and repurposing of transposable elements (TEs) and retroviruses ([@evaa026-B41]; [@evaa026-B10]; [@evaa026-B51]). This hypothesis is supported by the observation that many of the genes expressed early in mouse placental development are evolutionarily much older than the mammalian radiation, whereas genes expressed later showed a higher propensity to be rodent-specific ([@evaa026-B38]). The ancient retroviruses, relics of whose genomes were exploited during the evolutionary origin of placental mammals, presumably also predate the eutherian divergence ([@evaa026-B47]). As a result, many genes and sequences that were ancestrally expressed in other organs were recruited to be expressed in the endometrium and gained new functions related to maternal--fetal interactions ([@evaa026-B52]). A novel type of cell, decidual stromal cells, also appeared along the stem branch of eutherian mammals ([@evaa026-B81]). Among the functions of these new cells is thought to be a role in preventing immune rejection and the inflammation caused by maternal--fetal conflict ([@evaa026-B8]). It even appears that part of this immune response has been repurposed to facilitate placenta formation: during trophoblast invasion, an innate immune response by natural killer cells is needed for proper placentation ([@evaa026-B24], [@evaa026-B25]).

The evolutionary origins of placentation are of interest for many reasons, including that of what it can tell us about the functions of the immune system and the nature of certain diseases. For instance, the invasion of trophoblast tissue into the maternal tissue incurs an innate immune response, whereas eutherian mammals have evolved mechanisms to repress maternal--fetal immune conflict to enable viviparity ([@evaa026-B24], [@evaa026-B25]). Similarly, our other example of using the evolutionary age of genes for biological exploration comes from the analogy between placentation and metastasis ([@evaa026-B45]). Because these two processes have similarities at the cellular level, placentation genes may also be overrepresented in genes that are apt to change expression due to cancer. Perhaps strikingly, malignant tumors are suppressed in mammalian lineages that have evolved noninvasive placentae ([@evaa026-B20]; [@evaa026-B72]). It is also known that phylogenetically older genes tend to be associated with human diseases more generally ([@evaa026-B18]) and that one pathway for cancer development can be mutations that disrupt the communication between new and old genes ([@evaa026-B76]). We have therefore identified a list of genes that are involved in both placenta invasiveness and in cancer metastasis and inferred their evolutionary ages on a rather finer timescale than the previous analyses ([@evaa026-B26]; [@evaa026-B31]; [@evaa026-B20]).

A reliable method for estimating gene age is to infer orthologous genes across species on a known phylogeny, that is, to find homologous genes that last shared a common ancestor at a speciation event ([@evaa026-B39]). Orthology relations are often used to infer evolutionary histories ([@evaa026-B42]) and to provide functional annotation for uncharacterized homologs ([@evaa026-B44]), as such genes tend to conserve their molecular and biological functions ([@evaa026-B40]). Many common orthology detection methods are based on sequence similarity. For example, OrthoDB ([@evaa026-B44]) first identifies a set of putative orthologs on the basis of reciprocal best hits ([@evaa026-B83]) which can rapidly identify pairs of orthologs between two species. It then uses hierarchical clustering of orthologs to refine the search ([@evaa026-B43]). We have previously developed a pipeline that adds the use of gene order information (e.g., synteny) to data on sequence similarity in order to infer orthology relationships between pairs of genomes ([@evaa026-B11]; [@evaa026-B2], [@evaa026-B3]). When these pairwise estimates are placed on a phylogeny, they allow us to estimate relative gene ages. Because this use of synteny distinguishes our pipeline from most existing orthology inference tools ([@evaa026-B96]) and pregnancy-related gene prediction tools ([@evaa026-B36]) and because it resolves more multigene families into one-to-one orthologous relations than these methods, we used a rewritten version of it to explore two questions: 1) Understanding the genomic changes underlying the genetic origins of mammalian pregnancy and tracing how the gene "toolkit" for placentation evolved and 2) understanding how placental-related genes and mammal orthologs more generally can be used to infer whether a gene has a propensity to be modified or exploited in cancer development.

Our results show that most genes expressed in placenta have evolutionary origins prior to the eutherian radiation, suggesting that they were repurposed or gained new functions in this evolutionary transition. We found that orthologs shared by eutherian ancestor but not any earlier ancestor are enriched in functions and protein classes such as transcription regulation by Krüppel-associated box (KRAB)-zinc-finger proteins (ZNFs), natural killer cell activity, keratinization, and the Melanoma Antigen Gene (MAGE) protein family, indicating that there were large-scale changes in transcriptional regulation and immune response at the origin of mammalian pregnancy. As expected from studies with coarser timescales, genes implicated in tumor development tend to be evolutionarily older, even within the relatively restricted timeframe of the mammalian radiation.

Materials and Methods
=====================

Genomic Data
------------

The complete set of coding sequences from the genomes of 44 amniotes (3 birds, 1 monotreme, 3 marsupials, and 37 eutherian mammals) was downloaded from the Ensembl database (release 84, [@evaa026-B96]). From these genomes, we extracted the longest transcript of each protein-coding gene.

Orthology Assignment Pipeline
-----------------------------

Using these data, we sequentially compared each other genome with the human genome using our pipeline for orthology inference: [OR]{.ul}thology [I]{.ul}nference using [S]{.ul}ynteny (ORIS; [@evaa026-B11]; [@evaa026-B2], [@evaa026-B3]). Since our previous analyses, we have improved ORIS both by dramatically expanding the set of genomes searched (44 versus 18) and by improving the sensitivity of the homology search with a multipass approach using GenomeHistory ([@evaa026-B12]; [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). We also included three bird genomes in our analysis to allow better basal representation at the root branch of the phylogenetic tree that we used for tracing the origins of each orthologous group. By varying the sequence similarity required for a match, as well as the seed parameters used when deciding to spawn alignments, we have improved the pipeline's performance when used with genomes more distantly related to humans. The parameter values used for the searches are given as [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. The nonsynonymous (*K*~a~) and synonymous (*K*~s~) divergences of the homologous pairs were estimated using maximum likelihood. With these divergence estimates, we first collapsed tandem duplicates in each genome. Then, for a given pair of genomes, the orthology inference procedure began by defining one-to-one pairs of homologs below a set *K*~a~ threshold as the initial orthologs ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). If either of the nearest syntenic neighbors to these anchors are also homologs, then these neighboring pairs are assigned to be orthologs. This process of identifying one-to-one homologs and extending the orthology relationships to their neighboring homologs continues until no further gene pairs meeting the criteria exist.

Here, we modified the pipeline by adding a second pass in the orthology detection process. Once the initial ortholog detection was completed, we went back and examined the ambiguous cases where one gene has multiple hits in the other genome, but where only one of those homologs gives a *K*~a~ distance to the human gene below a threshold ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). We then added these pairs as anchors, again running the search process until no further orthologs are found. We also modified how the algorithm handles those remaining cases where there are clusters of genes in each genome with similarity to genes in the other and where synteny cannot resolve the resulting ambiguity. We used a *K*~a~ threshold such that when one pair of genes from the two genomes is below threshold and no other pairs of genes are below two times that threshold, we defined this closest pair to be orthologs. In pairwise comparisons against human genome, the implementation of ORIS with these changes yielded more orthologs, sometimes many more, than our earlier approach, both for distant mammal species such as platypus and for species that are more closely related to human such as the chimpanzee. Our orthology inferences are provided as a [supplementary data](#sup1){ref-type="supplementary-material"} set, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, which comprises of 554,774 pairs of orthologs among these 44 genomes.

Ancestral State Inference
-------------------------

We considered the presence/absence of an ortholog of each human gene in each other taxon as a character and constructed the corresponding character matrix. In this matrix, the columns are the 19,826 human genes (after the merging of tandem duplicates), and each row represents the presence or absence of an ortholog of that gene in another aminote genome (43 in total). [Supplementary figure S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, gives a schematic of our ancestral state reconstruction inferences. To make these estimates, we inferred the most recent common ancestors (MRCAs) of orthologs by placing the origin of that gene at the most recent internal node in the tree such that every species possessing the gene is a descendant of that node. This approach corresponds to making the assumption that each ortholog can appear only once in the tree, an assumption we are comfortable with due to the robustness of our synteny-based approach for orthology inference.

Branch Length and Branch-Specific Selective Constraint Estimation
-----------------------------------------------------------------

Of the 19,826 human genes analyzed, 948 had orthologs in all 43 other genomes. We obtained the longest transcripts of these 948 orthologous genes from every genome and then translated them into protein sequences. Multiple sequence alignments of the protein sequences were performed using T-coffee with default settings ([@evaa026-B61]). We then made codon-preserving nucleotide alignments from the protein alignments. After removing gaps from those alignments, we concatenated them into one meta-alignment, which had 44 taxa and 29,293 codons. A branch-specific *ω* (d*N*/d*S*) was estimated by running two rounds of codeml analyses in PAML 4.9 ([@evaa026-B91]) with CodonFreq F3X4. In the first run, we estimated one universal *ω* value for the entire tree and obtained the estimated number of synonymous substitutions (d*S*) on each branch. Then, in the second round, we set the initial branch lengths in the guide tree to be the d*S* estimates from the first step and allowed *ω* to differ on each branch. Phylogenetic trees were visualized using the R packages ape ([@evaa026-B65]), phytools ([@evaa026-B67]), ggtree ([@evaa026-B93], [@evaa026-B92]), and ggplot2 ([@evaa026-B85]) in R 3.6.1 ([@evaa026-B66]).

The estimated d*S* for each branch was converted to a time inference using the published divergence time estimate between birds and mammals of 312 Ma, as retrieved from timetree.org ([@evaa026-B30]; [@evaa026-B46]). For branches leading to humans (A to N in [fig. 1](#evaa026-F1){ref-type="fig"}), the average rate of new orthologs appearing was estimated as the sum of number of orthologs appearing on the branch divided by the length of that branch (in millions of years). We can convert these values to an inferred rate of new ortholog occurrence on each branch, such that the number of new orthologs on each branch is given by the inferred rate multiplied by the branch length (in millions of years). We first calculated an average rate of ortholog appearance, which was given by the sum of number of orthologs on all 14 branches over the sum of their branch lengths ([fig. 2*a*](#evaa026-F2){ref-type="fig"}). This average rate of new ortholog appearance was 20.32 genes per million years, a value comparable with previous estimates ([@evaa026-B50]). Because the rate of new gene appearance is known to be variable ([@evaa026-B50]), we then estimated branch-specific rates, which are given by the number of orthologs appearing on a branch divided by its length ([fig. 2*b*](#evaa026-F2){ref-type="fig"}). We also compared our gene age inference results with those from a previous study ([@evaa026-B97]). Because of the difference in taxon sampling between the two studies, we estimated one rate for branches F, G, and H combined, one rate for I and J combined, one rate for L and M combined, and separate rates for the remainder of the individual branches, which allows for a comparison to the data of [@evaa026-B97].

![---The appearance of mammalian orthologous genes in a phylogenetic context. (*a*) Shown is a vertebrate orthology tree with the number of synonymous substitutions per synonymous site (d*S*) as branch lengths. Black slices in the pie charts show the proportions of genes in the nonhuman genomes that have an orthologous human gene. We inferred the first appearance of each ortholog to give monophyletic groups that possess that gene (Materials and Methods): The inferred ancestral ortholog percentages are thus shown on the internal nodes. d*S* values were estimated using codeml with an unrooted guide tree and alignments of 948 orthologs (see Materials and Methods). The topology used was adapted from [@evaa026-B54]. The letters A--N label the internal branches leading to humans. The box region is expanded to show the primate lineage. The species images were downloaded from PhyloPic.org. (*b*) Bar chart of the number of orthologs on the internal branches leading to human (A--N in *a*). The number of such new orthologs with identified homologs (but not orthologs) in birds is shown in red (see Materials and Methods).](evaa026f1){#evaa026-F1}

![---The number of the orthologs on internal branches leading to human and branch-specific selection estimates from the alignments of 948 orthologs. (*a*) Scatterplot of number of the orthologs versus the estimated branch length in million years (see Materials and Methods). The dashed line represents the average rate of new ortholog occurrence across the 14 branches. (*b*) Comparison of two estimates of the relative rate of new gene appearance on different phylogenetic branches in the mammalian tree. The blue dashed line shows the average rate of occurrence on each branch given by dividing the number of orthologs by the branch length. The yellow line is the rates estimated using new gene counts from [@evaa026-B97] divided by the same branch lengths. (*c*) Branch-specific d*N*/d*S* over time: the *x* axis gives time (in millions of year ago) estimated from cumulative d*S* and the *y* axis gives the estimated average d*N*/d*S* for the corresponding branch (see Materials and Methods).](evaa026f2){#evaa026-F2}

Functional Analysis of Orthologs
--------------------------------

Human orthologs were grouped into five bins based on their inferred first appearance on the phylogenetic tree ([fig. 2](#evaa026-F2){ref-type="fig"}). We will refer to the branch leading to the MRCA of birds and mammals as the "root branch (R)," the branch to mammalia as the "mammal branch (A)," the branch to the theria subclass (marsupials and eutherians) as the "therian branch (B)" and the branch to the stem group of eutherian mammals (i.e., the common ancestor of all placental mammals) as the "eutherian branch (C)." We merged orthologs that appeared on later branches within the eutherian mammal clade to make up a new "post-eutherian (D-N)" branch and its associated new genes. Functional annotations, including molecular functions, biological processes, pathways, protein domains, and tissue expression for each gene list, were obtained from the DAVID Bioinformatics Resources version 6.8 ([@evaa026-B33]). The enrichment threshold EASE score was set to be 0.1 ([@evaa026-B32]).

Gene Ontology (GO) analyses were also performed to identify overrepresented GO terms associated with each group of orthologs using Fisher's exact tests in PANTHER release 14.1 ([@evaa026-B55]). In [supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, we list GO terms that are overabundant among genes appearing on each internal branch of [figure 1](#evaa026-F1){ref-type="fig"} relative to genes present on the root branch (False Discovery Rate-corrected significance level of *P *\<* *0.05 \[[@evaa026-B4]\]). We then obtained the functional annotations for each gene list using the PANTHER Gene List Analysis option. For each branch, we counted the number of genes with molecular function GO annotations that contain any placenta-related terms, such as placental/placenta development (GO:0001890) and embryonic placenta morphogenesis (GO:0060669). Similarly, we counted the number of genes which encode "transposable-element-derived," "retrotransposon-like," or "transposase" proteins for each branch. The percentages of transposable-element-like genes and placenta-related genes relative to the total number of genes appearing on each branch are given in [supplementary figure S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online.

We also compared the semantic similarity scores of the molecular function ontology annotations for genes from all possible pairs of our five age bins (i.e., the root, mammal, therian, eutherian, and postplacental branches). For this analysis, we used the GO annotation file from the R package org.Hs.eg.db ([@evaa026-B7]). The GO semantic similarity itself was calculated using the package GOSemSim release 3.10 ([@evaa026-B94]) in R 3.6.1, with the graph-based Wang's method ([@evaa026-B82]) and the Best-Match Average method for combining the similarity scores for GO terms in each cluster ([@evaa026-B1]).

In order to find potential genes that are related to both tissue invasion and cancer, we uploaded each gene list to the Cancer Genome Atlas (TCGA) data portal and obtained the number of cancer census genes within each ortholog group ([@evaa026-B6]; [@evaa026-B71]). Cancer hallmarks for these cancer census genes were then downloaded from the Catalogue of Somatic Mutations in Cancer (COSMIC) database release v88 ([@evaa026-B75]). From the TCGA data portal, we also downloaded, for each gene in our data set, the number of simple somatic mutations known for it. Within each age group, the number of mutations per gene was normalized by the length of the gene's longest transcript, and the median and average mutation counts were calculated for each branch ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The median mutation count of each branch was compared with that of the previous branch using the Wilcoxon rank-sum test ([@evaa026-B53]). The mutation counts were visualized using the R package ggridges 0.5.1 ([@evaa026-B88]).

Results
=======

Orthology Assignment and Gene Age Inference
-------------------------------------------

We performed orthology inference for 43 vertebrate species' genomes by comparison to the human genome (Materials and Methods). [Supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, shows the parameters used for orthology inference and the number of orthologs identified in each pair of species. We used slightly relaxed homology search parameters for orthology estimates in species more deeply diverged from humans. The pie charts at the tip nodes in [figure 1](#evaa026-F1){ref-type="fig"} show the percentage of genes with human orthologs in each genome. We next inferred the latest point of origin for every human gene in the phylogeny of [figure 1](#evaa026-F1){ref-type="fig"}: in other words, the internal node of the tree that partitions all of the genomes possessing an ortholog of that gene into a monophyletic group. The number of orthologs so appearing on each internal branch is also illustrated in [figure 1](#evaa026-F1){ref-type="fig"}. The percentage of genes with human orthologs under our definition is thus 60.9% in chicken (*Gallus gallus*), 54.1% in platypus (*Ornithorhynchus anatinus*), 72.9% in opossum (*Monodelphis domestica*), 80.8% in mouse (*Mus musculus*), and 85.9% in gorilla (*Gorilla gorilla*).

Next, we grouped orthologs into bins based on their phylogenetic branch of first appearance: In other words, each ortholog was inferred to have appeared on one of the internal branches leading to human. We placed 13,487 human orthologs on the root branch, the MRCA of birds and mammals, 1,407 orthologs on the mammal branch (branch A in [fig. 1](#evaa026-F1){ref-type="fig"}), 2,005 on the therian branch, 1,334 on the eutherian branch, and 1,593 on branches within the placental mammals. It is natural to ask where these genes that appeared after the split with birds came from, and in particular if they are novel duplications of genes with ancient homologs in birds. We thus examined the subgroup of human genes that have homologs in the three bird genomes but were not assigned bird orthologs by our pipeline. We computed the proportion of all genes appearing on each later branch of the tree that fell into this set of genes with bird homologs ([fig. 1*b*](#evaa026-F1){ref-type="fig"}). These comparisons allow us to begin to assess the relative importance of three explanations for the appearance of an orthology relationship on a particular post-bird branch: that it results from 1) a new duplication of a gene with existing homologs in birds, 2) an ancient gene whose orthologs in birds were not identified due to high sequence divergence, domain rearrangement, or annotation artifacts, or 3) a gene that lacks homologs in birds because it was produced by *de novo* processes, TE-derived open reading frames, or horizontal gene transfer ([@evaa026-B37]; [@evaa026-B73]). Based on annotated genomes alone, we cannot distinguish explanations \#2 and \#3, but as [figure 1*b*](#evaa026-F1){ref-type="fig"} suggests, the genes appearing on the oldest post-bird branches show a relatively large fraction (∼74.3%) of genes with homologs in birds, suggesting that many of these genes are duplications of existing vertebrate genes, especially because the necessity of using stringent homology cutoffs leads us to believe that we have underestimated the proportion of human genes with homologs in birds. As we move to more recent branches of the tree, relatively few genes that appear have bird homologs, which is expected, because for a gene to appear on these latter branches, its ortholog must be absent not merely in birds but also in, for instance, monotremes, marsupials, and early diverging placental mammals such as elephants. Hence, it is reasonable to attribute these (relatively rare) newer genes to the other processes described in \#3 above.

At a higher level, these results confirm that our pipeline generally meets our design goal of generating ortholog sets with very few false positives (pairs of genes identified as orthologs that actually diverged through other events such as gene duplications) at the expense of some lack of sensitivity to detecting the most diverged ortholog pairs in distantly related genomes. The inclusion of the phylogenetic ancestral reconstruction further reinforces this conclusion, because the absence of a gene from a single genome will generally be "corrected for" by the presence of other genomes equally related to humans that possess that gene. For example, the three marsupial genomes (opossum, tammar wallaby, and Tasmanian devil) lack orthologs to 27.1%, 38.6%, and 31.6% of human genes, respectively. However, they collectively lack orthologs to only 18.0% of human genes ([fig. 1*a*](#evaa026-F1){ref-type="fig"}).

The Number of Orthologs Appearing on the Eutherian Branch Is Not Unusual
------------------------------------------------------------------------

The branch lengths in [figure 1*a*](#evaa026-F1){ref-type="fig"} are the mean number of synonymous substitutions per synonymous site (d*S*) across a set of 948 orthologs that are identified across all 44 species (Materials and Methods). These estimates are consistent with known species divergence time estimates ([@evaa026-B30]). We estimated branch-specific selective constraints (d*N*/d*S* or *ω*) using codeml ([@evaa026-B91]), finding that the d*N*/d*S* estimates for all of the internal branches indicate that these orthologs experience (on average) purifying selection ([fig. 2](#evaa026-F2){ref-type="fig"}).

We further asked if the therian or eutherian branches showed an unusually high number of new orthologs appearing along them, given their estimated length in terms of d*S*. We therefore compared the timespan of each branch (extrapolated from the d*S* estimates using the mammal--bird divergence time of 320 Ma as the single calibration point) with the number of orthologs appearing along that branch ([fig. 2*a*](#evaa026-F2){ref-type="fig"}). As shown in [figure 2*b*](#evaa026-F2){ref-type="fig"}, we also compared our relative estimates of branch-specific gene appearances with those from [@evaa026-B97], controlling for the differing total numbers of genes analyzed. Notably, we do not observe an unexpectedly large number of orthologs appearing along the eutherian branch. However, both our results and those of [@evaa026-B97] did show two trends: 1) the branches immediately before and immediately after the common placental mammal branch *do* show elevated rates of gene appearance and 2) a relatively large number of apparent new genes are seen along the human-specific branch, an observation that we would tend to attribute to the extensive annotation work done on the human genome rather than to any actual evolutionary trend ([fig. 2](#evaa026-F2){ref-type="fig"}). We do not see any unusual patterns of selective constraint on the common mammal branch, the therian branch, or on the eutherian branch ([fig. 2*c*](#evaa026-F2){ref-type="fig"}), indicating that the evolution of mammalian pregnancy was not associated with a genome-wide "burst" of gene duplication or a relaxation of selective constraints.

The Majority of Human Placenta-Related Orthologs Appeared before the Origin of Placental Mammals
------------------------------------------------------------------------------------------------

We named the five groups by the internal branch they were placed on, that is, orthologs on the root branch (R), mammal branch (A), therian branch (B), eutherian branch (C), and branches within the eutherian lineage (D to N). We obtained functional annotations for all human genes appearing on each branch from the PANTHER database (release 14.1).

We first compared the percentage of genes that are associated with placenta-related biological process GO terms on each branch ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). We found that most of the placenta-related orthologs appeared on the root branch (72 genes out of 91), indicating that ancient genes were likely being repurposed to contribute to mammalian pregnancy after the origin of mammals ([@evaa026-B72]). The proportion of all genes that are annotated as placenta related is similar for the root branch and the eutherian branch ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), again suggesting no particular process of genic innovation coincident with the evolution of placentation. Notably, only a handful of placenta-related genes were inferred to have appeared *after* the eutherian branch, reinforcing our confidence in our gene age estimates.

We also calculated the percentage of genes that encode retrotransposon-like, transposable-element-like, or transposable-element-derived proteins on each branch ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Interestingly, the orthologs appearing on the eutherian branch encompass most of the TE-related genes so identified. These retrotransposon-like elements could be survivors of ancient host--viral interactions that were retained in eutherian genomes. It is possible that these TE-like protein-coding genes also contributed to the origin of placental mammals through neofunctionalization ([@evaa026-B5]).

Orthologs Appearing on the Eutherian Branch Are Enriched for Functions in Transcriptional Regulation and Immune Responses
-------------------------------------------------------------------------------------------------------------------------

We next obtained functional annotations for orthologs that appeared on the internal branches of interest (A, B, C, and D--N). Our analysis used five functional annotation sets in DAVID: *UP_Keywords*, which provides functional summaries, *GOTERM_MF_DIRECT* (molecular functions), *GOTERM_BP_DIRECT* (biological processes) for GO information, *KEGG* for biochemical pathways, and *Interpro* for protein domains ([@evaa026-B33]). Functionally related annotation terms were clustered together, and an enrichment score was calculated for each cluster ([@evaa026-B34]). The top ten enriched functional clusters for each ortholog group are shown in [table 1](#evaa026-T1){ref-type="table"}.

###### 

Summary of Top Ten Functional Annotation Clusters for Each Group of Orthologs

                                                                                        Top 10 Annotation Clusters                                             Functional Description                     DAVID Enrichment Score
  ----------------------------------------- ------------------------------------------------------------------------------------------------------------------ ------------------------------------------ ------------------------
  Orthologs on the mammal branch                                                                    1                                                          SPRY domain, B-box ZNF, RING-type ZNF      5.45
  2                                                                                       Innate immune response                                               3.6                                        
  3                                                                          Cell adhesion molecule binding, cell recognition                                  3.01                                       
  4                                                                           Adaptive immune response, cytokine, interferon                                   2.95                                       
  5                                                                             Interleukin, interleukin receptor binding                                      2.67                                       
  6                                                                                           EF-hand domain                                                   2.52                                       
  7                                                                     Chemical carcinogenesis, drug metabolism---cytochrome P450                             2.41                                       
  8                                                                        Steroid metabolic process, sulfotransferase activity                                2.07                                       
  9                                                                            Phototransduction guanylate cyclase activity                                    2.02                                       
  10                                                                         Biomineralization, biomineral tissue development                                  1.93                                       
  Orthologs on the therian branch                                                                   1                                                          Mammalian taste receptor activity          10.15
  2                                                                                    Glycoprotein, disulfide bond                                            9.41                                       
  3                                                                         Olfactory and sensory transduction, GPCR activity                                  8.19                                       
  4                                                                                 Glycoprotein, transmembrane helix                                          7.1                                        
  5                                                                                Keratinization, peptide crosslinking                                        4.01                                       
  6                                                                              KRAB, C2H2-ZNF, transcription regulation                                      3.4                                        
  7                                                                                       Peptidase S1 activity                                                3.25                                       
  8                                                                                C-type lectin, carbohydrate binding                                         2.56                                       
  9                                                                                DNA binding HTH domain, endonuclease                                        2.11                                       
  10                                                                           Herpes, measles, influenza A related pathway                                    2.09                                       
  Orthologs on the eutherian branch                                                                 1                                                          KRAB, C2H2-ZNF, transcription regulation   28.51
  2                                                                                    MAGE protein, tumor antigen                                             20.78                                      
  3                                                                 Innate immune response, defense response to bacterium, β-defensin                          17.46                                      
  4                                                                      Olfactory receptor, sensory transduction, GPCR activity                               17.05                                      
  5                                                                               Keratin, intermediate filament protein                                       11.36                                      
  6                                                                                C-type lectin, carbohydrate binding                                         7.06                                       
  7                                                                            Immune response, cytokine, cytokine receptor                                    4.6                                        
  8                                                                                       Immunoglobulin domain                                                4.2                                        
  9                                                              MHC I/II like antigen recognition protein, natural killer cell activity                       4.16                                       
  10                                                                               Keratinization, peptide crosslinking                                        3.66                                       
  Orthologs on the branches postplacental                                                           1                                                          Olfactory transduction, GPCR activity      32.51
  2                                                                              Protein deubiquitination, peptidase C19                                       13.21                                      
  3                                          Histone, epigenetic regulation of gene expression, transcriptional misregulation in cancer, viral carcinogenesis  10.93                                      
  4                                                                              KRAB, C2H2-ZNF, transcription regulation                                      8.79                                       
  5                                                                       Innate immune response, antibacterial humoral response                               7.87                                       
  6                                                                             Defense response to bacterium, β-defensin                                      7.14                                       
  7                                                                                   Cadherin, cell--cell adhesion                                            3.92                                       
  8                                                                               Fungicide, defense response to fungus                                        3.04                                       
  9                                                                                   GRIP, protein targeting Golgi                                            3.03                                       
  10                                                                      Serotonin pathway, neurotransmitter receptor activity                                2.19                                       

[note]{.smallcaps}.---SPRY, domain in SPla and the RYanodine receptor; RING, really interesting new gene; EF-hand, helix-loop-helix domain; GPCR, G-protein coupled receptor; KRAB, Krüppel-associated box; C2H2-ZNF, C2H2 zinc finger; HTH, helix-turn-helix; MAGE, melanoma-associated antigen; GRIP, glutamate receptor-interacting protein.

The most significantly enriched functional cluster for the orthologs appearing on the eutherian branch contains annotation terms such as "KRAB domain," "DNA binding," and "transcriptional regulation" ([table 1](#evaa026-T1){ref-type="table"}). KRAB domains and KRAB-ZFPs are a family of transcription factors known to be controllers of sequences derived from transposable elements. They are involved in embryonic development and genomic imprinting ([@evaa026-B14]; [@evaa026-B69]). Moreover, multiple rounds of KRAB-ZFP turnovers have occurred during different phases of mammalian evolution ([@evaa026-B35]; [@evaa026-B22]; [@evaa026-B90]). In fact, KRAB-ZFP family expansion has occurred independently in every vertebrate lineage ([@evaa026-B49]). The enrichment of KRAB-ZFPs at the eutherian branch could therefore suggest that there were changes in transcription regulation that co-opted sequences from TEs, including retroviruses, for the regulation of internal embryo development, as has been previously argued ([@evaa026-B51], [@evaa026-B52]). We also detected enrichment of the MAGE protein family, which were dramatically expanded in the mammal lineage ([@evaa026-B26]; [@evaa026-B84]). Type I MAGE genes are often expressed in tumor cells and are involved with cancer metastasis, the cellular mechanism of which parallels placenta invasiveness ([@evaa026-B20]). MAGE-like proteins are also involved in embryo implantation ([@evaa026-B9]). Orthologs on the eutherian branch were also enriched in functions including innate immune responses, natural killer cell activity, and MHC-like antigen recognition. Other research suggests that these immune-related genes might play an important role in preventing maternal--fetal immune rejection during pregnancy ([@evaa026-B24], [@evaa026-B25]). Another enriched annotation cluster was related to olfactory transduction, probably due to an expansion of olfactory receptor genes during the evolution of eutherian mammals. This observation could indicate the occurrence of important physiological functional changes but is plausibly unrelated to the evolution of the placenta ([@evaa026-B60]).

Tests of enrichment of protein families comparing the ancient orthologs (those appearing on the root branch) to orthologs appearing on the mammal branch, therian branch, and eutherian branch yielded similar results ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). KRAB transcription factors were enriched upon the appearance of placentation on the eutherian branch. Overrepresentation of these zinc-finger transcription factors could indicate that KRAB-mediated rewiring of genetic regulatory networks played an important role during the evolution of mammalian pregnancy. Cell adhesion molecules and immunoglobulin proteins were enriched in the mammal lineage, suggesting changes in immune modulation to avoid maternal--fetal immune conflict even at the early stages of mammalian evolution. Curiously, when we looked quite broadly at the functional roles of these genes of differing ages using the semantic similarity of the molecular function GO terms, we found that orthologs appearing on the eutherian branch have the most distinct GO terms relative to the other four groups ([supplementary fig. S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Gene Age as a Predictor of Tissue Invasion
------------------------------------------

In human and other eutherian mammals with hemochorial placenta, the trophoblasts will invade the uterine lining with behavior similar to metastatic cancer cells (e.g., proliferative signaling, evasion of apoptosis, and tissue invasion; [@evaa026-B31]). [Figure 3](#evaa026-F3){ref-type="fig"} shows orthologs from each internal branch that the cancer census has annotated with cancer hallmarks from the COSMIC database ([@evaa026-B29]; [@evaa026-B75]). Cancer census genes with unknown hallmarks are listed in [supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. The eutherian branch is again not unusual in the prevalence of cancer-associated genes that appeared along it. Instead, there is a direct relationship between the age of an ortholog and its tendency to be associated with cancer, with older orthologs showing the higher prevalence ([table 2](#evaa026-T2){ref-type="table"}, Fisher's exact tests). However, we observed a trend of an increasing relative prevalence of somatic mutations when moving from the root branch to the eutherian one ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), where this average somatic mutation propensity peaked. In other words, the orthologs appearing around the time of the eutherian origin are more often observed to have somatic mutations in humans, though we lack a clear hypothesis as to why. These results also complement previous findings that mutations in the very ancient genes shared by all metazoans can disrupt necessary interactions between those old genes and younger genes, driving inherited diseases or cancer ([@evaa026-B18]; [@evaa026-B76]).

![---Cancer hallmark genes for orthologs appearing on the mammal, therian, eutherian, and postplacental branches of the phylogeny in [figure 1](#evaa026-F1){ref-type="fig"}. Four groups of genes are listed based on their evolutionary time of appearance. Branch A is the mammal branch, B is the therian branch, C is the eutherian branch, and D--N are the postplacental branches. The association of these genes with cancer hallmarks were shown. Green circles stand for "promotes," dark blue circles represent "suppresses," and the aqua circles stand for both. The cancer hallmark annotations for cancer census genes were obtained from COSMIC release v88 ([@evaa026-B29]; [@evaa026-B75]). Ten cancer hallmarks are 1, proliferative signaling; 2, suppression of growth; 3, escaping immunic response to cancer; 4, cell replicative immortality; 5, tumor promoting inflammation; 6, invasion and metastasis; 7, angiogenesis; 8, genome instability and mutations; 9, escaping programed cell death; 10, change of cellular energetics.](evaa026f3){#evaa026-F3}

###### 

Number of Orthologs That Are Expressed in Placenta (from DAVID) and Are Annotated as Census Cancer Genes in COSMIC

  Branch                   Number of Genes    Number of Genes Expressed in Placenta[^a^](#tblfn2){ref-type="table-fn"}   Number of Cancer Census Genes[^a^](#tblfn2){ref-type="table-fn"}
  ------------------------ ----------------- -------------------------------------------------------------------------- ------------------------------------------------------------------
  Root branch              13,472                            7,215[\*\*](#tblfn4){ref-type="table-fn"}                               501[\*\*](#tblfn4){ref-type="table-fn"}
  Mammal branch            1,403                              616[\*\*](#tblfn4){ref-type="table-fn"}                                  24[\*](#tblfn3){ref-type="table-fn"}
  Therian branch           1,996                              785[\*\*](#tblfn4){ref-type="table-fn"}                                 34[\*\*](#tblfn4){ref-type="table-fn"}
  Eutherian branch         1,320                              388[\*\*](#tblfn4){ref-type="table-fn"}                                                   9
  Postplacental branches   1,174                                                277                                                                     4

The gene count on each branch was compared with the sum of gene counts on later branches, for example, number of orthologs on eutherian branch was compared with the number of orthologs on all branches after the eutherian branch. Fisher's exact tests were performed to examine whether the true odds ratios are \>1.

*P *\<* *0.05.

*P *\<* *0.0001.

At the level of specific genes, we found multiple cancer census genes that appeared on the therian and eutherian branches that are most likely involved with placenta formation. For example, on the therian branch, *NOTCH1* is involved in the Notch signaling pathway, which is crucial for maternal--fetal communication ([@evaa026-B16]); the *ERBB2* (erb-b2 receptor tyrosine kinase) gene is a member of the epidermal growth factor receptor family, which was previously shown to be expressed in both normal trophoblasts and malignant breast/ovary tumors ([@evaa026-B26]). On the eutherian branch, the *FAS* gene is involved in apoptotic processes, and it has been found under significant positive selection in mammalian evolution ([@evaa026-B78]). These findings hint at shared circuits between placenta invasion and cancer metastasis.

Discussion
==========

We illustrate the power of high-resolution comparative genomics to illuminate functional patterns in mammalian genomes. The evolutionary origin of mammalian pregnancy was a complex process involving maternal--fetal tissue fusion, gene regulatory network rewiring, as well as the emergence of a new cell type (decidual stromal cells) to repress the immune conflict between mother and developing offspring ([@evaa026-B87]; [@evaa026-B52]; [@evaa026-B8]). The predictions from our orthology inferences tend to support previous findings and illustrate again the value of the comparative method in resolving the origins of complex traits ([@evaa026-B86]). Most importantly, it appears that genes that allow for the key molecular and cellular events in placentation were preexisting and co-opted in the evolution of placenta ([@evaa026-B38]; [@evaa026-B8]; [@evaa026-B28]). A novel observation of our analyses is that although these "structural" genes were likely preexisting, there was a radiation of regulatory genes coincident with the placental radiation, potentially allowing for the evolution of new expression patterns among those existing genes. This hypothesis would in some respects parallel the observation that transposon-related sequences were repurposed at this same time to help regulate placental development ([@evaa026-B51]; [@evaa026-B23]). Indeed, the eutherian branch shows the strongest pattern of the development of new genes from TE substrates, suggesting the importance of this TE repurposing in placental mammal evolution. As an aside, we note that all of these estimates have resolution no finer than the total length of the common eutherian branch: Our results do not argue for an unusually rapid origin of such regulatory innovations ([fig. 1](#evaa026-F1){ref-type="fig"}).

This use of preexisting genes for new purposes is a very common, if perhaps insufficiently recognized, feature of evolution. Morphological examples include not only Gould's famous panda's thumb but also the skin glands responsible for milk production ([@evaa026-B27]; [@evaa026-B62]). At a molecular level, another part of the mammalian trait repertoire also shows strong evidence of such repurposing: a number of proteins seem to have been co-opted during the evolution of milk fat globule production, in one case perhaps through the employment of an enzyme in a nonenzymatic role ([@evaa026-B80]; [@evaa026-B63]; [@evaa026-B62]). Although none of the mammalian innovations discussed involve the duplication of the genes, there are clear parallels to that process: In both cases, the evolutionary innovation is (not surprisingly) based on the co-option of either primary or secondary protein activities ([@evaa026-B13]).

We have also explored the power of our comparative genomic pipeline for uncovering new features of the genomics of cancer by inferring a list of genes that are likely involved in tissue invasion in both placenta and tumors ([@evaa026-B15]). As seen in the analysis of placental genes, the eutherian branch is not unusual in the proportion of cancer-associated genes originating along it: Instead the genes so involved are generally evolutionarily ancient. In terms of somatic mutation counts, genes that appeared around the eutherian origin show slightly elevated average number of mutations, perhaps representing a tradeoff between the intolerance of older genes to such mutations and the chances of a mutation having a cell-lineage level selective effect. In this study, we only reported well-annotated cancer census genes originating on each branch. However, this criterion might be too stringent, with more genes related to cancer metastasis and placenta formation yet to be explored. Nonetheless, we did find some known associations: our pipeline finds that the placenta-specific protein 1 (*PLAC1*) gene arose at the origin of eutherian mammals, which agrees with previous findings ([@evaa026-B17]). Notably, *PLAC1* is a promising target for development of antibody-drug conjugates ([@evaa026-B59]; [@evaa026-B95]). A more striking finding is that even within the confined temporal window of the mammalian radiation, phylogenetically older genes are more likely to show associations with cancer.

More generally, as comparative genomics resources continue to improve, researchers in many areas of biological research will find it fruitful to refine their hypotheses regarding gene function in light of the evolutionary history of the genes in question: It makes little sense to explore old and evolutionarily conserved pathways by studying species-specific genes. On the other hand, the work above reminds us of the opportunistic and unexpected manner in which evolution can repurpose old genes to new functions. Detecting such innovations remains an important open problem in molecular evolution.

Supplementary Material
======================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at *Genome Biology and Evolution* online.

Supplementary Material
======================

###### 

Click here for additional data file.

We would like to thank V. Lynch, J.C. Pires, M. Lucy, T. Spencer, and J. Thorne for helpful conversations and H. Truong for computational advice. This work was supported by National Science Foundation Grants NSF-IOS-1339156 (Y.H. and G.C.C.), NSF-DBI-1358997 (M.B., H.J.L., K.B., T.M., and G.C.C.), NSF-CCF-1421765 (G.C.C.), European Union\'s Seventh Framework Programme Grant 613689 (Y.H. and G.C.C.) , and by the Triangle Comparative and Evolutionary Medicine Center (Y.H.).
