Utilisation of Methylome Data to Identify Stably Unmethylated Regions in Plant Genomes

Background

Groundbreaking research into epigenetics has opened up possibilities for its application to human diseases, modern agriculture, synthetic biology, and studies of evolution. Covalent attachment of a methyl group to the 5′ carbon of cytosine in DNA (5-methylcytosine) is generally known as DNA methylation, and 5-methylcytosine is sometimes referred to as the fifth base [1]. DNA methylation can provide a mechanism for epigenetic inheritance of phenotypes and, accordingly, is often referred to as an epigenetic modification. DNA methylation plays critical roles in transposon silencing, genome stability and organisation, heterochromatin formation, gene regulation, development, and imprinting [2]. More recently, we have demonstrated that there is great utility in identifying regions of the genome that lack DNA methylation, referred to as unmethylated regions (UMRs) [3]. In this protocol, we provide a step-by-step guide to identify UMRs from DNA methylation sequencing data.

There are a variety of technologies that can be used to detect and quantify the levels of DNA methylation at a particular locus, and many of these methods can also be scaled to genome-wide profiling of the entire methylome [4,5]. Technologies include bisulfite based, digestion based, or affinity based. Whole-genome techniques include whole-genome bisulfite sequencing (WGBS), enzymatic methyl-seq (EM-seq), reduced representation bisulfite sequencing (RRBS-seq), which uses bisulfite coupled with enzymatic digestions, and methyl-DNA immunoprecipitation (MeDIP), which attracts methylation through antibody enrichment [6–9]. In addition, over the last few years, nanopore sequencing has emerged as another alternative, allowing direct detection of DNA and RNA methylation on long reads [10]. Importantly, the ability to identify methylation at single-base-pair resolution has allowed specific understanding of how methylation—or the lack thereof—influences specific regions of the genome, including regulatory regions that may be important for transcriptional changes and phenotypic variation for traits of interest. In this protocol, we analyse WGBS sequence data; we routinely apply the same analysis pipeline to EM-seq data, to which it is equally applicable. Together, these methyl-seq whole-genome sequencing approaches (WGBS and EM-seq) are the current gold standard for DNA methylation profiling, as they provide a whole-genome analysis of regions that are methylated and the types of methylation present [11,12].

The key step in WGBS involves treating genomic DNA with sodium bisulfite, which converts unmethylated cytosines to uracils by deamination (which are converted to thiamine following PCR amplification), while 5’-methylcytosines remain unaffected [8]. After bisulfite treatment is complete, converted DNA can be prepared for sequencing, commonly on the Illumina platforms, allowing inference of methylation at single-base resolution by comparison to a reference genome. Due to the cytosine conversion step in WGBS (and in EM-seq), a key requirement in the bioinformatic analysis is to use bisulfite aware mapping software to identify the T-converted unmethylated cytosines, which will appear as polymorphisms compared to a reference genome. Commonly used mapping software includes BSMAP, BWA-meth, and Bismark [13,14].

In plants and animals, DNA methylation occurs in different sequence contexts including CG, CHG, and CHH (where H is A, T, or C). Studies have found that the type of methylation can correspond to different functions; for example, CG-only methylation has been correlated with actively transcribed gene bodies, while transcriptionally silenced regions are associated with high levels of methylation in all contexts (CHH, CHG, and CG) [15]. The output of this analysis pipeline includes files for visualisation of each type of methylation. Recently, the research community has demonstrated that regions that lack DNA methylation in all contexts are of particular significance [3,15–23]. Partitioning genomic regions into different categories of methylation types allows us to identify unmethylated regions (UMRs). This approach can be very useful for annotating a genome [24] because UMRs tend to align with regions containing functional genes and also cis-regulatory elements (CREs) [3,25]. CREs are non-coding elements and include enhancers, promoters, and silencers, which can influence gene expression. These have the potential to be important targets for selection and breeding and also as targets for genetic engineering. Researchers have yielded promising results from genetic modification of CREs for tailoring gene expression without causing developmental problems [26,27]. For example, improvements in yield in rice [28,29], maize [30], and tomato [31–33] have been produced by gene editing non-coding cis-regulatory promoter alleles. However, efficiently identifying functional regions of regulatory importance is challenging [3]; identification of UMRs can narrow the search for genetic targets. As genomic sequencing and tools for genome annotation improve, our ability to understand the functionality of the genome is made increasingly powerful when combining knowledge of transcription binding sites, conservation of sequences through evolution, epigenomic markers for transcriptional activity, and other emerging technologies.
In the case study outlined below, we have provided a subset of raw WGBS reads extracted from SRX5532987, a published study of DNA methylation in maize leaf tissue [3]. This subset of reads aligns to a small 20 Mb region of maize chromosome 1; we provide this subset of the genome as a reference sequence for mapping. These files enable processing this example data with minimal computational resources and could be completed on a basic laptop computer with the appropriate software installed. The pipeline was originally optimised for maize, but it is generally applicable to other species without modification [34]; however, users could consider modifying the thresholds in the UMR-calling step if applied to a plant genome that has particularly unusual levels or distribution of DNA methylation. The key steps of the workflow are outlined in Figure 1.

Overview of the workflow to identify unmethylated regions (UMRs) in a plant genome. The bioinformatic workflow presented in the protocol includes six major steps: 1) read trimming, 2) read mapping, 3) filtering, 4) extraction and quantification of methylation levels per cytosine, 5) data visualisation, and lastly 6) summarisation and identification of UMRs. The output file formats are indicated above each step and the software tools and their purpose are summarised.