Comparative restriction enzyme analysis of methylation (CREAM) reveals methylome variability within a clonal in vitro cannabis population

2.2. CREAM libraries preparation

Two sequencing libraries were prepared in parallel for each sample with the extracted DNA using the Comparative Restriction Enzyme Analysis of Methylation (CREAM) approach (Figure 1) at the Plateforme d’Analyses Génomiques² at the Institut de Biologie Intégrative et des Systèmes (IBIS) of Université Laval, Quebec, Canada. The CREAM method builds on the 3D-GBS approach (de Ronne et al., 2023). Briefly, in both libraries, DNA molecules were cleaved at one end either by the restriction enzymes NsiI (New England Biolabs) or PstI (New England Biolabs), which have distinct restriction sites (5’-ATGCA/T-3’ and 5’-CTGCA/G-3’, respectively) to anchor DNA fragments at specific locations in the genome. Since they recognize specific sequences of six nucleotides (six-cutters) and each nucleotide has one of the four nucleic bases, these enzymes cut at a theoretical frequency of 4⁶ or 4096 base pairs (bp). The combined use of the enzymes reduces this theoretical frequency by half, anchoring the DNA fragments at every 2048 bp and thus covering a larger part of the genome. The actual frequency of restriction sites varies between species and within a same species and is mostly influenced by the percentage of GC bases (Roberts and Murray, 1976; Lindsay and Bird, 1987). The other end of the DNA fragments was cleaved by either MspI (New England Biolabs) or HpaII (New England Biolabs), depending on the library. These enzymes are isoschizomers and share the same restriction site (5’-C/CGG-3’) but have different sensitivity to full DNA methylation on the second cytosine of the recognition sequence (at the CpG site), MspI being insensitive and HpaII being sensitive to methylation (Waalwijk and Flavell, 1978). As illustrated in Figure 1A, the variation in sensitivity of the restriction enzymes to DNA methylation results in the generation of fragments of different lengths within the same genomic region. All four possibilities expected from this comparative analysis of the restriction fragments in both libraries, based on the filtering of the fragments with the size selection step, are represented in Figure 1B. After digestion, the sample-specific barcodes and universal adapters were ligated. Since the cohesive ends for the MspI and HpaII restriction sites are identical, the same adapters were used in both libraries. A size selection step using a BluePippin apparatus (Sage Science, Beverley, MA, USA) was performed to capture digested fragments of 200-400 bp. Finally, DNA libraries were amplified by PCR and sequenced with an Illumina NovaSeq 6000 System at the Centre d’expertise et de services Génome Québec (Montreal, QC, Canada), generating 204 M and 173 M paired-end reads of 150 bp for the MspI and HpaII libraries, respectively.

2.3. Bioinformatics pipeline

2.4. Accumulation and distribution of methylated positions

To examine the accumulation of methylated loci in the population (all clones from both sister lines) as the number of subcultures increased, a Kruskal-Wallis test (one-way ANOVA on ranks) was conducted. This test was chosen since the assumptions for conducting an ANOVA were not met in this case. Then, a principal component analysis (PCA) was performed to assess the distribution of methylated positions within the population. The PCA aimed to determine clustering patterns in the samples based on the methylated loci using the R packages FactoMineR and factoextra (Lê et al., 2008; Kassambara and Mundt, 2016).

The distribution of methylated regions across the cannabis genome was visualized using the RIdeogram R package (Hao et al., 2020). This also included the gene density information obtained from the NCBI Gene table for the cs10 reference genome⁶ (accessed February 15, 2023). A Spearman’s rank correlation coefficient was calculated, using the cor command from the R package stats (R Core Team, 2022), to assess the monotonic relationship between the gene density and the methylation density across the genome since the data for both variables were skewed towards 0. Other visualizations were generated with ggplot2 in R (Wickham, 2011).

3.2. Identification of methylated positions

From the mapped paired-end reads, a bioinformatics analytical pipeline was developed to determine the methylated positions captured with the CREAM approach. Briefly, the pipeline includes quality checks for inserts (composed of a pair of reads) and calls loci that encompass all potentially overlapping inserts within a region (See Supplementary Figure 2 for details). This methylotyping pipeline successfully divided the inserts from the CREAM libraries into the four expected possibilities, resulting in a total of 5,235 loci (See Figure 1B for the four expected categories). Of these confidently called high-quality loci, 3,762 (71.86%) were identified as methylated regions, while 1,473 (28.14%) were identified as unmethylated regions with perfectly overlapping inserts in both libraries (possibility A). The most common type of captured methylated regions (2,949) were loci with an insert only in the MspI library (possibility B), accounting for 56.33% of all captured loci. Loci with an insert only in the HpaII library (possibility D) accounted for approximately 15.53% of all captured loci (813). Finally, a total of 186 loci were excluded from the analysis due to challenges in accurately categorizing them into the predetermined possibilities. These included methylated regions composed of loci with inserts of different length in both libraries (possibility C), which were far less common and loci that exhibited ambiguous characteristics or lacked clear patterns for classification, rendering them less reliable for further analysis.

3.3. Distribution of methylated positions in the population

The methylotype of each sample for the list of called loci was used to assess the variability in methylated positions across the clonal population. The frequency of the methylated loci revealed four main peaks (Figure 3A). The largest peak consisted of the 1,490 (39.61% of total methylated loci) monomorphic methylated loci, which are DNA methylation positions captured in all 68 remaining samples, indicating a shared methylotype across the population (all clones from both sister clonal lines). On the other hand, among polymorphic loci, 427 (11.35%) were unique to a single sample, indicating specific methylation patterns within individual samples. The two other peaks in the distribution of the number of methylated loci corresponded to the number of clones derived from the two sister lines. Specifically, there were 155 loci unique to the 13 samples in the AT-4 line and 117 loci unique to the 55 samples in the AC-150 line, indicating distinct methylotypes within the sister clonal lines.

In this study, we also investigated whether the number of subcultures that the clonal lines went through in tissue culture had an impact on the number of DNA methylation positions captured. Although a significant methylome variability within clonal lines was observed, we did not find a significant correlation between the number of subcultures and the total number of methylated loci (p-value = 0.19; one-way ANOVA on ranks (Kruskal-Wallis test)) (Figure 3B). This suggests that the variability in DNA methylation within clonal lines was not influenced by the number of subcultures the plants went through in this study. Finally, the principal component analysis (PCA) performed on the methylated loci (Figure 3C) revealed that the two principal components explained a significant portion of the variation, accounting for 81.94% and 10.99% of the total variation, respectively. The samples displayed clear clustering patterns, with two distinct groups corresponding to the clones derived from the two sister lines. Within each cluster, the variability was relatively limited, although some samples showed variation compared to overall trend of their respective group.

4.1. CREAM: a new efficient and cost-effective methylotyping method

The CREAM method was successful in identifying significant variation in methylotypes among cannabis clones, specifically distinguishing the two subpopulations derived from two sister lines. This success can be attributed to the generation of high-quality sequencing fragments, which facilitated the accurate identification of methylated and unmethylated loci using our newly developed methylotyping pipeline. Importantly, the CREAM method offers a cost-effective solution compared to the main approaches based on RRBS, like epiGBS (Van Gurp et al., 2016) and bsRADseq (Trucchi et al., 2016), with an estimated sequencing cost per sample of around 30$, representing a 30-80% decrease in costs (Werner et al., 2020). Additionally, the CREAM approach avoided the extensive DNA damage typically associated with the bisulfite conversion in other techniques (Tanaka and Okamoto, 2007), preserving DNA integrity. Other large-scale methylation detection approaches have been developed to investigate methylation patterns on a large scale, such as Nanopore sequencing (Tourancheau et al., 2021) and whole-genome bisulfite sequencing (Urich et al., 2015). These established methods have demonstrated enhanced success in identifying methylation patterns (Ni et al., 2021). Additionally, ongoing efforts to enhance methylation calling (Ni et al., 2021) have further improved detection rates. However, while these methods may offer a thorough and insightful evaluation of methylated regions across the genome, they are still generally considered expensive to conduct studies at a population level (Suzuki et al., 2018; Kong et al., 2023).

The CREAM method achieved a satisfactory average genome coverage of ~0.4% and a mean depth of coverage of around 100X. In comparison, well-regarded RRBS methods such as epiGBS (Van Gurp et al., 2016) and its modified version (Werner et al., 2020) covered 0.37% of the 135 Mb of the A. thaliana genome and 0.28% of the 246 Mb almond genome, respectively. When examining the genome coverage and mean depth of coverage specific to each library (Figure 2), it becomes evident that the genome coverage tends to be higher in the MspI library compared to the HpaII library. This observation is consistent with the higher proportion of loci with fragments found only in the MspI library (category B) among the captured regions. Therefore, the differences in genome coverage and the number of loci found between the two libraries can be attributed to the variations in the efficiency of the two restriction enzymes (methylation sensitive and methylation insensitive) used in the CREAM method.

The CREAM method builds on methylation-sensitive amplified polymorphism (MSAP) approaches coupled with sequencing technologies, which aim to provide a direct genomic localization of differentially cleaved DNA fragments using methylation-sensitive restriction enzymes. Typically, DNA fragments are cleaved at one end by the EcoRI enzyme and at the other end by one of the two isoschizomers MspI or HpaII (restriction enzymes sharing the same restriction site, 5’-CCGG-3’). The presence and absence of EcoRI/MspI (unmethylated) and EcoRI/HpaII (methylated) fragments allows methylation patterns to be determined in a comparative manner (Fulneček and Kovařík, 2014), and their sequencing allows for their genomic localization. Such methods, like MSAP-Seq (Chwialkowska et al., 2017) and others (Baránek et al., 2016) have achieved the direct analysis of DNA methylation in plants using only EcoRI as an anchoring enzyme in combination with MspI and HpaII. Using a permissively low coverage of 2 reads per sample, the MSAP-Seq was able to capture a total of 190,000 CCGG sites for DNA methylation analysis in barley, which has a genome of about 5.3 Gb (around 6 times the size of the cannabis genome) (Chwialkowska et al., 2017). The CREAM method focuses on higher quality CCGG sites, with a coverage of at least 20, to avoid artefacts in captured methylated sites. It also utilizes two anchoring enzymes (NsiI and PstI) to capture more diverse restriction fragments, therefore providing a better overview of methylation sites across the genome. While it provides interesting elements to the family of MSAP methods coupled with sequencing technologies, the CREAM approach also only focuses on CCGG sequences recognized by restriction enzymes, thus ignoring DNA methylation in CHG and CHH contexts, notably involved in the regulation of genes and transposons (Agius et al., 2023).

The choice of the restriction enzymes in the CREAM method has a direct impact on the number of loci captured and their distribution across the genome. In our findings, we observed a notable enrichment of fragments within genic regions. To explain the abundance of captured fragments, we propose the influence of GC content as a contributing factor. In the case of the cannabis reference genome (cs10), it exhibits an overall GC content of 33%, whereas the GC content differs among the restriction sites of the restriction enzymes. Specifically, the GC content of the NsiI and PstI restriction sites is 66% and 33%, respectively, while the MspI and HpaII restriction sites have a GC content of 100%. This indicates that the distribution of restriction fragments is not uniform across the genome and tends to be more concentrated in regions that are relatively richer in GC than AT. In plants, these GC-rich regions are predominantly associated with genic regions (Serres-Giardi et al., 2012; Glémin et al., 2014). This supports the observed correlation between the gene density and methylation density in the data (Figure 3D). Therefore, the CREAM approach might result in an overrepresentation of genic regions due to their higher accessibility to the restriction enzymes.

While our approach offers several advantages, it is important to acknowledge its limitations. First, its resolution is currently limited. It provides a binary methylation status (either 0 or 1) for a given cytosine in a locus, rather than providing a methylation quantitative score at a base level. Although this binary representation still provides valuable insights, a finer resolution would be desirable to gain a more comprehensive understanding of the methylation landscape across the genome. Second, we encountered challenges in classifying 186 captured loci (3.43% of all loci) into the expected categories (A, B, C or D). The comparative analysis of the libraries was challenging for these specific restriction fragments because of the complexity and diversity of their alignment. Lastly, it is important to note that our approach did not operate at a single-cell level. As a result, the methylation patterns obtained from a single sample represent a mixture of methylation patterns from various cells, introducing some level of heterogeneity. However, these limitations might be overshadowed by the compelling low-cost and high-throughput features of the approach.

4.2. Methylation variation in cannabis clones

Upon analyzing the DNA methylation frequency in the population, we observed contrasting patterns between the two subsets of clonal lines. These patterns can be used as methylation tags or methylation fingerprints, allowing us to clearly differentiate individual clones within their respective clonal line. Notably, the methylation patterns exhibited by each clone strongly reflect their clonal lineage, suggesting a unique and identifiable DNA methylation profile for each clonal line. This is consistent with various studies of methylation patterns in tissue culture. For example, in maize, stable changes in differentially methylated regions (DMRs) were observed between independently regenerated clonal lines, highlighting potential “epigenetic footprints” unique to certain lines and maintained over several generations (Stelpflug et al., 2014). A further study has shown consistent (monomorphic) and rare (polymorphic) DMRs maintained across generations that distinguished tissue culture-derived plants, providing more evidence for specific DNA methylation patterns acting as “tags” (Han et al., 2018). It is likely that these patterns originated from the sister seeds that were used to initiate clonal lines and were preserved to some extent throughout multiple generations of subculturing. In sexual reproduction, such as seed production, various mechanisms exist to limit defective and potentially problematic epigenetic variations from one generation to the next, including resetting of DNA methylation state (Quadrana and Colot, 2016). However, clonal propagation, being asexual, lacks these reset processes and is expected to affect DNA methylation stability across the genome (Ibañez and Quadrana, 2023), leading to potential slight variations in methylation patterns.

Different DNA methylation dynamics have been observed in the vegetative state. For instance, high levels of CHG methylation in Arabidopsis thaliana shoot apical meristematic stem cells (Gutzat et al., 2020) and low levels of CHH methylation in species with extensive clonal propagation histories (Niederhuth et al., 2016) have been reported. However, trends regarding CG methylation levels in tissue culture are less conclusive and differ among species, with increasing levels reported in gentian (Fiuk et al., 2010), bush lily (Wang et al., 2012), banana (Peraza-Echeverria et al., 2001) and tomato (Smulders et al., 1995), decreasing levels observed in triticale (Bednarek et al., 2017), barley (Li et al., 2007), grapevine (Baránek et al., 2010) and Freesia (Gao et al., 2010), and no significant changes detected in pea (Smýkal et al., 2007) or apples (Li et al., 2002). Since the cannabis methylome is poorly understood, our results represent an encouraging first report of methylation patterns for this economically important crop in a tissue culture context.

Tissue-specific differences in DNA methylation can also be significant (Lloyd and Lister, 2022). Cells from stem tissues are generally less studied and methylation patterns are less decisive, making comparisons between studies and species more intricate. While no methylome studies have been done in cannabis, it is still interesting to present the results of this current study in the context of the closest work achieved in terms of species of interest, DNA methylation capture approach and stem tissue studied in a tissue culture propagation system. In a study conducted in hops (Humulus lupulus L.), a close relative to C. sativa (Kovalchuk et al., 2020), micropropagated stem tissues (branches) were tested for methylation changes (Peredo et al., 2009). The study found that most the methylated loci (56.34%) were monomorphic, meaning they were shared by all 80 clones obtained from two cultivars, while 13.24% of the variation were unique to individual samples (referred to as singletons). These results are consistent with the distribution of methylated loci observed in the present study (Figure 3A), with 41.22% and 11.81% of the total loci being monomorphic and unique to a single sample, respectively. The difference in proportions between the methylation patterns observed in our study and those reported in the MSAP method used in the hops study (Peredo et al., 2009) could be attributed to the higher genomic resolution of our approach and the variability between the species. Further methylome studies in cannabis would also benefit from tracking methylation patterns in several tissues, such as female flowers that accumulate many compounds of interest.

The lack of significant differences in the number of methylated loci with the number of subcultures may be attributed to the dynamic coordination of processes involved in establishing, maintaining, and removing specific DNA methylation states (Zhang et al., 2018). It has been reported that while the total number of methylated positions may not change significantly, their specific identity is likely to vary (Zhang et al., 2018). This suggests that the stability of methylation patterns is established in the earlier stages of production. The stability of DNA methylation patterns during subculturing has been observed in other plant species as well. For example, in garlic, methylation patterns were found to stabilize after 6 months of micropropagation (Gimenez et al., 2016). However, it is important to note that the trends in DNA methylation variations are highly species-dependant. Some studies have reported increases in methylation levels during subculturing (Rival et al., 2013; Fraga et al., 2016), while others have observed decreases (Xu et al., 2004; Huang et al., 2012; Machczyńska et al., 2014) over different subculturing durations. There have also been cases where methylation levels initially decrease in early generations and then recover over longer periods of time, as shown in sweet orange after 30 years (Wang et al., 2022). Considering that our study is the first investigation of the cannabis methylome, it is plausible that methylation patterns stabilized in our clonal population over several subcultures. All clones were acclimatized to tissue culture and could have reached a steady state by the 6^th to 10^th subculture. Further research is needed to gain a better understanding of the dynamics and stability of DNA methylation patterns in cannabis and how they may influence the phenotypic characteristics of clonal plants. A focus should be made on the induction phase in tissue culture, where plants are transitioning from regular growth conditions to in vitro conditions, which could affect methylation patterns.

This work has previously appeared online as a preprint (Boissinot et al., 2023).