Cannabis labelling is associated with genetic variation in terpene synthase genes
grid.55602.340000 0004 1936 8200Department of Plant, Food and Environmental Sciences, Dalhousie University, Truro, Nova Scotia Canada
Bedrocan International, Veendam, the Netherlands
grid.4818.50000 0001 0791 5666Biosystematics Group, Wageningen University, Wageningen, the Netherlands
Abstract
Analysis of over 100 Cannabis samples quantified for terpene and cannabinoid content and genotyped for over 100,000 single nucleotide polymorphisms indicated that Sativa- and Indica-labelled samples were genetically indistinct on a genome-wide scale. Instead, we found that Cannabis labelling was associated with variation in a small number of terpenes whose concentrations are controlled by genetic variation at tandem arrays of terpene synthase genes.
Web Summary
By quantifying over 100 Cannabis samples for terpene and cannabinoid content and genotyping them for over 100,000 single nucleotide polymorphisms, this study finds that Cannabis labelling is associated with genetic variants in terpene synthase genes.
Main
Cannabis has been consumed for its psychoactive properties for over 2,500 years, and its estimated global market value is US$340 billion1–3. Because it is a widely used drug that is increasingly being legalized for medicinal and recreational use, it is critical that Cannabis’s genetic and chemical variation be accurately quantified and communicated. The vernacular labels Sativa and Indica (not to be confused with the taxonomic names C. sativa sativa L. and C. sativa indica Lam.) are routinely assigned to Cannabis cultivars by breeders, retailers and users to describe a cultivar’s morphology, aromas and/or psychoactive effects4. However, it is unclear whether these labels capture meaningful information about Cannabis genetic and chemical variation.
Cannabis genomics research has thus far largely focused on the characterization of genes underlying the production of the cannabinoids cannabidiol (CBD) and tetrahydrocannabinol (THC)5–8. However, Cannabis produces hundreds of aromatic terpenes that drive consumer preference and are frequently associated with Sativa and Indica labels4,9. In addition, there is evidence to suggest that a cultivar’s terpene profile affects its psychoactive properties10,11. To date, various terpene synthase genes have been identified in Cannabis; however, the genetic control of terpene variation across Cannabis cultivars remains largely unexplored12–15.
Here we re-analysed 297 samples of drug-type Cannabis that were previously quantified for 40 terpenes and cannabinoids using gas chromatography–mass spectrometry (GC–MS)16 (Supplementary Table 1 and Extended Data Fig. 1), and we paired these data with 116,296 newly generated single nucleotide polymorphisms (SNPs) from 137 of these samples from which sufficient high-quality DNA could be extracted. We determined the degree to which the genomic and GC–MS data corresponded to a five-point labelling scale ranging from 1 (100% Sativa) to 5 (100% Indica) as reported by sample sources.
Principal component analysis (PCA) of the genomic data showed no clear clustering according to sample labels (Fig. 1a). Even though PC1 and PC2 were significantly correlated with the Sativa–Indica scale, the variance explained by the primary PCs was low (PC1: R2 = 0.12, P = 2.1 × 10−5; PC2: R2 = 0.12, P = 1.8 × 10−5). Furthermore, the overall genetic structure (captured by including the first ten PCs of the genomic data in a linear model) explained only 37% of the variance in labelling (Fig. 1c). Sativa–Indica labels thus do not accurately reflect genetic relatedness, which is consistent with previous work17,18. In addition, we determined that pairs of samples with identical cultivar names (for example, OG Kush) were often as genetically and chemically distant from each other as pairs of samples with different names (Extended Data Fig. 2). This is consistent with previous studies indicating that cultivar names were not reliable indicators of a sample’s genetic or chemical identity17,19–21.
Similar to the PCA of the genome-wide SNP data, the PCA of the terpene and cannabinoid profiles provided poor separation of samples according to their Sativa–Indica labels (Fig. 1b). Nevertheless, we observed significant correlations between the first two PCs and the Sativa–Indica scale (PC1: R2 = 0.049, P = 7.5 × 10−5; PC2: R2 = 0.24, P = 3.7 × 10−19). Including the first ten PCs from the terpene and cannabinoid profiles in a linear model accounted for only 41% of the variance in labelling (Fig. 1c). The pairwise genetic and chemical relatedness matrices were correlated (Mantel r = 0.21, P = 1 × 10−3, Extended Data Fig. 3), and a linear model including the first ten PCs from both the genomic and chemical profiles captured only 41% (Fig. 1c; P = 3.1 × 10−10) of the variance in labelling. Since the overall patterns of genetic and chemical relatedness could not fully account for the labels applied to Cannabis samples, we aimed to determine which individual chemicals were the strongest predictors of Sativa–Indica labelling.
Of the 40 measured terpenes and cannabinoids, 12 (30%) were correlated with the Sativa–Indica scale at P < 0.01 (Fig. 2a and Supplementary Fig. 1). Sativa content was positively correlated with the concentrations of bergamotene (R2 = 0.12, P = 9.26 × 10−8) and farnesene (R2 = 0.11, P = 1.09 × 10−7), which impart tea-like and fruity aromas, respectively22,23. This is consistent with descriptions of Sativa cultivars as having a ‘sweet’ or ‘herbal’ aroma4,9. The strongest correlation was between Indica content and myrcene, whose concentration explained 21.2% of the variation in labelling (P = 2.29 × 10−15; Fig. 2a). The sedative effect and earthy aroma attributed to high myrcene content are often reported by recreational users to be characteristic of Indica cultivars10,24–26. We also observed significant positive correlations between Indica labelling and three sesquiterpenes: guaiol (R2 = 0.18, P = 7.7 × 10−13), γ-eudesmol (R2 = 0.11, P = 3.8 × 10−7) and β-eudesmol (R2 = 0.21, P = 8.2 × 10−15). Hillig27 found that these three sesquiterpenes were associated with plants from Afghanistan, which is considered the region of origin for Indica cultivars.
Previous chemical analyses of Cannabis have suggested that the distinction between Sativa and Indica is best explained by differences in the concentrations of specific monoterpenes and sesquiterpenes19,28–30. In addition, the contrasting aromas that have been associated with Sativa (that is, sweet) and Indica (that is, earthy) were key discriminators in a sensory evaluation of Cannabis cultivars and mediated customers’ perceptions of potency and quality9. As a previous study suggested31, we hypothesize that Cannabis growers and breeders have been assigning labels to cultivars primarily on the basis of aroma profiles and purported effects, rather than genetic ancestry or overall chemical similarity. The primary differences between cultivars labelled as Sativa and Indica may thus be driven by a small set of genomic regions controlling the concentrations of a small number of contrasting aromas. To examine this, we conducted a genome-wide association study (GWAS) of the 40 chemicals examined here (Supplementary Fig. 2 and Supplementary Table 2).
We identified three regions of the Cannabis genome associated with the four terpenes most strongly associated with Sativa–Indica labelling (Fig. 2). The optimal model from the multilocus mixed linear model (MLMM) GWAS for myrcene identified two significantly associated SNPs 1.2 megabases apart that tag independent blocks of linkage disequilibrium (LD) on the proximal end of chromosome 5 (Fig. 2b). The first SNP (chr5:1348048) is located 6.4 kilobases (kb) from a block of terpene synthase genes composed of four copies of TPS30, which is known to encode myrcene synthase12 (Supplementary Table 3). The second SNP (chr5:2576403) is 46.7 kb from another tandem array of terpene synthase genes spanning ~200 kb (Supplementary Table 3). Within this gene cluster are two sequences highly similar to the myrcene synthase gene, TPS3 (refs. 12,13). These observations suggest that myrcene synthesis is mediated by genetic variants at two independent terpene synthase gene clusters on chromosome 5. The other three sesquiterpenes (guaiol, β-eudesmol and γ-eudesmol) strongly associated with Sativa–Indica labelling are correlated with each other (Extended Data Fig. 4) and share a common GWAS hit on chromosome 6: the single SNP identified from the MLMM (chr6:76790611) is 51.9 kb from a gene cluster comprising sesquiterpene synthase genes related to TPS7FN (δ-selinene synthase), TPS8FN (γ-eudesmol/valencene synthase)12 and TPS20CT13 (hedycaryol synthase) (Fig. 2c and Supplementary Table 3).
Our results demonstrate that the Sativa–Indica scale currently used to label Cannabis poorly captures overall genomic and metabolomic variation. Cannabis labelling is instead probably driven primarily by a small number of key terpenes whose concentrations contribute to the characteristic aromas commonly associated with Sativa and Indica and whose variation we genetically mapped to tandem arrays of terpene synthase genes on chromosomes 5 and 6. While the vernacular labels ‘Sativa’ and ‘Indica’ are derived from taxonomic names that were originally used to categorize plants according to ancestry4, these terms have been co-opted by contemporary Cannabis culture and now probably reflect locus-specific genetic variation affecting terpene synthesis. Our results suggest that a practical and reliable classification system for Cannabis that is consistent with contemporary understanding of the terms ‘Sativa’ and ‘Indica’ may be achievable by quantifying a small number of terpenes and/or genotyping genetic markers associated with key Cannabis aromas.
Methods
Samples
The samples come from a previous study of 460 Cannabis chemotypes16. The samples were collected from Bedrocan International BV (n = 37), HempFlax (n = 205) and Dutch ‘coffee shops’ either directly or indirectly through the TRIMBOS Institute (n = 55). Samples labelled as ‘Hemp’ were excluded from the analysis. We retained and analysed 297 samples that were classified along a five-point scale according to ancestries reported by the sources: ‘Sativa’ (100% Sativa), ‘Hybrid-Sativa’ (75% Sativa, 25% Indica), ‘Hybrid’ (50% Sativa, 50% Indica), ‘Hybrid-Indica’ (25% Sativa, 75% Indica) and ‘Indica’ (100% Indica). These five groups were encoded as 1 (100% Sativa) to 5 (100% Indica) for the statistical analyses described below.
Gas chromatography
A total of 297 samples were previously quantified for terpene and cannabinoid content, and we conduct a re-analysis of these data here. The chemical analyses of the samples are described in detail in ref. 16. Briefly, for each sample, 500 mg of ground homogenized dried flower material was mixed with 40 ml of ethanol, agitated for 10 minutes and centrifuged. The supernatant was collected, and the process was repeated twice more on the pellet. An internal standard consisting of 200 μl of 1% solution of 1-octanol was added to the combined supernatant, the volume was adjusted to 100 ml with ethanol and the combined sample was centrifuged again. The combined sample was analysed using an Agilent GC 6890 series (Agilent Technologies) equipped with a 7683 autosampler and a flame ionizing detector. The instrument was equipped with a DB-5 column (length, 30 m; internal diameter, 0.25 mm; film thickness, 0.25 μm; J&W Scientific). Peaks from the sample chromatograms were manually integrated, and the peak area was recorded with correction for the internal standard peak area. Peak identification was conducted by analysing selected samples using GC–MS and then comparing compounds’ mass spectra and retention times with authentic standards and literature reports as described in ref. 16. Compounds without authentic standards are marked with an asterisk in the figures to indicate that they were tentative identifications. Peak areas of monoterpenes, sesquiterpenes and cannabinoids were quantified (in mg per g of plant material) using calibrated standards of β-pinene, α-humulene and CBD, respectively. We re-assessed the compound identifications in Hazekamp et al.16, and in certain cases we renamed compounds on the basis of the inability to distinguish stereoisomers using a DB-5 column. For example, in the case of the compound listed by Hazekamp et al.16 as ‘(−)-linalool’, we renamed this to ‘linalool’. There are also two compounds that could not be reliably identified; they are listed as ‘unidentified compounds’ (Supplementary Table 3). THC, δ-8-THC and CBN were combined into a single value, ‘Total THC’, because δ-8-THC and CBN are degradation products of THC. Peaks of R-limonene and β-phellandrene were indistinguishable and were therefore combined into a single value and reported as ‘limonene’. Thymoquinone, geraniol, thymol and carvacrol were removed because they were not present in any samples, and cineol was removed because it was present in only one sample. Pearson correlations were calculated between each pair of chemicals using the cor.test function in R v.3.5.132. According to previous work33, the samples analysed here were nearly all drug-type Cannabis (that is, type I) (Extended Data Fig. 1), except nine samples with THC > 0.3% and CBD > 0.5% (that is, type II).
Genomic analysis
Whole-genome DNA was extracted using a NucleoSpin 96 Plant II kit (Machery-Nagel) and quantified using the QuantiFluor dsDNA System and the GloMax-Multi + Microplate Multimode Reader with Instinct (Promega). Genotyping-by-sequencing libraries were prepared using the restriction enzyme ApeKI34, and the libraries were sequenced on two lanes of an Illumina Hi-Seq 4000 (Illumina). The DNA sequence data are available as NCBI BioProject PRJNA713792. Calling of SNPs was performed in TASSEL (v.5.0)35 by aligning to the CBDRx reference genome8. SNP calling was performed before the implementation of the new chromosome numbering of the CBDRx genome in April 2020. Chromosomes were recoded for analyses to reflect the new chromosome numbering system. We used VCFtools (v.0.1.15)36 to retain only bi-allelic SNPs and samples with <70% missing data, which resulted in 155 remaining samples and 284,988 SNPs. Genotype imputation was performed using LinkImputeR37 with a minor allele frequency threshold of 0.01, a minimum read depth for masking of 20 and the number of masked genotypes set to 5,000. We chose to impute with a minimum read count of 2 and a maximum missingness threshold of 70%, which resulted in an imputation accuracy of 92.88%. After imputation, 149 samples remained. An additional 12 samples were removed because they had no phenotype data. This resulted in a final set of 137 samples with both genetic and chemical data. The SNP data were filtered using PLINK (v.1.90)38 to exclude SNPs with a minor allele frequency less than 0.05 and SNPs with excess heterozygosity resulting in Hardy–Weinberg P values less than 1 × 10−5. The final SNP dataset used for GWAS consisted of 116,296 SNPs from 137 samples. For PCA, 1,257 unanchored SNPs were removed, and the remaining 115,039 SNPs were LD-pruned using PLINK (command: –indep-pairwise 10 3 0.5), resulting in 80,939 SNPs.
Genetic and chemical analysis
The chemical distance between cultivars was calculated as the Euclidean distance using the ‘dist’ function in R from the matrix of metabolomic data—that is, 40 terpenes and cannabinoids quantified across 297 samples. The genetic similarity between samples was calculated as an inverse identity-by-state matrix generated in PLINK. The correlations between the matrices were computed using a Mantel test in R32 by first reducing the chemical matrix to the 137 samples with both chemical and genetic datasets. PCA was performed on the scaled genetic and chemical data using the prcomp function in R. To calculate the variance in labelling explained by the chemical and genetic data, linear models including the top ten PCs from the genetic data, the chemical data and both the chemical and genetic datasets together were performed. Pearson correlations between chemical concentration and the 1-to-5 Sativa–Indica scale were performed with the cor.test function in R. A Bonferroni correction was applied to the P values from the correlation test between chemical concentration and the Sativa–Indica scale.
Genome-wide association
We performed GWAS for 40 terpene and cannabinoid phenotypes, using both normalized and non-normalized data. Normalizing was conducted to generate values for a chemical concentration in a sample relative to the total abundance of its chemical class (that is, monoterpene, sesquiterpene or cannabinoid) in that sample. Thus, a sample’s myrcene content was divided by the total concentration of all monoterpenes in that sample to generate a normalized value for myrcene. GWAS was performed using an MMLM39 accounting for relatedness using a kinship matrix created in TASSEL (v.5.0)35. The MLMM incorporates significant SNPs as cofactors using stepwise regression (maxsteps = 10), and the optimal model was chosen on the basis of the extended Bayesian information criterion. We also present the first step of the MLMM, which is equivalent to an MLM where relatedness is accounted for but no SNPs are included as cofactors. Using the simpleM40 package in R, the effective number of independent tests (Meff) was generated, and the threshold for significance was then calculated using −log10(α/Meff), where α = 0.05. Quantile–quantile and Manhattan plots were created using the qq function in R. Genomic regions with significant GWAS hits were explored, and the physical locations of genes within these regions were retrieved using annotations from the CBDRx reference genome8 in Geneious Prime (v.2020.1.2). The GWAS results and LD regions of interest were visualized using code adapted from ref. 41.
Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Supplementary information
Appendix Group
Appendix
Extended data
Untitled section
Extended data
is available for this paper at 10.1038/s41477-021-01003-y.
Supplementary information
The online version contains supplementary material available at 10.1038/s41477-021-01003-y.
Acknowledgements
We thank A. Hazekamp, M. E. Schranz and F. F. M. Becker for their contributions to this work. We thank C. F. Forney and T. Soomro for their assistance. This work was funded primarily by Bedrocan but was also supported by the National Science Foundation Plant Genome Research Programme grant no. 154686 to Z.M. and a Vanier Scholarship from the National Sciences and Engineering Research Council of Canada to S.W.
Data availability
The authors declare that the data supporting the findings are available within the paper. The sequence data are available in the NCBI Short Read Archive under BioProject No. PRJNA713792. The genotype files are available at 10.5061/dryad.gqnk98smm.
Code availability
All code used for the analyses is available through GitHub at https://github.com/MylesLab/cannabis-labelling.
Competing interests
R.v.V. and H.M. are employed by Bedrocan. Bedrocan funded this work, and R.v.V. played a role in the conceptualization, design, data collection, analysis, decision to publish and preparation of the manuscript. The remaining authors declare no competing interests.