{'fraction_non_alphanumeric': 0.02881835421394597, 'fraction_numerical': 0.02835203780271707, 'mean_word_length': 5.003732736095558, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genome-wide enhancer - gene regulatory maps in two vertebrate genomes

Gene regulation by enhancer sequences controls the spatiotemporal expression of target genes, often over long genomic distances. Current methods used to identify these interactions in the human genome require complex experimental setups and are restricted to specific cell types. Here we use PEGASUS, an approach that captures evolutionary signals of synteny conservation to predict such interaction potentially in any biological context. We apply it to the human and zebrafish genomes and exploit the 1.3 million human and 55,000 zebrafish enhancers that we predicted to uncover fundamental principles of gene enhancer function in vertebrates. We show that the number of enhancers linked to a gene positively correlates with expression breadth in space and time and that the enhancer-target distance is evolutionarily neutral. We uncover almost 2000 regulatory interactions ancestral to vertebrates, which are strongly enriched in core developmental processes. Using PEGASUS and its evolutionary view of enhancer-gene interactions, we provide a highly complementary resource to functional assays which uncovers key principles of enhancer biology.


15
Enhancers are short DNA sequences that bind transcription factor and contact promoters in cis to activate or repress the transcription of protein coding genes into messenger RNA 1 . This control -or regulation -of gene expression by enhancers ensures the fine tuning of protein abundance in cells, in the course of their differentiation during development but also in response to various stimuli in specialised cell types [1][2][3] . Understanding how enhancers 20 regulate gene expression is therefore essential to understand how organisms develop and function throughout life. In fact, disruption of enhancer function has been shown to lead to abnormal gene expression and thus to diseases [4][5][6][7][8] . In addition, the majority of variants identified in Genome Wide Association Studies (GWAS) are found outside coding sequences 9,10 . Together with the observation that many patients remain undiagnosed after 25 genome sequencing because no plausible coding variant can be incriminated 11 , these results underscore the importance of identifying enhancers and their target genes to better understand genome biology. Numerous methods have been developed to identify enhancers. Some are based on the analysis of the conservation of non-coding sequences [12][13][14][15] . Others analyse the binding of 30 enhancer-associated proteins on the genome 16,17 , or identify chromatin signatures (e.g. histone modifications [18][19][20] or open chromatin 21,22 ). However, identifying the target genes of a given enhancer remains difficult because enhancers can regulate genes located at long distances, typically several hundreds of kilobases (kb) away 23,24 . Methods that have been developed to address this problem fall in three broad categories. The first uses 35 chromosomal conformation capture techniques to identify physical interaction between two loci in the genome [25][26][27][28][29][30] , but it is not clear which of these interactions are linked to a regulatory role. The second measures the correlation of transcription activity between noncoding sequences and nearby genes 31,32 , assuming the two are signatures of a coordinated regulatory function. Finally, single nucleotide variants (SNVs) can be associated to significant 40 differences in gene expression, thus qualifying as expression QTL (eQTL) that presumably reside within or close to enhancers 33 . These methods all rely on a specific experimental set up, i.e. a particular cell line or tissue upon which the experiment can be carried out, therefore investigating regulatory interactions in a specific but narrow biological context. We previously developed PEGASUS (Predicting Enhancer Gene ASsociations Using Synteny), 45 a method to identify the target genes of predicted enhancers using an evolutionary signal instead of an experimental assay 34 , therefore potentially identifying all functional interaction that are under negative selection and hence conserved. This is of particular interest for regulatory functions important during human development, where experimental assays are limited. We tested this method on the human X chromosome followed by experimental 50 validations of more than 1,000 predicted interactions using transgenic assays 34 .
Here, we applied this method on the entire human and zebrafish genomes in two independent analyses. We predict the interactions between 1,311,643 human enhancers and 18,339 target genes, and between 55,515 zebrafish enhancers and 17,363 target genes. By comparing human and zebrafish predictions, we outline a set of genes with conserved cis-55 regulation in vertebrates enriched in brain and development functions. We find evidence for a direct link between the expression breadth of human genes and the number of enhancers predicted to target them, explaining at least in part the broader expression patterns of genes with promoters overlapping CpG islands. We find that the average distance separating enhancer and their target genes scales with genome size, demonstrating that little selective 60 pressure acts to preserve this distance. In addition, highly conserved enhancers are more likely to reside within the same Topologically Associating Domain (TAD) as their target gene, irrespective of the distance separating them. This evolutionary signal provides a basis to explain the observed conservation of TADs in vertebrates 35,36 and sheds light on the complex relationship between 3D genome organisation and the regulation of gene 65 transcription.

Identification of enhancers & target genes in two vertebrate genomes
We applied the PEGASUS method 34 to identify predicted enhancers and their most likely target genes in the human and zebrafish genomes. We first analysed the UCSC 100-way 70 multiple genome alignment restricted to Sarcopterygii (all vertebrates except fish but including Coelacanth) to identify human conserved non-coding elements (CNEs). We then computed an evolutionary linkage score between each CNE and neighbouring genes within a 1 Mb radius, thus identifying the linkage that is the most conserved. This resulted in predicting 1,311,643 CNEs linked to 18,339 human genes. We replicated the analysis in 75 zebrafish using a multiple alignment of 7 teleost fish genomes, leading to the identification of 55,515 CNEs linked to 17,363 genes. Of note, the zebrafish was compared to 6 genomes and shares a 300 million years old ancestor with the closest species of this set, while in contrast human was compared to 35 genomes, including many placental mammals who share a ~90 million years old common ancestor (Supplementary Figure 1). By reducing the 80 phylogenetic sampling used to analysing the human genome to a situation equivalent to that of zebrafish, we show that this contrast is sufficient to explain the vastly different numbers of CNEs identified in zebrafish (Supplementary material). Therefore, the number of zebrafish CNEs identified here are likely to be substantially underestimated. CNEs are linked to a median number of 3 genes in human and 4 genes in zebrafish, while genes are linked to a 85 median number of 301 enhancers in human and 5 enhancers in zebrafish. These CNEs are short, their size ranging from 10 bp to 1,709 bp in human (with a median of 27 bp), and from 10 bp to 1,628 bp in zebrafish (with a median of 90 bp; Table 1). CNEs are relatively close to their linked genes: the median distance to TSS was 353,488 bp in human and 288,643 bp in zebrafish (Supplementary Figure 2), with 29% of enhancer -TSS distances in human and 37% 90 in zebrafish smaller than 200 kbp.
Because PEGASUS can predict CNE -gene interactions that skip over unrelated genes, we defined three categories of CNEs: intronic (located on an intron of their linked gene), flanking (intergenic CNE flanking their linked gene, and jumping (CNEs skipping over intervening genes to connect to their linked genes). We found that 20% of predicted 95 enhancers are intronic, 53% flanking and 27% jumping in human, while in zebrafish, 15%, 36% and 49% of predicted enhancers are respectively intronic, flanking and jumping. Interestingly 34% of intronic human CNEs (37% in zebrafish) are linked to a gene that is not their "host" gene. These results suggest that choosing the nearest neighbour as potential target gene, as if often done short of a better alternative, is likely to fail in more than 27% of cases (49% in 100 zebrafish).  CNEs are assigned a linkage score as long as a gene lies within a 1 Mb radius. However, not all CNEs are regulatory enhancers and we wished to confirm that the linkage score is correlated with evidence for a regulatory role in the human and zebrafish genomes, as noted previously for human chromosome X 34  PEGASUS target gene predictions overlap experimentally defined regulatory interactions 125 We sought to measure the agreement between PEGASUS target gene predictions and those from tissue or cell-type dependent experimental approaches, restricting our analysis to the 30% of CNEs linked to a single gene. While 28% of our target gene predictions agree with Capture Hi-C interactions 27,30 , the level of agreement depends on cell type: from 22% for Embryonic Stem Cells (ESC) to 29% for GM12878 (a lymphoblastic cell line). The level of 130 agreement with eQTLs-based predictions 33 , is higher, with a 57% overlap. Interestingly, in all dataset that we analysed, the agreement decreased with distance to the TSS (Supplementary with dataset obtained on a limited number of cell or tissue types. The level of agreement between PEGASUS predictions and methods based on functional assays is within the same range as when the latter are compared to each other (between 10 and 48% of eQTLs-based predictions agreed with Capture Hi-C predictions, depending on the cell type), and here also, the overlap decreases with the distance to the TSS (Supplementary Figure 6). We 140 therefore conclude that PEGASUS predictions enhancer-TSS predictions are consistent with results from genome wide functional assays.

Long-range regulatory interactions within TADs
Topologically Associating Domains (TADs 35 ) have been shown to coincide well with the 155 regulatory landscape governing gene expression 6,41 . Here, 57% and 66% of CNEs linked to a single gene indeed reside within a TAD in hESCs and IMR90 cells respectively, compared to . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018; an average of 34% and 40% respectively when we shuffle TAD intervals. CNEs are linked to their target gene with a higher score when inside the same TAD (mean scores, 0.72 inside vs 0.67 outside for hESCs, 0.71 vs 0.68 for IMR90, Figure 2), overlapped more with functional 160 marks (13.9% vs 9.7% for H3K4me1, 14.9% vs 11.1% for H3K27ac in hESCs) and were also closer to each other (median distance to TSS, 332 kb inside vs 522 kb outside for hESCs, Figure 2). Finally, we see a striking link between predicted enhancer conservation and their localisation within a TAD. First, enhancers-target gene pairs located within the same TAD have more often orthologous enhancer-targets conserved in zebrafish than other pairs 165 (0.23% vs. 0.18%). More importantly, we see a positive link between the conservation depth of enhancers and the co-localisation of enhancers and target genes within the same TAD (Supplementary Figure 7). The association between TAD co-localisation and both distance to TSS and conservation depth cannot be explained by chance alone (Supplementary Figure 7). In line with previous observations 42 , our results therefore support the existence of a link 170 between the control of gene expression and 3D chromatin organisation.

Figure 2: Overlap between topologically associating domains and regulatory interactions.
Percentage of predicted enhancer -target gene interactions located in the same TAD as a function 175 of distance (left panel) or number of species in which the predicted enhancer is conserved (right panel). TADs were obtained from two cell types (IMR90 and hESCs 35 ).

Function of regulatory interactions conserved in vertebrates
The human and zebrafish CNE-gene maps described here were built independently using 180 distinct sets of genome, allowing us to use their intersection as increased evidence for functional and highly conserved enhancer -target gene interactions. In addition, these interactions were probably inherited from their last common Euteleostomi ancestor. We identified 1,986 human CNEs and 1,949 zebrafish CNEs conserved between human and zebrafish (methods) and linked to 567 human genes and 608 zebrafish orthologous genes 185 respectively. Gene Ontology term analysis (methods) shows that these ancestral regulatory circuits are overwhelmingly enriched in neuronal, skeleton, muscle, vascular and renal development ( Table 2). To illustrate this, we used the CRISPR-Cas9 technology to delete a 236 bp zebrafish CNE predicted to target the irx1b gene. Its orthologous CNE targets IRX1 in the human genome. This gene plays multiple roles during pattern formation of vertebrate 190 embryos 43 , and we expect its expression pattern to be tightly regulated by a complex array of enhancers. The deletion of the CNE greatly decreases the expression of the endogenous gene in several structures of the zebrafish embryo ( Figure 3) establishing it as a bona fide . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018; developmental enhancer. The strong enrichment in core developmental functions observed with orthologous PEGASUS predictions is consistent with earlier observations, as enhancers 195 identified through sequence conservation are often found to be active during development, especially in the nervous system 12   CNE -gene distances scale with genome size The "action range" of enhancers is known to encompass a wide span, from within the target 215 gene itself to more than 1 Mb away. It is therefore reasonable to postulate that little constraints exist on maintaining a specific genomic distance between enhancers and their target genes, especially as they are thought to physically contact each other through DNA looping in the nucleus. However, to our knowledge, this has never been formally tested at a . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018; genomic scale. The genomes of human and zebrafish are very different in size, 3.1 Gb and 1.5 220 Gb respectively, allowing us to examine this question using the ~ 2000 orthologous humanzebrafish CNE -gene linkage: do their respective genomic span scales with genome size or instead remain conserved to some degree? We used the sizes of orthologous introns to estimate the relative neutral evolution of genomic distances, since introns are thought to be under negligible size constraint above a small minimal threshold. Results show that distances 225 between orthologous CNEs and their orthologous target genes scale precisely with intron size (Figure 4). There is therefore no evidence for a constraint on CNE -gene distance, suggesting that 3D folding principles governing enhancer -promoter contacts are not sensitive to this distance.

Figure 4: Comparison of enhancer -target genes interaction distances.
Pairwise comparison of interactions distances between human -zebrafish orthologous genes with conserved predicted enhancers, and of intron length in the same set of genes. Interaction distances are computed as enhancer-TSS distances. Comparisons are computed as log2 values of human/zebrafish ratios, with each pair of orthologous genes being studied independently.

235
Regulation complexity is linked to expression breadth Enhancers are largely responsible for the spatio-temporal specificity of their target gene's expression (e.g. SOX2 47 , SOX9 48 or HOX genes 49 ). It naturally follows that the complexity of a gene's expression pattern might be correlated with the number of enhancers that regulates 240 its expression, although this has never been demonstrated. We used our predicted enhancers to investigate this question. We observe a positive link between the number of enhancers targeting a gene and the number of adult tissues where this gene is expressed (ρ = 0.23, p-value < 10 -15 , Figure 5). Genes cover a broad range of tissue specificities, from 'housekeeping genes' required for generic cellular functions and expressed in most tissues to 245 tightly regulated developmental genes sometimes expressed in just a few cells in a short window of time. The presence of a high frequency of CpG dinucleotides in the promoter of genes (so called CpG islands) is a useful quantitative marker of broad expression 50-53 . We split our target genes between those which TSS overlap a CpG island (referred to as CpG genes) and other genes (referred to as non-CpG genes). First, we confirm that CpG genes 250 are more broadly expressed than non-CpG genes (one-sided Kolmogorov-Smirnov test pvalue < 10 -15 , Supplementary Figure 8). Second, we observe that CpG genes are targeted by a higher number of predicted enhancers than non-CpG genes (365 vs 314 on average, Kolmogorov-Smirnov test p-value < 10 -15 , Supplementary Figure 8), thus providing an immediate mechanism for their broader expression. The link between predicted enhancer 255 number and expression breadth is stronger for non-CpG genes than for CpG genes (ρ = 0.29 & 0.16, respectively, p-values < 10 -15 ). But this is explained by the narrower (but higher) . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018; range of tissue number where CpG genes are expressed ( Figure 5). Indeed, reducing the range of expression for non-CpG genes to a level similar to CpG genes decreases correlations to similar levels as CpG genes (Supplementary Figures 9). Our results show that 260 expression patterns are likely to be driven by regulation complexity, and that the effect is stronger for non CpG genes. We asked if, in addition to expression complexity, a similar influence of predicted enhancers could be observed for expression levels of human genes. However, the number of predicted enhancers only weakly correlates with expression levels in 17 out of 21 tissues (p-value <= 0.05, spearman's ρ from -0.08 to 0.15). Regulation 265 complexity does not seem to be correlated with expression levels in human tissues.

Discussion
We applied the PEGASUS method to predict ~1,300,000 human and ~55,000 zebrafish enhancers targeting the majority of the genes in their respective genomes. After finding evidence for a regulatory role of these interactions, we used this dataset to make a number of fundamental observations related to enhancer function. We show that the number of 280 enhancers that regulate a gene positively correlates with its breadth of expression, that the . CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018; distance between predicted enhancers and their target gene evolves neutrally, that chromatin organisation and regulation of gene expression are linked and that regulatory interaction ancestral to vertebrates concentrate on core functions necessary to build an organism (skeleton, muscles, neurons, vascular system). 285 The challenges of predicting long-distance regulatory interactions Most methods currently employed to predict long-range regulatory interactions in the human genome rely on specific cell lines 26-31, 33 . These methods usually differ in their approach and the cell types or tissues they study, which might limit expectations to observe overlap in 290 their predictions, especially given that most enhancers are tissue-specific 29,[54][55][56] . It is therefore interesting to observe that when comparing these methods with one another, the level of agreement between eQTLs-based predictions ranges from 10% to 48% (between 10% and 16% for ESCs and NECs cells, between 38% and 48% for CD34 and GM12878 cells). In contrast, PEGASUS is agnostic to cell-type or tissue context. The sole rationale underlying 295 the predictions is that the interactions are functional, therefore under sufficient evolutionary conservation to be picked up by comparisons with other genomes. This property introduces its own bias, in that conserved interactions are more likely to relate to developmental processes 16,57 , which are typically harder to identify in differentiated cell lines. We show that our predictions agree with functional assays in a manner consistent with the intersection 300 between the latter studies, suggesting that an evolutionary-based approach can capture more than 50% of interactions observed in experimental assays. Remaining, non-overlapping interactions likely include false positives in either type of methods, which are difficult to quantify at present. Another source of discrepancy may be that PEGASUS misses cell-type specific enhancers that do not show enough sequence conservation to be identified by 305 evolutionary linkage: given the rapid evolutionary turnover of enhancer regions during evolution [58][59][60] , it is likely that a fraction of cell-type specific enhancers have had little time to leave detectable footprints of selection in a genome. We observed a general trend for predictions, whether experimental or evolutionary based, to overlap consistently less with increasing distance between the predicted enhancer and the TSS of the target gene. Our 310 data suggests that this could be explained by long distance enhancers being more tissuespecific than short distance enhancers. This has important implications for regulatory interactions predictions.
No evidence for natural selection acting on enhancer -target gene distances Our results suggest that there is no selective pressure acting on the enhancer -target gene 315 interaction distance to ensure proper regulation. These results can easily be explained by how enhancer regulation is mediated through the 3D organisation of the genome. Metazoan genomes are organised in TADs (for topologically associating domains), large units of chromosomal interactions 61,62 that are mostly conserved between cell types and species 35,36,63 . Enhancer -target genes interactions occur mostly within TADs 64 , via DNA looping 65 . 320 Interestingly, it has been recently shown that the position of an enhancer relative to its target within a TAD has no effect on its regulatory potential 66 . These results are consistent with our finding that the relative position of an enhancer relative to its target gene is effectively neutral. As a result, the same forces that affect genome sizes will affect the interaction distances. A recent analysis of genomic regulatory blocks (or GRBs) in metazoans 325 based on the analysis of clusters of conserved non-coding elements showed that a) these blocks correlated well with known TADs and b) their sizes seem to correlate well with genome size 67 , providing further evidence that interaction distances between enhancers and target genes are under the same forces that affect genome size in metazoans.
. CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018;

Conclusions 330
This study provides a unique view of the conservation and evolution of enhancers in vertebrate genomes. Our results support the model where gene expression complexity is driven by regulatory region diversity and highlight the biological functions which regulation has been conserved since the vertebrate ancestor. Moreover, the PEGASUS method provides a robust tissue and life stage agnostic target gene prediction method that opens 335 research possibilities in the study of gene regulation in a wide number of species.

Materials & Methods
Defining conserved non-coding elements and their most likely target genes We used a previously published method to identify predicted enhancers and their most likely target genes 34 . This method first identifies predicted enhancers as conserved non-340 coding elements (or CNEs for short) in multiple alignments, and second identifies a predicted enhancer's most likely target gene as the gene in its vicinity with the most conserved synteny, through the computation of a co-segregation score measuring this conservation.
We identified CNEs in the human genome (hg19 version) using the UCSC 100-way multiple 345 alignments that we restricted to 35 Sarcopterygii species (a full list is available in Supplementary Table 3). CNEs had to be identified in at least 6 species (including human) with one non-primate species to be considered. We identified CNEs in the zebrafish genome (danRer7 version) in custom-made multiple alignments that include 6 other Neopterygii species (a full list is available in Supplementary Table 4). Multiple alignments 350 were built first by making pairwise alignments between zebrafish and other species using LastZ 68 , then by using these alignments to build multiple alignments using Multiz 69 . Alignment parameters can be found in Supplementary Materials. In these multiple alignments, CNEs had to be identified in at least 3 species (including zebrafish) to be considered. We used a seed of 10 bp (with 90% conservation) and a maximum fraction of mismatches for a column to be 355 considered as pseudo-conserved of 40% for both analyses.
Most likely target genes were identified with the same method in both genomes. We first identified all protein-coding genes in a 1Mb radius around CNEs in the reference genome and computed a co-segregation score for each gene as a measure of the conservation of synteny between a CNE and a gene. Considering only species in which the CNE has been 360 identified, the score sums over all species and adds a term when the gene is identified, on the same chromosome and in the radius, and subtracts a term when the gene is either absent from the genome, on a different chromosome or outside the radius. Terms are normalized by pairwise rearrangements rates and genome sequencing coverage (see 34 for more details). For a given CNE, the most likely target gene is the gene with the highest co-365 segregation score. Adjacent CNEs targeting identical gene, present in the same species, having identical co-segregation scores and distant by less than 100bp were merged together. Predicted enhancers located at 100bp or less away from an exon were discarded.

Overlap with functional marks & enhancer predictions
We analysed the overlap of our predicted enhancers with functional marks and previous 370 enhancer predictions. We first looked at the overlap with ChIP-seq peaks of histone modifications (namely H3K27ac, H3K4me1, H3K4me3 and H3K27me3) in embryonic stem cells in human 20 and across various developmental stages in zebrafish 37 . We also looked at the overlap with peaks of ATAC-seq in zebrafish 70 . We also looked at the overlap of our predicted enhancers with other enhancer predictions. In human, we looked at enhancers 375 from the FANTOM project 31 or from the Vista database 39 . In zebrafish we looked at enhancer predictions from differentially methylated regions 40 .

Gene expression data
Gene expression values and calls were downloaded from the Bgee database 71 . For each gene, we computed the number of adult tissues for human or developmental stage for zebrafish in 380 which a gene is called as expressed. We filtered out terms that had daughter terms in our lists of tissues or stages. For human, we also retrieved expression levels in 21 tissues.
Comparing target gene predictions with chromosome conformation capture predictions We compared our enhancer -target gene predictions in human with other predictions 385 based on correlation of transcriptional activity between predicted enhancers and target genes 31 , on eQTLs in various tissues 33 or on Capture Hi-C data in CD34 and GM12878 cells 27 and in embryonic stem cell (ESCs) and ESC-derived neurodectodermal cells (NECs) 30 . The latter contains interactions present in either cell types or in both simultaneously (labelled "Both"). We compared predictions between our datasets and others firstly by 390 identifying overlapping enhancers between our predictions and other datasets, then by identifying identical interactions predictions in these overlapping regions.

Defining orthologous enhancers & target genes between human & zebrafish
We downloaded human-zebrafish and zebrafish-human pairwise chain alignments from UCSC. We defined orthologous predicted enhancers as human and zebrafish enhancers that 395 overlapped by at least 10bp on either pairwise alignment. We next downloaded humanzebrafish orthologous genes from the Ensembl database (version 75) to identify orthologous enhancers targeting orthologous genes.
Because of the evolutionary distance between human and zebrafish, some orthologous regions are difficult to align and are thus impossible to detect. To circumvent this problem, 400 we used the spotted gar genome 72 to identify orthologous predicted enhancers. We downloaded human-spotted gar pairwise chain alignments and used our custom-made zebrafish-spotted gar pairwise chain alignments to respectively map human and zebrafish predicted enhancers on to the spotted gar genome. We considered human and zebrafish enhancers as orthologous if they overlapped by at least 10bp on the spotted gar genome.

405
We identified orthologous targets by looking at the orthologous genes set used above. Both direct and via spotted gar sets of orthologous enhancers -orthologous target genes were merged together.

Gene expression enrichment analysis
We performed GO terms enrichment analyses using the PantherDB website 73,74 . The 410 website performs a Fisher's exact test with FDR multiple test corrections.

Distance to transcription start sites
We downloaded transcription start sites (or TSS) locations from the Ensembl database (version 75). For each gene, we considered only the transcript annotated in RefSeq giving the longest protein.
We computed for each enhancer -target gene the distance to the TSS 415 as the shortest distance from enhancer boundary to the target's TSS.
. CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was . http://dx.doi.org/10.1101/244475 doi: bioRxiv preprint first posted online Jan. 7, 2018;

Topologically associating domains
We downloaded topologically associating domains (or TADs) coordinates for two cell types, human embryonic stem cells (hESCs) and IMR90 fibroblasts 35 . We converted these coordinates from hg18 to hg19 using the liftOver webtool available at the UCSC genome 420 browser 75 . For enhancers targeting only one gene, we computed for each cell type whether both an enhancer and its target gene were located in the same TAD. We also generated shuffled TAD set by shuffling the TADs coordinates along the human genome.
In situ hybridization: 435 In situ hybridization were performed as described 78 , using an Irx1b probe corresponding to exon 2.