Development and validation of a minimal SNP genotyping panel for the differentiation of Cannabis sativa cultivars

Plant material and growing conditions

Fresh leaf material from 31 clonally propagated C. sativa cultivars was provided by Organigram Inc., a Health Canada registered licensed producer of cannabis (Moncton, New Brunswick, Canada). Fresh leaf material from a hemp variety, ‘Anka’, was obtained from plants grown at Université de Moncton (‘Anka’ is a variety from UniSeeds, and seeds were obtained from Céréla, Saint-Hughes, Québec, Canada).

Seeds were germinated by submerging them in distilled water at room temperature overnight, then leaving them on a damp paper towel enclosed in a sealed plastic bag until germination. Germinated seeds were sown in Pro-mix Mycorrhizae BX (PremierTech Horticulture, Rivière-du-Loup, QC). The plants were grown in a Conviron PGR15 (Controlled Environments Limited, Winnipeg, MB) at 60% relative humidity under an 18:6 light-dark photoperiod cycle, at 23^oC and 22^oC during the day and night respectively. Lighting was a combination of fluorescent and incandescent bulbs, providing approximately 360 µmol/m²/s at the canopy. Plants were grown under these conditions for six weeks. Leaves from two individuals per cultivar were used to gather leaf punches from fresh samples for DNA extraction. Leaf punches were pooled by cultivar before being weighed for extraction. In fresh leaf samples provided by Organigram Inc., both individuals for each cultivar were clones propagated from a mother plant. In the fresh leaf samples harvested from ‘Anka’, both plants were grown from seed.

Nine dried flower samples of three commercially available cultivars from three different LPs each (Organigram Inc. and two other LPs per cultivar, nine samples total) were obtained through Cannabis NB, the crown corporation responsible for distributing cannabis products in New Brunswick, Canada. Dried flower samples were composed of dried flowers from one dried flower sale unit (i.e., 1-gram or 3.5-gram container) per cultivar. It was not possible to confirm whether the dried flowers samples were composed of flower from more than one individual per sale unit. Table 1 summarizes the sample origin, test condition, primary cannabinoid type, and source of cannabinoid data.

DNA extraction

DNA was extracted from both fresh leaf tissues and dried flower samples with the DNeasy ^® Plant Mini Kit (Qiagen, Hilden, Germany). For fresh leaf tissues, the standard protocol was followed using samples with masses ranging from 90 to 100 mg but with a doubling of the recommended RNase A volume (8µL instead of 4 µL). Three extractions were performed using three samples per cultivar. For dried flower samples, the lyophilised material protocol was used with a sample mass of 19–25 mg but with a doubling of the recommended RNase A volume (8µL instead of 4 µL). Three extractions were performed using three samples per dried flower sale unit. DNA quality was assessed using a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific), following manufacturer specifications. DNA yield and quality is reported in Supplementary Table 1.

Library preparation and sequencing

Whole Genome Sequencing DNA libraries were prepared and sequenced on 11 lanes, using the Illumina HiSeq X v4 technology (PE 150 bp) (Illumina, San Diego, CA, USA) at the Centre d’expertise et de services Génome Québec (CESGC; Montreal, QC, Canada). The CESGC generated the libraries with the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, USA) using the manufacturer’s recommendations. Adapters and PCR primers were acquired from Integrated DNA Technology (IDT; Coralville, IA, USA). Size selection was performed using SparQ beads (Qiagen), and libraries were quantified using the KAPA Library Quantification Kits– Complete Kit (Universal; Kapa Biosystems, Wilmington, MA, USA). The LabChip GX (PerkinElmer, Sheldon, CT, USA) was used to assess average fragment size. Libraries were normalized, pooled, then denatured in 0.05 N NaOH and neutralized using HT1 buffer (Illumina). ExAMP (Illumina) was added to the mix following the manufacturer’s instructions. The pool was loaded at 360 pM on a cBot (Illumina) and the flowcell was run on a HiSeq X for 2 × 151 cycles (paired-end mode), using a PhiX library (Illumina) mixed with libraries at 1% as a control. The Illumina control software was HCS HD v. 3.4.0.38 and the real-time analysis program was RTA v. 2.7.7 (Illumina). BCL2FASTQ v2.20 (Illumina) was then used to demultiplex samples and generate fastq reads.

Marker identification pipeline

Scripts to identify the markers used in this study are available at https://github.com/alexcull/SNPGenotyping (doi.10.5281/zenodo.10933656). Filtering, alignment and variant calling was completed using the Digital Research Alliance of Canada’s computing resources.

Initial filtering of the raw reads provided by CESGC was done using FASTP version 0.20.0 [51]. Filtered sequencing data for the 32 C. sativa cultivars were aligned to reference genome cs10, built using ‘CBDRx’, a CBD dominant cultivar, using Bowtie 2 version 2.3.5.1 [38, 52, 53]. Variants were called using SAMTOOLS version 1.10 and BCFTOOLS version 1.10.2 [54]. Variant Call Format (VCF) files were then merged, normalized, multiallelic sites were split to biallelic, and were then sorted with BCFTOOLS version 1.10.2 [55]. Allele frequencies, SNP_IDs and tags were generated using BCFTOOLS version 1.10.2. The VCF was then filtered using BCFTOOLS version 1.10.2, VCFTOOLS version 0.1.16 and custom scripts [55]. Python package pandas was used to sort the highly filtered SNPs by Minor Allele Frequency (MAF) [56, 57]. MAF is a value between 0 and 0.5 that denotes the proportion of the population carrying the second most common allele for a given SNP in a population. SNPs with high MAF values were extracted using VCFTOOLS version 0.1.16 to generate subsets of SNPs. Genotyping data for the SNPs of each subset was then used to conduct Identity-by-state (IBS) analyses using the SNPRelate version 1.18.0 package in R version 4.0.0 [58, 59]. Genetic distance matrices were created from the IBS matrices and were then tested using custom Python scripts to determine the discriminatory power of the subsets. The set of 4,318 high-quality SNPs were split into 10⁶ bp bins according to their position along the cs10 reference genome. Counts per bin were charted in Microsoft Excel to obtain a distribution of SNP density. To identify features, 10⁵ bp regions with the highest density of SNPs were explored in the National Center for Biotechnology Information (NCBI) genome data viewer (cs10 genome) browser [38, 60].

rhAmp assay design

A subset of 20 SNPs was retained to test genotyping results in vitro. 200 bp flanking sequences for these SNPs were generated using the getfasta function from BEDTOOLS version 2.29.2. Primers targeting these sequences were tested for specificity using NCBI Primer Basic Local Alignment Search Tool (BLAST) on the cs10 genome [38, 61, 62]. IDT’s rhAmp SNP genotyping system was used to validate the custom assays. IDT’s genotyping design tool was used to generate rhAmp assays, and outputs were compared in terms of chemistry and thermodynamics (% GC, length, T_m, primer dimers) using IDT’s OligoAnalyzer Tool [63, 64]. IDT synthesized the rhAmp genotyping assays. Primer details are provided in an additional file (see Additional File 1).

Validation

Validation of SNP genotyping assays was performed using rhAmp SNP Assays and rhAmp Genotyping Master Mix (www.idtdna.com/rhAmp-Genotyping). 10 µL reactions using 2 µL of 5ng/µL gDNA diluted in sterile H₂O and 8 µL of reagent mix were used for genotyping. PCR amplification and data capture were performed using CFX Connect Real-Time PCR Detection System (BioRad, Hercules, CA, USA). The DNA used for validation was extracted from the same samples and at the same time as the DNA submitted for sequencing. PCR amplification was conducted using the IDT protocol: An initial 10 min denaturing step at 95˚C was followed by 40 amplification cycles with a 10 s denaturing at 95˚C, a 30 s annealing at 60˚C for 30 s, and a 20 s extension at 68˚C. Each assay was tested three times per sample. Genotyping for dried flower samples proceeded as described above.

Allelic discrimination data was exported from CFX Connect software (BioRad) for analysis in Excel (Microsoft, Redmond, WA, USA). Call data was encoded based on the relationship between the two sets of relative fluorescent unit values. Calls were encoded as “heterozygous”, “homozygous for reference allele”, “homozygous for alternate allele” or “no call”. Call consistency and accuracy was determined in Excel. Call consistency represents the rate at which calls match between replicates. Call accuracy is the rate at which the majority of replicates (i.e., 2 of 3 or 3 of 3 replicates) match the genotype as predicted by the sequencing data.

Encoded genotype data produced in Excel was exported as a comma separated value (.csv) file, which was used for analysis in R versions 4.0.0 and 4.3.0 [59]. The dist.gene() function from the R package ape versions 5.4 and 5.8 was used to generate genetic distance matrices between individuals [65]. The discriminatory power of the subsets was tested using the custom Python script as described in the Marker identification pipeline section [66]. The most parsimonious marker set, or minimal marker set, was determined by finding the smallest group of SNPs that could still produce a genetic distance matrix with a non-null genetic distance between each cultivar pair. Neighbour-joining trees were produced using the genetic distance matrices and the BIONJ algorithm from the R package ape versions 5.4 and 5.8, then compared using R packages SNPRelate version 1.18.0, gdsfmt version 1.24.1, ape versions 5.4 and 5.8, adegenet versions 2.1.3 and 2.1.10, dendextend versions 1.13.4 and 1.18.0, ade4 version 1.7–18 and RColorBrewer version 1.1-3 [58, 65–72]. As the genetic distance matrices were generated using limited subsets of markers, the resulting dendrograms may not provide representative portrayals of genetic relationships between cultivars. The mantel.rtest() function from the R package ade4 version 1.7–18 and the mantel.test() function from the R package ape version 5.4 were used to compare genetic distance matrices produced by the various subsets of SNPs with the full set of 4,318 high-quality SNPs to determine correlations and similarity between matrices [65, 72, 73]. Mantel’s tests allow the comparison of symmetrical matrices with identical dimensions to determine whether the matrices are similar or correlated. Mantel’s test for similarity between two matrices generates a z value and a p value. Mantel’s z-statistic (z) is the sum of the products of corresponding elements of the compared matrices. A z-statistic is generated for the data and for a permutated version of the data. These statistics are compared to generate a p value, which indicates statistical significance. The threshold for significance set for the Mantel’s tests for similarity was p < 0.01. Mantel’s test for correlation generates r, which is the correlation coefficient, and a p value. The correlation coefficient represents the strength of the relationship between the genetic distance matrices, and the p value, its statistical significance. The threshold for significance of the correlation tests was p < 0.01. Only correlation results are presented in this study. Further analysis of rhAmp assay results was conducted using Excel. 95% binomial confidence intervals were generated using Agresti-Coull method of the binom.confint() function from the R package binom version 1.1–1.1 [74]. Heatmaps presented in Fig. 4a-b were generated using the heatmap() function in R version 4.3.0 [59]. Significance and variance testing were conducted using the functions t.test() and var.test() in R version 4.3.0 [59].

Sequencing results, coverage and read counts

Sequenced cultivars were categorized into three types (1– THC-dominant, 2– balanced THC: CBD, 3– CBD-dominant) using THC and CBD contents provided by OrganiGram. THC and CBD contents of the cultivars could not be published due to their proprietary nature. When data on THC and CBD content was not provided by the LP, cultivar THC and CBD profile data was obtained from Leafly [50]. In total, 901 Gb of raw sequencing data was generated, which represents 5.08 billion 2 × 150 bp paired-end sequences. Pre- and post-processing read counts, coverage, alignment rate and variant counts are presented in Table 2. An average of 1.64 × 10⁸ (SD = 1.91 × 10⁷) reads were produced per cultivar. After filtering using FASTP, with filtering parameters for read quality (average quality score filter = 20), base correction based on overlap, and automatic adapter trimming, 70.3% of reads were retained across all cultivars, corresponding to an average of 1.15 × 10⁸ (SD = 1.70 × 10⁷) reads per cultivar [51]. ‘CBD Shark’ (CBSh1_Seq), ‘Acadia’ (Acad1_Seq), and ‘Armageddon’ (Arma1_Seq) had the highest coverage, while ‘Nukem’ (Nuke1_Seq) and ‘BC GodBud’ (BCGB1_Seq) had the lowest (35.7x, 34.8x, 33.6x, 22.4x, and 20.8x, respectively). Average coverage pre-QC was 28.2x (SD = 3.3), dropping to 19.8x (SD = 2.92) after QC. The variance between pre- and post-QC read counts and between pre- and post-coverage was not significantly different, though mean values were significantly lower (p < 0.01). Filtering and alignment were reproduced using the ‘Pink Pepper’, ‘CB2’, ‘Purple Kush’, and ‘White Widow B’ reference genomes, which are Type 3, Type 2, Type 1, and Type 1 cultivars, respectively [39, 75–77]. Supplementary Table 2 presents this additional coverage and alignment data.

SNP detection and distribution across chromosomes

SNPs detected per cultivar are reported in Table 2. 15,509,562 SNPs were called in the pre-processed merged VCF file. Cultivars belonging to Type I (n = 25) had a mean variant count of 3,856,816 (SD = 3.24 × 10⁵), Type II had a mean variant count of 3,566,463 (SD = 1.32 × 10⁵), and Type III had a mean variant count of 3,267,781 (SD = 1.38 × 10⁶). No significant difference was detected between the means of any Type or between any Type and the full sample. After several rounds of stringent filtering based on mean depth (33x), linkage disequilibrium (0.5), MAF scores (excl. < 5%), heterozygosity (excl. > 60%), missingness (excl. > 20%), and base quality (bq = 30), 4,318 SNPs were retained.

The distribution of these SNPs is presented in Fig. 1. There was a significant correlation between chromosome length and SNP density (R² = 0.737, f < 0.01), but no significant correlation between the chromosome length and retained SNP count. Beyond the higher density of variants in certain chromosomes, some chromosomes host a disproportionate number of SNPs with high MAF values. Chromosomes 2, 5, 6, 8, and X account for roughly half the bp of the genome (excluding contigs), and 52% of the 4,318 high-quality SNPs, but 17 of the 20 SNPs targeted by the genotyping assays are from these chromosomes. The average number of SNPs per Mbp was not significantly different between the chromosome 2, 5, 6, 8, and X group and the group consisting of the remaining chromosomes.

In silico marker selection

After identifying high-quality SNPs with high MAFs, a variety of subsets were tested iteratively to find the minimal marker set. The selection of top markers by MAF was sufficient to differentiate most cultivars. Certain cultivars, notably ‘Nukem’ and ‘Armageddon’, required a unique SNP to be differentiated due to their apparent genetic similarity. After accounting for the SNPs required to differentiate highly similar cultivars, a total of 6 SNPs were found to be sufficient to differentiate all 32 cultivars in silico. A further 14 SNPs were retained for in vitro validation to provide redundancy. 75% of SNP targets (15 of 20) were in intergenic regions, which are non-coding areas between genes. The five remaining SNP targets were located in coding regions, with two of the five being in C. sativa-specific genes. Table 3 summarizes location and mutation types of the 20 SNP targets.

Genetic distance matrices produced with SNP subsets were compared using Mantel’s t-test to test similarity and correlation with either the full set of 4,318 SNPs or the validation subset of 20 SNPs. Figure 2 presents a hierarchical clustering dendrogram using all retained SNPs, with clustering and Type information.

Table 4 reports the results of Mantel’s correlation test between genetic distance matrices created using progressively fewer SNPs and the genetic distance matrix generated using the full set of SNPs. Table 5 reports results for the same test while using the 20 SNPs retained for the genotyping assay as the reference matrix.

There was a weak (r = 0.135) but non-significant (p = 0.041) correlation between the genetic distance matrices of the full set of 4,318 SNPs and the validation subset of 20 SNPs. The genetic distance matrix generated by the 20-SNP validation subset was significantly correlated to all the genetic distance matrices produced by the further reduced SNP subsets (p < 0.01). The correlation between the 20-SNP validation subset and the reduced subsets decreased steadily from the 10-SNP subset (r = 0.738) to the 3-SNP subset (r = 0.433). Figure 3a-i illustrate the differentiation of the cultivars using progressively fewer SNP markers. Figure 3h illustrates the first failure in the complete discrimination of cultivars.

Figure 6 presents a tanglegram comparing genetic relationships generated using in silico and in vitro data. Cultivars ‘TB1004’ (TB141_Seq, TB141_FG) and ‘Ultrasour’ (UlSo1_Seq, UlSo1_FG), ‘Nukem’ (Nuke1_Seq, Nuke1_FG) and ‘Armageddon’ (Arma1_Seq, Arma1_FG), and ‘Critical Kali Mist’ (CrKM1_Seq, CrKM1_FG) and ‘BCGB’ (BCGB1_Seq, BCGB1_FG) remained proximal across both datasets.

In vitro assay validation

In vitro assay validation was conducted using rhAmp assays. Each SNP marker was tested in triplicate on DNA from each of the 32 cultivars. 4.6% of assays failed to amplify a target (no call rate). We compared results for all replicates (32 cultivars, 20 SNP assays, three replicates per assay– 1920 total calls) and for consensus calls (32 cultivars, 20 SNP assays, one consensus call per assay triplicate– 640 calls). Call consensus rate was defined as the proportion of cultivar x SNP observations where at least two of three triplicates reported the same genotype. Call accuracy was defined as the proportion of assays triplicates where at least two of three replicates called the same genotype as was predicted by the sequencing data. Table 6 reports consensus rates and accuracy across all assays during the fresh leaf validation testing and further testing on DNA extracted from dried flower samples of commercial cultivars. Across all cultivar x SNP observations, 81.1% (0.95CI– [0.783, 0.836]) of called genotypes were supported by three identical calls, and 98.9% (0.95CI– [0.979, 0.995]) were supported by at least two identical calls.

No significant difference was found between the consensus rates of the validation (FG) and commercial data sets (DG), or with either of the prior subgroups and the combined group. Figure 4 present call accuracy and consensus using heatmaps. The consensus rate is consistently very high, whereas the accurate call rate varies between assays, and is overall middling.

Call accuracy was significantly higher in the validation set than in the commercial sample set (45.2%, 31.7% respectively, p < 0.05). However, when commercial sample genotypes were compared to the validation genotypes, call accuracy was not significantly different (47.4% and 45.2% respectively, p > 0.05). For the validation data set, call accuracy rate was dispersed across assays (mean = 45.2%, SD = 16.8%). For the commercial data set, the dispersion was greater (mean = 47.4%, SD = 30.1%). No significant difference was found when comparing the overall accuracy of calls on a per replicate versus a per triplicate basis. No significant difference was found between assays targeting intergenic vs. non-intergenic SNPs or between transition- and transversion-type SNPs.

To further assess the utility of the SNP assays, three sequenced cultivars (‘Acadia’ a.k.a. ‘Blue Dream’, ‘Ultrasour’, and ‘Wabanaki’ a.k.a. ‘Jack Herer’) were chosen to be compared to putatively identical (i.e., same cultivar name) cultivars available from the Canadian recreational market. Three products per cultivar were tested, one of which being a product from the LP which had provided the initial fresh leaf samples. Yield and quality of DNA extracted from the dried flower samples were generally poorer than those from fresh leaf extractions (see Supplementary Table 1). Each cultivar’s sequencing and genotyping data was compared to its dried flower equivalent from the LP which provided the sequenced sample, and two other dried flower samples from other LPs. We re-evaluated the accurate call rate from sequenced data to validation calls using this reduced data set which only included the three cultivars of interest and found the accurate call rate to be significantly higher than in the full data set (reduced set − 58.3% (95% CI [0.510, 0.653], full set– 45.2% (95% CI [0.429, 0.474])). 47.4% (95% CI [0.432, 0.516]) of genotypes called by assays deployed on the dried flower samples matched the genotype calls expected based on the in vitro validation data. When using the sequencing data as a reference point, the accurate call rate was significantly lower at 31.7% (95% CI [0.279, 0.357]).

Using the validation data as the reference calls (i.e., expected calls), accuracy for SNP assays tested on dried flower samples varied across SNPs. The lowest accurate call rate was 0% (SNP15 (95% CI [-0.023, 0.148]), and the highest was 100% (SNP12 (95% CI [-0.852, 1.023]). Certain assays consistently failed to amplify in tests on dried flower samples despite being successful in the previous in vitro validation. As the threshold for filtering by missingness was not set to 0 (i.e., some retained SNPs were detected in most, but not all cultivars), the number of No Calls predicted based on post-processing sequencing data was not 0. Based on the sequencing data, 3.9% of assays were not expected to detect any targets due to the absence of a complementary target (i.e., neither targeted allele was present). In the validation tests, 4.6% (95% CI [0.037, 0.056]) of assays failed to amplify any target, but only 7.4% of those No Calls matched the predicted No Calls. Using the reduced validation dataset, 1.7% of calls were predicted to be No Calls based on both sequencing data and validation data calls. In the tests on dried flower samples, 30.1% (95% CI [0.272, 0.350]) of assays failed to amplify a target. To identify factors that may influence assay success, we analyzed a range of variables to suss out any correlation with call data and assay success. As all assays were designed following a close study of assay chemistry and specificity, and SNPs were chosen using existing criteria from a highly filtered pool, no meaningful differences were found.

In vitro differentiation

Table 7 reports the number of cultivars differentiated with each subset of SNPs. In silico, all 32 cultivars could be differentiated with as few as 6 SNP markers. 30 of 32 cultivars could be differentiated in vitro with as few as 6 SNPs. As with the in silico differentiation, ‘Nukem’ (Nuke1_Seq, Nuke1_FG) and ‘Armageddon’ (Arma1_Seq, Arma1_FG) were particularly challenging to tease apart. In vitro differentiation of commercial cultivars was possible, with no samples sharing identical genotyping results.

Further analysis was conducted to determine whether in vitro genotyping data generated similar genetic relationships between cultivars as those generated using sequencing data. Mantel’s t-test was used to compare the genetic distance matrices produced by in silico genotyping data and in vitro genotyping data, also summarized in Table 7. There was a non-significant negative correlation (r = -0.058, p = 0.823) between these matrices, and a significant correlation between the reference matrix generated using genotyping data and all other matrices generated using genotyping data (p < 0.01). Figure 5a-f presents dendrograms illustrating the differentiation of the cultivars using progressively fewer SNP markers. Figure 5b and d illustrate sequential failures in the complete discrimination of cultivar.

Figure 7 presents a hierarchical clustering dendrogram using in vitro genotyping data from fresh leaf samples and dried flower samples from putatively identical cultivars produced by the same LP in addition to genotyping data from the dried flower samples from putatively identical cultivars produced by different LPs. Genotyping results indicate that the putatively identical cultivars were not genetically similar, nor did they cluster together. Genotyping results for fresh leaf samples (Acad1_FG, UlSo1_FG, and Waba1_FG) were not identical to those of their dried flower sample counterparts (Acad1_DG, UlSo1_DG, and Waba1_DG), despite being grown by the same LP.

Sampling strategy

Sampling strategies in genomic research vary broadly depending on experimental objectives. Research on C. sativa genomics has used a range of sampling densities, ranging from single individuals to ten or more per cultivar [25, 27, 78, 79]. Varying levels of heterozygosity within seed lines, an indicator of genetic instability in seed-propagated varieties, have been found in hemp and cannabis cultivars [80–82]. As such, significantly higher genetic diversity is expected between seed-grown plants of a given cultivar as compared to clonally propagated plants of a given cultivar, requiring a commensurately greater number of individuals to be sampled in order to adequately characterize this increased genetic diversity [82, 83]. In the context of our study, 31 of 32 samples were composed of clonally propagated cannabis cultivars. Two individuals per cultivar were sampled on the presumption that there would be minimal intra-cultivar variation between these clonally propagated plants, which were generated from the same mother plant [81]. The sole hemp variety that was sequenced, ‘Anka’, was grown from seed. As such, the sampling of only two individuals may have captured an insufficiently broad range of genetic diversity to fully account for intra-varietal variability in ‘Anka’. Users of our publicly available ‘Anka’ sequencing data (NCBI Sequence Read Archive SRR14857087) should therefore be mindful of the limitations of this data in terms representativeness of the variety’s genetic diversity. However, filters applied to our variants data filtered out SNPs that were not well-represented across the sample population, thus severely curtailing the data on unique elements of ‘Anka’s genetic composition (discussed further in Variant calling and SNP selection) and limiting the utility of a broader sample of ‘Anka’ in achieving our experimental objectives.

With regards to the sequencing data for the other 31 cultivars, the genetic variability between clonally propagated plants is expected to be limited, thus sampling from a limited number of plants per cultivar is appropriate for the experimental goals of this study, though secondary users of these datasets should be mindful that representativity of the genetic composition of cultivars within one LP does not equate to representativity of putatively identical cultivars across LPs.

This latter limitation is expected to contribute to the performance of assays when tested on samples from putatively identical cultivars (i.e., cultivars with the same name but grown by different LPs). In short, significant research on both hemp and cannabis cultivars has demonstrated intra-cultivar variability between individuals, particularly when grown by different cultivators [27, 80, 82]. Further, it was not possible for us to determine whether flowers in the dried flower samples were harvested from one or many individuals. Assuming other LPs also use clonal propagation, the inclusion of flowers from many individuals may result in some intra-cultivar variability being captured by the genotyping assays, though the effect of this factor on assay results is questionable due to the relative stability in genetic composition between clonally propagated individuals [80, 81]. The impact of these sources of variability is discussed further in Assay performance.

Sequencing

As our experimental goal was not genome assembly, many of the issues inherent to short-read sequencing, particularly in C. sativa, pose lesser risks to our experimental design [37, 38, 40]. As such, the unparalleled sequencing depth offered by short-read sequencing and scalability at significantly lower costs than the technically superior hybrid-sequencing approaches (i.e., using short- and long-read technologies) makes short-read sequencing considerably more appealing. As financial and computational constraints are relevant factors for most labs weighing sequencing approaches, prioritizing sequencing depth and sampling density (i.e., individuals/cultivars sequenced) is recommended in projects depending on SNP discovery [48, 84]. Sufficient depth is particularly important when seeking to discover SNPs, as it minimizes the risk of both false positives and negatives while increasing the likelihood of identifying rare SNPs [85–87]. In our samples, coverage of 25x was reached or surpassed for all but three cultivars. Further filters were applied once the SNPs were called, setting a 33x threshold for coverage for each retained variant. While there is a case to be made for lower coverage thresholds, 20-30x coverage is generally sufficient for SNP discovery [86, 88]. As sequence repeats (e.g., Short-Tandem Repeats) and Copy Number Variation (CNV) are common in C. sativa, their impact on mapping quality and false-positive heterozygous calls are worth mitigating [89]. Longer sequencing reads provide more base pair comparisons per read to grant greater overlap, and paired-end reads can contribute to reducing the risk of confusing very similar regions of the genome [89, 90]. While our sequencing approach did use paired-end reads, using longer short read technology, such as 300 bp paired-end sequencing, rather than 150 bp paired-end sequencing, might have enhanced mapping results without resorting to the gold standard of hybrid sequencing. Regarding issues stemming from CNVs, approaches mining SNPs for information and applying new bioinformatic techniques may help to rule in/out CNVs [91].

Mapping and choice of reference genome

Filtered reads were aligned to the ‘CBDRx’/cs10 reference genome. Comparisons of key quality metrics for existing assemblies continue to support and generally favour the use of the cs10 reference genome and at the time of the development of the panel, the ‘CBDRx’ assembly was the reference genome of choice [39, 40, 92]. Since the development of our SNP panel, the ‘Pink Pepper’ assembly has supplanted the ‘CBDRx’ assembly as NCBI’s reference genome of choice for C. sativa [39]. Notably, both ‘Pink Pepper’ and ‘CBDRx’ are Type III cultivars. As the majority (25 of 32) of cultivars in our sample were of Type I, it is possible that the difference in genetic composition underlying the differences in cannabinoid production may have generated suboptimal results when aligning Type I (or Type II) cultivars to a Type III reference genome [93]. While a robust comparison of reference genomes falls outside the scope of this work, we repeated the filtering and mapping process of the pipeline using a variety of other reference genomes (see Supplementary Table 2). Compared to the alignment rates of all sequenced cultivars to the ‘CBDRx’ assembly, we found that alignments rates were significantly higher when aligned to the ‘CB2’ (Type II) and ‘White Widow B’ (Type I), assemblies (p < 0.01) [75, 77]. No significant differences were found when compared to alignment rates to the ‘Pink Pepper’ (Type III) and ‘Purple Kush’ (Type I) assemblies (p > 0.01) [39, 76]. While alignment rates were sometimes similar between reference genomes of different Types, variability between reference genomes has been discussed elsewhere [76, 77, 94]. This variability will influence the variants detected [94].

When seeking to differentiate between cultivars as we have, the ubiquity of SNPs and the ability to filter detected variants to ensure their suitability for the experimental outcome limit the impact of considerations beyond assembly quality when selecting a reference genome. As detected variants were first filtered for quality, then were subjected to filters that targeted prevalence of the SNPs and their constituent alleles in the sample population (i.e., the 32 sequenced cultivars), the final set of SNPs is culled from variants with high mapping scores that are well represented across the sample population. Considering the available controls and the abundance of SNPs, it is not clear that one reference genome should be selected over another of comparable quality if seeking to differentiate between cultivars without consideration of trait-specific loci.

However, it is conceivable that differing experimental goals may be better served by choosing an alternate assembly. As noted, if one’s experimental objectives include trait or locus-specific targets, the selection of a reference genome may require that some additional consideration be given to the properties of the cultivar used to generate the assembly. Currently, there is a wide variety of chromosome-level assemblies, allowing for careful selection of the assembly most similar (or dissimilar) to a given sample population [37–39, 75–77, 95]. In our case, the use of an assembly derived from a Type I cultivar may have allowed for screening of genes coding for THCA synthase, which is not fully present in the ‘CBDRx’ assembly [38]. The selection of an alternate reference genome (e.g., the ‘White Widow’ assembly instead of the ‘CBDRx’ assembly) would likely result in a completely different set of SNPs being detected due to the genetic differences between the two cultivars, and it is expected that an entirely different set of filtered SNPs would have been considered for the 20-SNP genotyping panel [94]. The increasing number and diversity of C. sativa assemblies offer a great deal of choice, thus careful consideration can and should be taken in selecting among the available reference genomes if experimental outcomes are expected to be significantly affected by the properties of the cultivar used to build the assembly.

Variant calling and SNP selection

No significant differences were found in the number of variants called between Types. A lack of chemotypic diversity in our sample should temper any broad inferences based on intergroup comparisons; however, it is not shocking that cultivars with the most and the least variants were assigned to Type III. ‘ACDC’, a Type III cultivar known for its high CBD and low THC content, was bred for similar purposes as ‘CBDRx’ [38, 50]. ‘Anka’, also a Type III cultivar, is prized for its grain and fibre yield, and was bred for distinctly different purposes than the rest of the cultivars included in our sample [96]. This vast intra-Type difference exemplifies the vast genetic diversity that is ignored when the focus is uniquely on THC and CBD contents and supports the practical categorization of cultivars by their intended use (i.e., hemp and cannabis, in keeping with Canadian regulatory definitions).

SNPs gain several advantages as genetic markers due in part to their vast numbers across genomes. Despite potential genetic bottlenecks in the recreational market, the ubiquity of SNPs found in C. sativa still allow for the development of powerful genotyping arrays [81, 97]. In our study, the abundance of variants allowed for the application of more stringent filters across multiple stages of analysis. In short, filters excluded variants that exhibited high heterozygosity (> 60%), low MAF (< 0.05%), and poor-quality bases (depth < 33x, BQ < 30). Variants that were absent in over 20% of samples were also excluded. Excluding low-frequency variants (i.e., due to missingness and MAF filters) is standard practice, though the loss of information can cause inaccurate reports of similarity, as was the case in the aggregation of ‘Anka’ with other cultivars in dendrograms produced using our in silico data, as it ought to have been a clear outlier (see Fig. 2). Based on Fig. 2, which was generated using filtered SNPs, there would seem to be a similar level of genetic diversity across the primarily Type I clusters (clusters 1, 2, and 3, represented by red, blue and green, respectively) as there is across the Type II and III clusters (clusters 4, 5 represented by purple and orange, respectively, on Fig. 2). We expect that this is a consequence of the filtering parameters which favour more stable loci (i.e., low heterozygosity) and which eliminated variants and alleles that are uncommon in our primarily Type I sample population due to the missingness and MAF filters. This tailoring of filtered SNPs to Type I cultivars may be a benefit when seeking to differentiate cultivars with relatively low inter-cultivar genetic diversity but may also reduce the applicability of the retained SNPs to a more genetically diverse sample population.

MAF was chosen as a key selection criterion following the example of in silico research on hops (Humulus lupulus), a close relative of C. sativa [49]. The focus on MAF arose from early research into the efficient selection of tagging SNPs [98]. As our goal is to develop and test a genotyping assay capable of distinguishing between C. sativa cultivars, we are, in effect, searching for a minimal set of tagging SNPs. The development of minimal marker sets and tagging SNP algorithms has been the focus of a great deal of research, and multiple metrics were developed and contrasted over time to aid in the selection of informative tagging SNPs [99–104]. There are multiple technical considerations when deciding which SNPs to retain for in vitro validation, such as targeting assay chemistry, target specificity, and SNP calling quality score. Further, experimental objectives may require the consideration of additional metrics, such as allele prevalence and whether the SNP is in a coding or non-coding region. We retained 20 SNPs as targets based on the aforementioned technical considerations, but did not favour intragenic or intergenic targets. 15 of the 20 SNPs targeted by our assays were in intergenic regions, and we found no significant difference in the performance of assays targeting intragenic and intergenic SNPs.

From a technical perspective, we used a standard suite of tools (Bowtie2, SAMTOOLS, BCFTOOLS, and VCFTOOLS) for the alignment, variant calling, filtering, and analysis of our data. While our approach was in line with standards established in contemporaneous work, inventories of available pipelines and comparisons of their performance highlight potential improvements [89, 90, 105]. Notably, the use of hybrid sequencing approaches, which integrate long reads and short reads, can mitigate many mapping issues, and significant research has been conducted to develop better mapping and haplotype resolution algorithms [106]. Further, post-processing using tools such as GATK VariantRecalibrator, quality benchmarking, fragment length thresholds to prevent false-positive variant detection caused by paired-end overlap, and the use of novel variant calling pipelines (e.g., Dynamic Read Analysis for Genomics (DRAGEN)) could all contribute to improving the stock from which we selected our SNP markers for in vitro validation [89, 105, 107, 108]. Inventories and subsequent evaluation of variant calling pipelines have noted that approaches should be tailored to experimental goals, but a 2020 study found that the use of Novoalign combined with the Genome Analysis ToolKit (GATK) for alignment and variant calling provided generally superior results [105]. Finally, new medium-density genotyping platforms such as HASCH and CoreSNP may provide streamlined and sophisticated processes to facilitate initial reductions in SNP targets [109, 110].

Assay performance

We tested our 20-SNP genotyping panel’s ability to differentiate between the 32 sequenced cultivars in vitro using rhAmp assays. Differentiation of all 32 cultivars was possible using as few as seven SNPs in vitro (SNPs 1 through 6 and SNP10), whereas differentiation of all 32 cultivars required only six SNPs in silico (SNP1 through 6). While other studies have used SNPs to group C. sativa cultivars by chemotype, by crop type and geographical origin, and a recent study using SSRs was able to cluster five cultivars in highly cultivar-specific groups, our findings demonstrating the feasibility of leveraging minimal marker approaches to differentiate between C. sativa cultivars using SNPs are novel [23, 34, 111].

The proportion of our assays which failed to amplify a target (no call rate) was 4.6%, resulting in a conversion rate of 95.4%. While the no call rate of 4.6% was higher than expected based on research evaluating rhAmp assay performance using validated SNPs, research which validated medium and low-density SNP assays in other crops found conversion rates ranging from 70 to 100% across assays [44, 112, 113]. Alleles detected by the SNP assays were highly consistent, with 81% of genotype calls being supported by all three replicates, whereas 98.9% of genotype calls were supported by two or more replicates. However, while differentiation with seven SNPs demonstrated the feasibility of the approach, seven is not the smallest number of SNPs technically required to encode 32 unique states. As biallelic SNPs can be detected in one of three states (homozygous for the reference or the alternate allele and heterozygous), the maximum number of cultivars that could be differentiated by n SNPs is 3ⁿ. Tools such as Minimal Marker may facilitate the selection of parsimonious SNP subsets from both in silico data and in refining subsets after in vitro testing [99].

While cultivars were differentiated, the genotypes detected by the SNP assays did not consistently concord with the genotypes detected by sequencing. While some research compares the concordance of genotypes detected in silico using a variety of bioinformatic approaches, concordance between various in vitro methods of genotyping for cultivar differentiation in C. sativa remains poorly studied [109]. We found that 45% of genotypes detected via SNP assays matched the genotypes detected via sequencing. Figure 6 provides a clear representation of the impact of these discrepancies on the genetic relationship established using both sets of genotypes. As both sets of genotypes were captured from the same sample population, this level of discrepancy was not expected. While refinements to the sequencing, mapping and variant calling approach have already been discussed, an additional source of uncertainty, paralogous sequence variants (PSV), is worthy of consideration. PSVs are caused by gene duplication, a common event in plants that has also been studied in C. sativa, and can result in the erroneous alignment of similar sequences to the same position [93, 114–116]. Any variation to the PSV having occurred after the duplication event may then result in the erroneous detection of a SNP, which generally results in erroneous heterozygous calls [117, 118]. The use of novel bioinformatic tools, such as the Putative Paralogs Detection pipeline, may help reduce the risk of erroneous SNP detection caused by the presence of PSVs [115]. In addition to these refinements, higher density sequencing (i.e., sampling more individuals per cultivar) may improve resolution of genotypes in C. sativa [34].

Even when provided with good targets, genotyping assays may amplify non-target loci. In short, non-target amplification, or non-specific amplification, is the result of an assay binding to an unexpected target, which may generate unexpected results. Genotyping technology (i.e., rhAmp assays) such as that used in our work is designed to reduce non-target amplification, and specificity was tested using BLAST [62, 112]. DNA concentration was in accordance with the assay manufacturer’s specifications (5 ng/µL), thus limiting the risk of allelic drop out due insufficient template DNA. However, allelic dropout can also be the result of other causes, such as the presence of variants in or proximal to binding sites, leading to the amplification of only one allele, which then results in the detection of unexpected homozygosity [119, 120]. Despite these issues, improvements to assay performance and further reduction of the number of SNPs used to discriminate between cultivars could be achieved through the integration of assay performance data as a SNP retention criteria. An iterative approach has previously been applied to successfully differentiate between 94 varieties of cabbage (Brassica oleracea var. capitata) using 10 SNPs [44]. In this study, researchers selected 348 SNPs with polymorphic informative values of 0.11 to 0.37, MAFs ranging from 1 to 30%, and call rates ranging from 72 to 100%. These 384 markers successfully discriminated between 94 varieties. Further successful discrimination was achieved through subsequent genotyping with a subset of 87 markers, then reduced sets of 24 and ten markers were proposed [44].

Assay performance was poorer when tested using DNA extracted from the dried flower samples of commercial cultivars. Notably, the no call rate was significantly higher than in the validation set (30.1% and 4.6%, respectively, p < 0.01). While DNA extracted from dried flowers was generally less abundant and of lower quality (see Supplementary Table 1), the concentration met the assay manufacturer’s specifications for all samples. As such, it is unlikely in this case as well that alleles were not detected due to insufficient template DNA. A more likely explanation is that the targeted alleles did not exist in the putatively identical cultivars, resulting in failed amplification. Considering the findings of studies comparing genetic differences between cultivars with the same label or name, this divergence between the commercial samples is not wholly unexpected [27, 80–82]. In Fig. 7, we note that there are differences at a minimum of five loci between the most similar dried flower samples (Waba2_DG and Waba3_DG, UlSo2_DG and UlSo3_DG). Further, there was limited similarity between genotyping results for dried flower and fresh leaf samples from the same cultivar grown by the same LP, with only UlSo1_FG and UlSo1_DG having fewer differences than the least dissimilar dried flower samples. As there were only 20 loci considered, this is a large divergence between cultivars sharing the same name. Another potential source of genetic variability is the genetic composition of the dried flower samples. It was not possible to determine whether each dried flower sample was composed of genetic material from a single individual. As such, there may be some level of intra-cultivar variability between the individuals composing each sample, in addition to the intra-cultivar variability that is expected between samples from the same cultivar grown by different LPs [27, 80]. Higher density sampling, including sampling of individuals of the same cultivar grown by multiple LPs, has been used by others to improve resolution of cultivars with distinct genetic identities but that share the same name [80, 81]. Increased sampling density may have resulted in the detection and selection of SNPs that accounted for the diversity within individuals of putatively identical cultivars, thus resulting in improved aggregation of genetically similar individuals. Due to our limited sampling, uncertainty regarding genetic composition of dried flower samples, and the limited number of loci successfully assayed in our genotyping panel, we cannot confidently assert that our genotyping results demonstrate significant inconsistencies in genetic identity between putatively identical cultivars. Considering the challenges in validating the genetic composition in dried flower samples, modification of the bio-informatic pipeline to incorporate tri-allelic SNPs might improve discriminatory power and mitigate some of the uncertainty by identifying mixed samples in diploid species, though this latter benefit has only been well-studied in humans [121, 122]. As C. sativa is diploid, the presence of a third allele in a given sample might indicate the presence of more than one individual, or more than one cultivar, in the sample, though recent research studying naturally occurring polyploidy in C. sativa has introduced an interesting caveat to this generalization [123].

Finally, DNA of dried flower samples may have been exposed to radiation as part of post-harvest sterilization by irradiation procedures, potentially resulting in DNA damage. There is significant research on the effects of sterilization by irradiation on a variety of traits of agronomic interest in plants, including the production of C. sativa metabolites, a good deal of research on the effects of radiation on plant DNA, but research on the effects of radiation on post-mortem and ex vivo DNA is primarily focused on human forensics [124–134]. Due to the paucity of research at the intersection of these fields and because our methodology did not incorporate appropriate controls to adequately characterize the potential effects of irradiation on C. sativa DNA, the influence of these sterilization procedures cannot be discussed without resorting to speculation. Should further research elucidate the potential impact of post-harvest sterilization procedures on C. sativa DNA, research on tri-allelic SNPs in forensic applications has demonstrated their potential in accurately genotyping degraded DNA samples [121].

Limitations of minimal marker approaches in determination of genetic relationships

While cultivars were successfully differentiated from one another, genetic relationships between cultivars were not maintained between sequencing and SNP assay results (see Fig. 6). Additionally, the reduced number of SNPs limited our ability to establish genetic identity thresholds, as have been proposed in other eukaryotes [135–137]. Further, genetic relationships between cultivars were generally not maintained as the number of SNPs used to differentiated between cultivars were reduced, as demonstrated by the changes to cluster assignments and neighboring cultivars for any given cultivar in Figs. 3a-i and 5a-f. Despite these differences, there remained a significant, though progressively weaker, correlation between the genetic distance matrix generated using 4,318 SNPs and those generated using 1,000 (r = 0.507), 100 (r = 0.424), 50 (r = 0.248), and 25 SNPs (r = 0.151) (p < 0.01) (see Table 4). This non-linear decrease in correlation is counterintuitive to the impression gleaned when studying Fig. 3a-i. Interestingly, correlation values were comparable between the 4,318 SNP and 1,000 SNP matrices (r = 0.507) and the 20 SNP and 5 SNP matrices (r = 0.554). In both comparisons, eliminating roughly three-quarters of the SNPs results in a correlation of roughly 0.5. These significant correlations demonstrate that while representations of genetic relationships and clustering results are altered, the reduction of SNPs does not result in random genetic relationships between cultivars. However, the fact that genetic relationships established by smaller numbers of SNPs are correlated to the genetic relationships established by larger numbers of SNPs does not mean the relationships established by the former group are meaningfully representative of the true genetic relationships between cultivars. Furthermore, the assumption that the genetic relationships between cultivars established by the 4,318 SNPs (Fig. 2) are representative of the true genetic relationships between cultivars is itself worth questioning. Notably, cultivars belonging to Type III did not cluster together, and a cultivar from Type II clustered with Type I cultivars. As these unexpected cluster assignments may be a result of the overrepresentation of Type I cultivars in the sample population and filtering parameters, the development of representative genetic relationships between cultivars of different Types may require increased chemotypic diversity in the sample population. While accurate representation of genetic relatedness is not strictly necessary for effective cultivar differentiation, our findings highlight the limitations of applying low-density genotyping to the study of genetic relatedness. However, the inclusion of genetic relatedness analysis (e.g., STRUCTURE) to identify closely related cultivars might be a worthwhile consideration when working with large sample populations to ensure that SNPs differentiating these cultivars are not culled early in the reduction process [138].

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.