Eutheria

TBA [19] (link) alignments of the human genome (hg18) to 43 other vertebrate species were obtained from the UCSC genome browser [20] (link), [21] (link) together with a phylogenetic tree with the generally accepted topology (Fig S1) and neutral branch lengths estimated from 4-fold degenerate sites. Both the tree and alignments were projected to the 34 mammalian species. The alignment was compressed to remove gaps in the human sequence, and GERP++ scores were computed for every position with at least 3 ungapped species present, or approximately 88.9% of the 3.08 billion positions on the 22 autosomes and X/Y chromosomes. We used the HKY85 [13] (link) model of evolution with the transition/transversion ratio set to 2.0 and nucleotide frequencies estimated from the multiple alignment.
To limit memory requirements and allow parallelization of the constrained element computation, each chromosome was broken up into regions of approximately 2 megabases, with long segments where no RS score was computed chosen as boundaries. These boundary segments contain no information usable by GERP++ and because the algorithm never annotates constrained elements spanning them, excluding such segments did not sacrifice any predictive ability. These boundary regions made up approximately 6.8% of the human genome, including a 30.2 megabase region that made up more than half of chromosome Y. Constrained element predictions were generated using default parameters and a 5% false positive cutoff measured in terms of number of predictions; the estimated nucleotide-level false positive rate was under 1%. As additional validation, we computed overlap between our predictions and a set of ancestral repeats (L2) annotated by RepeatMasker. We found the overlap to be in line with what we expected given our estimated false positive rates: about 5% of the repeats overlap a predicted CE, with around 1.6% nucleotide-level overlap.
Gene, noncoding RNA, and PhastCons conserved element annotations were obtained from the UCSC genome browser's [20] (link), [21] (link) Known Genes [22] (link), RNA Genes, and Conservation [4] (link) tracks respectively. To avoid skewed statistics due to alternative splicing, gene annotations were resolved to a consistent nonoverlapping set where any segment belonging to multiple conflicting annotations was assigned a single annotation in the following order of priority: coding exon, 5′ UTR, 3′ UTR, intron. For meaningful comparison against phastCons, separate GERP++ scores and constrained elements were generated according to the same procedure as above but using only placental mammal data (ignoring platypus and opossum in the alignment and projecting them out of the phylogenetic tree).
PolII binding regions were defined as 50 bp upstream and downstream of PolII binding ‘peaks’ as identified from ChIP-seq experiments performed by the ENCODE Consortium [3] (link). A 100 bp window allows capture of the likely PolII binding site and its flanking sequence. We obtained data from nine ChIP-seq experiments conducted in two labs (the Snyder lab at Yale and the Myers lab at Hudson Alpha) on six cell types. Data was downloaded through the DCC at UCSC (ftp://encodeftp.cse.ucsc.edu). All data have passed publication embargo periods. Overlap statistics were calculated as described above for other annotation sets and averaged across all nine experiments.

The algorithm for mapping protein binding sites on the RNA sequences is based on our
WR approach (24 (link)), previously exploited in the
SFmap web server for mapping SF binding sites (23 (link)). The mapping algorithm considers the clustering propensity of the
binding sites and the overall tendency of regulatory regions to be conserved (24 (link)). In RBPmap we have improved the algorithm
by adding new features including the ability to map PSSM motifs, a
conservation-based filtering to reduce the rate of false-positive predictions and a
new background model which is specific to different genomic regions, namely intronic
regions flanking the splice sites, internal exons, exons in 5’ and
3’ UTR regions, non-coding RNAs and mid-intron/intergenic regions (a
detailed description of RBPmap algorithm is given in Supplementary file 1). A
pipeline summarizing RBPmap algorithm is shown in Figure 1. Briefly, given an experimentally defined motif (provided as
either a consensus sequence or a PSSM) and a query sequence (Figure 1A), RBPmap computes the match score for the
motif per each position in the sequence in overlapping windows (Figure 1B). The match score is then compared to a
background that is calculated specifically per each motif, filtering out all matches
below a significant threshold (default P-value<0.005)
(Figure 1C). At the next step, the WR function
is employed to calculate the multiplicity score which reflects the propensity of
suboptimal motifs (default P-value<0.01) to cluster around
the significant motif in a window of 50 nts, weighted by their match to the motif of
interest (24 (link)) (Figure 1D). Further, to reduce false-positive predictions, the final WR
scores are compared to a background model that is calculated independently per each
motif for the relevant genomic region. A Z-score is calculated for each WR score and
coupled to a P-value, which represents the probability of obtaining
a specific Z-score, considering a normal one-tailed distribution. RBPmap requires
that the final WR score of a site will be significantly greater (with
P-value<0.05) than the mean score calculated for the
background, in order to consider this site as a predicted binding site (Figure 1E). The new position-specific background model
provides more accurate and specific thresholds for the different regulatory regions
on the RNA (see above). For sequences from genomes other than human, mouse or
Drosophila, the WR scores are compared to a theoretical
threshold instead of the genome-specific background model which cannot be obtained
(see Supplementary file 1). This threshold is calculated for each motif separately,
according to its length and complexity (23 (link)).
At the last stage, we have added to the WR approach a conservation-based filtering,
which exploits the tendency of regulatory regions to be evolutionary conserved. The
conservation filter is optional and is applied only to sites that are mapped to
mid-intron/intergenic regions on the query sequence. These positions are removed
from the results if the mean conservation score of their environment is lower than
the mean conservation score calculated for intronic regulatory regions (Figure 1F). For sequences from human and mouse, the
conservation information is retrieved from the UCSC phyloP conservation table (28 (link)), based on the conservation of all placental
mammals. For Drosophila sequences we use the phastCons insect
conservation table (28 (link)). Both the
position-specific background model and the conservation filtering are applied only
for motifs which are searched in human, mouse or Drosophila sequences.

We compared AGORA’s v.92 eutherian reconstructions to DESCHRAMBLER’s^{33 (link)} (300 kb resolution: APCF_hg19_merged.map from http://bioinfo.konkuk.ac.kr/DESCHRAMBLER/). Because DESCHRAMBLER uses segments of the human genome as units of the reconstruction and was based on the hg19 genome assembly, we converted those regions to their protein-coding gene content and selected the genes still found in Ensembl v.92 and descendants of ancestral boreoeutherian genes. The Oxford grid plot was generated with the AGORA src/misc.compareGenomes.py script in ‘matrix’ mode.

The phylogeny was built based on two sets of UCEs: 5,472 baits for 5,060 UCEs in tetrapods^{57 (link)} and 2,628 baits for 1,314 UCEs in acanthomorphs^{69 (link)}. We used the Phyluce software^{70 (link)} to locate the probes in the reference genomes of our 68 species with 6 additional species contained in our original dataset. We extracted a flanking region of ±1,000 bp for each probe and aligned them with Mafft aligner version 7.470 (ref. ^{71 (link)}). We then created a 75% completion matrix, that is, each alignment contains at least 75% of the taxa (55 species), resulting in 63 alignments from the acanthomorph set and 2,742 probes from the tetrapod set (all alignments are available on Figshare). A phylogenetic tree was built using IQ-TREE version 2.0.3 (ref. ^{72 (link)}), with the appropriate substitution model inferred for each of the 2,805 alignments, a maximum likelihood tree search and 1,000 bootstrap replicates. To validate our tree, we also estimated a second tree based on a MultiZ alignment to the human genome and obtained similar results (Extended Data Fig. 9). The phylogenetic tree was calibrated to absolute time using the chronos function of the ‘ape’ package in R, with a smoothing parameter lambda of 0 and a ‘relaxed’ model^{73 (link),74 (link)}. Fourteen nodes were calibrated following previously published calibrations. The robustness of the tree was assessed by removing each node independently (see Extended Data Fig. 3).

Actinopterygii/Sarcopterygii: divergence time 416 million years ago (Ma), upper bound 425.4 Ma^{75 (link)}

The first node in the Actinopterygii group: divergence time 378.2 Ma^{76 (link)}

Sauropsida (birds and reptiles)/Synapsida (mammals): divergence time 313.4 Ma^{77 (link)}

Archosauria (birds)/Testudines: divergence time 260 Ma^{78 (link)}

The basal nodes of the Lepidosauria: divergence time 222.8 Ma^{79 (link)}

First mammalian node, Eutheria/Metatheria: divergence time 160.7 Ma^{75 (link)}

Galloanserae/Neoaves: divergence time 66 Ma^{77 (link)}

Glire/Primates: divergence time 61.7 Ma^{77 (link)}

Basal gekkotan node: divergence time 54 Ma^{80 (link)}

Passeriformes/Psittaciformes: divergence time 51.81 Ma^{81 (link)}

Cynoglossidae/Paralichthyidae: divergence time 50 Ma^{76 (link)}

Sus scrofa/other Cetartiodactyla: divergence time 48.5 Ma^{77 (link)}

Canidae/Arctoidea: divergence time 37.1 Ma^{75 (link)}

Hominoidea/Cercopithecoidea: divergence time 23.5 Ma^{77 (link)}

We estimated the BLAST-based [232 (link)] PAIs for a given human gene whose NCBI Entrez gene number served as input for our Orthoscape plug-in [229 (link),230 (link)] within the Cytoscape software suite [231 (link)]. The output was the most recent common ancestor of all the animal species whose DNA sequence of this gene is already known. The following evolutionary rank scale was used: 0, Cellular organisms; 1, Eukaryota; 2, Opisthokonta; 3, Metazoa; 4, Eumetazoa; 5, Bilateria; 6, Deuterostomia; 7, Chordata; 8, Craniata; 9, Vertebrata; 10, Gnathostomata; 11, Teleostomi; 12, Euteleostomi; 13, Sarcopterygii; 14, Dipnotetrapodomorpha; 15, Tetrapoda; 16, Amniota; 17, Mammalia; 18, Theria; 19, Eutheria; 20, Euarchontoglires; 21, Primates; 22, Haplorrhini; 23, Simiiformes; 24, Catarrhini; 25, Hominoidea; 26, Hominidae; 27, Homininae; and 28, Homo.