Evaluation of cytosine conversion methods for whole-genome DNA methylation profiling

Plant growth and DNA isolation

Seeds of soybean genotype Williams 82 were sown in a Miracle-Gro Moisture Control potting medium and grown at 25°C, 16 h d⁻¹ light at 230 to 365 μM m⁻² s⁻¹, and 90% relative humidity in a growth chamber. At the V2 developmental stage (Fehr and Caviness, 1977), a trifoliate leaf was collected for DNA isolation (Figure 2).

Figure 2. Schematic illustration of experimental design and project approach for WGBS and EM-seq library preparation.

Created with BioRender.com. WGBS, whole-genome bisulfite sequencing; EM-seq, enzymatic methylation sequencing; NEB, New England BioLabs.

A detailed protocol can be found on protocols.io. Genomic DNA (gDNA) was isolated using the Zymo Research Quick-DNA Plant/Seed Miniprep kit (Cat #D6020; Irvine, CA, USA). DNA concentration was measured using a Qubit fluorometer (Thermo Fisher Scientific; Waltham, MA, USA) coupled with a Thermo Fisher Scientific dsDNA High Sensitivity Assay kit (Cat #Q32851). DNA purity was assessed from 260/230 and 260/280 nm absorbance ratios using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). Furthermore, the quality and size of obtained DNA was determined using gel electrophoresis (Fisher Scientific gel casting and electrophoresis apparatus (Cat #FB-SB-1316); 1% agarose gel with 1X Tris acetate-EDTA buffer (ThermoFisher; Cat # B49). Nucleic Acid was stained with Biotium GelRed Nucleic Acid Gel Stain (Fisher Scientific; Cat # NC9594719) and gel imaging was performed using LiCOR ODYSSEY imager (Serial # OFC1321) under a 520 nm channel. As a non-methylated control, the soybean gDNA was spiked with Escherichia coli gDNA (Zymo Research, Cat #D5016) representing 1% of the total gDNA input. The spiked sample was then separated into three aliquots (each 65 μL), which were sheared using Qsonica sonicators (Newtown, CT, USA). Sonication settings entailed 15 sec on/90 sec off for 8 cycles at a 20 kHz frequency, resulting in DNA fragments of approximately 350 bp. Sheared aliquots were used for library preparation, as described below.

Whole-genome Illumina sequencing library preparation

Sequencing libraries were prepared using the NEBNext^® Ultra™ II DNA Library Prep Kit for Illumina^® (New England Biolabs; Cat #E7645S; Ipswich, MA, USA). A total of 200 ng (for NEB EM-seq kit) or 1 μg (Zymo and QIAGEN Conversion kits) of starting DNA was used based on the manufacturer’s instructions (Figure 2). The sheared DNA fragments were repaired, and Illumina Methylated Adaptors were ligated using NEBNext^® Multiplex Oligos for Illumina (NEB Cat #E7535S/L). Adapter-ligated fragments were then purified using Solid Phase Reversible Immobilization (SPRI) magnetic beads (Beckman Coulter Inc.; Cat #B23317; Brea, CA, USA). For the cytosine conversion, the adapter-ligated, purified DNA samples were processed with either a Zymo Conversion kit (EZ DNA Methylation™ Kit, Zymo Research; Cat #D5001), QIAGEN Conversion kit (EpiTect Bisulfite Kit, QIAGEN; Cat #59104), or NEB EM-seq kit (NEBNext^® Enzymatic Methyl-seq Conversion Module, New England Biolab; Cat #E7125S/L) (Figure 2).

For SB conversion using the Zymo Conversion kit, DNA was incubated with CT Conversion Reagent (Zymo Research; Cat #D5001-1) (Table 3 and Figure 1) in a thermocycler for approximately 16 hrs. The converted DNA samples were then desulphonated using M-Desulphonation Buffer (Zymo Research; Cat #D5001-5). For SB conversion using the QIAGEN Conversion kit, the sample was incubated in a thermocycler for approximately 5 hrs with the Bisulfite mix (Table 2 and Figure 1) and the DNA protection buffer provided in the kit. Converted DNA samples were then treated with Buffer BD from the kit for desulphonation. For EM conversion, the sample was first treated with TET2 (NEB Cat #E7130AVIAL; part of the NEBNext^® Enzymatic Methyl-seq Conversion Module) and APOBEC (NEB Cat #E7133AVIAL) (Table 1 and Figure 1) enzymes in a series of steps, which ultimately converted cytosines to uracil, allowing for the detection of 5 mC and 5 hMC.

Figure 3. Comparison of sequencing data quality between the cytosine conversion kits.

Quality was determined by library recovery (%), number of unprocessed and preprocessed raw sequences, average length of preprocessed sequences (bp), duplication rate (%), GC content (%), mapping efficiency (%), total coverage (X), and cytosine conversion efficiency (%). NEB, New England BioLabs.

Converted samples were purified using Beckman Coulter SPRIselect magnetic beads (Cat #B23317). The barcoding of converted DNA fragments was carried out with NEBNext^® Multiplex Oligos for Illumina (Methylated Adaptor, Index Primers Set 1) (NEB Cat #E7535S/L) using PCR (16 cycles) with EpiMark^® Hot Start Taq DNA Polymerase (NEB Cat #M0490S) following manufacturer instructions (more details can be found on protocols.io).

Quality control and sequencing of WGBS and EM-seq libraries

After amplification, prepared libraries were purified using SPRIselect magnetic beads (Cat #B23317). The concentration of prepared libraries was measured using the Thermo Fisher Scientific dsDNA High Sensitivity Assay kit (Cat #Q32851) with an Invitrogen Qubit Fluorometer as well as a NanoDrop ND-1000 spectrophotometer (Cat# E1123552; Thermo Fisher Scientific). Prepared libraries were then sent to Novogene Corporation Inc. (Sacramento, CA, USA) for downstream quality assessment and whole-genome sequencing. Fragment size and concentration were validated with an Agilent Bioanalyzer (Agilent Technologies; Santa Clara, CA, USA). Libraries were sequenced on an Illumina Hi-Seq sequencer to obtain 150 bp paired-end reads.

Data analysis

The quality of raw reads obtained from the three samples was assessed using FastQC (RRID:SCR_014583) (v.0.11.8) (Andrews, 2010) and MultiQC (RRID:SCR_014982) (v1.9) (Ewels et al., 2016). Adapter sequences and poor-quality bases (Phred <20) were removed using Trim Galore (RRID:SCR_011847) (v0.4.2) (Bolger et al., 2014). The pre-processed reads were then aligned to the soybean reference genome (Gmax_275_v2.0) obtained from Phytozome (Goodstein et al., 2012) using Bismark (RRID:SCR_005604) (v.0.22.3) aligner with default parameters (Krueger and Andrews, 2011). After removing duplicate reads, Bismark (v.0.22.3) (Krueger and Andrews, 2011) was used to detect cytosine methylation at a single-base resolution. Cytosine methylation levels at each region of the sequence were calculated as the number of methylated C vs. total C and T present. To estimate and evaluate false-positive methylation levels and efficiency of cytosine deamination, reads were aligned to the non-methylated controls, soybean chloroplast genome (NC_007942.1), and E. coli non-methylated genomic DNA sequences. Data were visualized with RStudio (RRID:SCR_000432) (v1.1.463) (R Core Team, 2020) implementing the ggplot2 package (RRID:SCR_014601) (v.3.3.5) (Wickham, 2016). The code is available from GitHub and is archived with Zenodo (Wijeratne, 2022).

The QIAGEN kit outperformed other kits in time requirement and DNA recovery

Essential factors to consider when selecting a kit are the final yield of libraries and the ease of the protocol. Based on the time required for cytosine conversion steps, the QIAGEN Conversion kit was the fastest method to prepare the whole library (≈ 24 hrs), followed by the NEB EM-seq kit (≈ 26 hrs) and Zymo Conversion kit (≈ 35 hrs). When evaluating the cytosine conversion step, the Zymo Conversion kit needed the longest time of ≈ 16 hrs. However, the QIAGEN Conversion kit and the NEB EM-seq kit required a total of ≈ 5 hrs, which included the post-conversion cleanup step (Tables 1-3).

Following cytosine conversion, the highest recovery of DNA was observed with the QIAGEN conversion kit (100%), followed by the Zymo Conversion (20%) and NEB EM-seq (7.40%) kits (Figure 3). Previously, researchers noted factors that resulted in the reduction of libraries while using SB-based cytosine conversion. In a study conducted by Tanaka and Okamoto (2007), a real-time PCR experiment indicated that DNA degradation during SB conversion occurs primarily due to the hydrolytic reaction or the prolonged exposure to the SB. Additionally, prolonged exposure to SB results in the formation of abasic sites and subsequent DNA strand damage (Tanaka and Okamoto, 2007). An abasic site is a region of DNA that lacks both purine and pyrimidine bases due to DNA damage formed by spontaneous hydrolysis of the N-glycosidic bond (Boiteux and Guillet, 2004). According to the manufacturer's manuals for SB-based conversion kits, DNA is denatured and fragmented before or during cytosine conversion in the thermocycler to chemically convert non-methylated cytosines into uracil. Thus, the prolonged exposure to SB during denaturation, as well as DNA fragmentation, likely contributed to the loss of small DNA fragments during the downstream purification steps, reducing library yield. Conversely, high recovery of DNA when using the QIAGEN conversion kit may be attributed to the presence of DNA protection buffer in the kit (Izzi et al., 2014; Leontiou et al., 2015; Tierling et al., 2018). The drastic loss of DNA with the NEB EM-seq kit may be due to comparatively more purification steps, increasing the chances of DNA loss due to pipetting.

All cytosine conversion kits yielded quality sequencing data

Prepared samples were sequenced to see if the cytosine conversion method had any effect on the sequencing. The initial quality of sequencing data was analyzed based on the following criteria: the total number of sequencing reads obtained, total coverage (X), duplication rate (%), Phred quality score, Per-base sequence, and GC content (%) in the sample. These parameters indicated no major differences among the three samples as discussed below. The number of raw sequencing reads per sample ranged from 33–50 million, with coverage ranging from 20X to 30X (Figure 3). The Phred quality score of each data set was found to be excellent, ranging between 31 to 40 (Figure 4). Based on the per-base sequence plot, thymine was the most abundant base in each library with 50% of all four bases, whereas cytosine was the least abundant (Figure 5). GC content was similar between samples, with the Zymo Research library having the highest level, followed by the QIAGEN and NEB kits (Figure 3). As expected, the per-sequence GC content for each sample was biased due to cytosine conversion into thiamine during PCR amplification. Thus, our data suggest that the cytosine conversion method has little effect on the initial sequence quality.

Figure 4. Mean quality scores of WGBS and EM-seq libraries.

The result was obtained using MultiQC. WGBS, whole-genome bisulfite sequencing; EM-seq, enzymatic methylation sequencing.

Figure 5. Per Base Sequence Content of the (a) Zymo-research kit, (b) Qiagen kit, and (c) NEB kit.

Results were obtained using MultiQC. NEB, New England BioLabs.

Nevertheless, initial data quality may not reflect data utility. Therefore, we mapped reads to the soybean genome to see the mapping efficiency (defined as how many reads can be aligned to the genome proportional to the total reads). Reads generated from all three methods had similar mapping efficiencies. The NEB kit library had the highest mapping efficiency (74.4%), followed by the QIAGEN kit (72.5%) and the Zymo Research kit (70.8%) (Figure 3). The observed mapping efficiency (>70%) is comparable with previously reported efficiencies for cytosine-converted samples (≈ 60–75%) (Hari and Parthasarathy, 2019). These results implied that sequences obtained from all three methods could be aligned correctly to the soybean genome.

We also calculated the duplication rate for each sample, as a higher rate of duplication distorts the true sequence proportion in a given sample. Duplication rates between 13 and 15% were detected for each library (Figure 3), which is lower than previously reported duplication rates during soybean methylome profiling (Rambani et al., 2015). Duplication can occur for several reasons. For instance, biased PCR enrichment and overamplification of DNA fragments can result in an overrepresentation of library fragments. Duplication can also occur when an identical template binds to numerous clusters on a flow cell. Duplicates of PCR products distort the actual proportion of sequences in the sample (please see more details on the Babraham Institute website).

Cytosine conversion efficiency

It is crucial to analyze the cytosine conversion efficiency during DNA methylation analysis. If a library has low conversion efficiency, it can result in false methylation detection (Singer, 2019). Adding methylated, non-methylated, or both types of genomic controls in a sample can help identify these limitations presented by bisulfite-converted DNA. Here, the cytosine conversion efficiency was calculated using reads aligned to non-methylated E. coli and chloroplast gDNA for each sample. According to these alignments, the bisulfite conversion efficiency of samples ranged between 90 and 99.8% (Figure 3). QIAGEN and Zymo Conversion kits had a similar conversion efficiency of roughly 99% (Figure 3). Similar results were reported for both bisulfite kits in a study conducted on human chromosomal DNA (Tierling et al., 2018). For the NEB EM-seq kit, the cytosine conversion efficiency was found to be 90.9%, around 9% lower than the data obtained from a published study (Figure 3) (Feng et al., 2020). The efficiency and accuracy of a conversion kit are dependent primarily on the PCR cycling procedure, conditions, the length of the desulphonation process, and the addition of reagents to prevent DNA degradation (Izzi et al., 2014; Tierling et al., 2018). In this case, the QIAGEN kit included a DNA protection buffer containing a pH indicator dye to ensure the correct pH for cytosine conversion. Alternatively, the Zymo Research kit protocol distinguishes the alkaline or denaturation step from the conversion step (Izzi et al., 2014). However, it is unclear what contributed to poor cytosine conversion efficiency in the NEB EM-seq kit. It is possible that further optimization of the input DNA quantity is necessary for this kit when it is used with the NEB library preparation kit.

Furthermore, the data were aligned to the soybean reference genome to detect the number of methylated cytosines in each library. The highest number of methylated cytosines in all three contexts (i.e., CpG, CHG, and CHH) was found in the library prepared with the NEB EM-seq kit, with around 1.295 × 10⁹ total methylated cytosines, followed by 1.169 × 10⁹ in the QIAGEN library and 703 million in the Zymo Research library. Moreover, the QIAGEN and Zymo Research libraries displayed a similar distribution of methylated cytosine across contexts. The NEB library showed a distinct distribution compared to the other kits, with an increase in methylated cytosines at CHH context (Figure 6). This can be attributed to the comparatively poor cytosine conversion efficiency, which may have resulted in a high number of erroneously methylated cytosines (Simpson et al., 2017). As expected, the majority of cytosines in each library were methylated in CpG context, followed by CHG and CHH (Figure 6). These data are congruent with a past soybean DNA methylation study, which revealed the greatest levels of methylation in CpG context, followed by CHG and CHH contexts, under normal developmental conditions (Song et al., 2013).

Figure 6. Percentage of methylated cytosines detected in CpG, CHG, and CHH methylation contexts in sequencing data obtained from libraries prepared using the Zymo-research kit, Qiagen kit, and NEB kit.

NEB, New England BioLabs.