Genetic insights into agronomic and morphological traits of drug-type cannabis revealed by genome-wide association studies

Sequencing and genotyping

Genetic analysis

Population structure analysis

Genome-wide association analysis

Marker-trait association analysis was performed using the method BLINK³⁹ in GAPIT v3²¹, using the 282 K high-quality SNPs and the phenotyping data for nine different traits. The identification of false positive was minimized by incorporating population structure (i.e., P matrix generated with fastStructure for K = 3) and kinship (i.e., K* matrix generated with TASSEL5) for the analysis. The threshold of significance for marker-trait associations in both methods was set to ensure a false discovery rate < 0.05, adjusted with a Benjamini–Hochberg correction. Markers with a phenotypic variance explained (PVE) less than 3% were excluded from the analysis as they were considered uninformative and of limited interest. Manhattan plots showing –log₁₀(p) distribution of markers by chromosome were generated with rMVP²² using ‘plot.type = ”m”’ and quantile–quantile (QQ) plots were created with GAPIT v3²¹ (Supplementary Fig. ^S6). Boxplot of the allelic classes of significant markers were generated with ‘ggplot2’ in R (Supplementary Fig. ^S7).

Genetic diversity in the GWAS-panel revealed by dense genotyping

To achieve comprehensive marker coverage across the Cannabis genome, an HD-GBS approach was used. Sequencing of HD-GBS libraries generated 486 M reads, averaging 2.8 M reads per sample. This extensive sequencing effort resulted in an average per-sample coverage of 7.7% of the cs10 v2 assembly, achieving a cumulative coverage of 34.1% across the entire genome for the entire population. The analysis of variant calling from our sequencing data initially yielded a substantial dataset of 2.7 M raw variants that met the quality criteria. Following filtering for missing data and minor allele frequency (MAF of 1%), we successfully identified ~ 800 K polymorphic variants, with an overall proportion of missing data reaching 61% before imputation step. This SNP catalog meets the criteria required to perform a relevant missing data imputation⁵⁶. Subsequently, we performed a secondary round of filtering, primarily aimed at retaining common variants, as defined by a MAF of 6%, retaining approximately 39% of the raw data. While this filtering step may exclude rare variants that could potentially influence complex traits, it is essential to reduce the risk of false-positive associations and ensure that a minimum of 10 accessions carries the significant allele, thereby preventing overfitting in GWAS models²³. The HD-GBS approach and the filtering procedures resulted in a catalog of 282 K high-quality SNPs (all details in Supplemental Table ^S2). Within this catalog, 25.5% of the genotypes were found to be heterozygous and the SNPs exhibited an average MAF of 21.7%. For a detailed overview of filtering steps and the number of variants retained at each stage, refer to Supplementary Table ^S3. Overall, this SNP catalog represents an extensive genetic resource for the subsequent GWAS and underscores the robustness of the genotyping strategy used in this study.

Markers were exceptionally well distributed across the genome, ensuring coverage of gene-rich regions. On average, there was one marker per every ~ 3 kb of the genome, which significantly enhances the likelihood of identifying markers in strong LD with putative candidate genes or regions (Fig. 2a). Across the entire physical map, only 12 gaps exceeding 1 Mb, with the largest being 1.2 Mb, were identified. Comparing our dataset with the RAD-Seq method used in the study of Petit et al.^32,33, by employing comparable filtration criteria, the HD-GBS approach yielded a comparable number of markers while utilizing only one-tenth of the sequencing efforts (averaging 2.8 M vs. 29.7 M reads per sample). Therefore, the density and genomic distribution of SNPs provided by the HD-GBS approach make it a cost-effective option for conducting GWAS on large Cannabis panel. Furthermore, this approach is compatible with the miniaturization of sequencing libraries using the NanoGBS procedure, which further contributes to substantial cost reduction in genotyping⁵⁷.

The average extent of LD decay to its half ranged from 22.6 to 89.0 kb across different chromosomes (Fig. 2b). It is important to note that LD decay is a relative value and does not precisely reflect to reality recombination rates throughout the entire genome, particularly between heterochromatic and euchromatic regions⁵⁸. However, this measure proved valuable for comparing the impact of domestication and selection on recombination rates among different populations. In this context, the LD observed in the GWAS-panel showed rapid decay compared to modern cultivars of comparable genome size, such as soybean (where LD may extend over 100 kb⁵⁹) and tomato (where LD can extend over 1 Mb⁶⁰). Nevertheless, LD decayed to its half more slowly compared to a recent study of 110 domesticated and landrace Cannabis accessions from various worldwide origins, where LD decayed over approximately 10 kb⁶¹. This resulted in a large number of small HBs with an average size of ~ 4 kb (Supplemental Table ^S3). It is worth noting that the LD decay on the sex chromosome was almost twice slower (p < 0.001) compared to autosomes. These observations were consistent with the recent history of Cannabis cultivation in Canada, characterized by extensive hybridization efforts by breeders with a particular focus on sexual characteristics, such as the production of female flowers².

Low level of population structure

The population structure within the GWAS-panel was assessed using the 282 K high-quality SNPs. Initially, the degree of admixture of individuals and clustering inference was estimated by fastStructure (Supplemental Table ^S4). While the model maximizing the marginal likelihood suggested a K value of 6, the optimal number of principal components (PCs) to explain the structure of the population was determined to be 3. The K value of 3 revealed two clusters (clusters 1 and 3, Fig. 2c) with low admixture compared to a K value of 6 (Supplemental Fig. ^S1), indicating a more robust assignments with more homogeneous individuals within each cluster. Using the BIC criterion, DAPC inferred three clusters (Fig. 2d, Supplemental Fig. ^S2ab, Supplemental Table ^S4). The minimal BIC values were obtained with K values ranging from 3 to 6, consistent with the optimal number of clusters determined by fastStructure, where 3 represents the minimum value. Thus, a K of 3 was chosen to explain the structure of the GWAS panel. Comparing both methods, 94.3% and 90.0% concordant assignment were observed for K values of 3 and 6, respectively (Supplemental Fig. ^S4).

Nucleotide diversity (θ_π ) across the three clusters varied from 8.44 × 10⁻⁴ to 1.20 × 10⁻³. A lower level of genome-wide genetic diversity was observed here in drug-type cannabis (mean θ_π = 1.05 × 10⁻³) compared to broader cannabis populations worldwide (θ_π = 3.0 × 10⁻³)⁶¹. This level of diversity is also lower than that found in other major crops such as soybean (mean θ_π = 1.36 × 10⁻³)⁶², rice (θ_π = 4.0 × 10⁻³)⁶³ and corn (θ_π = 6.6 × 10⁻³)⁶⁴. Relatedness analysis among individuals revealed low intra- and inter-cluster genetic diversity, with accessions appearing neither significantly similar nor significantly distant (Supplemental Fig. ^S3). This is consistent with the cumulative variance explaining genetic variation in the population, showing gradual increase with number of retained PCs up to 176 PCs (number of accessions in the GWAS panel) rather than reaching a plateau (Supplemental Fig. ^S2a). Despite the overall genetic homogeneity, significative differences were observed between clusters for traits such as SM, HM, height, NC and ILI (Supplemental Fig. ^S5). In particular, cluster K3 exhibited significant differences from the cluster K1 for these five traits, while the cluster K2 displayed intermediate trait values between K1 and K3. In different studies, similar clustering patterns related to drug-type and hemp-type accessions^61,65–67 or geographic origins⁶⁸ were documented, where each clusters grouped independently, albeit with low intra- and inter-cluster genetic diversity. Due to limited information on the pedigree of the GWAS panel, no correlation was observed between cluster assignment and geographic or germplasm origins. Additionally, no correlation was observed between the clustering and cannabinoid composition (data not shown) of these accessions.

The limited genetic diversity observed in cultivated drug-type Cannabis has historically been attributed to intensive clandestine breeding practices since the 1970s², coupled with the impact of the war on drugs, which led to the destruction of many plants and seeds, effectively reducing the gene pool^5,69. Despite the limited genetic diversity, Cannabis exhibits a remarkable phenotypic variation that are highly desirable for breeding programs. Hence, it could be hypothesized that a portion of the observed phenotypic variations in Cannabis may be attributed to transcriptional variations, along with potential contributions from epigenetic factors. In both plants and animals, factors such as variation in the number of gene copies (CNVs)⁷⁰, epigenetic elements⁷¹, and the insertion/deletion of transposable elements (TEs) in gene control regions⁷², impact phenotypic diversity, especially those crucial in domestication and breeding⁷³. Therefore, an associated SNP may be in strong LD with either a candidate gene, where an allelic variant alters the phenotype, or with a regulatory region that either enhances or suppresses the expression of the phenotype⁷⁴.

The constrained availability of germplasm resources and low genetic diversity observed in Cannabis pose significant limitations for breeding, which, in turn, hinder innovation and the long-term sustainability of the crop⁷. In contrast to other crops where wild-type or landrace varieties are promising genetic pools to enrich genetic diversity in breeding programs⁷⁵, the situation in Cannabis is more complex. Although hemp-type and drug-type Cannabis genetically diverged⁷⁶, they still share a considerable common pool of genetic variation, limiting the ability to mine rare alleles⁶⁵. Given the growing demand for cannabis products, there is a critical necessity to pinpoint suitable genetic resources that can not only support production but also serve as a source of genetic diversity to help ongoing breeding efforts⁷.

Identification of genomic regions controlling key agronomic and morphological traits

The GWAS analysis was performed using the method BLINK with the incorporation of population structure (P) and cryptic relatedness (K*) as covariates to minimize the risk of false-positive associations. In total, 18 markers associated with the nine traits were identified (Fig. 3, Table ​Table1).1). For all significant markers identified, the three genotypes were observed (Supplemental Fig. ^S7). Six of these SNPs (SNP_1, 4, 7, 8, 9 and 11; Table ​Table1)1) demonstrated significant phenotypic impact, with the proportion of phenotypic variance explained (PVE) ranging from 18 to 45% while the remaining identified markers have a modest influence on the phenotype (PVE < 10%). Interestingly, several SNPs associated with different traits were located in close proximity to each other. For instance, SNP_9and _17 were situated within a region of about 38 kb on chromosome 1 (Chr01: 87456694–87494979) and were associated with ILI and height. The identification of 2 SNPs associated with correlated traits is consistent and suggests that this region of chromosome 1 plays a crucial role in modulating plant size in the GWAS panel. These markers are associated with key characteristics for Cannabis cultivation and are therefore of particular interest to breeders and growers. For instance, markers associated with smaller size can be advantageous for maximizing indoor cultivation, where smaller plants are preferred. Regarding the markers associated with a shorter flowering or maturation, they are advantageous for cultivators aiming for a quicker crop turnover. Similarly, the allele T at Chr09:59690286 (SNP_4) is associated with reduces canopy size and slightly the height, which can help maximize plant density in cultivation.

Given that the traits under study appear to be governed by a complex genetic control involving multiple QTL, BLINK appeared to be the most suitable method as it can capture intricate interactions among several loci through multi-locus analysis³⁹. Furthermore, this method has proven its effectiveness with large catalog of SNPs^77–79 and was ranked as the most statistically powerful method for multi-locus analyzes for GWAS in plants^21,39, 80. As the identification of high-value markers for Cannabis is in its early stages, the practical implementation of these markers by breeding programs will nevertheless require preliminary cross-validation. This can be achieved through meta-GWAS⁸¹, QTL mapping with biparental population and BSA. Additionally, comprehensive functional analyses of the candidate genes will be crucial. .

The authors wish to thank Fuga Group Inc. for supporting this project. The authors also extend their sincere appreciation to Justine Richard-Giroux, Rosemarie Boulanger and Sean Kyne for their valuable contributions to the tedious phenotyping of the GWAS-panel.