Materials and Reagents

Pipette tips
Note: Low retention tips are recommended.
1.5 ml microcentrifuge tubes (Eppendorf, catalog number: 0030108051 )
PCR tubes (Thermo Fisher Scientific, Applied Biosystems^TM, catalog number: A30588 )
Clean razor blade
Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, Invitrogen^TM, catalog number: Q32850 ) (used to quantify double strand DNA in combination with Qubit Fluorometer)
MseI (New England Biolabs, catalog number: R0525S ) (used to digest maize genomic DNA, CutSmart Buffer is included)
Note: The enzyme could be replaced according to the in silico enzyme selection result in Procedure A.
Agarose
Note: For the gel size selection, we recommend low melting point agarose.
TBE buffer (Thermo Fisher Scientific, Invitrogen^TM, catalog number: 15581044 )
AMPure XP (Beckman Coulter, catalog number: A63881 )
Absolute ethanol
Note: To dilute, use DNase/RNase-free distilled water.
1 M Tris-HCl (Thermo Fisher Scientific, Invitrogen^TM, catalog number: 15568025 )
Note: To dilute, use DNase/RNase-free distilled water.
NEBuffer 2 (New England Biolabs, catalog number: B7002S )
Klenow Fragment (3’→5’ exo-) (New England Biolabs, catalog number: M0212S ) (used to end-repair and A-tail DNA)
T4 DNA Ligase (New England Biolabs, catalog number: M0202S , T4 DNA ligation buffer is included) (used to ligate sequencing adaptors)
100-bp DNA ladder
QIAquick Gel Extraction Kit (QIAGEN, catalog numbers: 28704 , 28706 )
EpiTect Fast Bisulfite Conversion Kit (QIAGEN, catalog number: 59824 )
Note: Keep this kit at an appropriate temperature according to the manufacturer’s guide.
TruSeq DNA LT Sample Prep Kit (Illumina, catalog number: FC-121-2001 )
PfuTurbo DNA Polymerase (Agilent Technologies, catalog number: 600250 )
DNase/RNase-free distilled water (Thermo Fisher Scientific, Invitrogen^TM, catalog number: 10977015 )
Deoxynucleotide (dNTP) Solution Set (New England Biolabs, catalog number: N0446S )
MinElute PCR Purification Kit (QIAGEN, catalog number: 28004 )
10 mM dNTP mix (see Recipes)
5x RRBS dNTP mix (see Recipes)

Equipment

Pipettes:
1. Thermo Fisher Scientific, Thermo Scientific^TM, model: P5000, catalog number: 4641110N
2. Thermo Fisher Scientific, Thermo Scientific^TM, model: P1000, catalog number: 4641100N
3. Thermo Fisher Scientific, Thermo Scientific^TM, model: P300, catalog number: 4641090N
4. Thermo Fisher Scientific, Thermo Scientific^TM, model: P200, catalog number: 4641080N
5. Thermo Fisher Scientific, Thermo Scientific^TM, model: P100, catalog number: 4641070N
6. Thermo Fisher Scientific, Thermo Scientific^TM, model: P20, catalog number: 4641060N
7. Thermo Fisher Scientific, Thermo Scientific^TM, model: P10, catalog number: 4641030N
8. Thermo Fisher Scientific, Thermo Scientific^TM, model: P2, catalog number: 4641010N
Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Invitrogen, model: Qubit^® 2.0 Fluorometer , catalog number: Q32866)
Heating block (Eppendorf, model: ThermoMixer comfort , catalog number: 5355000011)
Incubator (Only Science, Firstek, model: S300 )
Thermocycler
1. Thermo Fisher Scientific, Applied Biosystems^TM, model: GeneAmp^TM PCR System 9700 , catalog number: 4413750
2. Bio-Rad Laboratories, model: T100^TM Thermal Cycler, catalog number: 1861096
Centrifuge
1. Eppendorf, model: Centrifuge MiniSpin^® plus , catalog number: 5453000011
2. Eppendorf, model: Centrifuge 5424 , catalog number: 5424000410
3. GYROZEN, model: 1730R
Magnetic stand
1. Thermo Fisher Scientific, DynaMag^TM-2 Magnet, catalog number: 12321D
2. Thermo Fisher Scientific, Applied Biosystems^TM, DynaMag^TM-PCR Magnet, catalog number: 492025

Software

System requirements: Linux/Unix or Mac OS, python 2.7+, R 3.2.0
Python and R scripts and a tutorial are available at: https://gitlab.com/fmhsu0114/maize_RRBS

Procedure

The complete pipeline contains 1. in slico enzyme selection, followed by 2. RRBS library preparation. The resulted library is to be sequenced to generate DNA methylation profiles at genomic regions of interest.

In silico enzyme selection
1. Clone the GitLab project from: https://gitlab.com/fmhsu0114/maize_RRBS.
2. Download the genome FASTA file of your interest, i.e., maize_genome.fa, and use the following command to build the genome database:
  $ python Build_Gnome.py AGPv3.fa
  2 output files (genome.shelve and log-genome.txt) will be generated.
3. Download the annotation refFlat file, i.e., refFlat.txt, and use the following command to build the gene annotation database:
  $ python Build_refgene_shelve.py refFlat.txt log-genome.txt
  1 output file (refgenes.shelve) will be generated.
4. Download the annotation gff3 file and grep ‘repeats’ into an individual file, i.e., repeats.gff3, and use the following command to build the repeat annotation database:
  $ python Build_repeats_Gnome.py repeats.gff3, log-genome.txt
  1 output file (repeat.shelve) will be generated.
5. Build the genome CG, CHG and CHH sites database with the following command:
  $ python Build_CG_CHG_CHH_Gnome.py genome.shelve log-genome.txt
  1 output file (CG_CHG_CHH.shelve) will be generated.
6. Perform in silico digestion with the provided enzyme list in the GitLab project (enzyme_cutting_sites.txt) with the following command:
  $ python Enzyme_digestione.py
  fragment_info (per MSRE and per size range) for each enzyme and a Summary-enzyme-coverage.txt file will be generated.
7. Place the output files from 1.6 into one directory ./fragment/ and run the enrichment analysis with the following command:
  $ python RRBS_fragments_enrichment.py repeats.shelve, refgenes.shelve
  log₂_fold_enrichment (per MSRE and per size range) files will be generated (see Table 1).
  
  Table 1. Fold enrichment values in log₂scale of the MseI 40-300-bp fragments
8. Move all of the log₂_fold_enrichment files into one directory ./output/ and plot the enrichment bar charts with the following command:
  $ Rscript Barchart_enrichment.R
  See Figure 2 to Figure 4 as examples.
  
  Figure 2. Bar chart of the MseI 40-300-bp fragments shows enrichment in promoters
  
  Figure 3. Bar chart of the CviQI 40-280-bp fragments shows enrichment in exons, introns, UTRs and splice sites
  
  Figure 4. Bar chart of the BfaI 50-300-bp fragments shows enrichment in promoters
9. Use the following command to calculate the covered promoter and genebody region:
  $ python RRBS_fragments_enrichment_gene_level_per_MRE_Size.py -i Enzyme -s size_range -c 1
  This script has 3 parameters:
  -i: restriction enzyme
  -s: size selection range, e.g., 40-280
  -c: number of cytosines per promoter/genebody to be included
  An output file, Fragment-Enrichment-genebody-pmt-Enzyme-size_range-1.txt, will be generated (Table 2 and Table 3). This file contains the percentages of the promoters and genebody covered in the reduced representation (RR) genome.
  Note: In Table 2 and Table 3, Pct refers to percentage; Pmt refers to promoter.
  
  Table 2. Promoter and genebody covered in the MseI 40-300-bp fragments
  
  Table 3. Promoter and genebody covered in the BfaI 50-300-bp fragments
RRBS library preparation
1. Use Qubit Fluorometric Quantification with Qubit dsDNA BR assay to quantify maize genomic DNA according to the manufacturer’s instructions.
  Note: To acquire high quality double strand genomic DNA, we recommend spin column-based DNA extraction kit, e.g., QIAGEN DNeasy Plant Mini kit. If Pheno-Chloroform method is used, we suggest another round of spin-column based purification.
2. Aliquot 1 μg of maize genomic DNA into a 1.5 ml tube, add 5 μl of 10x CutSmart Buffer, 1 μl of MseI and ddH₂O to a final volume of 50 μl. Mix carefully with a pipette and incubate at 37 °C on the heating block for 4 h.
  Note: MseI used here is designed for promoter-enriched RRBS in maize. For different genomes (organisms), researchers should conduct ‘in silico enzyme selection’ and replace MseI with the predicted enzyme, and modify the digestion time, buffer and temperature according to the manufacturer’s instructions.
3. [Optional] Take 5 μl for gel electrophoresis (100 ng, 1.5%, 0.5x TBE) and add 5 μl ddH₂O to the enzyme digestion product.
  Note: This step is to confirm the enzyme digestion efficiency.
4. Post-digestion AMPure beads purification
  1. Add 150 μl of well-inverted AMPure XP beads solution to 50 μl of DNA, vortex for 10 sec, briefly centrifuge and incubate at RT for 5 min.
    Note: When vortexing solutions with AMPure XP beads, avoid high vortexing strength causing bubbles. After vortexing, use microcentrifuge to spin down quickly to remove liquids from the lid and the wall of the tube.
  2. Place the tube on the magnetic stand for 3 min.
  3. Remove the supernatant.
  4. Add 200 μl of 80% EtOH without disturbing the beads and incubate at RT for 30 sec. Remove the EtOH.
    Note: Do not remove the tube from the magnetic stand.
  5. Repeat the washing Step B4d.
  6. Air-dry the beads for 5 min.
  7. Remove the tube from the magnetic stand.
  8. Add 32 μl of 10 mM Tris-HCl and pipet up and down to mix thoroughly. Incubate at RT for 2 min.
  9. Place the tube on the magnetic stand for 3 min.
  10. Transfer the supernatant to a PCR tube.
5. End repair and A-tailing
  1. Add 6 μl of 10x NEBuffer 2, 12 μl of 5x RRBS dNTP mix and 3 μl of 3’→5’ Klenow exo- to the purified DNA and add ddH₂O to a final volume of 60 μl.
  2. Incubate using thermocycler at 30 °C for 20 min and 37 °C for 20 min, then place the tube on ice.
6. DNA purification with AMPure XP beads
  1. Add 180 μl of well-inverted AMPure XP beads to the 60 μl of DNA, vortex for 10 sec, briefly centrifuge and incubate at RT for 5 min.
  2. Place the tube on the magnetic stand for 3 min.
  3. Remove the supernatant.
  4. Add 200 μl of 80% EtOH without disturbing the beads and incubate at RT for 30 sec. Remove the EtOH.
    Note: Do not remove the tube from the magnetic stand.
  5. Repeat the washing Step B6d.
  6. Air-dry the beads for 5 min.
  7. Remove the tube from the magnetic stand.
  8. Add 24 μl of 10 mM Tris-HCl and pipet up and down to mix thoroughly. Incubate at RT for 2 min.
  9. Place the tube on the magnetic stand for 3 min.
  10. Transfer the supernatant to a PCR tube.
7. Ligate with a TruSeq methylated Illumina adaptor
  1. Add 3 μl of 10x T4 DNA ligation buffer, 2 μl of adaptor and 1 μl of T4 DNA ligase to the 24 μl end-repaired and A-tailed DNA to a final volume of 30 μl.
  2. Incubate at 16 °C overnight in a thermocycler.
8. Gel size selection (Figure 5)
  1. Load all the adapter-ligated DNA to 1.5% 0.5x TBE gel (15 μl per well, in total 2 wells) for gel electrophoresis with a 100-bp DNA ladder.
  2. Use a clean razor blade for each sample to cut 40-300-bp fragments and purify the DNA using the QIAGEN Qiaquick Gel Extraction Kit following the manufacturer’s instructions. Elute DNA twice with a total volume of 34 μl of EB buffer. Transfer the solution to a PCR tube. Take 1 μl for Qubit quantification (optional).
  Note: One should avoid UV light during gel size selection. Since the adaptor sequence is about 120 bp, we cut 160-420-bp fragments from the gel. For each well in the previous step, use one QIAGEN Qiaquick Gel Extraction reaction. For the first round of elution, add 10 μl of EB buffer to each MinElute spin column, and collect the 20 μl of the solution in a new tube. Then use another 14 μl of EB buffer to elute the remaining DNA from each spin column one by one.
  
  Figure 5. Gel size selection
9. Bisulfite conversion
  Note: Confirm that the reagents in the EpiTect Fast Bisulfite Conversion Kit are ready to use.
  1. Add RNase-free water to the DNA solution to a final volume of 40 μl, mix with a pipette and aliquot 20 μl to another PCR tube.
  2. Add 35 μl of DNA Protect Buffer and 85 μl of a Bisulfite Solution to each tube to a final volume of 140 μl and mix thoroughly with gentle pipetting.
  3. Run in a thermocycler with the following program:
    
    Then hold at 20 °C. Lid temperature: 100 °C.
10. Briefly centrifuge the PCR tubes containing the bisulfite reaction with microcentrifuge to remove the liquid from the tube lid and transfer to clean 1.5 ml microcentrifuge tubes. For each tube, add 310 μl of freshly prepared Buffer BL containing 10 μg/ml carrier RNA, mix thoroughly with pipetting and then centrifuge briefly.
11. Add 250 μl of 96-100% EtOH, gently pipet for 10 sec and then centrifuge briefly.
12. Transfer all of the solution to MinElute DNA spin columns with collection tubes. Centrifuge at maximum speed for 1 min and discard the flow-through.
13. Washing
  Add 500 μl Buffer BW to each spin column, centrifuge at maximum speed for 1 min and discard the flow-through.
14. Desulfonation
  1. Add 500 μl of Buffer BD to each spin column and incubate for 15 min at RT.
  2. Centrifuge at maximum speed for 1 min and discard the flow-through.
  3. Add 500 μl Buffer BW to each spin column, centrifuge at maximum speed for 1 min and discard the flow-through.
  4. Repeat the washing Step B14c.
  5. Add 250 μl of 96-100% EtOH to each spin column and centrifuge at maximum speed for 1 min.
  6. Place each spin column into a clean 1.5 ml microcentrifuge tube, open the lid and incubate using a 60 °C heating block to evaporate the liquid.
  7. Elute the DNA by adding 15 μl of Buffer EB onto the center of the spin column membrane, close the lid and incubate at RT for 1 min.
  8. Centrifuge at 15,000 x g for 1 min.
15. PCR amplification
  1. Add 2.5 μl Illumina TruSeq PCR primer cocktail, 0.5 μl of 10 mM dNTP, 5 μl of 10x Pfu Turbo Buffer, 0.5 μl of Pfu Turbo DNA polymerase and 26.5 μl of ddH₂O to the eluted DNA solution to a final volume of 50 μl. Mix thoroughly and centrifuge briefly.
  2. Run in a thermocycler with the following program:
16. DNA purification with AMPure XP beads
  1. Add 150 μl of well-inverted AMPure XP beads to the 50 μl of DNA, vortex for 10 sec, briefly centrifuge and incubate at RT for 5 min.
  2. Place the tube on the magnetic stand for 3 min.
  3. Remove the supernatant.
  4. Add 200 μl of 80% EtOH without disturbing the beads and incubate at RT for 30 sec. Remove the EtOH.
    Note: Do not remove the tube from the magnetic stand.
  5. Repeat the washing Step B16d.
  6. Air-dry the beads for 5 min.
  7. Remove the tube from the magnetic stand.
  8. Add 22 μl of 10 mM Tris-HCl and pipet up and down to mix thoroughly. Incubate at RT for 2 min.
  9. Place the tube on the magnetic stand for 3 min.
  10. Transfer the supernatant to a 1.5 ml microcentrifuge tube.
17. Quantify the final library concentration with the Qubit assay.
18. Sequence the RRBS library with Illumina paired-end sequencing. Any regular length of pair-end reads, e.g., 100-250 bp, is fine.

Data analysis

Our in silico pipeline creates a table and bar chart for each enzyme with its associated range of size selection. Figure 1 is a flowchart of this pipeline, and users could see input and output files of each step. Figure 2 is a bar chart showing the enrichment of each potential ROI from approximately 40-300 bp size-selected MseI-digested DNA. The X-axis is in log₂ scale, so values > 0 can be considered to be an enrichment of these MseI-digested DNA fragments. Table 1 includes the values used to plot Figure 2. From Figure 2 and Table 1, we can see that the MseI-digested 40-300-bp fragments are enriched in promoters. As an alternative, Figure 3 shows the result from another enzyme (CviQI) 40-280 bp, with fragments more enriched in exons, splice sites, introns and UTRs.
Different enzymes produce similar enrichment values. Figure 2 and Figure 4 suggest that MseI 40-300-bp and BfaI 50-300-bp fragments are enriched in promoters. Table 2 and Table 3 show the statistics of the promoters and genebody percentages covered. For the MseI RR genome, by sequencing 566 Mbp of the 2.1 Gb genome, 84.3% of promoters are included. As for the BfaI RR genome, sequencing 23.0% of the genome includes 74.8% of the promoters. Users need to determine the tradeoff between the cost of 4.4% of the genome and 8.5% of the promoters.

Conclusion

We present an in silico pipeline to perform ROI-directed RRBS in maize with an experimentally verifiable prediction. The most critical steps for users to run our pipeline are: 1) Preparing the input files; 2) Choosing an optimal enzyme to show clear enrichment of the selected ROIs (see Figure 2).
Input files can often be downloaded from public databases, such as Ensembl (Zerbino et al., 2017). It is important to note that our pipeline requires an individual repeat annotation file in .gff3 format. In addition to the bar charts and table (Figures 2 and 3, Table 1), when judging which size-selection region is optimal, the total RR genome size and read length of the sequencing platform should also be considered. For example, if two RR genomes show 2x enrichment in promoters and one covers 50% of the whole genome while the other covers 25%, the tradeoff between sequencing cost and quantity of information should be considered.
For experimental validation, size selection and bisulfite conversion are the most critical steps. In this protocol, we used gel size selection to ensure that the fragment length could be visualized. As an alternative, AMPure XP beads can also be used. To ensure efficient bisulfite conversion, 0.01% sheared lambda phage DNA could be pooled into original genomic DNA prior to enzyme digestion. Since lambda phage DNA lacks methylated cytosine, non-conversion rate can be estimated by calculating the methylation level of lambda phage DNA to evaluate bisulfite conversion efficiency.

Recipes

10 mM dNTP mix
10 mM dATP
10 mM dTTP
10 mM dCTP
10 mM dGTP
5x RRBS dNTP mix
200 μM dATP
200 μM dTTP

Acknowledgments

This study was supported by funding from the Institute of Plant and Microbial Biology, Academia Sinica and grants obtained from Taiwan Ministry of Science and Technology (106-2633-B-001-001) to P.-Y.C. We thank Yi-Jing Lee for providing the equipment information. The authors declare that they have no competing financial interests.

References

Cokus, S. J., Feng, S., Zhang, X., Chen, Z., Merriman, B., Haudenschild, C. D., Pradhan, S., Nelson, S. F., Pellegrini, M. and Jacobsen, S. E. (2008). Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature 452(7184): 215-219.
Gent, J. I., Ellis, N. A., Guo, L., Harkess, A. E., Yao, Y., Zhang, X. and Dawe, R. K. (2013). CHH islands: de novo DNA methylation in near-gene chromatin regulation in maize. Genome Res 23(4): 628-637.
Hsu, F. M., Yen, M. R., Wang, C. T., Lin, C. Y., Wang, C. R. and Chen, P. Y. (2017). Optimized reduced representation bisulfite sequencing reveals tissue-specific mCHH islands in maize. Epigenetics Chromatin 10(1): 42.
Law, J. A. and Jacobsen, S. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet 11(3): 204-220.
Li, Q., Gent, J. I., Zynda, G., Song, J., Makarevitch, I., Hirsch, C. D., Hirsch, C. N., Dawe, R. K., Madzima, T. F., McGinnis, K. M., Lisch, D., Schmitz, R. J., Vaughn, M. W. and Springer, N. M. (2015). RNA-directed DNA methylation enforces boundaries between heterochromatin and euchromatin in the maize genome. Proc Natl Acad Sci U S A 112(47): 14728-14733.
Schnable, P. S., Ware, D., Fulton, R. S., Stein, J. C., Wei, F., Pasternak, S., Liang, C., Zhang, J., Fulton, L., Graves, T. A., Minx, P., Reily, A. D., Courtney, L., Kruchowski, S. S., Tomlinson, C., Strong, C., Delehaunty, K., Fronick, C., Courtney, B., Rock, S. M., Belter, E., Du, F., Kim, K., Abbott, R. M., Cotton, M., Levy, A., Marchetto, P., Ochoa, K., Jackson, S. M., Gillam, B., Chen, W., Yan, L., Higginbotham, J., Cardenas, M., Waligorski, J., Applebaum, E., Phelps, L., Falcone, J., Kanchi, K., Thane, T., Scimone, A., Thane, N., Henke, J., Wang, T., Ruppert, J., Shah, N., Rotter, K., Hodges, J., Ingenthron, E., Cordes, M., Kohlberg, S., Sgro, J., Delgado, B., Mead, K., Chinwalla, A., Leonard, S., Crouse, K., Collura, K., Kudrna, D., Currie, J., He, R., Angelova, A., Rajasekar, S., Mueller, T., Lomeli, R., Scara, G., Ko, A., Delaney, K., Wissotski, M., Lopez, G., Campos, D., Braidotti, M., Ashley, E., Golser, W., Kim, H., Lee, S., Lin, J., Dujmic, Z., Kim, W., Talag, J., Zuccolo, A., Fan, C., Sebastian, A., Kramer, M., Spiegel, L., Nascimento, L., Zutavern, T., Miller, B., Ambroise, C., Muller, S., Spooner, W., Narechania, A., Ren, L., Wei, S., Kumari, S., Faga, B., Levy, M. J., McMahan, L., Van Buren, P., Vaughn, M. W., Ying, K., Yeh, C. T., Emrich, S. J., Jia, Y., Kalyanaraman, A., Hsia, A. P., Barbazuk, W. B., Baucom, R. S., Brutnell, T. P., Carpita, N. C., Chaparro, C., Chia, J. M., Deragon, J. M., Estill, J. C., Fu, Y., Jeddeloh, J. A., Han, Y., Lee, H., Li, P., Lisch, D. R., Liu, S., Liu, Z., Nagel, D. H., McCann, M. C., SanMiguel, P., Myers, A. M., Nettleton, D., Nguyen, J., Penning, B. W., Ponnala, L., Schneider, K. L., Schwartz, D. C., Sharma, A., Soderlund, C., Springer, N. M., Sun, Q., Wang, H., Waterman, M., Westerman, R., Wolfgruber, T. K., Yang, L., Yu, Y., Zhang, L., Zhou, S., Zhu, Q., Bennetzen, J. L., Dawe, R. K., Jiang, J., Jiang, N., Presting, G. G., Wessler, S. R., Aluru, S., Martienssen, R. A., Clifton, S. W., McCombie, W. R., Wing, R. A. and Wilson, R. K. (2009). The B73 maize genome: complexity, diversity, and dynamics. Science 326(5956): 1112-1115.
Smith, Z. D., Gu, H., Bock, C., Gnirke, A. and Meissner, A. (2009). High-throughput bisulfite sequencing in mammalian genomes. Methods 48(3): 226-232.
Zerbino, D. R., Achuthan, P., Akanni, W., Amode, M. R., Barrell, D., Bhai, J., Billis, K., Cummins, C., Gall, A., Giron, C. G., Gil, L., Gordon, L., Haggerty, L., Haskell, E., Hourlier, T., Izuogu, O. G., Janacek, S. H., Juettemann, T., To, J. K., Laird, M. R., Lavidas, I., Liu, Z., Loveland, J. E., Maurel, T., McLaren, W., Moore, B., Mudge, J., Murphy, D. N., Newman, V., Nuhn, M., Ogeh, D., Ong, C. K., Parker, A., Patricio, M., Riat, H. S., Schuilenburg, H., Sheppard, D., Sparrow, H., Taylor, K., Thormann, A., Vullo, A., Walts, B., Zadissa, A., Frankish, A., Hunt, S. E., Kostadima, M., Langridge, N., Martin, F. J., Muffato, M., Perry, E., Ruffier, M., Staines, D. M., Trevanion, S. J., Aken, B. L., Cunningham, F., Yates, A. and Flicek, P. (2017). Ensembl 2018. Nucleic Acids Res.