# On the importance of evolutionary constraint for regulatory sequence identification

François Giudicelli, Hugues Roest Crollius

► To cite this version:

François Giudicelli, Hugues Roest Crollius. On the importance of evolutionary constraint for regulatory sequence identification. Briefings in Functional Genomics and Proteomics, 2021, 20 (6), pp.361-369. 10.1093/bfgp/elab015. hal-03455906

HAL Id: hal-03455906

https://hal.science/hal-03455906v1

Submitted on 27 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the importance of evolutionary constraint for regulatory sequence identification

François Giudicelli and Hugues Roest Crollius

Institut de Biologie de l'École Normale Supérieure (IBENS), INSERM U 1024 / CNRS UMR 8197 : 46 rue d'Ulm, 75005 Paris

## Abstract

Regulation of gene expression relies on the activity of specialized genomic elements, enhancers or silencers, distributed over sometimes large distance from their target gene promoters. A significant part of vertebrate genomes consists in such regulatory elements, but their identification and that of their target genes remains challenging, due to the lack of clear signature at the nucleotide level. For many years the main hallmark used for identifying functional elements has been their sequence conservation between genomes of distant species, indicative of purifying selection. More recently, genome-wide biochemical assays have opened new avenues for detecting regulatory regions, shifting attention away from evolutionary constraints. Here, we review the respective contributions of comparative genomics and biochemical assays for the definition of regulatory elements and their targets and advocate that both sequence conservation and preserved synteny, taken as signature of functional constraint, remain essential tools in this task.

**Key words:** comparative genomics; gene regulation; enhancer; vertebrate evolution

## Introduction

Understanding how genetic information concealed in DNA sequence is turned into biological function in live cells, organisms or ecosystems remains one of the main goals driving current biological research. In this endeavour, it is of prime importance to properly identify the functional elements, which orchestrate the expression of genomes into phenotypes. Several types of genomic elements affect gene expression: proximal promoters, which define transcription start sites, distal enhancers and silencers, which bring regulatory complexes to the promoters, and insulators, which define the broad genomic domains where regulatory interactions take place. Here we will use the term 'regulatory elements' exclusively to mean distal cis-regulatory elements such as enhancers and silencers. It should be noted that behind the apparent simplicity of their conceptual definition, the operational definition of enhancers/silencers, i.e. the type of experimental evidence required for their validation, has fluctuated depending on research fields, times or even people. For example, the original description of enhancers in cells transfected with episomal vectors demanded that they

function independently of their orientation and distance to the target [1, 2], requirements that are seldom tested nowadays, at least in genomic contexts. The difficulties related to the definition of enhancers are discussed in detail in a recent review [3].

Indeed, while it is now relatively straightforward to identify the coding sequences of genes, recognizing regulatory elements that control their expression remains a much harder task. Two main reasons account for this contrast. First, coding sequences reside in exons that need to be transcribed and processed in characteristic ways (e.g. polyadenylation), which allow for their specific identification by biochemical isolation (e.g. polyA + RNA sequencing). In addition, their nucleotide sequence obeys the universal genetic code, which induces characteristic constraints on nucleotidic arrangements at small (codon triplets) as well as large (Open Reading Frame syntax) scales. In contrast, regulatory elements are more difficult to isolate biochemically and do not seem to obey a recognizable code at the sequence level.

However, it is generally admitted that the function of regulatory elements is dependent on their nucleotide sequence,
for example via the presence of transcription factor binding sites, offering a unique avenue for comparative genomics to fill this gap. Indeed, sequences driving evolutionary conserved functions are expected to be themselves evolutionary conserved. Conversely, since genomes result from hundreds of million years of evolution along which virtually every base has had an opportunity to vary, genomic elements that have resisted variation (i.e. conserved elements) have a high probability to be functional elements. More precisely, functional elements need not only be conserved, but they also need to be more conserved than expected under a neutral mutation rate, which can be estimated for example by counting synonymous mutations in coding sequences, assumed to be necessarily neutral, or by measuring divergence of ancestral repeats, viewed as archetypal neutral genomic elements devoid of organismal function. This equivalence between function and conservation led to the concept of 'phylogenetic footprint', which bases the identification of functional DNA elements on their evolutionary conservation in genomic comparisons [4, 5].

## Phylogenetic footprints define regulatory elements

In the pre-genomic era, the quest for phylogenetic footprints was necessarily restricted to small regions around genes of interest, yielding only scarce and spatially biased knowledge. For example, a considerable number of studies were dedicated to the globin loci, paradigmatic cases for the cis-regulatory control of gene expression by conserved elements [5-7]. With longer spans of genomes being sequenced in multiple species, the use of phylogenetic footprinting allowed the identification of far-acting enhancers located up to 1 Mb from their target [8-11]. When the 1st complete sequences of vertebrate genomes became available two decades ago, the field of comparative genomics rapidly expanded its discoveries. Comparing human and mouse genomes with the 1st fully sequenced fish genome, the pufferfish *Fugu rupripes*, led to the identification of several thousands highly conserved non-coding elements, named CNEs, which survived 450 M years of diverging evolution [12]. Similar findings were reported upon comparison of the human and elephant shark genomes [13], albeit without experimental confirmation in this case. Here, evolutionary constraints could be observed since Gnathostomata, the last common ancestor of humans and elephant sharks, which lived 530 M years ago. The initial focus was indeed placed on conservation over such extreme evolutionary distances as the most reliable indicator of regulatory function. This focus on ancient conservation helped install the zebrafish as one of the key species used in vertebrate enhancer studies, mainly because of its status as widespread experimental model for developmental biology. Not only could it be used to infer regulatory regions, but it also allowed experimental validation of these inferences using a wide range of transgenesis techniques [14]. It should be noted however that teleost fishes like zebrafish, fugu or medaka, which make one of the largest vertebrate groups with ~25 000 species, may not be the most appropriate basal group to study the gene regulatory networks underlying the developmental programme of vertebrates. According to [15], the accelerated evolution of their genomes that followed a whole-genome duplication event at the root of teleosts led to loss of many ancient regulatory elements, which may explain why Venkatesh et al. [13] found more CNEs conserved between human and elephant shark (4782) than between human and the less distant

fugu (2107) or zebrafish (2938) with the same sequence similarity thresholds.

The reason why conservation signals are more easily detected between distantly related species is because longer time spans since divergence mean higher accumulation of mutations at neutral sites, yielding better contrast against the lack of such mutations at constrained sites. However, high phylogenetic distance also translates in phenotypic divergence, itself underlain by the turnover of some ancient functional elements and the emergence of novel ones. Hence, although they offer more evolutionary contrast, comparisons between distant species are of course restricted to functions at least as ancient as the time of their divergence, inherently missing all subsequent regulatory innovations or exaptations. A kind of methodological trade-off ensues, where embracing a larger range of shared functions by comparing closer species comes with the price of lowering the signal-to-noise contrast. Comparisons between less distant genomes nevertheless proved to be also informative as long as more stringent scoring alignment criteria were used, such as perfect conservation between human, mouse and rat [12, 16]. Subsequent studies thus revealed thousands of CNEs recruited at various stages of mammalian evolution [17-19]. More recently, genome-wide characterization of zebrafish regulatory elements showed that while only a small proportion (~10%) had identifiable orthologues in humans, 88.6% of the remaining part showed a significant conservation among fish species [20]. This confirms that the amount of conserved regulatory sequence highly depends on phylogenetic distance. Attempts at extending that distance further back beyond the origin of vertebrates have so far yielded very few results. Sequencing of the Amphioxus genome, once thought to be the closest sister species to Vertebrata, allowed the alignment of only 56 CNEs susceptible to fulfil enhancer function [21]. Accordingly, most regulatory sequences later identified in Amphioxus had no counterparts in vertebrates [22]. Looking even further back, Clarke et al. [23] found a total of five vertebrate CNEs that were also shared with at least one protostome species. Two of these 'Bicores' (standing for Bilateria conserved regulatory elements) were found in the vicinity of the same genes in vertebrates and invertebrates. Interestingly, both vertebrate and invertebrate orthologues of these two Bicores were able to drive specific expression in the nervous system of zebrafish embryos, suggesting very deep functional conservation. Besides these spectacular archaisms, it seems nevertheless that turnover and innovation in cis-regulation have been so extensive in the lineages preceding vertebrates as to leave very little identifiable sequence homology, in contrast to the situation of coding sequences, which have been much more conserved and have allowed major advances in the evo-devo field along the past decades.

With the number of sequenced genomes expanding, it became possible to estimate the extent of constraint that natural selection has impressed on their evolving sequence. Hence, it appeared that while 40% of the human genome could be aligned with the mouse genome, only 5% was under purifying selection, i.e. had evolved more slowly than the neutral substitution rate [24]. While coding exons represented around one-third of these 5%, the remaining 3.5% (~100 million nucleotides) could be taken as an estimate of the amount of functional non-coding genome, most of which is thought to play a role in gene regulation. In support of this view, in a large-scale mouse transgenesis test of 31 000 ultra-CNEs, half of the candidates exhibited tissue-specific enhancer activity during mouse development [25]. Considering that this approach inherently misses enhancers active only in
adult, small numbers of cells or during very short developmental
timespans, as well as purely repressive elements, this suggests
indeed that a majority of CNEs function as regulatory elements.

Probing genomes for function

Large-scale projects were started to question these initial find-
ings by systematically probing genome-wide biochemical activ-
ities. The ENCODE consortium set out to annotate the human
genome by systematic molecular assays in various cell lines
and tissues. Among these assays, several were associated with
enhancer function, including DNase hypersensitivity. The initial
release of ENCODE, characterizing 1% of the human genome,
confirmed that 5% of the human genome was under positive
selection and reported in addition that only 60% of these con-
strained regions could be assigned a biochemical activity by any
of the assays [26]. These findings initiated a rupture between
the concepts of molecular function as defined by biochemical
assays and conservation as detected by sequence alignment:
one could occur without the other. Soon, the development of
next-generation sequencing (NGS) allowed the generalization of
genome-wide biochemical approaches, such as ChIP-seq, at an
unprecedented scale. Consequently, many studies came to rely
primarily on biochemical data as direct evidence of function,
often disregarding evolutionary conservation as a perhaps less
direct and less reliable proxy. Comforting this tendency, a com-
parative ChIP-seq study performed in the liver of five Amniota
species (mouse, human, dog, chicken and opossum) observed
that while the recognition motifs of the transcription factors
CEBP and HNF4 were the same between species, very few ChIP-
seq peaks for these factors could be aligned at orthologous
positions among the five genomes. Moreover, there seemed
to be no correlation between the computed affinity score of
the transcription factor binding motif found in a given peak
and the evolutionary conservation of that peak [27]. The weak
orthology between the position of putative enhancers, identified
by chromatin marks characteristic of active enhancers, among
20 mammals was further evidenced, at least for those active
in liver [28].

This departure from what we may tentatively call the
'Dobzhansky dogma' (i.e. 'Nothing in biology makes sense
except in the light of evolution') was somewhat substantiated
by the 2nd release of ENCODE, which used NGS techniques to
extend its initial analysis to the entire human genome [29]. Here,
biochemical annotations covered up to 80% of the total genome,
far surpassing the evolution-constrained fraction, of which
estimates range between 3% and 8%. Of course, a large part of
these 80%, half of which being broad chromatin state domains,
transcribed RNA of elusive function and decaying transposable
elements, may not correspond to what we think of as 'functional
elements'. However, even considering only transcription factor
binding elements and DNase hypersensitivity sites (DHS), these
made for ~20% of the total sequence, still much higher than the
largest estimates of constraint. Consistent with this, systematic
analysis of conservation of transcription factor binding sites
showed a weak signal, with sites bound by their cognate TF in
ChIP-seq experiments only slightly more conserved than sites
containing the same motif but with no ChIP-seq signal [29].

Other biochemical techniques were harnessed in massive
large-scale projects in order to identify regulatory regions. The
FANTOM consortium extensively characterized all RNA products
resulting from genomic transcription, including the so-called
eRNAs: transient, unprocessed transcripts that are produced by
RNA pol II at active enhancers [30]. Even though it is still not

known whether eRNAs have a function by themselves or if they
are contingent by-products of enhancer activity, they are now
considered as reliable markers of this activity and can therefore
be used to highlight enhancers [31]. The 2014 FANTOM release,
FANTOM5 [32], thus identifies 43 000 enhancers active in various
tissues and cell lines. Meanwhile, the Roadmap Epigenomics
project prolonged the analyses of ENCODE, this time with a focus
on tissues rather than on immortalized cell lines [33].

A common conclusion of these genome-wide collections of biochemical interactions is that an unexpectedly high fraction of the human genome is linked to at least one kind of annotation, and therefore candidate to a functional role, whereas a large part of these candidate functional elements do not show signs of purifying selection. This is illustrated in Figure 1, which shows that only 10–25% of the candidate regulatory elements or transcription factor binding sites annotated in the FOCS [34], ReMap [35] or UniBind [36] datasets exhibit conservation scores exceeding what can be observed for neutrally evolving elements such as ancient repeats, most of which are remnants of mobile elements with presumably no effect on organism fitness.

Resolution of this discrepancy is still subject to debate but
probably resides in a combination of explanations conciliating
two extreme standpoints [37]: i) that many functional sequences
may not be conserved ii) that many biochemical annotations
may not indicate biological functions.

Several non-exclusive reasons can explain why functional
regions may not be conserved at the nucleotide level. For
example, if a large number of regulations were species specific,
then the expectation of conservation upon evolution would be
clearly refuted. This would, however, conflict with the view
that most known biological mechanisms seem to be shared
among diverse species, sometimes up to large phylogenetic
depths. It seems unlikely that conservation be paramount
in the evolution of coding sequences yet simultaneously be
of so little importance for their regulation. Species-specific
functional regulatory sequences devoid of any measurable
evolutionary conservation cannot be ruled out but are likely
rare and restricted to precise biological functions prone to
accelerated evolution, for example those related to immunity or
fertilization. There are also functions that, while not restricted
to single species, are recent innovations of a specific lineage. In
these cases, neutral substitutions measured among the closely
related species may not be sufficiently numerous to leave a
signal of negative selection on the new regulatory sequences. In
fact, while current resources (computational tools and genomic
sequences) make it difficult to identify purifying selection
acting on recent innovations, these might on the contrary be
detected by signals of positive (adaptive) selection, indicating
that they evolved faster than the neutral rate. For example, it
has been shown that the recent evolution of the human lineage
has involved numerous cases of adaptations in regulatory
elements active in the brain [38]. Compared with purifying
selection, however, adaptive selection is much more difficult
to characterize, not only because it is much less pervasive,
but also because many unrelated genomic mechanisms (e.g.
GC-biased gene conversion) tend to produce confounding
signals [39].

It is also possible that regulatory elements are subject
to higher-order conservation rules, which escape current
nucleotide-based alignment methods. For example, their
function may reside on flexible arrangements of transcription
factor binding sites (TFBS) in modules [40, 41] or rely of DNA
secondary structure properties. Furthermore, TFBS themselves
may be too labile to be reliably detected by current alignment
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of mean phyloP scores for six sets of human putative genomic regulatory elements. The red area indicates negative (purifying) selection as the excess of conservation over expectation assuming neutral evolution. The total number of elements as well as the constrained fraction (in %) are indicated. A corpus of neutrally evolving sequences was compiled using ancestral repeats of the human genome, which present >30% divergence (retrieved at https://hgdownload.soe.u csc.edu/goldenPath/hg19/bigZips/hg19.fa.out.gz). The mean PhyloP score was computed on fixed-length windows (20, 50, 150, 400, 1000 or 2500 nucleotides) over these repeats. For each set of functional elements, a neutrally evolving control set was generated by picking repeat windows with length most closely matching the lengths of the elements. The distribution of mean PhyloP scores of these null sets is plotted in dashed lines. (A–D) FOCS enhancers defined from DHS (ENCODE, Roadmap), eRNA (FANTOM5) or global run-on (GRO-seq) assays [34]. (E) ChIP-seq peaks for various transcription factors, compiled in the ReMap 2018 database. Source data: http://remap.univ-amu.fr/storage/remap2018/hg19/MACS/remap2018_nr_macs2_hg19_v1_2-bed.gz). (F) Transcription factor binding sites defined by UniBind, based on the ReMap database combined with transcription factor binding matrices. Source data: https://unibind2018.uio.no/static/data/20190717/bulk_Permissive/Homo_sapiens/pwm_hg38_TFBS_per_TF.tar.gz. Binding sites positions were translated to the hg19 genome assembly using liftOver. For each transcription factor, multiple occurrences of identical sites were merged together so as to remove redundancy between assays.</figcaption></figure>

algorithms [42]. Several recent studies thus identified orthologous enhancers with no remaining sequence similarity, but with a conserved syntax of TFBS eliciting similar tissue-specific activity in distant species [43, 44].

A 3rd possible cause of low conservation of regulatory sequences may reside in their high redundancy. Indeed, if most gene regulation networks involve transcription factors binding multiple sites in multiple nearly redundant enhancers, as suggested both by experimental evidence and by theoretical work [45], then a significant portion of regulatory elements may experience rapid turnover in separate lineages, resulting in lineage-specific combinations of distinct enhancers regulating the target gene in each lineage [42]. Such widespread redundancy is supported by numerous cases where deletion of single putative enhancers does not alter the regulation of the target gene in vivo, whereas combined deletions of several such enhancers may have massive effects [46].

On the other hand, a large part of ENCODE annotations may correspond to 'molecular functions' (e.g. binding of a protein to DNA), which do not translate to phenotypes at the organism level, meaning they are not exposed to selection and therefore evolve neutrally. Some authors have argued that these would be more accurately described as evolutionary fortuitous biochemical properties rather than as genomic functions [47–50]. While this interpretation has to be true to a certain extent, simply because biochemical assays do capture such coincidental interactions, this extent remains extremely difficult to estimate as long as there is no simple test of high-level functionality, such as high-throughput genome editing coupled to systematic high-level functional assays.

Finally, an often-discussed cause for the gap between sequence conservation and regulatory function may lie in the reconciliation between approaches used to identify them. Thanks to deep multigenome alignments, statistical models
such as PhastCons [51] or GERP [19] can identify signals of
negative selection down to single nucleotides. In contrast,
functional assays are less precise and typically identify regions
spanning tens to hundreds of nucleotides where the truly
functional bases may represent only a small fraction. Finer-
grained methods are still required to tease them apart and would
likely resolve part of the current discordance.

Finding targets

Beyond the above discrepancy between function and conservation, a common challenge in many approaches to enhancer identification is to relate the putative enhancers to their cognate target genes. Neither conservation alone, nor molecular marks such as histone modification, open chromatin or production of eRNA provide direct insight into which genes are regulated by a candidate element. It is often assumed that the nearest gene is the most likely target, but there is no solid ground supporting that assumption. Indeed, while methodological biases have likely favoured the discovery of enhancers that are close to their targets, still a fair number of well-established regulatory elements are located at large distances, often closer to other, unrelated genes. Recent estimates indeed suggest that up to 50% of enhancers control a gene more distal than their closest promoter [32].

Several approaches have been used to determine which
genes are targeted by regulatory elements.

Physical gene–enhancer association

Chromosome conformation capture methods, 3C [52], measure
physical association between genomic regions. Coupled to NGS
techniques, Hi-C [53], they offer an opportunity for detecting reg-
ulatory interactions between genes and enhancers, since these
interactions involve the formation of large DNA loops that bring
the enhancer element, together with its bound chromatin/trans-
scription factors complex, in close contact to the gene promoter.
Capture Hi-C, an adaptation of Hi-C to genome-wide detection of
enhancer-promoter associations [54], identified over 1.6 million
interactions of genomic elements with 22 000 promoters in two
blood cell types.

Since Hi-C has no intrinsic limitation regarding the genomic
distance between interacting elements, the largest distances
observed between promoters and putative regulatory elements
can be taken as an upper limit for such interactions. For 99%
of 1.6 million long-range (>20 kb) interactions in Mifsud and
colleagues, the distance from the gene promoter is <1 Mb, sug-
gesting that most regulatory elements of the human genome lie
within 1 Mb from their target.

Physical proximity as detected by 3C technologies cannot
precisely delineate functional fragments within the loop-
forming regions: resolution of the putative regulatory elements
is limited to a few kilobases, which is often too large to conclude
about their conservation and constraint. Therefore, it is still
unknown how well-conserved the interacting sequences are.

Correlation between gene expression and enhancer activity

Since enhancers control expression of their target genes, one
may confront genome-wide gene expression data (e.g. measured
by RNA-seq) with enhancer activity data such as chromatin
accessibility (ATAC-seq, DHS) or characteristic histone marks,
the combination of which can be interpreted in chromatin states

[55], obtained jointly in large numbers of tissues/conditions.
Genes targeted by enhancers will be those whose expression
levels correlate convincingly with enhancer activity.

This has been done by the ENCODE consortium on DNase hypersensitivity data [56], by the FANTOM consortium on eRNA data [32] and by the Roadmap Epigenomics consortium on DNA accessibility data [33]. More recently, an effort has been made to improve the correlation methods and provide a unified procedure with a statistical validation scheme to produce high confidence inference of enhancer-gene regulations based on correlated activity patterns across many samples from heterogeneous sources, FOCS [34]. Compiling ENCODE, Roadmap Epigenomics, FANTOM and newly produced global run-on sequencing data (GRO-seq), FOCS thus identifies ~300 000 regulatory interactions affecting ~16 000 genes.

While this method is extremely powerful in that it addresses
directly the regulatory relationship between enhancers and
genes, it has several limitations, the main one being its depen-
dence on enormous amounts of experimental data, meaning
high costs, strong technical constraints and complex quality
control. It is also inherently limited to the study of cell-/tissue-
specific regulations, remaining blind to broader regulation types
such as those of housekeeping genes. Furthermore, is has been
so far restricted to the detection of positive regulations, leaving
aside the question of silencer elements, which might be of as
much interest as the better investigated enhancers.

So far, there seems to be little overlap between the regulatory
regions identified by these various approaches: for example,
only 11% of bases annotated as predicted enhancers by FOCS-
ENCODE intersect with FOCS-Roadmap Epigenomics data, 2% of
bases from the latter intersect with FOCS-FANTOM enhancers
and altogether, the intersection between ENCODE, FANTOM,
Roadmap Epigenomics and GRO-seq predicted enhancers in the
FOCS dataset covers only 17 000 bases, representing a mere 0.03%
of the combined annotated bases. This suggests either that they
are far from completion or that they comprise a large number
of false positives and perhaps both. Efforts to compile data from
various sources into unified collections aim at filling this gap,
such as GeneHancer [57], the original version of which reports
285 000 human enhancers out of which 94 000 derive from more
than one source.

Importantly, all the above methods based on experimen-
tal assays are of course dependent on the availability of tis-
sues and cell lines of interest and resulting interactions are
restricted to these biological contexts. In order to gain insight
into functional data independent on these technical contingen-
cies, one has to turn to the trace that functional interactions
have imprinted onto genomes along their evolution. Indeed,
while genomic functions associated to very scarce cell types
or very transient processes might remain beyond the reach of
large-scale experimental models, they still have to be some-
how encoded in genomes and this trace has to be shaped by
evolution.

Evolutionary conservation of synteny between
enhancer and target gene

It relies on the idea that if an element has an important func-
tion in regulating a target gene, then selection will constrain
not only the evolution of its sequence, but also of its physical
proximity to the target upon the rearrangement of chromosomes
accompanying speciation.

Application of this idea led to the concept of genomic regu-
latory blocks (GRBs). A GRB is a large genomic domain around a
gene, which includes arrays of ultraconserved elements, possibly interspersed with other genes. The term was coined by Lenhard and Becker in 2007 in a landmark study [58] where they exploited the recently completed zebrafish genome assembly to detect syntenies conserved across the large divergence time separating teleosts from humans. They noted that genes encoding developmental transcriptional regulators tended to be surrounded by larger regions of synteny than other functional categories of genes, suggesting that a selective constraint exerted to exclude breakpoints from these areas. They further correlated this constraint to the presence of ultraconserved elements, often closer to other genes.

When Irimia et al. [59] systematically identified pairs of unrelated genes with conserved physical proximity upon 600 M year evolution (Bilateria ancestor), they found a surprisingly high number of such cases (600), suggesting again that selection pressure acted to preserve microsyntenies. They estimated this pressure resulted from two main mechanisms. First, due to co-regulation of the two syntenic genes, presumably driven by common regulatory elements. Second, due to negative selection of breakpoints between the genes to preserve the integrity of putative GRBs.

These studies confirmed that gene regulatory interactions tend to prevent the fixation of chromosomal rearrangement breakpoints and highlighted microsynteny as a hallmark of such functionally important interactions.

With more and more assembled genomes, microsyntenies could be more systematically characterized. A two-step method was developed by Mongin, Dewar and Blanchette to detect putative regulatory relationships. First, whole-genome alignment of multiple species (UCSC Vertebrate 28-way) was screened to identify a large set of CNEs in the human genome. Then, human CNEs and neighbouring genes within 1 Mb radius were mapped to their orthologous counterparts in 16 selected genomes, chosen for the quality of their assembly and their evolutionary distance from human. Each orthologous CNE-gene pair was labelled either as associated, separated or incomplete in the case either the CNE or the gene could not be mapped. Ancestral association states were then inferred, and probability of functional association was computed for each pair based on models of synteny conservation/loss upon constraint [60].

In total, 168 000 human CNE-gene pairs with high probability of functional association were thus identified, involving nearly 7000 genes. More than half of these CNE-gene pairs were located >100 kb apart, and a large proportion of them jumped over genes. The properties of these associations, such as the correlation between complexity of the regulatory landscape (number of CNEs linked to a gene) on the one hand and gene function and promoter composition on the other hand, appeared consistent with known gene regulation features [60].

The Mongin method was expanded in the PEGASUS method, 1st developed in [61] where it was applied to the human X chromosome, then extended to the whole human genome and to the zebrafish genome by [62]. PEGASUS uses whole-genome alignments of a large number of species (46 vertebrate species in the 1st version) to delineate all CNEs. Then, for each CNE, the most likely target is determined by considering all genes in its vicinity (within 1 Mb for human, scaled according to genome size for other species) and score them according to the number of species in which this proximity is conserved. The scoring system favours ancient microsyntenies and explicitly penalizes cases where synteny has clearly been broken, with CNE and gene ending either on separate chromosomes or on very distant loci of the same chromosome.

High synteny scores are interpreted as high probability that the dual conservation (in sequence and in synteny) reflects conserved regulatory interaction between the element and its target gene.

The large number of species used for human PEGASUS allowed for an unprecedented number of CNE-gene associations to be called: >6 million in total for 1.3 million CNEs and 18 000 target genes. In parallel, applying the same method to the zebrafish genome, using whole-genome alignment of seven Actinopterygii fishes, yielded 55 000 CNEs associated to 17 000 genes. The large difference in the numbers underlies the importance of phylogenetic sampling: while human was compared with 34 other species at well-spaced evolutionary distances, zebrafish was compared with a smaller number of species (6), all of which had diverged from it >300 million years ago. When human PEGASUS analysis was reperformed on a subsampling of vertebrates mimicking the phylogenetic relationships of the zebrafish tree, the number of human CNEs detected dropped to 200 000, confirming that the density of the genome tree used for the analysis is a limiting step and highlighting the need to include as many good quality genomes as possible in such studies.

Currently used whole-genome alignments (e.g. the 100-way vertebrate Multiz alignment computed at UCSC; http://genome.u csc.edu/) probably still fall short of the optimal statistical power required to reliably align small conserved sequences, such as certain transcription factor binding sites, when they are isolated within large neutrally evolving regions. As a consequence, substantial amounts of constrained nucleotides may still remain undetected by comparative genomics approaches. Conversely, insufficient accumulation of neutral mutations along short time scales hampers proper identification of recently emerged functional sites and leads to higher rates of false positives when looking for lineage-specific regulatory elements.

Together with the sequencing of more genomes, especially in undersampled branches of the vertebrate tree, the development of new multiple genome alignment methods, such as the reference-free tool Progressive Cactus [63], should help improving both the false negative and false positive rates. In support of this perspective, a recent publication shows that genomic alignment of 240 eutherian mammals results in a doubling of the number of bases presumably under selection constraint in the human genome, compared with previous Multiz alignments involving three times less species [64].

In addition to predicting the target genes controlled by putative regulatory elements, ancient conservation of microsynteny between genes has also been used to identify conserved gene regulatory networks for which conservation of the enhancer sequence is too low to be detected by global alignment methods. Thus, a recent study identified enhancers in microsyntenic regions of the sponge Amphimedon queenslandica using ChIP-seq against histone marks characteristic of enhancers (high H3K4me1 and H3K27ac, low H3K4me3) [44]. Although BLAST queries could find no orthologous counterparts of these segments in vertebrate genomes, they were nevertheless able to activate expression of reporter constructs in both mouse and zebrafish along tissue-specific patterns corresponding to the genes involved in the microsynteny. Thus, combinations of comparative genomics with biochemical assays can help deciphering the primary syntax that governs the activity of regulatory elements.

Finally, scalable methods for directly assessing regulatory relationships between enhancers and target genes are beginning to emerge. Hints are given by eQTLs, collections of naturally
occurring variants associated with quantitative expression data [65]. These can detect subtle regulatory effects but are limited to sites that do vary in populations, therefore presumably missing the most essential sites, which do not vary. Moreover, eQTLs often come in groups of genetically linked variants, from which it is very difficult to identify the one that causes variation in gene expression. More recently, Crispr-based enhancer deletions and Crispr interference have been scaled up to near-genomic scales, targeting several thousands genomic elements and simultaneously measuring expression of whole transcriptomes in single cells in order to detect precise sites that alter gene expression when they are silenced [66]. In the future, further scaling of such techniques may become a powerful, albeit costly, tool for the discovery of new enhancer-target pairs and help progress towards establishing the genomic atlas of gene regulation.

## Conclusion

The last 15 years have seen a burst of new methods for probing the biochemical properties of chromatin at whole-genome scale. These technologies have been applied mostly to the functional exploration of the human genome, in various cell types, tissues, individuals or conditions. As a consequence, recent discoveries in functional genomics often tend to remain human centred or, at best, extended to a handful of model species, by lack of equivalent material resource for the vast majority of other species. Yet, sequencing technologies now permit assembly of complete genomes more rapidly than ever, with dozens of new species deciphered every month and massive projects undertaken in order to sequence a significant portion of biodiversity. These genomes are the repository of all important biological functions selected along hundreds of million years of evolution, a large part of which is played by the regulation of gene expression. It is thus likely that new understandings of common biological functions as well as the discovery of unsuspected lineage-specific innovations driven by regulatory changes will come from comparative genomics methods.

## Key Points

* Evolutionary sequence conservation has been used to successfully identify regulatory sequences and to estimate their contribution to the genome.

* Genome-wide assays that detect features linked to regulatory function have collectively uncovered millions of potential regulatory sequences.

* The overlap between these putative enhancers and sequences under evolutionary constraint remains poor.

* Possible explanations for this discrepancy are discussed.

## Acknowledgements

The authors wish to thank Camille Berthelot for insightful comments on the manuscript.

## Funding

This work has received support under the program «Investissements d'Avenir» launched by the French Government and implemented by ANR with the references ANR-10-LABX-54 MEMOLIFE and ANR-10-IDEX-0001-02 PSL* Université Paris.

## References

1. Banerji J, Rusconi S, Schaffner W. Expression of a beta-globin gene is enhanced by remote SV40 DNA sequences. *Cell* 1981;27:299-308.

2. Moreau P, Hen R, Wasylyk B, et al. The SV40 72 base repair repeat has a striking effect on gene expression both in SV40 and other chimeric recombinants. *Nucleic Acids Res* 1981;9:6047-68.

3. Gasperini M, Tome JM, Shendure J. Towards a comprehensive catalogue of validated and target-linked human enhancers. *Nat Rev Genet* 2020;21:292-310.

4. Wasserman WW, Palumbo M, Thompson W, et al. Human-mouse genome comparisons to locate regulatory sites. *Nat Genet* 2000;26:225-8.

5. Tagle DA, Koop BF, Goodman M, et al. Embryonic epsilon and gamma globin genes of a prosimian primate (*Galago crassicaudatus*). Nucleotide and amino acid sequences, developmental regulation and phylogenetic footprints. *J Mol Biol* 1988;203:439-55.

6. Hardison R, Miller W. Use of long sequence alignments to study the evolution and regulation of mammalian globin gene clusters. *Mol Biol Evol* 1993;10:73-102.

7. Moon AM, Ley TJ. Conservation of the primary structure, organization, and function of the human and mouse beta-globin locus-activating regions. *Proc Natl Acad Sci U S A* 1990;87:7693-7.

8. Nobrega MA. Scanning human gene deserts for long-range enhancers. *Science* 2003;302:413-3.

9. Goode DK, Snell P, Smith SF, et al. Highly conserved regulatory elements around the SHH gene may contribute to the maintenance of conserved synteny across human chromosome 7q36.3. *Genomics* 2005;86:172-81.

10. Poulin F, Nobrega MA, Plajzer-Frick I, et al. In vivo characterization of a vertebrate ultraconserved enhancer. *Genomics* 2005;85:774-81.

11. la Calle-Mustienes d E, Feijóo CG, Manzanares M, et al. A functional survey of the enhancer activity of conserved non-coding sequences from vertebrate Iroquois cluster gene deserts. *Genome Res* 2005;15:1061-72.

12. Sandelin A, Bailey P, Bruce S, et al. Arrays of ultraconserved non-coding regions span the loci of key developmental genes in vertebrate genomes. *BMC Genomics* 2004;5:99.

13. Venkatesh B, Kirkness EF, Loh Y-H, et al. Ancient non-coding elements conserved in the human genome. *Science* 2006;314:1892.

14. Müller F, Blader P, Strähle U. Search for enhancers: teleost models in comparative genomic and transgenic analysis of cis regulatory elements. *Bioessays* 2002;24:564-72.

15. Lee AP, Kerk SY, Tan YY, et al. Ancient vertebrate conserved noncoding elements have been evolving rapidly in teleost fishes. *Mol Biol Evol* 2011;28:1205-15.

16. Bejerano G, Pheasant M, Makunin I, et al. Ultraconserved elements in the human genome. *Science* 2004;304:1321-5.

17. Mikkelsen TS, Wakefield MJ, Aken B, et al. Genome of the marsupial *Monodelphis domestica* reveals innovation in non-coding sequences. *Nature* 2007;447:167-77.
18. Prabhakar S, Poulin F, Shoukry M, et al. Close sequence comparisons are sufficient to identify human cis-regulatory elements. *Genome Res* 2006;16:855-63.

19. Cooper GM, Stone EA, Asimenos G, et al. Distribution and intensity of constraint in mammalian genomic sequence. *Genome Res* 2005;15:901-13.

20. Yang H, Luan Y, Liu T, et al. A map of cis-regulatory elements and 3D genome structures in zebrafish. *Nature* 2020;588:337-43.

21. Putnam NH, Butts T, Ferrier DEK, et al. The amphioxus genome and the evolution of the chordate karyotype. *Nature* 2008;453:1064-71.

22. Marlétaz F, Firbas PN, Maeso I, et al. Amphioxus functional genomics and the origins of vertebrate gene regulation. *Nature* 2018;564:64-70.

23. Clarke SL, VanderMeer JE, Wenger AM, et al. Human developmental enhancers conserved between deuterostomes and protostomes. *PLoS Genet* 2012;8:e1002852.

24. Mouse Genome Sequencing Consortium, Waterston RH, Lindblad-Toh K, et al. Initial sequencing and comparative analysis of the mouse genome. *Nature* 2002;420:520-62.

25. Pennacchio LA, Ahituv N, Moses AM, et al. In vivo enhancer analysis of human conserved non-coding sequences. *Nature* 2006;444:499-502.

26. ENCODE Project Consortium, Birney E, Stamatoyannopoulos JA, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. *Nature* 2007;447:799-816.

27. Schmidt D, Wilson MD, Ballester B, et al. Five-vertebrate ChIP-seq reveals the evolutionary dynamics of transcription factor binding. *Science* 2010;328:1036-40.

28. Villar D, Berthelot C, Aldridge S, et al. Enhancer evolution across 20 mammalian species. *Cell* 2015;160:554-66.

29. ENCODE Project Consortium, Bernstein BE, Birney E, et al. An integrated encyclopedia of DNA elements in the human genome. *Nature* 2012;489:57-74.

30. Kim T-K, Hemberg M, Gray JM, et al. Widespread transcription at neuronal activity-regulated enhancers. *Nature* 2010;465:182-7.

31. Natoli G, Andrau J-C. Noncoding transcription at enhancers: general principles and functional models. *Annu Rev Genet* 2012;46:1-19.

32. Andersson R, Gebhard C, Miguel-Escalada I, et al. An atlas of active enhancers across human cell types and tissues. *Nature* 2014;507:455-61.

33. Roadmap Epigenomics Consortium, Kundaje A, Meuleman W, et al. Integrative analysis of 111 reference human epigenomes. *Nature* 2015;518:317-30.

34. Hait TA, Amar D, Shamir R, et al. FOCS: a novel method for analyzing enhancer and gene activity patterns infers an extensive enhancer-promoter map. *Genome Biol* 2018;19:56.

35. Chêneby J, Gheorghe M, Artufel M, et al. ReMap 2018: an updated atlas of regulatory regions from an integrative analysis of DNA-binding ChIP-seq experiments. *Nucleic Acids Res* 2018;46:D267-75.

36. Gheorghe M, Sandve GK, Khan A, et al. A map of direct TF-DNA interactions in the human genome. *Nucleic Acids Res* 2019;47:e21.

37. Kellis M, Wold B, Snyder MP, et al. Defining functional DNA elements in the human genome. *Proc Natl Acad Sci* 2014;111:6131-8.

38. Liu J, Robinson-Rechavi M. Robust inference of positive selection on regulatory sequences in the human brain. *Sci Adv* 2020;6:eabc9863.

39. Galtier N, Duret L. Adaptation or biased gene conversion? Extending the null hypothesis of molecular evolution. *Trends Genet* 2007;23:273-7.

40. Romano LA, Wray GA. Conservation of Endo16 expression in sea urchins despite evolutionary divergence in both cis and trans-acting components of transcriptional regulation. *Development* 2003;130:4187-99.

41. Hare EE, Peterson BK, Iyer VN, et al. Sepsid even-skipped enhancers are functionally conserved in Drosophila despite lack of sequence conservation. *PLoS Genet* 2008;4:e1000106.

42. Dermitzakis ET, Clark AG. Evolution of transcription factor binding sites in mammalian gene regulatory regions: conservation and turnover. *Mol Biol Evol* 2002;19:1114-21.

43. Minnoye L, Taskiran II, Mauduit D, et al. Cross-species analysis of enhancer logic using deep learning. *Genome Res* 2020;30:1815-34.

44. Wong ES, Zheng D, Tan SZ, et al. Deep conservation of the enhancer regulatory code in animals. *Science* 2020;370:eaax8137.

45. Paixão T, Azevedo RBR. Redundancy and the evolution of cis-regulatory element multiplicity. *PLoS Comput Biol* 2010;6:e1000848.

46. Osterwalder M, Barozzi I, Tissières V, et al. Enhancer redundancy provides phenotypic robustness in mammalian development. *Nature* 2018;554:239-43.

47. Cooper GM, Brown CD. Qualifying the relationship between sequence conservation and molecular function. *Genome Res* 2008;18:201-5.

48. Graur D, Zheng Y, Price N, et al. On the immortality of television sets: 'function' in the human genome according to the evolution-free gospel of ENCODE. *Genome Biol Evol* 2013;5:578-90.

49. Doolittle WF. Is junk DNA bunk? A critique of ENCODE. *Proc Natl Acad Sci U S A* 2013;110:5294-300.

50. Linquist S, Doolittle WF, Palazzo AF. Getting clear about the F-word in genomics. *PLoS Genet* 2020;16:e1008702.

51. Pollard KS, Hubisz MJ, Rosenbloom KR, et al. Detection of nonneutral substitution rates on mammalian phylogenies. *Genome Res* 2010;20:110-21.

52. Dekker J, Rippe K, Dekker M, et al. Capturing chromosome conformation. *Science* 2002;295:1306-11.

53. Lieberman-Aiden E, van Berkum NL, Williams L, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. *Science* 2009;326:289-93.

54. Mifsud B, Tavares-Cadete F, Young AN, et al. Mapping long-range promoter contacts in human cells with high-resolution capture Hi-C. *Nat Genet* 2015;47:598-606.

55. Ernst J, Kellis M. Discovery and characterization of chromatin states for systematic annotation of the human genome. *Nat Biotechnol* 2010;28:817-25.

56. Thurman RE, Rynes E, Humbert R, et al. The accessible chromatin landscape of the human genome. *Nature* 2012;489:75-82.

57. Fishilevich S, Nudel R, Rappaport N, et al. GeneHancer: genome-wide integration of enhancers and target genes in GeneCards. *Database* (Oxford) 2017;2017:bax028.

58. Kikuta H, Laplante M, Navratilova P, et al. Genomic regulatory blocks encompass multiple neighboring genes and maintain conserved synteny in vertebrates. *Genome Res* 2007;17:545-55.

59. Irimia M, Tena JJ, Alexis MS, et al. Extensive conservation of ancient microsynteny across metazoans due
to cis-regulatory constraints. *Genome Res* 2012;**22**: 2356-67.

60. Mongin E, Dewar K, Blanchette M. Mapping association between long-range cis-regulatory regions and their target genes using synteny. *J Comput Biol* 2011;**18**:1115-30.

61. Naville M, Ishibashi M, Ferg M, et al. Long-range evolutionary constraints reveal cis-regulatory interactions on the human X chromosome. *Nat Commun* 2015;**6**:6904.

62. Clément Y, Torbey P, Gilardi-Hebenstreit P, et al. Enhancer-gene maps in the human and zebrafish genomes using evolutionary linkage conservation. *Nucleic Acids Res* 2020;**48**:2357-71.

63. Paten B, Earl D, Nguyen N, et al. Cactus: algorithms for genome multiple sequence alignment. *Genome Res* 2011;**21**:1512-28.

64. Zoonomia Consortium. A comparative genomics multitool for scientific discovery and conservation. *Nature* 2020;**587**:240-5.

65. GTEx Consortium. The GTEx consortium atlas of genetic regulatory effects across human tissues. *Science* 2020;**369**:1318-30.

66. Gasperini M, Hill AJ, McFaline-Figueroa JL, et al. A genome-wide framework for mapping gene regulation via cellular genetic screens. *Cell* 2019;**176**:1516.