Genetic architecture of phenological, morphological, and phytochemical traits in Cannabis landraces

Authors: Mehdi Babaei, Davoud Torkamaneh

Abstract

Despite its long history of cultivation and diverse applications, Cannabis sativa remains underexplored at the genomic level, particularly in landrace populations that harbor untapped genetic diversity. In this study, we investigated the genetic architecture of 145 Iranian cannabis landrace accessions, including both male and female plants, using 233K common SNPs and genome‐wide association studies. Our analysis revealed three genetically distinct subpopulations shaped by geography, climate, and traditional cultivation practices. We identified 91 significant genomic regions associated with 40 phenological, morphological, and phytochemical traits, including 15 key loci with pleiotropic effects linked to multiple traits, including flowering time, plant architecture, biomass accumulation, and cannabinoid biosynthesis. These findings highlight the complex interplay between developmental and metabolic pathways in cannabis. The high heritability of most traits and rapid linkage disequilibrium decay underscore the potential of these landraces for high‐resolution mapping and genetic improvement. This work provides a valuable genomic resource for marker‐assisted selection, supporting the development of improved cultivars with tailored cannabinoid profiles and agronomic traits.

INTRODUCTION

Cannabis (Cannabis sativa L.), an annual flowering plant belonging to the Cannabaceae family, stands as one of the earliest domesticated crops in human history, with its origins tracing back approximately 3000–8000 years to East/Central Asia (S. A. Ahmed et al., 2008; Babaei et al., 2022; Babaei, Boissinot, et al., 2025; Crocq, 2020; G. Ren et al., 2021; Schilling et al., 2020; Small, 2017; Srinivasababu, 2014). Throughout millennia, cannabis has been extensively cultivated and utilized for diverse applications, including fiber production, oil extraction, and its distinct medicinal and psychoactive properties (Anwar et al., 2006; Barcaccia et al., 2020; Russo et al., 2008; Sun, 2023; Warf, 2014). Additionally, early Iranian medical texts authored by physicians, such as Rhazes (854–925) and Avicenna (980–1037), further attest to the therapeutic use of cannabis (Bachir et al., 2022; Mahdizadeh et al., 2015). This rich history has led to a wide array of genetic and phenotypic diversity across its indigenous populations, or landraces, which have adapted to various environments and human selection pressures (Babaei & Ajdanian, 2020; Babaei et al., 2024; Babaei, Boissinot, et al., 2025; Babaei, Nemati, et al., 2025).

Despite its historical significance and versatile applications, scientific research into cannabis genetics and trait inheritance has lagged compared to other major crops, largely due to decades of prohibition and clandestine breeding practices (Halpin‐McCormick et al., 2024; Mudge et al., 2018; Peng & Shahidi, 2021; Torkamaneh & Jones, 2021). While sampling bias in earlier studies often overrepresented drug‐type varieties, recent genomic analyses have revealed substantial genetic diversity across cannabis populations (Aina et al., 2025; Carlson et al., 2021; Lynch et al., 2016, 2025; G. Ren et al., 2021). Whether prohibition led to genetic diversity reduction or its preservation through limited commercial consolidation of germplasm remains debated (de Ronne & Torkamaneh, 2025; de Ronne et al., 2024; Lapierre, de Ronne, et al., 2023; Naim‐Feil et al., 2021, 2022, 2023). However, with shifting legislation and increasing legalization globally, the cannabis market is experiencing unprecedented growth, driving a pressing need for a deeper understanding of its fundamental biology and genetic architecture (Hammond, 2021; Statista, 2024).

Cannabis is predominantly dioecious (male XY, female XX), with a diploid genome (2n = 20), although monoecious forms and hermaphroditism can occur (Braich et al., 2020; Carpentier et al., 2012; Monthony et al., 2024; Razumova et al., 2016). Its taxonomy remains a subject of ongoing debate, with perspectives ranging from a single species (C. sativa) to multiple distinct species (C. sativa, Cannabis indica, and Cannabis ruderalis) (Anwar et al., 2006; Flores‐Sanchez & Verpoorte, 2008; Lapierre, Monthony, et al., 2023; McPartland, 2018; Small & Beckstead, 1973). Modern classifications often extend beyond these taxonomic definitions, incorporating legal status (hemp vs. drug‐type based on Δ^9^‐tetrahydrocannabinol (THC) concentration < 0.3%), phytochemical profiles (chemotypes based on allelic status of cannabinoid synthase Type I THC‐dominant, Type II intermediate THC:CBD, Type III cannabidiol [CBD]‐dominant, Type IV cannabigerol [CBG]‐predominant, and Type V cannabinoid‐free), ecological adaptation (ecotypes), and biological characteristics (biotypes) (Babaei & Ajdanian, 2020; Cherney & Small, 2016; De Meijer & Hammond, 2005; Jang et al., 2020; Lapierre, Monthony, et al., 2023; McPartland & Small, 2020; Sawler et al., 2015; Small, 2015).

Recent policy shifts have expanded access to genomic and transcriptomic data and streamlined regulations for cannabis research (Grassa et al., 2021; Hesami et al., 2020, 2024; Monthony et al., 2024). This has catalyzed the adoption of advanced sequencing technologies, such as next‐generation sequencing, which offer cost‐effective, high‐throughput data generation. Coupled with sophisticated bioinformatic tools, these innovations have greatly expanded the scope of genotype–phenotype association studies across diverse crop species (Abdi et al., 2024; A. Ahmed, 2024; de Ronne et al., 2023; Satam et al., 2023; Torkamaneh et al., 2016, 2018). A critical challenge in cannabis breeding lies in accurately evaluating the genetic contributions from diverse or exotic and indigenous germplasms. This difficulty is compounded by the quantitative nature of many desirable traits and the significant environmental influence on plant performance, factors that have historically limited the widespread incorporation of such valuable genetic resources (Barcaccia et al., 2020; Ingvardsen & Brinch‐Pedersen, 2023). Over the past decade, genome‐wide association studies (GWASs) have emerged as a powerful approach proving highly effective in dissecting the genetic basis of variation for intricate phenotypes, including plant physiological and agronomic characteristics (Belzile & Torkamaneh, 2022). GWAS is now recognized as a premier tool for pinpointing genetic markers linked to traits of interest, particularly excelling where conventional methods fall short for complex traits (Alqudah et al., 2020; Belzile & Torkamaneh, 2022; Torkamaneh & Belzile, 2022).

In cannabis, flowering time is a critical agronomic trait, influenced by photoperiod (short‐day plant) and strong genetic control (Amaducci et al., 2008; Babaei et al., 2022, 2024; Lisson et al., 2000; M. Zhang et al., 2021). Flowering directly impacts biomass accumulation, fiber quality, and cannabinoid production, with cannabinoids rapidly accumulating during early flowering stages (Amaducci et al., 2005, 2008; Petit, Salentijn, Paulo, Denneboom, Trindade, et al., 2020; Petit, Salentijn, Paulo, Thouminot, et al, 2020; Salentijn et al., 2015). Recent genetic studies have advanced our understanding of cannabis flowering, identifying key loci like Autoflower1 and Early1 (Toth et al., 2022), and Autoflower2 (FT1) (Dowling et al., 2024) alongside multiple quantitative trait loci (QTL) associated with photoperiod sensing, circadian rhythms, and hormone signaling (Petit, Salentijn, Paulo, Denneboom, Trindade, et al., 2020). Beyond flowering, association mapping analysis (e.g., GWAS and QTL mapping) has been instrumental in elucidating genetic markers linked to a wide spectrum of cannabis traits, encompassing sex determination, cannabinoid biosynthesis pathways, fiber characteristics, key morphological and agronomic features, and stress resilience (de Ronne et al., 2024; Petit, Salentijn, Paulo, Denneboom, Van Loo, et al., 2020, 2020b; Sun et al., 2023; Toth et al., 2022; Welling et al., 2020).

However, a key limitation in cannabis GWAS to date, particularly those using commercial cultivars, is their restricted genetic diversity. While valuable for specific traits, such studies often miss the broader genetic variation in landraces crucial for comprehensively dissecting complex traits. For instance, recent studies by de Ronne et al. (2024) on morphological traits and de Ronne and Torkamaneh (2025) on cannabinoid profiles, both conducted on Canadian commercial cultivars, identified extensive genetic variants within cultivated lines but were inherently limited by the narrower diversity of cultivated lines compared to landraces.

This study aims to address these gaps by providing a comprehensive characterization of 145 Iranian cannabis landrace accessions. We employed high‐density genotyping‐by‐sequencing (HD‐GBS) to generate a catalog of 233,624 high‐quality single nucleotide polymorphisms (SNPs) and GWAS to identify significant genetic loci and putative candidate genes associated with 42 key phenological, morphological, and phytochemical traits. Our findings will contribute significantly to the fundamental understanding of cannabis biology and provide valuable genetic resources and tools for accelerating molecular breeding efforts.

MATERIALS AND METHODS

Plant materials

Cannabis seed samples from 25 native populations (Table S1) of Iran were sourced from local markets and native farmers in various regions, ensuring genetic diversity and regional adaptation. The regions from which the seeds were collected encompassed five defined climatic zones based on the Köppen–Geiger climate classification method, with geographic latitudes between 25° and 40° north and longitudes 45° and 65° east.

Experimental design and growing conditions

Twenty seeds from each collected population were planted (Babaei et al., 2024). Thirty days after sowing (DAS), the seedlings were transferred to growth bags with a mixture consisting of garden soil, leaf mold, sand, and perlite in a 1:1:1 proportion. Based on the study by Amaducci et al. (2008), the planting was done 60 days before the photoperiod switch‐off in Mashhad, Khorasan, Iran, when the days start to shorten (the day length at sowing was 13 h and 17 min) to ensure that the plants would have an adequate vegetative period before entering the reproductive phase. The plants were arranged in a randomized complete block design with three blocks and five observations, totaling 375 plants (50% male and 50% female), in the greenhouse complex of Ferdowsi University of Mashhad, Iran (36^°^16ʹ N and 59^°^36ʹ E with an altitude of 985 m). According to spatial analysis, each unit, measuring 60 m^2^, was divided into 25 rows and 15 columns, where each block included five columns. The aisles within each unit were arranged in order to place 16 plants m^−2^. At the beginning of the reproductive stage (appearance of solitary flowers), male and female plants were separated and transferred to a unit with similar conditions, reducing the number of plants per square meter by half. The growing conditions and all cultural practices were previously described in Babaei et al. (2024).

Phenotypic data collection

Phenological and morphological data for this study were primarily derived from a previously published comprehensive phenotypic characterization of these Iranian cannabis landraces (Babaei et al., 2024). In that prior work, 12 phenological traits (covering both vegetative and reproductive stages) and 14 morphological traits were recorded on an individual basis for all 375 plants. For the current GWAS analysis, we specifically utilized six reproductive phenological stages and 12 morphological traits from this extensive dataset. A complete list of all phenological, morphological, and phytochemical (cannabinoid) traits, along with their defined categories, is provided in Table S2.

Phenological traits

Based on the detailed phenological descriptor developed in our previous study (Babaei et al., 2024), phenological traits were recorded and defined for both male and female plants. For the present study, the focus was specifically on six reproductive stages. These traits, measured in DAS, include (1) GV point (GVP), transition of bud phyllotaxis from opposite to alternate on main stem, marking entry into the reproductive phase (minimum 0.5 cm distance between alternate leaf petioles); (2) start flower formation time in individuals (SFFI), first bud (solitary flower) appearance on individual plants (bell‐shaped or closed sepals for male; symmetrical calyx with two styles for female); (3) start flower formation time in 50% population (SFFP), first flower appearance in 50% of the population; (4) start 10% flowering time in individuals (SF10I), with a minimum 10% of main inflorescence formed on individual plants; (5) start 10% flowering time in 50% population (SF10P), 10% of inflorescence formed in 50% of the population; 6) flowering time 50% in individuals (FT50I), 50% of the main inflorescence formed on individual plants.

Morphological traits

A total of 21 morphological traits were assessed in this study, categorized as node and branching architecture, growth and structural dimension, and biomass yield. These measurements were conducted at the end of the cultivation period, coinciding with the harvest of both male and female plants (Babaei et al., 2024). A schematic illustration of these key morphological traits is provided in Figure 1a. Node and branching architecture traits included (1) number of nodes to the main inflorescence (NTMI); (2) number of nodes to the first lateral shoot (NTFIS); (3) number of nodes on the main stem on harvest day (NNH); (4) number of nodes on the main inflorescence (NMI); and (5) number of lateral shoot (NLS).

The NNH was calculated as the sum of nodes on the main inflorescence and nodes to the main inflorescence (Equation 1). The NLS was determined based on the NTMI and the first lateral shoot (Equation 2). (1)NNHn=NMI+NTMI (2)NLSn=NTMI−NTFIS×2

Growth and structural dimension traits included (1) length of shortest lateral shoot (LSLS); (2) length of main inflorescence (LMI); (3) length of longest lateral shoot (LLLS); (4) length of internode in the middle third of the main stem on harvest day (LIMTH); (5) height on harvest day (HH); (6) height to GV point (HGV); (7) relative growth rate (RGR); and (8) stem diameter on harvest day (SDH).

Traits such as LMI, HH, LIMTH, and HGV were measured using a tape measure with 0.01 m precision. SDH was measured at 5 cm above the soil surface using calipers with an accuracy of 0.01 mm. The RGR was also calculated using the following equation (Equation 3), and according to mg g^−1^ day^−1^, (3)RGR=lnTDW2−lnTDW1t2−t1

In this equation, the biomass increases from TDW1 to TDW2, on average during a time interval of t1‐t2. TDW2 and TDW1 are the dry weights of plants at the time of harvest and in the seedling transplantation stage, respectively. Consequently, t 2 and t 1 represent the time of harvest and seedling transfer, respectively. 1n signifies the natural logarithm (Babaei et al., 2024; Hoffmann & Poorter, 2002).

Biomass yield traits included (1) dry weight of flowers (DWF); (2) fresh weight of flowers (FWF); (3) fresh weight of leaves (FWL); (4) dry weight of leaves (DWL); (5) fresh weight of stems (FWS); (6) dry weight of stems (DWS); (7) total fresh weight (TFW); (8) total dry weight (TDW). A digital scale (AND‐GF3000) with 0.001 gr accuracy was used for all fresh weight measurements of individual plant parts. For dry weight determination, samples were air‐dried at room temperature (25°C) for 15 days before being re‐weighed using the same digital scale. TFW was calculated as the sum of fresh weights of stems, leaves, and flowers (Equation 4). Similarly, TFW was calculated as the sum of fresh weights of stems, leaves, and flowers (Equation 5), (4)TFWg=FWS+FWL+FWF (5)TDWg=DWS+DWL+DWF

Cannabinoid profiles

To prevent pollination, which is known to reduce cannabinoid content, male plants were separated from female plants prior to the onset of flowering (Lipson Feder et al., 2021). Sampling for cannabinoid analysis was conducted on female plants at the stage of flower maturity. This maturity was determined by daily visual inspection, specifically observing the browning of stigmas and the change in trichome color from clear to amber/brown, indicating peak cannabinoid production (Punja et al., 2023; Tran et al., 2025). Individual samples were collected from the main inflorescence of each plant as it reached this mature stage. Harvested samples were immediately placed in separate paper bags and air‐dried in a dark environment at 25°C until completely dry. Subsequently, dried samples from each population within each experimental block were homogenized using a ceramic mortar.

Extraction procedure

First, 0.05 g of dried plant tissue was extracted with 2 mL of methanol/chloroform (9:1 ratio). This mixture was then sonicated for 40 min in an ultrasonic bath, followed by centrifugation at 10,000 rpm for 15 min at 10°C. The supernatant was carefully separated and thoroughly dried by air (De Backer et al., 2009). The dried extracts were then solubilized in 2 mL of 80% methanol and sonicated for 10 min in an ultrasonic bath. After a final centrifugation for 5 min at 18,000 g at room temperature, the supernatants were filtered through a 0.22‐µm nylon filter and subsequently diluted by factors of 2 and 40 for analysis, following a modified method described by Mudge et al. (2018).

Ultra performance liquid chromatography analysis

Biochemical analysis of these extracts was conducted at the Metabolomics Platform in the Institute of Nutrition and Functional Foods, Université Laval, Québec, QC, Canada. Cannabinoid analysis was performed using an ultra performance liquid chromatography Acquity H‐Class system (Waters Corporation) coupled to an Acquity TUV detector (Waters Corporation). Compound separation was achieved on a Cortecs 1.6 µm, 2.1 mm × 150 mm (Waters Corporation) column, maintained at 30°C. The mobile phase consisted of (A) 20 mM ammonium formate at pH 2.92 and (B) 100% acetonitrile. The gradient program was as 0–6.4 min, 76% B; 6.5–8 min, 99% B. Conditions were reinitialized to 76% B from 8.1 to 10 min. The flow rate was 0.45 mL min^−1^, and the injection volume was 2 µL. Detection was set at a wavelength of 228 nm. Cannabinoids were quantified using 5‐point calibration curves prepared in the 1–100 mg L^−1^ range for all standards.

The cannabinoid profiles quantified in this study, expressed as concentration (% w/w), included Δ^9^‐tetrahydrocannabinolic acid (THCA), Δ^9^‐THC, cannabidiolic acid (CBDA), CBD, cannabigerolic acid (CBGA), CBG, cannabichromene (CBC), cannabinol (CBN), tetrahydrocannabivarin (THCV), cannabidivarin (CBDV), and Δ^8^‐tetrahydrocannabinolic (Δ^8^‐THC). It is important to note that Δ^8^‐THC was consistently below the limit of detection in all samples and was therefore excluded from subsequent analyses. To determine the total THC potential and total CBD potential, which represent the total decarboxylated forms, were calculated using the following equations (Equations 6 and 7, respectively): (6)TotalTHCpotential%w/w=Δ9THC+0.877×THCA (7)TotalCBDpotential%w/w=CBD+0.877×CBDA

The coefficient 0.877, used in both equations, is derived from the molar mass ratio between the decarboxylated cannabinoids (Δ^9^‐THC and CBD; 314.5 g mol^−1^) and their acidic precursors (THCA and CBDA; 358.5 g mol^−1^), accounting for the loss of the carboxyl group (CO2) during decarboxylation.

Additionally, THC:CBD and CBD:(THC + CBN) ratios were calculated to provide insights into the chemotype classification. Populations were categorized into four classes using dual‐ratio THC‐dominant (Type I), balanced (Type II), moderate CBD (Type II), and CBD‐dominant (Type III). This approach maintains the standard chemotype framework (De Meijer et al., 2003; De Meijer & Hammond, 2005), with “moderate CBD” representing a CBD‐leaning subcategory within chemotype II (heterozygous CBDAS/THCAS) to better capture variation relevant for germplasm characterization and breeding applications. A decision tree illustrating this hierarchical classification based on predefined thresholds for these cannabinoid ratios is presented in Figure 1b.

Genotyping and SNP calling

Genotyping

Sampling was performed from healthy and young leaves of all individuals. Samples were then dried using the Freeze Dry System Alpha 2–4 LD plus for 24 h. Individuals of the same sex within each block were pooled prior to DNA extraction. A total of 145 pooled samples were ground with metallic beads in a RETSCH MM 400 mixer mill (Fisher Scientific). DNA extraction was carried out using the Qiagen DNeasy Plant Mini Kit. The extracted DNA quantity was assessed by a Qubit fluorometer with the dsDNA HS assay kit (Thermo Fisher Scientific), and their quality was randomly evaluated using agarose gel electrophoresis. DNA concentrations were adjusted to 10 ng µL^−1^ for all samples. Final DNA samples were used to prepare HD‐GBS libraries with BfaI as described in Torkamaneh et al. (2021) at the Institut de biologie intégrative et des systèmes, Université Laval, QC, Canada. Sequencing was conducted on an Illumina NovaSeq 6000 with 150 paired‐end reads at the Genome Quebec Service and Expertise Center (CESGQ), Montreal, QC, Canada.

SNP calling and filtration

Sequencing data were processed with the Fast‐GBS v2.0 using the C. sativa cs10 v2 reference genome (GenBank acc. no. GCA_900626175.2) (Grassa et al., 2018; Torkamaneh, Laroche, et al., 2020). For variant calling, a prerequisite of a minimum of six reads to call an SNP was opted. After mapping against the cs10 v2 reference genome, an initial dataset of over 4.5 million variants was obtained. Raw SNP data were then filtered with VCFtools to remove low‐quality SNPs (QUAL < 10 and mapping quality < 30) and variants with proportion of missing data exceeding 80% (Danecek et al., 2011). Following the first filtering step, which included quality filtering for minor allele frequency (MAF > 1%) and missing data > 80%, about 763 K variants were retained. A subsequent filtering step applying a minimum and maximum allele count of two reduced the dataset to ∼584 K variants. A second round of filtration was then applied, retaining only biallelic variants with heterozygosity ˂50% and a MAF of >5%. Additionally, variants residing on unassembled scaffolds were removed. After these stringent filtering steps, the dataset was reduced to 233,624 high‐quality SNPs on the scaffold, which were retained for downstream analysis (Table S3).

Genetic diversity analysis

SNP‐based genomic summary

Genetic diversity parameters were calculated across the entire panel (by chromosome, sex, and cluster/subpopulation). Read count and coverage were quantified from aligned reads using SAMtools (Danecek et al., 2021). MAF and heterozygosity level were determined using TASSEL v5.2.94 (Bradbury et al., 2007). Nucleotide diversity (θπ) was measured in 1 kb sliding windows across the genome using the—window‐pi option of VCFtools (Danecek et al., 2011). The proportion of SNPs located on annotated genes was determined by intersecting the filtered VCF file with a gene annotation BED file using BEDTools v2.31.1 (Quinlan & Hall, 2010). Genome‐wide gaps larger than 1 Mb were identified by subtracting the genomic regions covered by SNPs from the total chromosome lengths using BEDTools v2.31.1 (Quinlan & Hall, 2010). To visualize the distribution of SNP density, a plot was produced with rMVP using the plot.type = “d” parameter, in combination with the gene density distribution (Yin et al., 2021). Additionally, variant types (e.g., multi‐nucleotide polymorphisms [MNPs], insertions, and deletions) and their predicted functional impacts (e.g., intergenic, intronic, synonymous, and missense), along with the overall transition‐to‐transversion (Ts/Tv) ratio, were assessed using SnpEff v5.2e based on the Cannabis sativa reference genome (cs10 v2; GenBank accession, GCA_900626175.2) (Cingolani et al., 2012; Grassa et al., 2018).

Linkage disequilibrium decay and haplotype blocks

Linkage disequilibrium (LD) decay was assessed by calculating the squared allele frequency correlation (r ^2^) using PopLDdecay v3.41 (C. Zhang et al., 2019). LD decay was calculated for the entire panel, for each chromosome, for each sex and for each genetic subpopulation, with a maximum physical distance of 500 kb (‐MaxDist 500). Haplotype blocks (HBs) were identified using PLINK v1.90b5.3 by computing pairwise LD (r ^2^), considering a window of 999 kb SNPs –ld‐window‐kb 999 (Purcell et al., 2007).

Population structure analysis

Population structure inference

Population structure and admixture were determined using the variational Bayesian inference algorithm implemented in fastStructure v1.0 for K (number of subpopulations) values from 1 to 10 (Raj et al., 2014). The complete set of 233,624 high‐quality SNPs was used for this analysis. The optimal K was estimated using the ChooseK tool, and admixture proportions were visualized with Distruct v2.3. Discriminant analysis of principal component (DAPC) was performed using the R package adegenet version 2.1.10 (Jombart et al., 2010). The optimal number of clusters (K) was estimated using the find.clusters function, exploring up to 40 clusters, and determined by the minimal Bayesian information criterion (BIC). To visualize the DAPC using the “scatter” function, the optimal number of principal components (PCs) (22 PCs) was estimated with two cross‐validation procedures using optim.a.score and xvalDapc (exploring PCs up to a maximum of 144).

Phylogenetic and kinship analysis

Phylogenetic relationships among accessions were inferred using the neighbor‐joining method in MEGA version 11 (Jin & Nei, 1990; Kumar et al., 2018). This was based on genetic distances calculated from SNP data, using three different nucleotide substitution models, Jukes–Cantor (TC), Tajima–Nei (TN), and maximum composite likelihood. A bootstrap consensus tree was formed from 1000 replicates. The resulting phylogenetic tree was edited and beautified using iTOL version 7.2 (Letunic & Bork, 2021). Genetic distances were also estimated in MEGA version 11 using the p‐distance model with Gamma distributed (G) rates among sites (Gamma parameter = 1) (Kumar et al., 2018). Pairwise distances, within‐group average distances, between‐group average distances, and overall mean distances were calculated. For all distance calculations, a partial deletion approach was used for gaps/missing data with a site coverage cutoff of 95%. Variance estimation for these distances was performed using 1000 bootstrap replicates. The pairwise genetic distances were subsequently visualized as a heatmap in R using the pheatmap package (Kolde & Kolde, 2015; Team, 2020). Kinship analysis was performed using TASSEL v5.2.94 with the Centered_IBS method, and the kinship matrix was plotted with GAPIT v3 (Bradbury et al., 2007; J. Wang & Zhang, 2021).

Genome‐wide association analysis

Genetic associations between markers and traits were investigated using the Bayesian‐information and linkage‐disequilibrium iteratively nested keyway (BLINK) method implemented in GAPIT (Huang et al., 2019; J. Wang et al., 2022). This analysis focused on 233,624 high‐quality SNPs and phenotypic data collected for 42 distinct traits. Phenological and morphological trait values from Babaei et al. (2024) were averaged across pooled individuals to correspond with genotypic data, while cannabinoid profiles were measured directly on the pooled samples. To mitigate the occurrence of false positive associations, the analysis incorporated both population structure (represented by a P matrix derived from fastStructure with K = 3) and kinship (a K* matrix generated using TASSEL v5.2.94) as statistical covariates (Bradbury et al., 2007; Raj et al., 2014). A stringent significance threshold for marker‐trait associations was established at a false discovery rate (FDR) < 0.05, which was achieved through adjustment using the Benjamini–Hochberg correction (Benjamini & Hochberg, 1995). Furthermore, markers explaining ˂3% of the proportion of phenotypic variance explained (PVE) were excluded from the final analysis, as they were deemed to provide limited informative value. Visual representations of the GWAS results included Manhattan plots, which depicted the −log10(p) distribution of markers across chromosomes (generated using rMVP with plot.type = “m”), and quantile–quantile plots, created with GAPIT, to assess model fit (J. Wang & Zhang, 2021; Yin et al., 2021). Additionally, boxplots illustrating the phenotypic effects of different allelic classes for significant markers were generated using the ggplot2 package in R (Team, 2020; Wickham et al., 2016). A comprehensive circular genomic map illustrating the genome‐wide distribution of significant SNP markers and their associated traits was generated using shinyCircos‐V2.0 (Y. Wang et al., 2023).

Candidate gene identification and annotation

For the identification of putative candidate genes, HBs were defined based on significant SNP markers identified through GWAS. Given the genetic diversity in cannabis, only markers exhibiting high LD (r ^2^ ≥ 0.75) with the significant SNPs were retained to delineate these HBs (de Ronne et al., 2024). Genes situated within these defined HBs, delimited by the 5′‐most and 3′‐most SNP positions of each block, were considered as putative candidate genes. Haploview v4.1 was used for visual inspection of marker positions within HBs, with the required input formats generated using PLINK's –recode HV function (Barrett et al., 2005). Functional annotation and interpretation of these candidate regions were performed using SnpEff v5.2e, leveraging the Cannabis sativa reference genome (cs10 v2; GenBank accession, GCA_900626175.2) (Cingolani et al., 2012; Grassa et al., 2018). Additionally, gene ontology terms were extracted from the National Center for Biotechnology Information (NCBI) Cannabis sativa Annotation Release 100 to characterize the potential biological roles and pathways associated with these candidate genes.

Statistical analysis

The data were first evaluated using the Outlier Grubbs approach to detect and remove any outliers (Adikaram et al., 2015). Basic descriptive statistics for all traits were calculated to summarize their distribution and variability. One‐way analysis of variance (ANOVA) was performed to assess the significance of trait variation among accessions and subpopulations, followed by post hoc comparisons between clades at the 95% confidence level according to Tukey's honestly significant difference test. Skewness and kurtosis analyses, assessing data normality and residual errors based on the distribution type, were conducted using Minitab Statistical Software 21 (Cain et al., 2017; Minitab, 2021). Missing data were imputed using the Classification And Regression Tree method implemented in the AllInOne package in R (Team, 2020; Yoosefzadeh Najafabadi et al., 2023). For phenological and morphological data, spatial analysis (AR1⊗AR1) was applied to account for potential spatial effects and minimize environmental impact, considering the experimental layout (15 columns and 25 rows). The data‐fitting approach considered population as a fixed factor and block as a random factor. The model's performance was evaluated using four primary metrics, root mean square error, mean squared error, normalized root mean squared error, and coefficient of variation (CV) (Babaei et al., 2024). Both broad‐sense heritability (H ^2^) and SNP‐based heritability (h ^2^) were estimated for all traits to quantify the genetic contribution to phenotypic variation. H ^2^ was determined by the AllInOne package, and h ^2^ was estimated using GCTA software with the restricted maximum likelihood method, based on the genomic relationship matrix (Yang et al., 2011; Yoosefzadeh Najafabadi et al., 2023). Pearson's correlation coefficients were calculated to assess the interrelationships between all phenotypic traits, utilizing the package corrplot in R (Wei et al., 2017). A heatmap of cannabinoid profiles was generated using the pheatmap package in R (Kolde & Kolde, 2015). Principal component analysis (PCA) was conducted using the Factoextra package in R (Kassambara & Mundt, 2017). Boxplots, violin plots, and the heatmap of the confusion matrix were generated using Python 3.12 with the matplotlib and seaborn packages (Bisong, 2019).

RESULTS

Phenotypic characterization of cannabis landraces

A comprehensive phenotypic analysis was conducted on 145 cannabis landrace accessions, comprising 72 females and 73 males, encompassing a wide range of 42 phenotypic and phytochemical traits. Statistical differences among accessions were assessed by ANOVA. A comprehensive list of evaluated traits is provided in Table S2, with overall distributions (Figure 2) and analysis of trait interrelationships through Pearson's correlation coefficients (Figure S1).

Phenological traits

A comprehensive analysis of six phenological stages revealed significant variation in reproductive phase traits (Table S4; Figure 2). The average timing for phenological stages ranged from 69.2 ± 9.3 days for SFFP to 117.1 ± 24.6 days for FT50I. Traits exhibited significant differences, with GVP showing significance at p ≤ 0.01, and the remaining five stages (SFFI, SFFP, SF10I, SF10P, and FT50I) at p ≤ 0.001. Skewness and kurtosis analyses indicated non‐normal distributions, particularly for SFFI and SFFP, suggesting tendencies towards earlier flowering and presence of extreme phenotypes. Furthermore, the H ^2^ for these traits ranged from 0.41 (GVP) to 0.60 (FT50I).

Morphological traits

The analysis of twenty‐one morphological traits demonstrated considerable variation across the accessions (Table S5; Figure 2). All morphological traits exhibited significant differences among the accessions (p ≤ 0.001), except FWF and DWF, which were not statistically significant. Key structural traits showed substantial diversity, such as plant height (HH) with a mean of 153.79 ± 34.62 cm and stem diameter (SDH) at 13 ± 2.83 mm. High variability was observed in traits like LSLS (CV% 37.54) and TFW (ranging 359.95 g; mean 215.88 ± 86.6 g). The H ^2^ for morphological traits varied widely, ranging from 0.41 (LIMTH) to 0.80 (TFW), suggesting a broad spectrum of genetic influence.

Cannabinoid profiles

Analysis of fifteen cannabinoid profiles revealed substantial variability in accumulation among the landrace accessions (Table S6; Figure 2). Significant variation was observed for Δ^9^‐THC (p ≤ 0.01), CBN (p ≤ 0.05), CBDV (p ≤ 0.01), total THC potential (p ≤ 0.05), and the CBD:(THC + CBN) ratio (p ≤ 0.001). Δ^8^‐THC was undetectable in all samples. Individual cannabinoid concentrations and ratios showed high CV% values (e.g., CBG 99.98%, CBDV 164.5%), reflecting diverse chemotypes. Skewness and kurtosis analysis indicated non‐normal distributions, often with high kurtosis, suggesting a deviation from normality characterized by a greater concentration of values in the tails. Heritability was generally high, with H ^2^ ranging from 0.46 to 0.91 and h ^2^ from 0.38 to 1.00. Total THC Potential (H ^2 ^= 0.90, h ^2 ^= 1.00) and Total CBD Potential (H ^2 ^= 0.74, h ^2 ^= 0.64) exhibited strong genetic control over primary cannabinoid biosynthesis.

Genotyping and population structure

High‐Density SNP coverage of the GWAS panel

Sequencing libraries generated ∼400 M reads, averaging approximately 2.7 M reads per sample. The final catalog of 233,624 (211,621 SNPs; 7500 MNPs, 6588 insertions, and 7915 deletions) common variants (MAF > 5%) displayed a density of 273.5 markers per Mb (one marker every ∼3.6 kb), with 20.8% of the genotypes being heterozygous and an average MAF of 12.6%. Notably, 21.2% of variants were located on annotated genes. Across the genome, fifteen gaps larger than 1 Mb were identified, with the largest being 2.3 Mb (Figure S2a; Table S7). Additionally, a sex‐based analysis revealed that males had a higher average MAF (13.61% vs. 11.48%) and heterozygosity (22.58% vs. 18.94%) compared to females (Table S8). Analysis of 233,624 variants revealed that 95.1% were modifiers with minimal functional impact, primarily in non‐coding regions [intergenic (40.4%), intronic (11.4%), upstream (13.8%), and downstream (16.2%)]. Among coding variants, 60.8% were synonymous and 38.7% missense, yielding a missense‐to‐silent ratio of 0.64. The overall Ts/Tv ratio was 1.76, indicating a balanced mutation spectrum and overall data quality.

LD decay across the analyzed chromosomes was compact, with the maximum r ^2^ values ranging from 0.42 (chr09) to 0.51 (chrX) (Figure S2b; Table S7). On average, r ^2^ dropped to half its maximum within 300 bp. The results from the LD and HB analysis were congruent for male and female plants, with both sexes showing rapid LD decay and similar values for maximum r ^2^ (0.44 for females and 0.45 for males) (Figure S2c; Table S8). However, slight differences were noted in the HB analysis. Male plants exhibited larger HB spans (mean span of 0.51 kb and maximum span of 164.06 kb) compared to females (mean span of 0.43 kb and maximum span of 87.64 kb). Additionally, the proportion of SNPs in HBs was slightly higher in males (36.85%) compared to females (32.53%).

Population structure, genetic diversity, and evolutionary dynamics

The analysis of the population structure using the full catalog of SNPs revealed three main subpopulations (K = 3), supported by DAPC and model‐based clustering by fastStructure (Figure 3). Three main subpopulations were identified as K1 (clade I = 64 accessions), K2 (clade II = 64 accessions), and K3 (clade III = 17 accessions) (Table S9; Figure S3). DAPC and fastStructure analyses consistently supported K = 3 as the optimal number of clusters, with high agreement in cluster assignments (Figure 3a(ii),b; Table S9; Figure S4a,b). BIC comparisons revealed minimal differences between K = 2 (1359.038) and K = 3 (1359.258), with fastStructure identifying K = 3 as the optimal population structure (Figure S4b). Maximizing the marginal likelihood confirmed K = 3 as optimal, revealing subpopulation diversity and admixture patterns (Table S9; Figure 3a(iii) and Figure S5).

Of the initial 233,624 variants, 192,465 sites were utilized in evolutionary and distance analyses, comprising 52,894 conserved sites and 139,571 variable sites. The phylogenetic tree [Figure 3a(i)] aligned with DAPC [Figure 3a(ii)] and fastStructure clustering [Figure 3a(iii)], though with a different pattern, while clade I was identified as dominant subgroup in admixture analysis, and some of its accessions were positioned between clades II and III in the phylogenetic tree [Figure 3a(i)]. Genetic distance analysis further supported the observed population structure. Within‐group distances were lowest in clade I (0.0414 ± 0.002) and highest in clade III (0.1631 ± 0.0106), with clade II showing intermediate diversity (0.0864 ± 0.0035). Between‐group distances were smallest between clades I and II (0.0655), and highest between clades II and III (0.1513), followed by clades I and III (0.1232) (Figure 3d and Figure S6).

Consistent with clustering results, genetic diversity metrics (Table S10; Figure 3c) and geographic distribution patterns (Figure 3a(iv) and Figure S7) revealed distinct characteristics among clades. Clade II, found in warmer regions (BWh, BSh, and Cfa climate zones), exhibited the highest genetic diversity ([*θπ *]; 8.68 × 10^−4^) with MAF of 17.3% and heterozygosity of 29.8%, significantly exceeding values in clade I (*θπ *, 3.98 × 10^−4^; MAF = 7%; heterozygosity = 11.2%), which was present across all climate zones. In contrast, clade III, mainly distributed in colder regions (BSk and BWk climate zones), showed intermediate genetic diversity (*θπ *, 7.66 × 10^−4^; MAF = 14.5% and heterozygosity = 22.6%). Clade II exhibited the highest HB density (34.4 HB/Mb) compared to clades I (15.4 HB/Mb) and III (6.8 HB/Mb). However, clade II's HBs were shorter (mean span 0.34 kb) than those in clade I (1.38 kb) and clade III (0.81 kb) (Table S10). Despite this, the LD analysis revealed similar patterns across clades (Table S10; Figure 3c). Furthermore, kinship analysis revealed notable relationships within clade I, suggesting shared ancestral backgrounds within this subgroup (Figure S8).

Phenotypic insights from population structure

Phenotypic variation among the three genetic clades revealed distinct patterns across phenological, morphological, and phytochemical traits. Phenological analysis (Table S11) showed notable differences in reproductive phase initiation among clades (Figure S9). Specifically, traits like GVP (p ≤ 0.01), SFFI (p ≤ 0.05), and SFFP (p ≤ 0.05) exhibited significant variation among clades. Clade I (e.g., SFFI, mean = 70.87 ± 9.75 days, range = 53.2 days) and clade III (e.g., SFFI, mean = 65.95 ± 18.29 days, range = 51.02 days) displayed wide flowering time ranges, encompassing both early and late accessions. In contrast, clade II (e.g., SFFI, mean = 72.96 ± 4.88 days, range = 25.21 days) consistently showed a narrower range (Table S11). This substantial phenological diversity was paralleled by the observed morphological variations within and between clades (Table S12; Figure S10 and Figure 3e).

Notably, clades I and III displayed a higher phenotypic variation range for traits related to node and branching architecture (e.g., NNH, NMI, NTFIS, and NTMI all showing p ≤ 0.001) as well as growth and structural dimension (e.g., LMI [p ≤ 0.001], LSLS [p ≤ 0.001], RGR [p ≤ 0.001], and SDH [p ≤ 0.01]) (Table S12). Specifically, clade III exhibited the lowest mean for NNH (23.04 ± 2.99 nodes), significantly fewer nodes than clades I (28.55 ± 3.56 nodes) and II (29.46 ± 1.54 nodes), aligning with the observation of longer internode lengths. For plant HH, clade III had the lowest mean (133.09 ± 58.84 cm) but the widest range (165.05 cm), indicating it indeed encompasses both the shortest and tallest plants (p ≤ 0.05). Interestingly, despite these significant structural and architectural differences, biomass yield traits (e.g., FWS, DWS, TFW, and TDW) generally showed no significant difference (Table S12). However, while FWF did not show significant differences, DWF exhibited significant differences (p ≤ 0.05) among clades (Figure S10; Table S12).

Most notably, the phytochemical profiles of these landraces aligned strongly with the established population structure and geographical origins (Table S13). Analysis of chemotype classes (THC‐dominant, balanced, moderate CBD, and CBD‐dominant) distribution across clusters revealed a clear pattern, clade I exhibited a mixed chemotype profile, predominantly characterized by THC‐dominant accessions, showing a mean Total THC potential of 1.832 ± 1.079% and a lower mean total CBD potential of 1.172 ± 0.924%. Clade II comprised accessions showing either THC‐dominant or moderate‐CBD profiles, often displaying the highest mean concentrations for several cannabinoids, such as THCA (2.554 ± 2.01%) and Δ^9^‐THC (0.269 ± 0.159%). Significantly, clade III was almost exclusively composed of CBD‐dominant chemotypes (Table S13; Figure 3a(v) and Figure S11). This clade was characterized by significantly higher levels of CBD‐related cannabinoids (e.g., CBDA, CBD, and CBDV), with a markedly higher CBD:(THC + CBN) ratio (7.27 ± 7.124, p ≤ 0.001) compared to clades I and II. Conversely, clade III showed significantly lower levels of THC‐related cannabinoids (e.g., THCA, Δ^9^‐THC, CBN, THCV, CBC), with total THC potential averaging 0.659 ± 0.652% (p ≤ 0.01) (Table S13). The heatmap and PCA of cannabinoid profiles consistently supported chemotype‐based clustering, further reinforcing the observed distinctions among accessions (Figures S12 and S13).

GWAS in cannabis landraces

A GWAS using the BLINK model, which accounts for population structure and kinship, was conducted on 145 cannabis landrace accessions using 233,624 variants, identifying 91 loci associated with 40 phenological, morphological, and phytochemical traits under a stringent threshold (p‐value = 1.77 × 10^−7^, FDR < 0.05). Since some loci were associated with multiple traits across or within trait categories, the counts per category reflect unique significant markers.

Phenological traits

GWAS for phenological traits revealed 10 significant markers associated with six distinct phenological stages, each exerting a measurable phenotypic effect, with PVE values ranging from 5.3% (Cs3 for GVP) to 44.7% (Cs8 for SFFP) (Figure 4; Table 1; Figure S14a). Key loci included a novel auto‐flowering locus, designated AutoFlower3 (CsFTL3), which encompasses two closely linked markers, Cs1 (Chr8:7039513; PVE, 19.9% for GVP) and Cs5 (Chr8:7038258, PVE, 37% for SSFI, 39.3% for SFFP, 25.6% for SF10I, 24.6% for SF10P, and 20.9% for FT50I). Another important locus, designated FloweringTime_Maturity4 (CsFTL4), was identified on Chr7:7258159 (Cs9) with PVE values of 13.3% for SF10I and 30.6% for FT50I. Additionally, CircadianFloweringLocus1 (CsCFL1), located on Chr9:30838575 (CS2), showed strong associations with multiple traits (PVE, 33.2% for GVP, 41.4% for SF10I, 41.3% for SF10P, and 24.3% for FT50I). The allelic impact of CsCFL1 was particularly pronounced. For example, in the case of FT50I, individuals carrying the “CC” allele exhibited a significantly higher mean flowering time (121 days) compared to those with the “AA” allele (79 days), demonstrating clear allelic differentiation (Figure S15). Furthermore, several significant markers were located on the X chromosome (e.g., Cs4, Cs8, and Cs10), indicating its involvement in the regulation of phenological traits. Overall, these markers are of particular interest to cultivators aiming for faster crop turnover, as they are associated with shorter flowering or maturation times.

Morphological characteristics

A total of 52 significant genetic markers were identified across 20 morpho‐agronomic traits in the GWAS panel. The results for these traits are grouped into three (i) node and branching architecture, (ii) growth and structural dimensions, and (iii) biomass yield. (i) Node and branching architecture

In this subcategory, 11 significant genetic markers were associated with five traits, demonstrating phenotypic contributions with PVE values spanning from 3.7% (Cs18 associated with NTFIS) up to an impressive 87.1% (Cs16 associated with NMI) (Figure 5a; Table 2; Figure S14b). Key loci identified include NodeNumberRegulator1 (CsNNR1), represented by marker Cs11 (Chr3:90476582; PVE, 25.2% for NNH and 74.2% for NTFIS), NodeNumberRegulator2 (CsNNR2) (Cs19 on Chr1:12025015; PVE, 51.4% for NTMI), TopNodesDensity (CsTND), which includes two complementary loci, CsTND1 (Cs12 on Chr4:79659045; PVE, 23.6% for NNH) and CsTND2 (Cs16 on Chr4:79296939; PVE, 87.1% for NMI) as well as CsCFL1 (Cs2 on Chr9:30838575; PVE, 47% for NLS). The CsFTL3 locus, previously represented in flowering time traits by Cs1 and Cs5, is characterized by Cs5 (Chr8:7038258; PVE, 7.8% for NLS) and Cs13 (Chr8:7035841; PVE, 35% for NNH and 34.3% for NTMI), highlighting its extended influence on plant architecture. The allelic impact of Cs16 (CsTND2) was particularly pronounced, where individuals carrying the “CC” allele exhibited a significantly higher node count on the main inflorescence (12.23 nodes) compared to those with the “AA” allele (6.8 nodes) (Figure S16). Together, these loci coordinate both general and region‐specific control over phyllotactic development, including basal‐to‐apical node and branch patterning. (ii) Growth and structural dimension

GWAS analysis identified 25 significant markers associated with eight growth and structural dimension traits, with phenotypic effects spanning 3% (Cs40 for SDH) to 71.2% (Cs37 for RGR) (Figure 5b; Table 2; Figure S14c). Notable loci included RelativeGrowthRateLocus1 (CsRGR1), represented by marker Cs37 (Chr7:29874731; PVE, 71.2% for RGR), which demonstrated the strongest individual marker effect in this trait category, closely followed by CsTND1 (Cs12 on Chr4:79659045; PVE, 70% for LMI). Together, CsTND1 and CsTND2 underline a coordinated genetic control over inflorescence architecture through the modulation of both node number and internodal length. Additionally, CsCFL1 (Cs2 on Chr9:30838575) showed substantial pleiotropic effects with PVE values of 52.7% for LLLS, 54.4% for LSLS, 47.9% for HH, and 57.8% for SDH. Additionally, CsFTL3 locus was represented by Cs5 (Chr8:7038258; PVE, 19.6% for SDH) and Cs13 (Chr8:7035841; PVE, 28.5% for HH), further supporting its role in plant architecture beyond flowering regulation. Several markers associated with different growth and structural dimension traits were located in close proximity, suggesting shared genomic regions influencing these correlated characteristics. The allelic impact of these loci was particularly pronounced. For CsRGR1, individuals carrying the “CC” allele exhibited a significantly higher RGR (104.17 mg g^−1^ day^−1^) compared to those with the “TT” allele (69.7 mg g^−1^ day^−1^). Similarly, for CsCFL1, the “CC” allele consistently promoted increased structural dimensions across multiple traits (HH, 161.38 vs. 85.29 cm; LLLS, 53.72 vs. 28.4 cm; LSLS, 14.87 vs. 6.85 cm; and SDH, 13.6 vs. 7.67 mm), demonstrating its broad architectural influence compared to the “AA” allele. Overall, these markers coordinate fundamental aspects of plant morphology and growth dynamics, offering valuable targets for architectural improvement in breeding programs (Figure S17). (iii) Biomass yield

Biomass yield analysis identified 22 significant markers spanning seven yield components, with effects ranging from 3.3% (Cs44 for DWS) to 46.9% (Cs1 for DWS) (Figure 5c; Table 2; Figure S14d). The genetic architecture was largely shaped by the CsFTL3 region, where marker Cs1 (Chr8:7039513) exerted the strongest effect, particularly on stem dry weight (DWS; PVE = 46.9%) and also contributed to TFW (13.9%). Meanwhile, its closely linked counterpart, Cs5 (Chr8:7038258), displayed a partially complementary effect pattern, with major influence on stem fresh weight (FWS; PVE = 24.1%) alongside a notable contribution to TDW (11.4%). The CsFTL4 locus emerged as a major yield regulator through Cs9 (Chr7:7258159), controlling multiple biomass components with remarkable consistency, 32.4% PVE for DWS, 36.0% for TDW, 24.3% for TFW, and 13% for FWS. Additionally, CsCFL1 (Cs2 on Chr9:30838575; PVE, 20.4% for FWS, 30.5% for TFW, and 13.1% for TDW) extended its pleiotropic influence into yield determination through enhanced fresh biomass accumulation and total plant weight. The CsFTL4 locus demonstrated exceptional yield enhancement, where presence of the “AA” allele nearly doubled TDW (114.92 vs. 64 g) and increased TFW by 79% (276.47 vs. 154.41 g) compared to the deletion variant. The CsFTL3 region also showed similarly dramatic effects, with the “AA” allele of Cs1 producing a sixfold higher DWS (51.73 vs. 8.41 g) compared to the “CC” allele. Consistent with its broad developmental control, the CsCFL1 locus (Cs2) with the “CC” allele showed substantial 2.7‐fold gain in TDW (94.94 vs. 35.58 g), a 2.5‐fold increase in TFW (229.81 vs. 91.87 g), and a 3.3‐fold gain in FWS (112.3 vs. 34.16 g) compared to the heterozygous “AC” allele (Figure S18). It is noteworthy that while the DWL trait was included in the GWAS for biomass yield, no significant markers were identified for this specific trait. These yield‐controlling loci represent prime candidates for marker‐assisted selection targeting biomass optimization.

Cannabinoid profiles

GWAS of cannabinoid profiles identified 34 significant genetic markers associated with 14 phytochemical (cannabinoid) traits. The results for these traits are categorized by subgroup. (i) THC‐related

THC‐related traits analysis identified 20 significant markers associated with nine major biochemical traits, with PVE ranging from 3% (Cs69 for CBC) to 94% (Cs59 for total THC potential) (Figure 6a; Table 3; Figure S14e). The cannabinoid biosynthetic pathway was primarily regulated by THCALocus1 (CsTHCAL1), where the marker Cs59 (Chr5:29925798) exerted dominant effects across multiple THC‐related metabolites, explaining 94% of the variance in total THC potential, 91.4% in THCA, 89.9% in CBN, and 72.9% in Δ^9^‐THC. Beyond CsTHCAL1, other major loci shaped specific branches of the biosynthetic network, CBGALocus1 (CsCBGAL1) was captured by marker Cs67 (Chr3:13775640; PVE, 90.6% for CBGA and 7% for CBG), CBGRegulator1 (CsCBGR1) by Cs68 (Chr5:24473236; PVE, 91.3% for CBG), and THC:CBDRatioLocus1 (CsTCBR1) by Cs66 (Chr1:24225734; PVE, 79.2% for THC:CBD ratio), highlighting distinct control over precursor and end‐product partitioning. The THCV biosynthesis appeared to be regulated by Cs64 (Chr5:2418037; PVE, 57.9%), designated as THCVRegulator1 (CsTHCVR1). In contrast to these high‐effect loci, CBC biosynthesis was governed by a more complex multi‐locus architecture, involving 10 markers across chromosomes 1, 3, 4, 6, 9, and X, with the strongest individual contribution from Cs78 (ChrX:48786650; PVE, 27.6%), suggesting a polygenic basis for this trait.

The CsTHCAL1 locus demonstrated exceptional potency enhancement, where the “TT” allele produced nearly sixfold higher total THC potential (9.4% vs. 1.65% in “CC” allele) and similar magnifications in THCA content (9.84% vs. 1.66%). For the precursor pathway, CsCBGAL1 showed the “AA” allele producing substantially elevated CBGA levels (0.178% vs. 0.04% in “GG” genotypes), representing a 4.5‐fold increase in this critical biosynthetic intermediate. The CsTCBR1 exhibited striking chemotype differentiation. The “TT” allele shifts THC:CBD ratios from 7.32 in “GG” backgrounds to 59.73 in “TT,” indicating an eightfold change in cannabinoid profile direction (Figure S19). (ii) CBD‐related

CBD‐related traits through GWAS revealed 14 significant markers controlling five CBD‐related traits, with genetic contributions spanning 3.1% (Cs86 for CBD:(THC + CBN) to 82.1% (Cs1 (CsFTL3) for CBDV) (Figure 6b; Table 3; Figure S14f). The CBDALocus1 (CsCBDAL1), anchored by Cs79 (Chr3:82815664), emerged as a pivotal regulator with the highest individual effect on CBDA production (67.2% PVE) and substantial influence on CBD:(THC + CBN) ratios (16.5% PVE) and total CBD potential (32.6% PVE). The X‐linked CannabinoidRatioLocus1 (CsCRL1), represented by Cs91 (ChrX:65977078), demonstrated exceptional control over cannabinoid ratios with 62.4% PVE for CBD:(THC + CBN), establishing sex‐linked inheritance patterns in cannabinoid composition. Remarkably, the flowering‐time regulator CsFTL3 extended its pleiotropic influence into cannabinoid metabolism, where Cs1 (Chr8:7039513) achieved the strongest single marker effect in this category (82.1% PVE for CBDV). For CsCBDAL1, the “TT” allele dramatically enhanced CBDA production (2.69% vs. 1.11% w/w in “CC”) and substantially increased CBD:(THC + CBN) ratios (7.9 vs. 1.09 in “CC”), establishing its role as a major cannabinoid biosynthesis enhancer. Most remarkably, CsFTL3 revealed an unexpected connection between flowering regulation and cannabinoid metabolism, where the “CC” allele produced eightfold higher CBDV concentrations (0.0097% vs. 0.0012% w/w in “AA”), highlighting novel pleiotropic pathways linking developmental timing to secondary metabolite production (Figure S20).

HB and pleiotropic loci characterization

To understand the genetic basis of observed phenotypic variation, candidate genes within HB regions encompassing 91 unique genome‐wide significant markers associated with 40 morphological, phenological, and cannabinoid traits were investigated (Table S14). These HBs, defined by variants with an r^2^ ≥ 0.75, varied considerably in size, ranging from 5 bp (e.g., Cs7) to 37.6 kb (e.g., Cs45). Notably, several significant markers were found within the same haploblock, such as Cs1, Cs5, and Cs13 on chromosome 8 (HB size, 3 kb). This analysis identified 34 putative candidate genes within these regions, although six were classified as having uncharacterized or unknown functions, providing clear avenues for future investigation (Table S14). Notable cannabinoid‐related candidates include 2‐acylphloroglucinol 4‐prenyltransferase (Cs80, ChrX) and plastid‐lipid‐associated protein 13 (Cs81, Chr8) associated with CBD, and ketol‐acid reductoisomerase (Cs34, Chr8) linked to internode length. Major cannabinoid loci CsTHCAL1 (Cs59, Chr5) and CsCBDAL1 (Cs79, Chr3) lack annotated candidates in the reference genome, suggesting potential unannotated regulatory regions. For morphological traits, candidates include protein FATTY ACID EXPORT 1 and glycoside hydrolase (Cs11, Chr3), and pectinesterase (Cs90, Chr9). Additional candidates include disease resistance proteins and long non‐coding RNAs at various loci.

A key finding of this study is the extensive pleiotropy observed for 15 major loci, which are represented by 17 key markers (Figure 7). These loci exert broad regulatory control across multiple, seemingly distinct trait categories, highlighting genomic hotspots with far‐reaching effects. Specifically, the CsFTL3 complex, located within the 3 kb haploblock on chromosome 8, serves as a prime example of pleiotropy. The markers within this locus (Cs1, Cs5, and Cs13) significantly influenced multiple phenological stages (GVP, SFFI, SFFP, SF10I, SF10P, and FT50I), morphological traits (NNH, NTMI, HH, NLS, and SDH), and biomass traits (FWS, TDW, DWS, and TFW), with marker Cs1 also showing a strong association with CBDV concentration (82.1% PVE). Furthermore, marker Cs2 (CsCFL1) on chromosome 9 displayed the broadest pleiotropic impact, associated with 12 distinct traits across phenological, morphological, and biomass categories, and its alleles contribute to the differentiation of early, medium, and late‐flowering accessions. CsFTL4 (Cs9) on chromosome 7 was identified as a homolog of CsFTL3, influencing later phenological stages and several biomass components. Similarly, the loci CsNNR1 (Cs11) on chromosome 3 and its homolog CsNNR2 (Cs19) on chromosome 1 were found to regulate node number and branching patterns. The loci CsTND1 (Cs12) and CsTND2 (Cs16) on chromosome 4 demonstrated specialized pleiotropy, with CsTND1 controlling node number and internode length and CsTND2 specifically regulating the NMI, thereby affecting main inflorescence size and density. The extensive pleiotropy observed for several key markers emphasizes the complexity of trait regulation and offer valuable resources for molecular breeding strategies in cannabis.

DISCUSSION

This study provides a comprehensive genetic and phenotypic characterization of cannabis landrace accessions, comprising both females and males, offering valuable insights into the genetic architecture underlying key phenological, morphological, and phytochemical traits. Utilizing high‐density GBS data (Torkamaneh et al., 2021) and GWAS, we identified 91 unique genome‐wide significant markers associated with 40 traits, providing valuable insights into genetic architecture underlying trait variation in this crop. Our most significant finding is the discovery of 15 highly pleiotropic loci with exceptional effects on multiple trait categories, particularly CsFTL3 and CsCFL1, which provide molecular control points for coordinated developmental processes (flowering time, plant architecture, and biomass production) and, notably, extend into secondary metabolism (cannabinoid biosynthesis), revealing unexpected cross‐system genetic architecture. Notably, CsTHCAL1 and CsCBDAL1 exhibited exceptionally high phenotypic variance explained (94% for multiple THC‐related and 67.2% for CBDA, respectively), underscoring their relevance in regulating cannabinoid pathways. These discoveries provide unprecedented insights into the genetic mechanisms governing cannabis diversity and offer powerful tools for molecular breeding programs.

The substantial allelic diversity observed in our cannabis landrace panel is consistent with patterns reported in landrace populations of other crops (e.g., maize, rice, and wheat), where naturally maintained variation has provided a foundation for the identification of adaptive alleles, rather than implying direct equivalence in magnitude or impact (Dias et al., 2024; Hurni et al., 2015; Z.‐H. Ren et al., 2005; Rubiales & Niks, 2000; Z.‐J. Zhang, 1995). Our Iranian cannabis landraces provided the genetic diversity essential for identifying causal variants that remain hidden in narrow commercial germplasm (de Ronne & Torkamaneh, 2025; de Ronne et al., 2024). LD analysis revealed a notably rapid decay rate, much faster than in previous cannabis studies, which range from 3.9 (basal) to 6 kb (drug‐type) (G. Ren et al., 2021), 22.6–89 kb in commercial Canadian drug‐type cannabis cultivars (de Ronne et al., 2024), and 6.7 kb in feral germplasm from the United States (Aina et al., 2025). This aligns with outcrossing or self‐incompatible species such as tea (100 bp–8 kb; Lei et al., 2023) and maize (500 bp–6.3 kb; Dinesh et al., 2016; Pavan et al., 2020; Remington et al., 2001; Yan et al., 2009) and contrasts with self‐pollinated crops like soybean (>100 kb), tomato (1 Mb), rice (150 kb), and wheat (8 Mb) (J. Liu, He, et al., 2017; X. Liu, Geng, et al., 2017; H. Liu et al., 2019; Torkamaneh, Chalifour, et al., 2020; Viana et al., 2022). This rapid LD decay, shaped by dioecy, wind pollination, and limited selection pressure, enables high‐resolution mapping and demands dense marker coverage (Belzile & Torkamaneh, 2022; Hyten, 2022; Mohammadi et al., 2020). Interestingly, within the subpopulations, the highest diversity value in our panel was nearly identical to the lowest value observed among subpopulations of commercial drug‐type genotypes (8.44 × 10^−4^) (de Ronne et al., 2024). This contrasts sharply with broader populations, where *θπ

• values can reach as high as 3.87 × 10^−3^, reflecting much greater genetic diversity (G. Ren et al., 2021). Consistent with our previous study (Babaei et al., 2024), this population also showed substantial phenotypic variation across key traits. Heritability estimates (H ^2^ and SNP‐based h ^2^) revealed strong additive genetic control, supporting high selection potential (Falconer & Mackay, 1996; Naim‐Feil et al., 2021; Wray & Visscher, 2008). For example, H ^2^ for RGR reached 0.95 in females and 0.99 in males, with striking sex‐specific differences in morphological traits like DWF (0.33 vs. 0.87) and FWF (0.23 vs. 0.92) (Babaei et al., 2024). The observed deviations from normal distribution, indicative of rare alleles contributing to extreme phenotypes, further emphasize the utility of landraces in capturing unique genetic variations often absent in elite breeding lines (Belzile & Torkamaneh, 2022; Gupta et al., 2019; Ingvarsson & Street, 2011; Mohammadi et al., 2020; Torkamaneh, Chalifour, et al., 2020; Zhu et al., 2008).

Comparing our HD‐GBS results with de Ronne et al. (2024) on commercial Canadian drug‐type cannabis shows similar SNP counts, densities, and gene‐annotated SNP proportions. However, the commercial genotypes exhibited higher MAF (21.7%) and heterozygosity (25.5%) due to intensive selective breeding and use of diverse parental lines, which stabilize traits and maintain allele frequencies. In contrast, landraces, shaped by prolonged natural selection and limited human intervention, exhibit lower MAF reflecting broader allelic diversity and rare alleles at low frequencies (Kovalchuk et al., 2020; Linck & Battey, 2019). Studies by Sawler et al. (2015) and Soler et al. (2017) confirm hemp populations have higher heterozygosity (16% and 40.5%) than drug‐types (12.5% and 28.2%), likely due to wider genetic bases and hybridization levels. Lynch et al. (2016) reported lower heterozygosity in European hemp (22%) versus drug‐types (31%), while Gao et al. (2014) found higher heterozygosity in Chinese hemp (35.5%–37%) compared to European hemp (18.2%), highlighting geographic and breeding history effects. Overall, these patterns reflect natural selection and limited hybridization in landraces versus targeted breeding in commercial cannabis to stabilize specific traits (Hurgobin et al., 2021). Thus, the extensive phenotypic and genetic diversity, coupled with strong heritability, and rapid LD decay provides a robust foundation for GWAS and targeted breeding (Alqudah et al., 2020; Belzile & Torkamaneh, 2022; Ingvarsson & Street, 2011; Torkamaneh & Belzile, 2022).

The genetic architecture of Iranian cannabis landraces, characterized by three distinct genetic clades, reflects the ecological and historical complexity of the Iranian Plateau. Its fragmented topography and climatic heterogeneity have fostered local adaptation, creating recombination hotspots and divergent subpopulations (Akhani et al., 2010; Djamali et al., 2012; Gurjazkaite et al., 2018; Manafzadeh et al., 2017; Shumilovskikh et al., 2016). Iran's role as a Bronze Age trade hub, notably along the Silk Road, likely facilitating the dispersal and hybridization of cannabis (Abdullaev, 2022; Kovalchuk et al., 2020; X. Liu & Brancaccio, 2022; McPartland, 2018; Ndlangamandla et al., 2024; Warf, 2014).

While G. Ren et al. (2021), proposed a single East Asian origin based on Chinese hemp, earlier hypotheses favored Central Asia (de Candolle, 1867; De Candolle, 1883; McPartland, 2017; McPartland, 2018; McPartland & Small, 2020). Notably, Ren's study excluded Iranian accessions, potentially missing key diversification events. Supporting this, Balant et al. (2025) identified Iranian samples as a genetically distinct subgroup, favoring a geography‐based classification over use‐type model (G. Ren et al., 2021). While previous studies on Iranian cannabis have utilized low‐density molecular markers for morphological and phytochemical characterization (Shams et al., 2020) or a limited set of ∼24 K SNPs for population structure analysis (Soorni et al., 2017), our study provides a significantly higher‐resolution genomic landscape. Our genomic analysis revealed three Iranian subpopulations shaped by geography, environment, and traditional selection (G. Ren et al., 2021). Similar to Chinese cannabis (Chen et al., 2022), Iranian landraces show ecotype‐based adaptation. Clade I's intermediate phylogenetic position suggests either ancestral ties to clades II and III or historical hybridization (Pérez‐Escobar et al., 2021, 2022; Ramos‐Madrigal et al., 2019). Clade III stands out genetically, nearing thresholds used in the Genetic Species Concept (Bradley & Baker (2001), though such metrics remain debated (Zachos, 2016). Despite phenotypic variation—from ruderalis‐like too tall, late‐flowering types—clade III forms a cohesive genetic group, consistently CBD‐dominant. This highlights evolutionary processes like hybridization and convergent evolution that blur species boundaries (Steenwyk et al., 2023).

The observed phenotypic diversity in cannabis is profoundly influenced by environmental factors, including geographical latitude, climate, and the selective pressures exerted by both natural and human forces (Babaei & Ajdanian, 2020; Babaei et al., 2022; Babaei, Boissinot, et al., 2025; Hazekamp & Fischedick, 2012; Lata et al., 2023; G. Ren et al., 2021). As a short‐day plant, cannabis synchronizes its biological activities, including the critical transition to flowering, with environmental rhythms through its internal circadian clock. This fundamental biological process acts as an upstream regulator, influencing not only flowering time but also cascading effects on plant architecture, growth, and biomass accumulation (Gottlieb, 2019; Harmer, 2010; Steed et al., 2021; Webb, 2003; Webb et al., 2019). For instance, key loci identified in this study, such as CsFTL3 on chromosome 8 (encompassing markers Cs1, Cs5, and Cs13) and CsCFL1 (marker Cs2 on chromosome 9), exemplify these broad regulatory impacts. These markers exhibit significant pleiotropic effects, influencing multiple traits across phenological and morphological (e.g., node and branching architecture, growth dimensions), and cannabinoid profiles, consistent with their roles in photoperiod and circadian rhythm regulation. Such pleiotropy, where a single genomic region affects multiple seemingly unrelated traits, is a common feature in plant development and adaptation, observed in other short‐day plants like soybean and rice (e.g., FLOWERING LOCUS T (FT) homologs influencing vegetative and reproductive branching and yield‐related traits), and long‐day plants like Arabidopsis (e.g., FLOWERING LOCUS T [FT], FRIGIDA [FRI], FLOWERING LOCUS C [FLC] affecting branching, seed germination, and water use efficiency) (Cao et al., 2022; Hiraoka et al., 2013; Maple et al., 2024; Mohamedikbal et al., 2024; Weng et al., 2022). The identification of CsFTL3 and CsCFL1 in our study, with their broad pleiotropic effects, strongly suggests their central role in this intricate regulatory network in cannabis landraces. We also identified specific loci with major effects on phytochemical traits. Major cannabinoid loci CsTHCAL1 on chromosome 5 and CsCBDAL1 on chromosome 3 show the strongest associations with cannabinoid traits yet lack clearly annotated candidate genes in the reference genome, consistent with observations that cannabinoid synthase regions are embedded within structurally variable genomic regions enriched for transposable elements and pseudogenized paralogs (Lynch et al., 2025). In contrast, CBD‐associated candidates at other loci include 2‐acylphloroglucinol 4‐prenyltransferase (Cs80) and plastid‐lipid‐associated protein (Cs81), which represent plausible functional candidates based on their annotations.

The genetic architecture of THC and CBD extends beyond canonical synthase genes, with contributing markers distributed across various chromosomes, and total cannabinoid content QTL often residing outside the synthase cluster (de Ronne & Torkamaneh, 2025; Hurgobin et al., 2021). Our phenotypic analysis, which categorized cannabinoids into THC‐related and CBD‐related groups based on strong correlations and PCA, aligns with their distinct biosynthetic pathways (Govindarajan et al., 2023; Gülck & Møller, 2020; Hurgobin et al., 2021; Ingvardsen & Brinch‐Pedersen, 2023; Welling et al., 2020).

Classical genetic studies described a Mendelian B locus model where codominant alleles at a single locus (Chr7, THCAS/CBDAS) determine chemotype. We detected only minor Chr7 signals, with our primary loci on Chr5, Chr3, and Chr1 being fixed in modern cultivars but polymorphic in landraces (de Ronne & Torkamaneh, 2025). For varin cannabinoids, we found no overlap with the BKR locus (Chr04; reported as chromosome 9 by Welling et al., 2020, based on Purple Kush assembly) or ALT haplotype (Chr7; Lynch et al., 2025). Note that Welling et al. used Purple Kush and Finola assemblies with different chromosome numbering compared to the cs10 v2 reference used here. Instead, our Chr5 (THCV) and Chr8 (CBDV) loci reveal landrace‐specific genetic architectures potentially shaped by geographic isolation and distinct selection pressures in Iranian germplasm.

CONCLUSION

In summary, this study provides a robust genetic framework for understanding trait variation in Cannabis landraces, identifying key markers linked to important traits. These findings offer valuable tools for marker‐assisted and gene‐editing approaches, supporting the development of improved cultivars. Preserving these landraces is crucial for maintaining genetic diversity and unlocking novel traits for future breeding innovations.

AUTHOR CONTRIBUTIONS

Mehdi Babaei: Conceptualization; data curation; formal analysis; investigation; methodology; project administration; resources; software; visualization; writing—original draft; writing—review and editing. Davoud Torkamaneh: Conceptualization; funding acquisition; project administration; resources; supervision; writing—review and editing.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

Supporting information