The cannabis clonal population used in this study was initiated in March 2019 at the University of Guelph (Ontario, Canada). It was developed using two sister lines (seedlings) derived from a cross of a cultivar exhibiting an indica-leaning growth habit and THC and cannabidiol (CBD) levels of approximately 13% and < 1%, respectively (Adamek et al., 2023). Nodal explants from the seedlings were subcultured in vitro (about every month) and maintained on DKW Basal Medium with Vitamins (Product ID D2470; Phytotechnology Laboratories, Lenexa, Kansas, USA), 1 mL/L plant preservative mixture (PPM; Plant Cell Technology, Washington, DC, USA), 0.6% agar (w/v) (A360-500; Fisher Scientific, Fair Lawn, New Jersey, USA) and pH adjusted to 5.7 to generate two clonal lines (Page et al., 2021). A total of 78 individual plants considered as clones (from this existing population of two sister clonal lines) of the same chronological age but maintained with different subculture frequencies (i.e., number of subcultures ranging from 6 to 11 every 45 days on average depending on the existing material) were selected (Supplementary Figure 1). The two sister lines that composed the population, named AC-150 and AT-4, contained 63 and 15 clones that were analyzed, respectively. The maternal material (the first clones) has been destroyed to produce new explants as part of the subculture process in the existing population. DNA samples from the 78 clones were extracted from stem tissue collected at the same time and from the same stem regions (middle of the explant) using a Qiagen DNeasy Plant Mini Kit according to the manufacturer’s protocol. The DNA concentration of each sample was quantified using a NanoDrop One spectrophotometer (Thermo Scientific, Waltham, MA) and then diluted to 10 ng/μl. A volume of 10 μl containing 100 ng of DNA was used to prepare each sequencing library for each sample.

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.

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).

A gene ontology (GO) analysis was performed to identify significant GO terms affected by DNA methylation captured with the CREAM approach. To overcome challenges in matching protein IDs and GO terms with the annotations of the cs10 cannabis reference genome, a combination of the GAWN v0.3.5⁷ and go_enrichment⁸ pipelines was used. The GAWN pipeline annotated the cs10 reference genome using the available transcriptome from NCBI and found all the methylated loci and the captured loci (unmethylated and methylated) within ± 1 kb of transcripts. The GO enrichment analysis was then performed with the go_enrichment pipeline using the list of methylated loci adjacent to transcripts as the target and considering all captured loci as the background gene set. GO terms with a significant adjusted p-value of p < 0.10 (Benjamini-Hochberg false discovery rate correction for multiple testing, FDR BH) were kept for further exploration of biological processes.

We have developed a low-cost and efficient method for assessing methylome variability at the population level. This method, CREAM, involves digesting DNA samples with restriction enzymes of different sensitivity to DNA methylation. Libraries were developed and sequenced for 78 cannabis clones produced in tissue culture for which no phenotypic data were measured. This has yielded an average of >188 M paired-end reads per sequencing library (Table 1). The demultiplexing, trimming and alignment of the reads to the reference genome led to an average of >2 M paired-end reads per sample.

Out of the initial 78 samples, 10 samples (8 and 2 samples from the AC-150 and AT-4 clonal lines, respectively) were excluded from the analysis as they yielded less than 100,000 reads on average per library. The genome coverage and mean depth of coverage were computed for the remaining samples and compiled for each library (Supplementary Table 1). On average, we captured ~0.4% of the cannabis genome with a mean depth of coverage of ~100X across the captured regions, indicating a sufficient depth to ensure reliable and accurate analysis of the captured regions. As shown in Figure 2, the genome coverage tends to be higher in the MspI library, while the HpaII library exhibits a higher mean depth of coverage. This disparity can be attributed to the fact that the HpaII library does not capture DNA fragments with DNA methylation. Consequently, the HpaII library contains fewer fragments for a comparable sequencing effort compared to the MspI library, covering smaller proportion of the genome with higher coverage.

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.

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.

We examined the distribution of the methylated positions across the genome to understand the patterns of DNA methylation. Precisely, we calculated the gene density, which represents the number of genes present in each bin of 500 kb in the cannabis genome, along with the methylation density, which indicates the number of captured methylated positions within the same 500 kb bin (Figure 3D). We then compared the gene density and methylation density to understand the relationship between DNA methylation and gene distribution in different genomic regions. The density of captured methylated loci across the genome showed a strong positive correlation with the density of genes across the same genomic bins (Spearman’s rank correlation coefficient = 0.6473), suggesting the positive monotonic relationship between the two variables. We also observed that regions of high genic density were also enriched in captured methylated loci. However, the observed correlation between gene density and methylation density could be influenced by the approach used, which tends to capture loci in genic regions (Supplementary Figure 3). Finally, the gene ontology (GO) analysis was performed to identify genes potentially affected by the DNA methylation patterns captured within a range of ± 1 kb of gene transcripts (Table 2). The analysis utilized all the loci captured by the CREAM approach as a background set for comparison. As a result, 11 GO terms were found to be significant. Most importantly, 9 out of the 11 significant biological process GO terms were related to metabolic processes, indicating a potential influence of DNA methylation on important biochemical pathways.

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.

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.

Upon examining the distribution of DNA methylation in the cannabis genome, we observed a correlation between methylation density and gene density along the chromosomes. As mentioned earlier, this suggests that the captured loci in this study may be biased towards genic regions due to their GC-rich content. To a lesser extent, non-genic regions were also captured with the CREAM approach. These genomic regions rich in transposable elements and repetitive DNA sequences have also been associated with higher DNA methylation rates in A. thaliana, both in heterochromatin and euchromatin (Zhang et al., 2006). It is therefore expected to observe methylated loci in these regions. With the current limited knowledge of the cannabis genome, loci found near genes still represent the most interesting regions to study for their immediate potential functional impact on the plants biological processes.

In this study, the GO analysis within a range of ± 1 kb around gene transcripts resulted in 11 significant GO terms. The effect of gene body methylation (gbM) on gene expression and phenotypic changes is still a subject of debate and not fully understood. However, the association between gbM and gene expression is recognized as important, and it represents a potential source of variation that can be subject to natural selection (Muyle et al., 2022), highlighting the importance of the significant GO terms found in this study. Of these 11 GO terms, 9 are related to metabolic processes, which encompass a range of biochemical reactions involved in the synthesis, breakdown, and transformation of various molecules within the plant. This suggests that DNA methylation may play a regulatory role in modulating the expression and activity of genes involved in key metabolic pathways in cannabis. To gain a deeper understanding of the cannabis methylome, it would be crucial to validate the biological effects of DNA methylation on the transcriptional activity of these genes. Such validation experiments would provide important insights into the functional implications of DNA methylation in cannabis and its potential impact on metabolic processes.

In conclusion, the present study addresses a critical gap regarding somaclonal variation and epigenetic factors in C. sativa. To the best of our knowledge, this is the first study to investigate the cannabis methylome and the first of its kind in the field of cannabis tissue culture. The primary objective of this study was to develop a novel and cost-effective approach to detect methylation variation in C. sativa clones derived from tissue culture. The results revealed significant variation in methylotypes among the population of cannabis clones, indicating the presence of methylation footprints between sister clonal lines and offering valuable insights into the epigenetic landscape of this important crop. Importantly, this approach overcomes the cost and technical challenges associated with existing methods and enables the high-throughput analysis of methylome at a population level. This knowledge has implications for crop improvement and the development of sustainable production systems in the cannabis industry. By highlighting the dynamics of methylome and its correlation with gene density, this study advocates for further investigations on the regulatory role of DNA methylation in key metabolic pathways of cannabis. Functional validation of the observed methylation patterns and their impact on gene expression will undoubtedly contribute to a more comprehensive understanding of the cannabis methylome.