How to find genomic regions relevant for gene regulation

Introduction

Completion of the human genome sequence [1], which contains roughly 3.3 billion nucleotides, quickly led to the recognition that only about 2 % of the human genome is protein coding, while the remaining nearly 98 % do not code for proteins. Functional high-throughput studies have indicated that a considerable fraction of non-coding sequences, with estimates on the order of 5–10 %, harbor key functional elements responsible for the regulation of complex temporal and tissue-specific gene expression in different cell types of the human body [2]. Genome-wide association studies (GWAS) have mapped numerous loci that are associated with complex phenotypic traits. More than 90 % of these GWAS variants are located in non-coding regions of the genome clustering in and around regulatory elements. How these variants act on the phenotype has thus fueled new interest into the functional interaction between these quantitative trait loci (QTLs) and their associated genes [3].

With the genome sequence available, the next pressing challenges were to annotate all gene regulatory regions and understand their function. Due to the rapid progress in DNA sequencing technologies over the course of the last decade, combined with the development of novel in silico analysis and modeling tools, our understanding of the non-coding human genome has finally advanced to the point where we can begin to address these challenges.

Figure 1

Model of promoter–enhancer interaction in a chromatin loop extruded by the cohesin complex. (A and B) The ring-shaped ATPase cohesin complex initiates the extrusion of the chromatin fiber and keeps translocating along the chromatin, forming a DNA loop in the process. (C) Translocating cohesin stops at a pair of convergent CTCF sites that form the boundary of the chromatin loop or the TAD. The CTCF sites have to be in a proper orientation to each other in order to stop cohesin. (D) Enhancers and promoters within a single loop are prone to interaction due to proximity, resulting in the recruitment of transcription factor(s) and of RNA polymerase II, thereby initiating gene transcription.

Transcriptional regulatory regions, also referred to as cis-regulatory regions or modules, are often grouped into promoters, enhancers, and insulators. Promoters are proximal regions around the transcription start sites (TSSs) of genes which contain DNA elements that recruit the basal transcriptional factors (TFs) and, ultimately, RNA Polymerase II (RNAP II), setting in motion the transcription machinery. Enhancers are distal cis-regulatory regions that recruit sequence-specific TFs along with chromatin remodeling co-factors, thereby exerting an influence on the transcription initiation at the promotor. A prominent feature of enhancers is their key role in tissue-specific gene expression. Highly regulated genes can be influenced by up to several dozen enhancers [4], [5]. Their coordinated activity with promoters dictates the dynamic temporal and spatial regulation of transcriptional programs in different cell types.

To this end, proximal and distal regions contain several short sequence elements called TF binding sites (TFBSs) that attract specific TFs (see Leitz et al., this issue). Distal regions with an effect on gene transcription can be located downstream or upstream of the promoter, as well as inside the gene or at a large distance ranging from hundreds to tens of thousands of nucleotides. Finally, insulators act as barriers to prevent the spread of heterochromatin and to block distal regions from acting on unrelated promoters. In vertebrate cells, this enhancer-restrictive function of insulators is mediated by binding of the CCCTC-binding factor (CTCF), which plays a crucial role in organizing the genome into topologically associating domains (TADs) [6]. Currently, the generation of TADs is best explained by the loop extrusion model (Figure 1A–C), where the ring-shaped cohesin complex pinches a loop out of the chromatin fiber and continues to propel chromatin via its ATPase activity until it reaches a pair of convergent CTCF sites. Promoters and enhancers within a loop are prone to interaction due to proximity (Figure 1D). Notably, while enhancers are the primary drivers for cell type-specific gene expression programs, insulator positions have largely been found to be similar in different cell types.

Regulatory regions of the genome are marked by epigenetic signals, which comprise DNA methylation and chemical modifications of the histones. Eukaryotic DNA is tightly packed into chromatin by means of nucleosomes, consisting of 147 bp of DNA wrapped around a histone octamer, separated by stretches of linker DNA. Modification of these histones by means of methylation or acetylation has an influence on the packing density and accessibility of the chromatin. Unlike the underlying static DNA sequence, these modifications can dynamically respond to and alter the cellular states. Specific epigenetic modifications and their particular combinations, the so-called “histone code” [7], are associated with distinct genomic regulatory features and their activity states [8], providing key spatial and temporal functional information.

In this review, we focus on genome-wide approaches for defining genomic regulatory regions, including high-throughput profiling of open chromatin, epigenomic marks, and long-range promoter–enhancer interactions. We also discuss advances in the genome and epigenome editing assays for functional analyses of the identified candidate regulatory regions, enabling their validation for potential causal roles in target gene regulation.

Annotation of candidate regions

Activity of candidate regions

DNase- and ATAC-seq assays reveal accessible genomic loci, but they do not distinguish between different kinds of regulatory regions and their activity, i. e., if they are engaged in regulating a gene in a particular context. The regulatory state of a region is however reflected in patterns of histone modifications. DNA is wrapped around a nucleosome core, which is a tetramer composed of different histone proteins with tails that can be modified, most frequently by methylation or acetylation. These modifications serve as road signs for gene regulation. Using Chromatin ImmunoPrecipitation (ChIP), Hebbes et al. (1988) first established a direct link between core histone acetylation and transcriptionally active chromatin [17]. As more and more types of histone modifications were discovered, scientists put forward the existence of a “histone code,” in contrast to the “genetic code,” which orchestrates the transcriptional program of the invariant genome. Table 1 lists the best-known and most robust histone marks. Due to the specificity of antibody–antigen recognition, ChIP has been successfully adapted to the genome-wide profiling of histone modifications, histone variants, and DNA methylation. This method is called ChIP sequencing (ChIP-seq), where all the DNA fragments are sequenced that are precipitated with the protein of interest, e. g., a transcription factor or a modified histone. This allows bioinformatics to remap the position in the genome from where the respective chromatin fragments were precipitated. Figure 2 shows a schematic genome browser view of the human MYOD1 locus with RNA sequencing (RNA-seq) tracks of plus and minus strands and ChIP-seq tracks of CTCF, H3K4me, H3K4me3, H3K27ac, H3K27me3, and H3K36me3 generated by the ENCODE Consortium from a human myotube culture derived from a skeletal muscle myoblast line.

Table 1

Examples for frequent histone marks. The table lists the names of five common histone marks, their preferred genomic locations, and their relations with respect to gene transcriptional activity; e. g., H3K4me3 = trimethylation at the fourth lysine residue of the histone H3 protein; H3K27ac = acetylation at the 27th lysine residue of the histone H3 protein.

Histone mark	Location	Activity
H3K4me3	around the TSS	active
H3K4me1/me2	at the enhancer	active or primed
H3K27ac	at the enhancer	active
H3K36me3	within the gene body	active
H3K27me3	at the promoter and within the gene body	repressive

Figure 2

Schematic depiction of the browser view of the human MYOD1 locus, overlaid with RNA-seq, CTCF ChIP-seq, and various histone ChIP-seq tracks. Sequencing data were generated from human myotubes in culture by the ENCODE Consortium. For the gene annotation track at the top, blue indicates protein coding genes, thick bars protein coding regions (exons), narrow bars untranslated regions (UTRs), lines introns, and arrows the direction of transcription. RNA-seq signals of plus and minus strands are shown, as are the most common and robust active and repressive histone marks. The positions of gene regulatory regions are depicted by vertical shades: promoters (blue), enhancers (green), TAD boundaries (yellow), actively transcribed gene bodies (purple), and repressed gene bodies (gray). In this example the gene MYOD1 is transcribed, while KCNC1 is repressed.

DNA methylation has long been associated with gene repression, but research over the course of the past decades has shown that DNA methylation participates in multiple cellular functions, some of which are still not fully understood. These include repression of transposons, inactivation of the X-chromosome, and genomic imprinting. Most DNA methylation happens at CpG dinucleotides. This short palindromic sequence instructs the methylation of daughter strands during DNA replication.

Methylation is classically investigated by bisulfite sequencing [18]. Treatment with bisulfite before sequencing converts cytosine residues to uracil with the exception of 5-methyl-cytosines. Therefore, only methylated cytosines are retained. The computational comparison of sequence reads from treated and untreated DNA samples allows the exact determination of CpG positions at single-nucleotide resolution. NGS allows genome-wide analysis of methylation patterns, a method called Bis-seq [19].

The low overall percentage of CpG dinucleotides in the mammalian genome stands in contrast to the high percentage found at CpG islands (CGIs), where CpG dinucleotides cluster in genomic regions of about 1 kb. Over two thirds of mammalian promoters, including almost all housekeeping genes and several developmental genes, coincide with CGIs. CGI promoters are rarely methylated; in these instances, gene silencing is achieved by H3K27 methylation [20]. While methylation of some TFBSs deters binding, many TFs in the extended homeodomain family did prefer methylated CpG in an in vitro assay [21]. Interestingly, the mammalian genome rewrites its methylation pattern twice during development, once after fertilization and once after germline specification.

Figure 3

Chromatin conformation capture methods profile long-range interactions between regulatory sequences, thereby identifying TADs. (A) TADs can be repressed by binding of chromatin remodelers such as Polycomb Repressive Complex 2 (PRC2), which compacts the regions and prevents transcription. (B) Depicted are two TADs, one active (green) and one repressed (red). In the active TAD, a long-range promoter–enhancer interaction leads to active transcription. In the repressed TAD, silencing histone marks such as H3K27me3 recruit epigenetic modifiers and chromatin remodelers such as PRC2 to compact chromatin, preventing transcription. TADs are identified by directional ligation events in the 3C-based assays.

Through ENCODE and other concerted efforts to systematically profile the epigenomic landscapes in human cell lines, primary cell types, and model systems, it became obvious that the number of possible combinations of histone combinations and other epigenetic states at genomic loci exceeded the capabilities of manual human exploration. This is therefore one area where computational analyses, and in particular the use of machine learning, has made a crucial impact. In 2010, Ernst and Kellis identified 51 “chromatin states” de novo with a simple combinatorial model of 38 histone modifications and other epigenetic information. These states fit with existing genomic annotations [8], although subsequent publications have proposed a less complex picture with fewer states [22], [23].

Validating candidate variants

The “genetic code” for translation of nucleotide triplets into amino acids has been known for decades. While the impact of mutations on the protein coding function can be linked to large numbers of clinical disease phenotypes, there appears to be no straightforward “regulation code” counterpart. Many enhancers can influence one gene, one enhancer can affect several genes, and their joint activity is encoded in the combination of functional sequence elements, their redundancy, and their relationship with each other. To understand what the function of a particular distal region might be, genetics has long relied on reporter assays. Applicable mostly in model systems and cell lines, one candidate variant at a time is placed next to a reporter gene that enables the quantification of its influence on the expression of the reporter, such as a fluorescent protein. While allowing for in vivo insights, interpretation of reporters is limited in a number of ways: candidate variants (i) are typically tested outside of their native sequence and chromatin context, (ii) may comprise one, or only a part of a larger, regulatory region; and (iii) in the case of plasmid reporter constructs, are studied in the absence of chromatin.

To address the limited throughput, recent developments include massively parallel reporter assays that streamline and scale up the cloning of functional candidate fragments next to the reporter gene, allowing for readouts of thousands of fragments [32]. Placing short sequence fragments into the same controlled context holds promise for eliminating unwanted noise, thereby facilitating an unbiased evaluation of the impact of sequence variants. One example is Self-Transcribing Active Regulatory Region sequencing (STARR-seq). This method is based on the knowledge that enhancers can work independently of their relative locations. Placing an enhancer candidate sequence downstream of a minimal promoter enables active enhancers to transcribe themselves, meaning that it can be read out and quantified by RNA-seq [33]. In spite of the limitations of these reporter assays, e. g., the differences observed when placing the same sequence within native chromatin or on plasmids [34], such data have laid the foundation for promising computational models that identify non-coding variants of clinical relevance [35].

To dissect the impact of a regulatory region within its native genomic context, one of the most promising directions has been opened by CRISPR genome editing tools. As is the case for protein coding genes, CRISPR editing can be targeted to specifically alter regulatory regions by means of complementary guide RNAs. As we often do not know which individual sequence features are relevant, saturation screens that introduce thousands of mutations help to comprehensively assay a larger genomic region, such as an entire TAD [36]. In pooled screens, the effect of a sequence or chromatin change is detected by phenotypes such as growth or proliferation [37]. In single-cell screens, the effect on gene expression is directly determined via single-cell RNA-seq [38].

Alternatively, “epigenome editing” perturbs the chromatin state via the introduction of inhibitory or activating marks (CRISPRi/a) [39]. In this approach, a dead Cas9 nuclease is fused to epigenome modifier domains, such as histone (de)acetylases, and is directed to a specific genomic location via a guide RNA [34]. As the whole regulatory region is affected instead of a small sequence, data acquisition scales better and interpretation is easier [40], [41].

Single-cell genomics

With the requirements on sample size and purity decreasing, large consortia have made progress, moving from profiling of immortalized cell lines (ENCODE [11]) to post mortem fetal tissues (Roadmap [12]) and post mortem adult samples of healthy organs (GTEx consortium [42]). Genomics protocols have traditionally required tens of thousands to millions of cells to generate high-quality data. This meant that small cell populations, or short-lived states during development, were either not well reflected or entirely obscured in the averaged pictures resulting from these approaches. This situation has changed rapidly with the rise of single-cell genomics where sequencing protocols on individual cells [43] using microfluidic devices, liquid handling robotics, and/or clever multiplexed experimental designs enable us to profile tens of thousands and even millions of cells in a single experiment [44]. In contrast to the traditional, single average dataset generated from millions of cells and encompassing millions of short sequencing reads, each single-cell dataset is composed of a limited number of reads in the range of 10,000 to 100,000 reads. However, obtaining tens of thousands single-cell datasets from one experiment enables us to uncover heterogeneous states and responses of neighboring cells in a tissue [45].

Due to the highly parallel nature of these experiments, protocols that could be implemented in a few steps were adopted first, starting with RNA-seq for expression profiling. Yet other protocols, such as single-cell ATAC-seq, quickly followed and are now available as standardized kits from commercial vendors. Meanwhile, large atlases of gene expression and open chromatin in multiple systems have been generated for multiple organs [46], [47], [48]. Bisulfite sequencing and ChIP-seq are also implemented, but are not yet as widely adopted [44], [45]. With the possibility of obtaining high-resolution information from small, complex, and primary samples, single-cell genomics will transform basic genomics research and its clinical applications.

Since the start of genome sequencing projects, computational biology algorithms have been indispensable for organizing and interpreting today’s massive heterogeneous sources of data from bulk and single-cell experiments. In particular, carefully designed and vetted machine learning methods hold great promise for interpreting the complex rules of gene regulation that are spread out over large genomic regions and involve multiple mechanisms [35], [49]. As with other applications, a critical aspect of successful adaptation of this technique will be the ability to explain physiology and provide meaningful interpretations.

Conclusion

High-throughput biology, and in particular deep sequencing, has made exploring the vast space of non-protein coding genomic regions possible. We now have a large catalog of candidate regulatory regions at our disposal, building the foundation that will enable us to include non-coding sequence variation into diagnostics and clinical practice in the context of the underlying causes of rare diseases.