Targeted DNA Methylation Analysis by Next-generation Sequencing

1. Nucleic Acid Isolation from Tissue

NOTE: Co-isolation of DNA and RNA from the same tissue sample offers the opportunity to collect DNA for epigenetic analysis and paired RNA for gene expression analysis. The example given here is a form of column purification from experimental tissue but alternative approaches for different sample types or other isolation methods can be employed. The primary requirement is for highly purified nucleic acid without contaminating protein or organic solvents.

1. Select a piece of fresh-frozen or fresh tissue for mechanical homogenization and place on ice. Store or transfer tissue to a 2.0 ml round-bottom microcentrifuge tube.
2. Add one 5 mm stainless steel bead to each tube containing tissue. Add an appropriate amount of lysis buffer, based on manufacturer’s suggestions of cell number or tissue weight.
3. Homogenize the tissue using a bead mill mechanical homogenization at 30 Hz for 30 sec.
   NOTE: The frequency and duration of mechanical homogenization is dependent on the tissue type. Softer tissues will completely homogenize using the recommended settings, while tougher tissues require setting optimization.
4. Carry out the DNA and RNA isolation using silica spin columns at room temperature following manufacturer’s protocol.
   1. Elute RNA in 50 µl of RNase-free water and place on ice. Use RNA for mRNA expression analysis by qPCR.
   2. Elute DNA in 100 µl of TE buffer. Repeat elution off DNA column using the eluent to increase DNA concentration and yield. Use DNA for targeted methylation quantitation.
5. Assess the quality of DNA and RNA using a standard spectrophotometer for absorbance at 230 nm, 260 nm, 280 nm, and 320 nm.
   NOTE: A260/A280 ratio should be ~1.8-2 for high quality DNA. A260/A280 ratio should be ~2 for high quality RNA. The presence of excess absorbance at 230 nm indicates contaminants from the isolation procedure.
6. Quality check the RNA further by using a capillary electrophoresis chip according to manufacturer’s instructions.
   NOTE: High quality RNA will have a RNA integrity number (RIN) >8. The presence of extra peaks in the electropherogram besides rRNA peaks indicate sample contamination.
7. Quantify DNA and RNA by fluorescent assay according to manufacturer’s protocol.
   NOTE: Using a fluorometric assay is more sensitive and specific for quantifying nucleic acids than spectrophotometry alone.
8. Store isolated DNA and RNA at -80 °C.

2. Target Identification and Primer Design

1. Determine the genomic region(s) of interest to be analyzed by BSAS based on the appropriate reference genome sequence. Gene promoters of differentially expressed genes are common targets of methylation quantitation.
2. Prepare an in silico bisulfite-converted reference genome for later alignment of sequencing reads by altering the .fasta sequence file in a text editor. In the 5’-3’ orientation, replace non-CpG cytosine’s with thymine’s (Figure 2).
3. Design primer sets to amplify regions of interest from bisulfite converted DNA.
   1. Select and copy the unconverted region of interest, in the 5’-3’ orientation, into a bisulfite-specific PCR primer design program (Figure 3).
      NOTE: An optimal bisulfite-PCR amplicon length is 250-400 bp per amplicon as bisulfite treatment fragments DNA and it is difficult to amplify large >400 bp regions. If, for example, a region of 1 kb is of interest, multiple primer pairs can be designed to cover this region. Design amplicons to be of equal bp size for easy pooling. Avoid amplicon lengths <250 bp because these may be of insufficient length for library generation. An optimal primer length for BSAS is > 20 bp per primer.
   2. Select a T_m range of 55-65 °C and a max T_m difference between forward and reverse primers of 1-2 °C.
   3. Select primer pairs that best cover the region of interest. Do not select primers that contain CpG sites or are directly adjacent to CpG sites, as this will cause bias in the PCR reactions. Use standard PCR primers that can be ordered from any number of academic or commercial providers.
   4. Reconstitute lyophilized primers in RNase/DNase free water at a working stock of 100 µM.
4. Store PCR primers at -20 °C. Dilute 100 µM primer stock to 10 µM working stock for PCR reactions.

3. Bisulfite Specific PCR Optimization

NOTE: For bisulfite conversion of genomic DNA, a number of different commercial kits for bisulfite conversion are available. Select the kit or protocol that best suits the planned experiment.

1. Use between 200 ng to 2 µg of genomic DNA. For optimization experiments, run multiple conversion reactions in order to have sufficient bisulfite converted DNA for multiple BS PCR reactions.
2. Use small (e.g., 10 µl) elution volumes of bisulfite converted DNA to maintain high DNA concentrations.
   NOTE: 1-2 µl of bisulfite converted DNA is sufficient for optimization PCR reactions if 1 µg genomic DNA was used for conversion input.
3. Amplify bisulfite converted DNA by bisulfite specific PCR. BS-PCR reactions require using a Taq polymerase capable of amplifying bisulfite converted DNA.
   1. Assemble the following reaction for target amplification optimization of a single amplicon. For multiple samples, assemble reactions in a 96-well PCR plate.
      25 µl 2x reaction buffer
      0.5 µl dNTP Mix
      5 µl 10 µM Forward Primer (Final 1 µM)
      5 µl 10 µM Reverse Primer (Final 1 µM)
      2 µl template bisulfite converted DNA
      0.4 µl (5 U/µl) DNA polymerase
      12.1 µl DNase free water (Final reaction volume 50 µl)
   2. Seal PCR plate with the appropriate adhesive or heat-sealed film.
   3. Place reaction in an appropriate thermal cycler and use the following cycling conditions with a heated lid:
      1. Perform an initial denaturation at 95 °C for 10 min.
      2. Denature at 95 °C for 30 sec.
      3. Anneal for 30 sec at the specific T_m of primers being used. Start at an annealing temperature a few °C below the T_m of the primer for optimization reactions.
         NOTE: Best annealing temperatures can also be determined using a gradient thermal cycler to test a range of temperatures.
      4. Perform an extension at 72 °C for 30 sec. Lengthen the extension time for longer amplicons.
      5. Repeat steps 3.3.3.2-3.3.3.4 for 35 cycles total for initial optimization. Higher cycle numbers might be necessary but should generally be avoided to prevent clonal amplifications that will create reaction artifacts.
      6. Perform a final extension at 72 °C for 7 min.
      7. Hold reactions at 4 °C.
   4. Visualize amplicons by PAGE electrophoresis. Alternatively, use a capillary electrophoresis DNA chip according to manufacturer’s protocol.
      1. Mix 24 µl of the PCR reaction with 6 µl of 5x loading dye on ice and vortex.
      2. Make 1 L of 1x TBE running buffer by mixing 200 ml of 5x TBE running buffer and 800 ml of ddH₂O. Dilute DNA ladder for a final concentration of 0.5 µg per well in ddH₂O and loading dye.
      3. Load wells with 30 µl of PCR reactions or ladder.
      4. Run the gel at 200 V for ~45 min. Stop the gel run when the first dye front reaches the end of the gel.
      5. Stain the gel with 5 µg/ml ethidium bromide by mixing 5 µl of 10 mg/ml in 100 ml of ddH₂O. Cover the entire gel and let incubate for ~5 min.
      6. Image the gel via a UV light box or multimodal imager at an excitation wavelength of 482 nm.
      7. In reference to the ladder, size the PCR amplicons and determine specificity of the reaction by looking for multiple bands other than the expected amplicon size (Figure 4A and C).
   5. Once optimal conditions for BS-PCR have been determined, perform BS-PCR on experimental samples.
      1. Purify amplicons from remaining primers and other reaction components. Perform silica column purification, SPRI bead, or gel-based clean-up and quantification.
      2. Elute amplicons in 30 µl TE buffer.
      3. Quantify amplicons using fluorometric assay using manufacturer’s protocols.
      4. If there are multiple amplicons per sample, pool them at equal weight per volume amounts and store at -20 °C.

4. NGS Library Preparation and Library QC

NOTE: BSAS NGS library preparation uses a simplified transposome-mediated protocol with dual-indexing (See Materials List for details on library preparation kit selection). This method provides an extremely rapid and high-throughput route for library construction that has been validated and shows little to no bias in bisulfite sequencing applications^16,19. Dual-indexing allows for a higher multiplexing of samples than single indexing. Other library preparation styles may be possible but have not been tested.

1. Prepare NGS libraries from BS PCR amplicons in a 96-well PCR plate. Dilute amplicons to 0.2 ng/µl for a final input mass of 1 ng, or 5 µl of the 0.2 ng/µl dilution. Perform a clean-up of generated libraries using SPRI beads. Use 50-90 µl SPRI beads for isolating libraries with this protocol. The volume of SPRI beads depends on target size of library and cycle number of desired sequencing reactions.
   NOTE: Use a microplate shaker or plate vortexer to resuspend beads.
   1. Elute the libraries off the beads and transfer to a new 96-well plate and seal with heat-sealing film for storage at -20 °C.
2. Determine the size and concentration of libraries (Figure 4B and D).
   1. Run libraries on a capillary electrophoresis DNA sizing chip according to manufacturer’s protocols.
      1. Use a high sensitivity capillary electrophoresis DNA assay to determine the average size of detected library peaks and an estimation of molarity according to manufacturer’s protocols.
   2. Quantify libraries by qPCR, using primers designed against the adapter sequences in the generated libraries and use standards with a known concentration.
      1. Run qPCR reactions on a real-time instrument using the absolute quantitation settings according to manufacturer’s protocols.
      2. Use the size of libraries and the molar quantitation from qPCR to calculate the molar concentrations of the libraries.
3. Dilute libraries to 4 nM using 10 mM Tris-Cl with 0.1% Tween-20 at pH 8.5. Store libraries in a 96-well PCR plate sealed with heat-sealing film at -20 °C.

5. Sequencing Libraries Using a Benchtop Sequencer

1. Thaw out and prepare reagents according to manufacturer’s instructions.
2. Thaw 4 nM libraries on ice.
3. Make 1 ml of 0.2 N NaOH in a 1.5 ml microcentrifuge tube.
4. Pool 4 nM libraries in equal volume amounts, if multiple libraries are to be used.
5. Dilute and denature pooled libraries according to manufacturer’s instructions to the desired final molar concentration.
   1. Keep diluted and denature pooled library on ice until ready to load on reagent cartridge.
6. Generate a sample sheet containing sample names and corresponding index information and specify to generate .fastq files only.
   1. Load sample sheet in the corresponding sample sheet folder on the sequencer.
7. Add a total of 600 µl diluted and denatured library to the reagent cartridge in the corresponding well.
8. Run the sequencing reactions to completion.

6. Methylation Quantitation Data Analysis

NOTE: There are multiple programs available for NGS data analysis including both commercial and open source. Details on running programs can be found with each specific package. General instructions are given at each step along with the specific commands for this software package. (See Materials List for details on software used in this protocol).

1. Import .fastq.gz files into appropriate sequencing analysis pipeline retaining quality scores (Qscores) for use in read trimming. (For this example program select ‘Import’ and select the sequencing method used to generate reads. Select .fastq.gz read files to import, and retain Qscores from .fastq files for read trimming applications by deselecting ‘Discard quality scores’ under ‘General Options’. Select a location for the imported reads and select ‘Finish’.)
   1. Trim and retain reads using the following settings in a NGS data analysis pipeline. (For this program select ‘Trim Sequences’ under ‘NGS Core Tools’, and select the reads to be trimmed.)
      1. Retain only reads containing ≥ Q30 scores. (For this program select ‘Trim using quality scores’ and set the limit to 0.001.)
      2. Retain only reads containing ≤ 1 ambiguous nucleotide. (For this program select ‘Trim ambiguous nucleotides’ and set the ‘maximum number of ambiguities’ to 1. Select a location to save trimmed reads and select ‘Finish’.)
   2. Map trimmed reads to the in silico bisulfite converted gene reference from 2.2. (For the example software select ‘Map Reads to Reference’ under ‘NGS Core Tools’. Select trimmed reads to be mapped and select ‘next’. Select the converted reference sequence and select ‘no masking’ under ‘reference masking’. Select ‘next’.)
      1. Assign a maximum penalty for mismatches, insertions, and deletions. Set parameters such that 100% of the read length is required to map with 90% sequence identity. For example software, assign a score of 3 for mismatches, insertions, and deletions. Assign a score of 1 for minimum fraction of read length mapped to reference, and assign a score of 0.9 for minimum fraction of identity between mapped read and reference.
      2. Discard or ignore non-specific reads mapped to reference. These reads can map equally well to multiple regions. For example software select ‘Ignore’ under ‘Non-specific match handling’ and select ‘next’.
      3. For alignment output generate each sample as an independent file. For example software, select ‘Create stand-alone read mappings’ under ‘Output options’ and ‘Save’ under ‘Result handling’. Select a location to save the read mapping and select ‘Finish’.
   3. Run a variant calling workflow on the mapped reads. For the example software, select ‘Probabilistic Variant Detection’ under ‘Resequencing Analysis’, select the read mapping to be analyzed and select ‘next’. Select ‘Ignore non-specific matches’ under ‘Read filters’ and select ‘next’.
      1. Ensure that the sequence is at a depth of ≥ 1000 X by setting the ‘minimum coverage’ under ‘Significance’ to 1,000.
      2. Ensure that the variant call is present in both forward and reverse reads. For example software, select ‘Require presence in both forward and reverse stands’ under ‘Variant filters’ and select ‘next’.
      3. Export aligned, variant called data as an annotated table. For example software, select ‘Create annotated table’ under ‘Output options’ and select a location to save the variant table file. Select ‘Finish’.
2. Calculate 5-mC percentage using the variant frequencies at annotated CpG sites in the region of interest from the variant table file. The frequency of cytosine at a reference C in CG is the 5-mC percentage.
   NOTE: When using variant calling algorithms for highly methylated CpG’s, replace CpG cytosines with thymines. This will identify the mapped cytosine’s as variants, resulting in cytosine frequency at CpG’s.
3. Determine bisulfite conversion efficiency (bsc) by first calculating the percentage of C’s not in CpG sites in the reference (C). Multiply the percentage of non CpG C’s by number of total T’s mapped (T_mp), and divide by total C’s mapped at the reference T’s (C_mpT). 100 minus the frequency calculated is the bisulfite conversion efficiency (bsc) of that library (Equation 1).

Equation 1. Bisulfite conversion efficiency.
NOTE: While mammalian genomes contain very low amounts of non-CpG cytosine methylation (with the exception potentially of stem cells²⁰), plant genomes contain a significant amount. In order to determine bisulfite conversion efficiency in these reference genomes, a suitable route is to sequence a portion of the mitochondrial genome or chloroplast genome or known unmethylated genes present in plant genomes²¹ and calculating conversion efficiency as described above. Bisulfite conversion efficiency should be ≥98% in most cases, with the unconverted C’s potentially arising from non CpG methylation.