Prediction and Validation of Gene Regulatory Elements Activated During Retinoic Acid Induced Embryonic Stem Cell Differentiation

Summary

In this work we provide an experimental workflow of how active enhancers can be identified and experimentally validated.

Abstract

Embryonic development is a multistep process involving activation and repression of many genes. Enhancer elements in the genome are known to contribute to tissue and cell-type specific regulation of gene expression during the cellular differentiation. Thus, their identification and further investigation is important in order to understand how cell fate is determined. Integration of gene expression data (e.g., microarray or RNA-seq) and results of chromatin immunoprecipitation (ChIP)-based genome-wide studies (ChIP-seq) allows large-scale identification of these regulatory regions. However, functional validation of cell-type specific enhancers requires further in vitro and in vivo experimental procedures. Here we describe how active enhancers can be identified and validated experimentally. This protocol provides a step-by-step workflow that includes: 1) identification of regulatory regions by ChIP-seq data analysis, 2) cloning and experimental validation of putative regulatory potential of the identified genomic sequences in a reporter assay, and 3) determination of enhancer activity in vivo by measuring enhancer RNA transcript level. The presented protocol is detailed enough to help anyone to set up this workflow in the lab. Importantly, the protocol can be easily adapted to and used in any cellular model system.

Introduction

Development of a multicellular organism requires precisely regulated expression of thousands of genes across developing tissues. Regulation of gene expression is accomplished in large part by enhancers. Enhancers are short non-coding DNA elements that can be bound with transcription factors (TFs) and act from a distance to activate transcription of a target gene¹. Enhancers are generally cis-acting and most frequently found just upstream of the transcription start site (TSS), but recent studies also described examples where enhancers were found much further upstream, on the 3' of the gene or even within the introns and exons².

There are hundreds of thousands of potential enhancers in the vertebrate genomes¹. Recent methods based on chromatin immunoprecipitation (ChIP) provide high-throughput data of the whole genome that can be used for enhancer analysis^3-9. Though data obtained by ChIP-seq experiments greatly increases the likelihood to identify cell and tissue-specific enhancers, it is important to keep in mind that detected binding sites do not necessarily identify direct DNA binding and/or functional enhancers. Thus, further functional analysis of newly identified enhancers is indispensable. In this work, we present a basic three-step process of putative active enhancer identification and validation. This includes: 1) selection of putative transcription factor binding sites by bioinformatics analysis of ChIP-seq data, 2) cloning and validation of these regulatory sequences in reporter constructs, and 3) measurement of enhancer RNA (eRNA).

Exposure of embryonic stem (ES) cells to retinoic acid (RA) is frequently used to promote neural differentiation of the pluripotent cells ¹⁰. RA exerts its effects by binding to RA receptors (RARα, β, γ) and retinoid X receptors (RXRα, β, γ). RARs and RXRs in a form of heterodimer bind to DNA motifs called RA-response elements, that is typically arranged as direct repeats of AGGTCA sequence (called as half site) and regulate transcription. Ligand-treatment experiments allowed the identification of several retinoic acid regulated genes in ES cells ^11,12. However, enhancer elements for many of these genes has not been described yet. To demonstrate how the here-described workflow can be used for enhancer identification and validation we show step-by-step the selection and characterization of two retinoic acid-dependent enhancers in embryonic stem cells.

Protocol

1. Enhancer Selection Based on Chip-seq Analysis

1. Download the RXR ChIP-seq raw data fastq file (mm_ES_RXR_24h_ATRA.fastq.gz) from http://ngsdebftp.med.unideb.hu/bioinformatics/
2. Download and extract the required BWA index file for the alignment (in our case: Mus_musculus_UCSC_mm10).(ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Mus_musculus/UCSC/mm10/Mus_musculus_UCSC_mm10.tar.gz
   NOTE: Visit to https://github.com/ahorvath/Bioinformatics_scripts for more information regarding the steps of bioinformatics analysis and to download the scripts used below.
3. Align the example fastq file to mm10 genome (use the script: perform_alignment.sh). This will create a folder with a .bam file and statistics of the alignment. Aligned RXR ChIP-seq data (mm_ES_RXR_24h_ATRA.bam), and the corresponding index file (.bai) are available here:
   http://ngsdebftp.med.unideb.hu/bioinformatics/
   NOTE: BWA is a software package that aligns relatively short sequences (e.g., ChIP-seq results) to a sequence database, such as the mouse reference genome (e.g., mm10) ¹³.
4. Run the script (callpeaks.sh) for peak calling and de novo motif analysis. Use the .bam file as the input. The script is based on Homer findPeaks ¹⁴.
   NOTE: The output files give us information about the total number of peaks, enrichment of the motifs, the percentage of the background, and the target sequences with motifs, etc. (Figure 1). Typically several motifs are listed in order of their significance.
5. Remap the #1 ranked motif (NR half `AGGTCA`) (use the script: remap_motif.sh)
   NOTE: As the result, the script generates a .bed file that will show which ChIP-peaks are covering genomic regions that contain the canonical binding site of the transcription factor of interest.
6. Download and visualize results of the aligned RXR ChIP-seq reads (mm_ES_RXR_24h_ATRA.bam (also download the .bai)) and the AGGTCA motif occurrences (mm_ES_RXR_24h_ATRA_homerpeaks_motif1_mm10s_200_remaped_mbed.bed) by Integrative Genomics Viewer (IGV)¹⁵ to identify putative RAR/RXR binding sites (Figure 2 and Figure 3).
   NOTE: Integration of various ChIP-seq data will help to more precisely predict good candidate enhancers^16,17. Recent studies revealed that active enhancers are typically enriched for P300 and correlate with H3K4me1 and H3K27ac, and are located in open chromatin regions, thereby displaying DNase-I hypersensitivity is also recommended ^18-21.
7. Use the "Define a region of interest" in IGV ¹⁵ to obtain 200 - 400 bp sequence of the selected genomic region and use these sequences for subsequent primer designs.

2. Reporter Assay

NOTE: The luciferin/luciferase system is used as a very sensitive reporter assay for transcriptional regulation. Depending on the enhancer activity luciferase enzyme is produced that will catalyze the oxidation of luciferin to oxyluciferin resulting in bioluminescense, which can be detected. As a first step, identified putative enhancer sequences should be subcloned into a reporter vector (e.g., TK-Luciferase ²², pGL3 or NanoLuc).

1. Primer Design
   1. Design PCR primers for amplification of 200 - 400 bp of putative enhancer regions by Primer3 (available here: http://bioinfo.ut.ee/primer3/)²³
   2. Copy-paste the genomic sequence from IGV that includes the putative binding site into Primer3 (see step 1.7). Set product size range to 250 - 300 bp in Primer3.
   3. Validate the PCR primers using the UCSC Genome Browser's In silico PCR tool (https://genome.ucsc.edu/cgi-bin/hgPcr)²⁴
      NOTE: This protocol in a later step uses HindIII and BamHI restriction enzymes for cloning. Check whether the PCR amplified sequence as predicted by the UCSC In-silico PCR tool will contain AAGCTT or GGATCC restriction sites (e.g.,  http://nc2.neb.com/NEBcutter2/).
   4. Obtain the required primers from a commercial vendor.
2. PCR Amplification of Enhancer Sequence from Genomic DNA
   1. Purify template genomic DNA from primary cells (e.g., primary embryonic fibroblasts). Use commercial kit for gDNA isolation and follow the manufacturer's instructions.
      NOTE: High quality of gDNA is essential for successful PCR amplification. Appropriate BAC clones can be used if PCR is not easy from gDNA.
   2. Prepare 50 µl PCR reaction containing 100 ng gDNA, 5 µl 10x buffer, 3 µl MgSO₄ (25 mM), 5 µl dNTP (2 mM each), 1.5 µl Fwd primer (10 µM), 1.5 µl Rev primer (10 µM), 1 µl DNA polymerase
      NOTE: Use high-fidelity DNA polymerase with a low error rate for the PCR reactions.
   3. Set the PCR cycle (2 min 95 °C, (20 sec 95 °C, 10 sec 58 - 65 °C, 18 sec 70 °C) x25 repeat, 2 min 70 °C, forever 4 °C). Set annealing temperature with the consideration of the primer's predicted Tm values.
   4. Purify the PCR product with a commercial PCR purification kit and follow the manufacturer's instructions.
   5. Introduce the restriction sites for cloning by repeating the PCR with primers as used above, but adding overhangs on their 5' ends (e.g., Fwd: 5'-ATATAAGCTTxxxxxxxxxxxxxxxxxxxx-3' (HindIII), Rev: 5'-TATAGGATCCxxxxxxxxxxxxxxxxxxxx-3' (BamHI)).
   6. Run 5 µl of the PCR product on agarose gel to check for unspecific PCR products.
      NOTE: If more bands are detected, run the entire PCR product and purify the appropriate band by a commercial gel extraction kit.
3. Restriction
   NOTE: To clone insert upstream of TK promoter²² the purified PCR product and the vector (e.g., TK-Luc) should be digested with HindIII and BamHI.
   1. Prepare a 50 µl reaction mixture containing 1 µg purified PCR product, 1 µl HindIII (20 U/µl), 1 µl BamHI (20 U/µl) and 5 µl 10x buffer.
   2. Prepare a 250 µl reaction mixture containing 5µg vector (TK Luc, or alternative reporter vector), 5 µl HindIII (20 U/µl), 5 µl BamHI (20 U/µl), 5 µl thermosensitive Alkaline Phosphatase (1 U/µl) and 25 µl 10x buffer.
      NOTE: AP treatment of the vector greatly decreases the self-ligation of single-digested vector and thus the background.
   3. Incubate the reaction mixtures for 4 hr at 37 °C, then purify the digested PCR products and the vector with a commercial PCR purification kit.
4. Ligation
   1. Set up the reaction for ligation in PCR tubes including 2 µl 10x T4 DNA ligase buffer, 10 ng TK-Luc vector, 50 ng insert (digested PCR product), 1 µl T4 DNA ligase (400 U/µl) and nuclease-free water up to 20 µl.
      NOTE: Prepare a negative control where the reaction does not contain insert.
   2. Gently mix the reaction by pipetting up and down, spin down briefly the PCR tubes using mini-centrifuge and then incubate at RT for 10 min.
   3. Heat inactivate at 65 °C for 10 min then chill on ice and use 5 µl of the product for transformation into 50 µl competent cells (e.g., DH5α).
   4. Use standard heat shock procedure to transform bacteria, pick up colonies and isolate plasmid DNA with a commercial kit. Validate all constructs by sequencing prior to the next steps.
5. Transfection of ES Cells
   1. Use 48-well plates. Add 200 µl 0.1% gelatin solution/well 30 min prior to cell plating.
   2. Plate embryonic stem (ES) cells a day before transfection. Use feeder-free ES cells and plate at a density of 3 x 10⁴ cells per well in 250 µl ES media/well.
   3. Prepare the plasmid mixes (each mix is calculated for 8 wells, 4 untreated and 4 retinoic acid-treated). Mix 1: 1,250 ng TK-Luc Empty (negative control), 950 ng β-Galactosidase coding vector (βGal), Mix 2: 1,250 ng TK-Luc Hoxa1 enhancer (positive control), 950 ng βGal, Mix 3: 1,250 ng TK-Luc PRMT8 enhancer (region of interest), 950 ng βGal.
      NOTE: ES cells express the desired transcription factors (RAR/RXR). Leave wells untransfected to determine base values later during the luciferase and β-galactosidase measurements.
   4. Add transfection quality reduced serum media to each plasmid mix up to 106 µl total volume (enough for 8 wells).
   5. Add 4.5 µl ES quality transfection reagent then mix carefully by pipetting 15 times up and down. Incubate the transfection mix for 10 min at RT.
   6. Add 13 µl of the transfection mix to the cells per well, mix thoroughly by pipetting. Place the cells back into the incubator (37 °C) O/N before ligand treatment.
   7. Remove media carefully by aspiration from cells and add 250 µl fresh media/well containing 0.01 - 1 µM all-trans retinoic acid (RA) as final concentration or DMSO as vehicle. Incubate cells for 24 - 48 hr.
      NOTE: Retinoic acid is light sensitive.
6. Preparation of Cell Lysates
   1. Dilute 5x Lysis Buffer (1.25 ml 0.5 M Tris pH = 7.8, 1 ml 1 M dithiothreitol (DTT) (in H₂O), 10 ml 0.1 M EDTA pH = 8.0, 50 ml glycerol, 5 ml Triton X-100, sterile water up to 100 ml) to 1x using sterile water, add 20 µl 1 M DTT/10 ml and equilibrate to RT before use. Prepare 10 ml for a 48-well plate on the day of the assay.
   2. Remove media carefully by aspiration from cells. Rinse cells once with 1x PBS. Add 200 µl 1x Lysis Buffer/well directly to the cells.
   3. Shake for 2 hr at RT (chemical lysis). Freeze the cell lysates on the plate at -80 °C (mechanical lysis). Lysates can be stored at -80 °C for several days.
   4. Thaw the lysates. After complete thawing, transfer cell lysates to a 96-well plate using an electronic multichannel pipette. Transfer 80 µl of the cell lysate to 96-well clear plate for β-galactosidase assay and 40 µl of the cell lysate to a white plate for luciferase measurement. Avoid forming any bubbles.
7. β-galactosidase Assay
   NOTE: Beta galactosidase or alternative second reporters should be used as an internal control to account for non-specific experimental variations.
   1. Prepare β-galactosidase substrate buffer: 80.0 g Na₂HPO₄•7H₂O, 3.8 g NaH₂PO₄•H₂O, 30.0 g KCl, 1.0 g MgSO₄•7H₂O, make up to 1,000 ml with sterile water. Adjust pH to 7.4. Filter Sterilize (0.2 micron) and store at 4 °C
   2. Mix 1 ml β-galactosidase substrate buffer with 4 mg ONPG (o-nitrophenyl-β-D-galactosidase). Add 3.5 µl of 2-Mercaptoethanol (βME) per 1 ml of buffer just before use. Add 100 µl β-galactosidase substrate buffer to 80 µl cell lysate.
   3. Incubate the reactions at 37 °C until a faint yellow color has developed. Read absorbance at OD 420 nm on a plate reader and export data.
      NOTE: The time varies depending on the transfection efficiency, usually no longer than 30 min.
8. Luciferase Assay
   1. Prepare 50 ml luciferin substrate reagent: 15.1 mg D-luciferin salt, 82.7 mg ATP salt, 185 mg MgSO₄•7H₂O, add to 50 ml polypropylene tube, then add 1.5 ml 1 M HEPES buffer, 48.5 ml ATP free water, vortex, make sure all compounds dissolve completely. Keep 5 - 10 ml aliquots at -80 °C.
   2. Warm 5 ml luciferin substrate reagent to RT before starting. Pipette 100 µl substrate reagent to each well of the white plate containing the 40 µl of the cell lysate and measure the signal immediately using a luminescent counter machine.
   3. Obtain activities from at least three independent experiments in triplicate transfections.
   4. Calculate first the base normalized luciferase activity and base normalized β-galactosidase values for each well by subtracting the average values measured for non-transfected cells from each value measured. Then calculate the ratio of luciferase/β-galactosidase activity.

3. Characterization of Enhancer RNA

NOTE: A more direct indicator of enhancer activity has emerged from recent genome-wide studies that identified many short non-coding RNAs, ranging in size from 50 to 2,000 nt, which are transcribed from enhancers, and are termed enhancer RNAs (eRNAs)^16,25,26 (Figure 4). eRNA induction highly correlates with the induction of adjacent exon-coding genes. Thus, signal-dependent enhancer activity can be quantified in vivo by comparing eRNA production between various conditions by RT-qPCR.

1. Guidelines for Primer Design for eRNA detection by RT-qPCR
   1. Select genomic regions for eRNA measurements that are at least 1.5 - 2 kb away from annotated transcription start sites.
      NOTE: Importantly, high levels of mRNA transcription across intragenic enhancers prevent accurate quantification of eRNA in the sense orientation, antisense eRNAs at intragenic enhancers are detectable and are similar in level to eRNAs at extragenic enhancers.
   2. As a general rule, use regions 200 - 1,000 bp from the center of the transcription factor binding site of interest for primer design either on the sense or anti-sense strand (Figure 4 and Figure 6).
      NOTE: If GRO-seq, DNase-seq, H3K4me1, H3K4me3 or H3K27ac ChIP-seq data are available from the desired cell types and conditions it can be used for more accurate primer design. eRNAs are transcribed from putative enhancer regions characterized by high levels of H3K4me1 and me3. Expression of eRNAs positively correlates with the enrichment of activated enhancer histone mark H3K27ac.
   3. Design PCR primers for fluorescent dye-based RT-qPCR measurement of eRNAs by Primer3 (http://bioinfo.ut.ee/primer3/)²³Generate anti-sense sequence (reverse complement of the sense strand) for primer design as per: http://www.bioinformatics.org/sms/rev_comp.html. Insert 200 - 300 bp of sense or anti-sense sequence for primer design.
2. Measurement of eRNA Transcription
   1. Isolate total RNA from treated cells with commercial acidic phenol/chloroform-based extraction reagent according to the manufacturer's recommendations.
   2. Use DNaseI to treat isolated RNA prior to using them in reverse transcription reaction. Inactivate DNase according to the manufacturer's recommendation prior to reverse transcription.
   3. Measure RNA concentration after DNase I treatment and use 1,000 ng per RT reaction. Use high quality reverse transcriptase.
      NOTE: As large number of eRNAs may not be polyadenylated, reverse transcription requires the use of random primers.
   4. Detect reverse transcribed eRNA by RT-qPCR using standard procedure. Measure a non-treatment-dependent mRNA for normalization (e.g., 36B4, Ppia, Gapdh, Actb).

Results

We used a pan-specific RXR antibody in order to identify genome-wide which RA-regulated genes have receptor enrichment in their close proximity. Bioinformatics analysis of RXR ChIP-seq data obtained from ES cells treated with retinoic acid revealed the enrichment of the nuclear receptor half site (AGGTCA) under the RXR occupied sites (Figure 1). Using a bioinformatics algorithm we mapped back the motif search result for the half site to the RXR ChIP-seq data (Figu...

Discussion

In recent years, advances in sequencing technology have allowed large-scale predictions of enhancers in many cell types and tissues ^7-9. The workflow described above allows one to perform primary characterization of candidate enhancers chosen based on ChIP-seq data. The detailed steps and notes will help anyone to set up a routine enhancer validation in the lab.

The most critical step in the luciferase reporter assay is the transfection efficiency. It is recommended to include a GFP...

Materials

Name	Company	Catalog Number	Comments
KOD DNA polymerase	Merck Millipore	71085-3	for PCR amplification of enhancer from gDNA
DNeasy Blood & Tissue kit	Qiagen	69504	for genomic DNA isolation
QIAquick PCR Purification kit	Qiagen	28106	for PCR product purification
Gel extraction kit	Qiagen	28706	for gel extraction if there are more PCR product
HindIII	NEB	R3104L	restriction enzyme
BamHI	NEB	R3136L	restriction enzyme
FastAP	Thermo Scientific	EF0651	release of 5'- and 3'-phosphate groups from DNA
T4 DNA ligase	NEB	M0202	for ligation
QIAprep Spin Miniprep kit	Qiagen	27106	for plasmid isolation
DMEM	Gibco	31966-021	ES media
FBS	Hyclone	SH30070.03	ES media
MEM Non-Essential Amino Acid	Sigma	M7145	ES media
Penicillin-Streptomycin	Sigma	P4333	ES media
Beta Mercaptoethanol	Sigma	M6250	ES media
FuGENE HD	Promega	E2311	transfection reagent
Opti-MEM® I Reduced Serum Medium	Life Technologies	31985-062	for transfection
All-trans retinoic acid	Sigma	R2625	ligand, for activation of RAR/RXR
96-well clear plate	Greiner	655101	for Beta galactosidase assay
96-well white plate	Greiner	655075	for Luciferase assay
D-luciferin, potassium salt	Goldbio.com	115144-35-9	for Luciferase assay
ATP salt	Sigma	A7699-1G	for Luciferase assay
MgSO4•7H2O	Sigma	230391-25G	for Luciferase assay
HEPES	Sigma	H3375-25G	for Luciferase assay
Na2HPO4 •7H2O	Sigma	431478-50G	for Beta galactosidase assay
NaH2PO4 •H2O	Sigma	S9638-25G	for Beta galactosidase assay
MgSO4•7H2O	Sigma	230391-25G	for Beta galactosidase assay
KCl	Sigma	P9541-500G	for Beta galactosidase assay
ONPG (o-nitrophenyl-β-D-galactosidase)	Sigma	N1127-1G	for Beta galactosidase assay
TRIzol®	Life Technologies	15596-026	RNA isolation
High-Capacity cDNA Reverse Transcription Kit	Life Technologies	4368814	reverse transcription of eRNA
Rnase-free Dnase	Promega	M6101	Dnase treatment
SsoFast Eva Green	BioRad	750000105	RT-qPCR mastermix
CFX384 Touch™ Real-Time PCR Detection System	BioRad		qPCR machine
BioTek Synergy 4 microplate reader	BioTek		luminescent counter