Cannabis sativa L. is both an ecologically and economically important crop that has been used to produce textiles, paper, building materials, food, and medicine. The first genomic draft of cultivated variety in C. sativa was published in 2011 [45], following genomic sequencing of wild-type varieties [46]. With the publication of genome sequencing, genome-wide identification of key gene families in this species has been carried out in recent years [47,48,49,50,51,52,53,54]. WAK gene family members trigger cell elongation during plant development by binding to cross-linked pectin in the cell wall and provide a response to pathogen or environmental stress factors by binding to OGs (oligo-galacturonides) formed by pathogens or mechanical damage [55]. Thus, in this study, a comprehensive genome-wide analysis of wall-associated kinases in C. sativa and their expression levels in different tissues and flower development was performed. A total of 23 CsWAK and 30 CsWAKL protein family members were identified. Previous studies showed that the number of WAK/WAKL genes differed according to species, regardless of genome size—for example, 26 in Arabidopsis [22], 29 in cotton [41] and 29 in tomato [8], 125 in rice [56], 91 in barley [41], 68 in rose [10], 115 in Brachypodium [7], and more than 100 in maize [14].

To clarify the evolutionary relationship among the members of WAK gene family, a phylogenetic tree was constructed using the WAK/WAKL protein sequences from Arabidopsis thaliana, Oryza sativa, Zea mays, and Gossypium hirsutum. We compiled a total of eight pathogen-defense-response-related WAK proteins including OsiWAK112, OsWAK91, OsWAK92, OsWAK14 from Oryza sativa, GhWAK5, GhWAK9, GhWAK26 from Gossypium hirsutum, and ZmWAK-RLK1 from Zea mays. Additionally, WAK proteins differentially expressed in tissues (GhWAK3, GhWAK4, and OsiWAK1) and defense response against excess copper (OsWAK11) were included. Phylogenetic analysis has shown that all CsWAKs/CsWAKLs were classified into four main groups (Figure 2). Arabidopsis WAKs were grouped together into Subgroup 1. The fact that all Arabidopsis WAKs were in the same subgroup and that other OsWAKs and GhWAKs were generally clustered close to each other was consistent with the findings of previous studies in cotton [40], in Junglas regia [12], in apple [57], and in barley [11]. However, the exception is that OsiWAK11 was not close to the others. The species–specific grouping of WAKs by phylogenetic analyses might be explained by the independent evolution of WAKs in each species. Moreover, the phylogenetic groups indicated that WAKs/WAKs have a high diversity of across species.

According to chromosomal localization, 23 CSWAK/CsWAKL genes were unevenly distributed on six of the cannabis chromosomes, but most of them (12 CSWAK/CsWAKL) were located on Chr 7 (Figure 4). Furthermore, two tandem duplications occurred on Chr 4 and Chr 7, including four genes (CsWAKL3-CsWAKL4 and CsWAKL7-CsWAK19). Additionally, the flanking region of CsWAK1 and CsWAKL24 genes contain reverse transcriptase-like and zinc-binding in reverse transcriptase sequences, respectively. This finding is consistent with a previous study [58] concerning Arabidopsis WAK1, WAK2, WAK4, WAK5, WAKL11, and WAKL16 genes carrying Hopscotch (Copia like retrotransposon). Kumar et al. [59] underlined that retrotransposons may carry all or part of the gene nearby, causing them to expand and spread throughout the genome. Similarly, the increase in the number of members of the WAK gene family by tandem or segmental duplications and retrotransposon and the acquisition of new functions in duplicated genes have been demonstrated in various studies [7,8,60].

To gain a better understanding of the CsWAK/CsWAKL gene family’s evolutionary relationship, we calculated the Ka/Ks ratios of 18 linked gene pairs (Table 2). We found that all but one of the Ka/Ks values was less than 1. This indicated that the pattern of evolution was mainly by purifying selection except for one with positive selection, hence, critical in sustaining the number of gene family members.

Classification of the members of a gene family is an essential step in exploring their structural and functional characteristics as well as evolutionary relationships [51]. In this respect, the protein motifs and their distribution were investigated to reveal their functional diversity. A total of 10 motifs, ranging from 17 to 50 amino acids length were found in 53 CsWAK/CsWAKLs. Half of the 10 predicted motifs were associated with WAK protein domains. Four motifs (motifs 1, 2, 5, and 6) corresponded to pkinase domain and one (motif 7) to the EGF domain. CsWAKs/CsWAKLs belonging to the same group/subgroup in the phylogenetic tree (Figure 6a) showed highly similar motif distribution. Also, the fact that each group/subgroup of CsWAK/CsWAKLs had different motifs from each other may be attributed to the diversity of their functional roles.

The presence of cis-elements in the promoter sequences suggests regulation of the transcription of these genes that may be responsible for effects on development and response to hormones and stress. Therefore, cis-acting elements were analyzed in this study. A total of 23 CsWAK/CsWAKL gene promoters contained light-response elements, especially CsWAKL11 and CsWAKL9 which had the most; additionally, 21 CsWAK/CsWAKL had Box 4 cis-elements. Eleven CsWAK/CsWAKL promoters had the cis-element associated with endosperm-specific negative expression. CsWAK3 contained only the cis-elements (as-1) for root-specific expression. The seed-specific regulation (RY-element) element was found only in the promoter of the CsWAKL11 gene. Therefore, the development-related cis-element analysis indicates tissue-specific expression of certain CsWAK/CsWAKL genes. The promoters of CsWAK7, CsWAK9, CsWAKL12, and CsWAKL10 genes had circadian-rhythm-responsive elements. This finding was consistent with previous studies pointing to the regulation of some members of receptor-like kinases by the circadian clock [61]. Considering hormone response, abscisic acid, Methyl jasmonic acid, auxin, gibberellin, salicylic acid, and ethylene hormone-responsive cis-elements were found in 13, 15, 12, 1, 19, and 21 CsWAK/CsWAKL genes, respectively. This means that the expression of these genes might be induced by treatment with the respective hormones. In addition, these CsWAKs/CsWAKLs that have hormone-related elements may play an important role in cell elongation and stress response, as hormones have important effects on plant growth, response to stress, and adaptation to adverse environmental conditions. Likewise, a previous study showed that the expression level of 13 and 10 GhWAK genes containing the gibberellin- and the auxin-responsiveness element, respectively, were increased after gibberellin and auxin treatment [41].

In this study, many genes closely located on the same chromosome and in the same subgroup of the phylogenetic tree constructed using 10 conserved motifs had similar exon–intron distribution, for example, CsWAKL2, CsWAKL3, and CsWAKL4 in Chr 4; CsWAK5, CsWAK6, CsWAKL9, and CsWAKL10 on Chr 7; and CsWAKL13-CsWAKL14 in Chr 9. These genes differed only in the length of introns and upstream regions. Similarly, Dou et al. (2021) [41] observed that genes with high sequence similarity in cotton are located close to each other on the same chromosome and reported that this would prevent loss of function during evolution.

Gene expression pattern analysis could predict systematic gene function [62]. Known WAK/WAKL genes in other plant species could give us an effective way to select candidate WAK/WAKL genes in cannabis [12,41]. Cannabis is one of the most important commercial crops and has a long history during the development of human civilization [63]. Cannabis flower tissue, especially trichomes, has a central place in industry and academic research [64]. We used homolog analysis to select four CsWAK/CsWAKL genes (predicted in flower development) and one gene (predicted in pathogen response). Interestingly, these four CsWAK/CsWAKL genes (CsWAK1, CsWAK4, CsWAK7, and CsWAKL7) were mainly expressed in leaf tissue, not flowers. This could give us new thoughts about the function of these genes. Therefore, our work suggests that further research to explore the function of these gene products in leaves is warranted.

Protein–protein interaction analysis was conducted to predict the potential function of CsWAK/CsWAKL proteins. WAK1 was seen to be the most interactive protein (Figure 9). This also confirms the strong interaction between Gly-rich extracellular protein AtGRP-3 [65] and a cytoplasmic type 2C protein phosphatase KAPP [66] with WAK1 to generate a 500-kD protein complex (Figure 9). GRP-3 expression is upregulated in response to ABA, salicylic acid, and ethylene but downregulated in response to desiccation. Furthermore, KAPP plays a role as a negative regulator of the CLV1 signaling in plant development. In addition, WAK1 protein is required during the plant’s response to pathogen infection and in plant defense against heavy-metal toxicity [5]. WAK1 had four CsWAK/CsWAKL proteins (Supplementary Table S2) corresponding to it, suggesting that four CsWAK/CsWAKLs may share in the fundamental roles of WAK1 and its complex with GRP-3 and KAPP. According to Wagner and Kohorn (2001) [23], WAK2 plays an important role in cell elongation. WAK2 corresponded to 24 CsWAK/CsWAKL proteins. Thus, these CsWAKs/CsWAKLs participate in cell elongation. Furthermore, interactions between RFO1(AtWAKL22) corresponding to CsWAKL12, CsWAKL18, and RLP3 could also be associated with their similarity in conferring disease resistance to a broad spectrum of Fusarium races, indicating that CsWAKL12 and CsWAKL18 have roles in response to abiotic stress. AT1G16260 matched with CsWAK22 and CsWAK9 functions as a signaling receptor of extracellular matrix components, and it interacted only with AT4G01860 and AT5G58412.

miRNAs are small, non-coding RNA molecules. They regulate gene expression at transcriptional and translational levels and complex biological processes [67]. Das et al. (2015) mined miRNAs in silico for the first time in C. sativa, followed by Hasan (2015) [44]. A total of 18 conserved miRNAs belonging to nine families were found by homology-based computational approach in Transcript Sequence Assemblies (TSA) and ESTs data in C. sativa [43,44]. The number of members in each of the miRNA family varied from one to seven. A total of 18 csa-miRNAs targeted 80 genes in the Arabidopsis genome, including transcription factors, transporters, kinase, and other enzymes as well as signalling and stress-responsive proteins. In the present study, csa-miR5021 and csa-miR172 were found to be related with regulation of three CsWAK genes (CsWAK5, CsWAK6, and CsWAK9) and four CsWAKL genes (CsWAKL4, CsWAKL5, CsWAKL13, and CsWAKL14) via the process of cleavage or translation (Supplementary Table S2). csa-miR172 targeted CsWAK5 transcript. It has been reported that the csa-miR172 family has seven members and csa-miR172-mediated gene regulation in C. sativa plays dominant roles. Moreover, csa-miR172 were found to act on transcription factors that regulate flowering time in C. sativa (Das et al. 2015) and also promoted flowering and adult leaf traits in Arabidopsis [68]. In addition, Das et. al. (2015) [43] concluded that the genes targeted by csa-miR5021 and csa-miR6034 were involved in ribosomal proteins, calcium-binding EF hand family protein, global transcription factor, disease resistance protein, leucine-rich receptor-like protein kinase, and so on. Thus, these csa-miRNAs participate in growth, development, and stress responses in C. sativa. Further studies are needed to elucidate CsWAK/CsWAKLs’ roles in plant growth and stress response through miRNA-mediated silencing of CsWAK/CsWAKLs.