In this study, we used the hemp variety Diku (DK), a female plant crossing Purple Kush with Dinamed Autoflowering CBD. This variety has a short growth cycle and is a variety with high CBD and low THC content. It is an excellent reference variety for cannabis breeding. The plants were cultivated in the experimental field of the Institute of Chinese Materia Medica of the Chinese Academy of Chinese Medical Sciences, China. Flowers, bracts, leaves, stems, seeds, and roots of hemp were collected to analyze the expression patterns of HDAC genes. Hemp seeds were cultured in MS medium for 3weeks. After the seedlings grew, they were transferred to an MS medium containing 0-μm (DMSO matrix), 0.25-μm, or 1.25-μm Trichostatin A (TSA) and cultured for 1week to analyze the effect of TSA on the HDAC gene of hemp. All sample materials collected in this study were frozen in liquid nitrogen immediately and then reserved at −80°C for further analysis.

The genomic data (GCA_900626175.1) of hemp variety female CRBRX, with a high-CBDA cultivar (15% CBDA and 0.3% THCA) used in this study, were downloaded from the NCBI¹ database. Arabidopsis thaliana and soybean HDAC genetic data were obtained from the PlantTFDB database.² We measured the transcriptome of five different tissues and organs of Diku. The transcriptome data of its flowers, bracts, stems, leaves, and seeds were available at NCBI (NCBI accession number: flowers: SAMN16122886~SAMN16122888; bracts: SAMN16122880~SAMN16122882; stems: SAMN16122883~SAMN16122885; leaves: SAMN16122889~SAMN1612281; and seeds: SAMN20474456~SAMN20474458).

We downloaded the full gene set and annotation files of C. sativa (GCA_900626175.1), Glycine max (GCF_000004515), and A. thaliana (GCA_000005425.2) from the NCBI Web site. We extracted the CDS sequence and protein sequence of hemp using the TBtools software (Chen et al., 2020). Using the protein sequence of AtHDAC and GmHDAC as query in BLASTp (score value of ≥100, e-value ≤ e⁻¹⁰), we ran the extracted hemp data against the AtHDAC sequence to identify potential HDACs in hemp. Then, we initially established identity-reliable domains by manually searching the NCBI reserved domain database.³ Finally, 14 HDAC family members were obtained. The basic characteristics, such as the protein molecular weight, isoelectric point, and amino acid number, were analyzed on the ExPASy Web site.⁴

CsHDAC, GmHDAC, and AtHDAC domain amino acid sequences were aligned using MEGA 7.0 software. The neighbor-joining (NJ) adjacency method was used to build a phylogenetic tree against the A. thaliana and G. max. The bootstrap test was repeated 1,000 times, and the BLASTp cutoff was set as value of E<1×10⁻¹⁰. The structure of the CsHDACs and their localization on the chromosomes were visualized using TBtools (Chen et al., 2020). The conserved motifs of CsHDAC proteins were carried out using MEME.⁵ The parameters are set as follows: The maximum number of themes is 20, and the best theme width is 10–200. Conserved structural domains were analyzed using NCBI-CDD.⁶ Subcellular localization was predicted using the wolf PSORT Web site.⁷ Sequences 2000bp upstream of the start codon of each CsHDAC gene were extracted, and the distribution of the cis-elements in the promoter regions was analyzed using PlantCARE software.⁸ We downloaded the chromosomal locations of 14 HDAC genes in the C. sativa genome from the NCBI.⁹ The HDAC gene locations were then represented using TBtools (Chen et al., 2020). Multiple collinear scanning toolkits (MCScanX) were used to detect gene duplication events (Wang et al., 2012).

Differential expression analysis of CsHDACs and all isoforms of CsHDACs were performed, and heat maps were created using TBtools software. In short, SpliceMap was used to detect the junctions of Illumina short reads obtained from samples, and IDP was then used to predict and detect isoforms by integrating both Illumina short reads and SMRT long reads. The results of IDP show the sites where isoforms are located on the genome. Integrative Genomics Viewer was used to visualize the isoform structure that may exist for each HDAC gene family member.¹⁰ Specific steps were referred to in our previous research (Gao et al., 2019; Xu et al., 2020). We took the original annotated genes as reference isoforms (ref isoforms). The other isoforms were divided into intron retained (IR), an exon skipping (ES), and alternative 3′ splice site (AA). Isoforms that did not belong to the aforementioned three types were defined as “others.”

According to the manufacturer’s instructions, we isolated total RNA from samples (Tiangen Biotech, Beijing, China). Using a Fast Quant RT Kit, the first cDNA strand was synthesized (Tiangen Biotech). The specific primers for qRT-PCR were designed using NCBI-Primer blast (Table 1). CFX96^™ real-time system (Bio-Rad, Hercules, CA, United States) was used for qRT-PCR. Program: 95°C for 2min, followed by 40cycles of 95°C for 15s, 56°C for 35s, and 72°C for 15s. A melting curve was generated and analyzed. The reference gene in this study was EF1-α (Nicot et al., 2005). Each sample was in triplicate experiments.

After adding liquid nitrogen in a mortar, leaves of hemp seedlings were crushed into a uniform powder. An ultrasonic ice water bath was used to dissolve approximately 0.1-g sample powder in 10.0-ml 95% aqueous methanol for 30min. Approximately 12,000rpm was applied for 10min at 4°C to centrifuge the solution; then, the supernatant was filtered by a 0.22-μm hydrophilic organic nylon microporous membrane.

The Agilent UPLC 1290II–G6400 triple quadrupole mass spectrometer (QQQ; Agilent Technologies, Santa Clara, CA, United States) was used to determine the relative quality of synthetic precursors. We used the peak area to express the relative mass under each condition.

A bidirectional solvent delivery system, an autosampler, and a column compartment were used in the UPLC. A column (2.1 ^* 100mm, 1.8μm) of C18 was used for chromatographic separations at a temperature of 40°C. Mobile phase A contained water containing 0.1% formic acid; phase B is 100% methanol. Elution was proceeded linearly at a flow rate of 0.3ml/min. In gradients, these settings were used as: 5% B→70%B (0–2min); 70% B→100% B (2–10min); 100% B (10–13min); 100% B→5% B(13–14min); and 5% B (14–15min). Autosampler was set to 4°C and 3-μl sample volume was injected. The experiment was performed in triplicate. The MRM diagrams of OA, GPP, and their standard products are shown in Supplementary Figure S1.

All the data were analyzed by the Prism 8 Statistics programs, and the means were compared by the least significant difference test at the 0.05 and 0.01 level of significance.

In our study, 14 HDAC genes were obtained, and RPD3/HDA1 subfamily contained 10 members, named CsHDA1–CsHDA10; two members of HD2 subfamily named CsHDT1–CsHDT2; and two members of SIR2 subfamily named CsSRT1–CsSRT2. Their sequences ranged in length from 432 to 1953 bp. The molecular weight of their translated protein ranged from 21702.17 to 56624.95Da, and they had between 203 and 504 amino acids. The other characteristics, such as theoretical IP and subcellular localization, were shown in Table 2.

To demonstrate the evolutionary relationships of CsHDAC in plants, we constructed a phylogenetic tree using amino acid sequences from model plant Arabidopsis, model food crop soybean, and hemp. The CsHDAC family was categorized into three subfamilies based on the difference in sequence features and the number of HDAC domains, SIR2, HD2, and RPD3/HDA1, with 2, 2, and 10 members, respectively (Zhang et al., 2020). Clades I–IV represented the RPD3/HDA1 subfamily, which is a type of Zn⁺-dependent histone deacetylase. Clade V represented the HD2 subfamily, which are plant-specific HDACs. Clade VI represented the SIR2 subfamily, a type of nicotinamide adenine dinucleotide (NAD⁺)-dependent HDACs. The 10 members of RPD3/HDA1 subfamilies were further divided into four subgroups based on the similarities of the HDAC domain amino acid sequences, with Subgroup I including CsHDA1 and CsHDA7, Subgroup II including CsHDA6 and CsHDA10, Subgroup III including CsHDA2, and Subgroup IV including CsHDA3-5 and CsHDA8-9 (Figure 2).

The conserved motifs of proteins may be related to their transcriptional activity, protein interactions, and nuclear localization (Ikeda and Ohme-Takagi, 2009; Causier et al., 2012). Thus, we studied the conserved motifs of all CsHDAC proteins combined with the evolutionary relationship of hemp (Figure 3A), and 20 conserved motifs were predicted (Figure 3B). The RPD3/HDA1 subfamily included numerous conserved motifs, while the SIR2 and HD2 subfamilies contained relatively fewer conserved motifs. According to the analysis of the conserved motifs, motifs 1 and 3 partly represented the distribution of the HDAC-conserved domain and were shared by all 10 members of the RPD3/HDA1 subfamily. Motifs 10, 13, and 14 were highly conserved and were included only in the Class II subgroup. The conserved motifs 16 and 17 were identified as a characteristic NPL domain, a structure that is unique to the HD2 subfamily. Motif 18 corresponded to the SIR2 conserved domain, which only existed in the SIR2 subfamily (Figure 3E). Multiple sequence comparison of important conserved motifs is shown in Supplementary Figure S2. Other conserved motifs were also found in CsHDACs, while the action mechanism of these motifs has not been analyzed. Overall, the conserved motif composition and gene structure within each family of CsHDAC members had a high similarity, and the results of the phylogenetic analysis advocated the validity of the population classification.

We investigated the exons and introns from their number and distribution to study the structural composition of the CsHDAC genes (Figure 3C). Two members of the SIR2 subfamily had eight exons, and two members of the HD2 subfamily contained 10 exons. The coding sequences of the entire RPD3/HDA1 subfamily were interrupted by introns, with the number of exons in the range of 1–17. The CsHDAC genes had three main domains, which were specifically distributed into three subfamilies (Figure 3B). The SIR2, RPD3/HDA1, and HD2 subfamilies consisted of the SIR2, HDAC, and NPL domains, respectively. Members of the HD2 subfamily of C. sativa had conserved MEFWG pentapeptide at the N-terminus and possibly had a zinc finger at the C-terminus. Generally, the number of exons within each subgroup was similar. The domain of each subfamily was conserved, although the structures of members within a particular subfamily were different. In addition, we analyzed the co-expression analysis of CsHDACs and key genes in the cannabinoid synthesis pathway (Figure 3D; Supplementary Table S1). We found that CsHDA1, CsHDA5, CsHDA7, and CsSRT1 have obvious correlations with key genes.

Chromosome mapping of the CsHDAC genes was performed using the latest C. sativa genome database. A total of 14 CsHDAC genes were unevenly distributed on six chromosomes (Figure 4A). Among them, chromosomes 1, 2, 3, 4, 6, and 8 harbored 5, 1, 4, 1, 3, and 1 genes, respectively. While chromosomes 5, 7, 9, and 10 did not contain any HDAC gene. We further studied the homologous genes of CsHDAC within and between species. No homologous genes are found on the chromosomes of C. sativa. The results of collinearity analysis between C. sativa and A. thaliana showed three orthologous genes (Figure 4B).

Cis-acting element analysis was performed on the CsHDAC gene family members, and the results were visualized (Figure 5). The HDAC gene family contained 13 cis-type elements, and each member of the family contained multiple response elements. Among them, HDAC9 contained many light-responsive elements, suggesting that HDAC9 may be highly sensitive to light stimuli. HDA10 contained the largest number of cis-elements, with six cis-elements, whereas HDA8 had the least number of cis-elements, with only two cis-elements.

Alternative splicing, which causes polymorphisms in the function and structure of proteins, is extensive in eukaryotic organisms (Chaudhary et al., 2019). Therefore, we performed hypomorphic detection and prediction on the basis of data obtained via second-generation sequencing and third-generation sequencing and identified 21 splicing events in eight HDAC genes (Figure 6A). Compared with the eight selected reference isoforms, splicing events included seven intron retentions (IR), three alternative 3′ splice site (A3SS; AA), one exon skipping (ES), and two other isoforms (“others”). One gene may generate different alternative splicing events. For instance, NC_044372.1: 65918457–65940183.9 (CsHDA5) contained two cutting types: IR and “others.” The expression of different transcripts of the same gene differed. In other words, some transcripts dominated the expression of the gene in all tissues. For instance, rna-XM_030622750.1 (CsSRT1) had higher expression in all tissues than the other transcripts. Additionally, different transcripts may be enriched in various tissues. For instance, the rna-XM_030643555.1 (CsSRT2) had the highest expression in roots, while another splice variant, rna-XM_030643556.1 (CsSRT2), had the highest expression in leaves.

To explore the expression patterns of the CsHDAC gene family in different organs, we drew heat maps based on the FPKM calculated from RNA-seq data of seeds, flowers, stems, leaves, and bracts of DK (Figure 6B). CsHDA1, 5, 7, and 13 were commonly expressed in all tissues. CsHDA1 and CsHDT2 were highly expressed in all tissues, especially in seeds. The expression of CsHDA6 was also relatively higher in seeds. CsHDA10 was highly expressed in the stem, while CsHDA4 was enriched in leaves. Genes in flowers and bracts had a similar expression pattern, and the majority was expressed at lower levels than in other tissues. Expression patterns of CsHDAC genes in different tissues can be used as the foundation for identifying functional genes in C. sativa. According to previously published research, MaDHAC1, AtHDA6, and AtHD2 were associated with plant secondary metabolism (Han et al., 2016). Therefore, we selected candidate homologous CsHDAC members by constructing NJ trees and performing motif analysis between these genes and CsHDACs (Supplementary Figure S3) and verified them via qRT-PCR analysis. CsHDA1 and CsHDA5 were evenly enriched in all tissues. CsHDA7, CsHDA8, and CsSRT1 were less strongly expressed in the bracts (Figure 6C).

TSA can broadly inhibit the HDAC subfamily (Li et al., 2014; Su et al., 2015). To study whether HDAC is related to the synthesis of cannabinoids, we compared the changes in gene expression in hemp seedlings before and after TSA treatment. HDAC gene expression decreased significantly with increasing concentration of TSA. TSA could effectively inhibit the expression of all the CsHDAC genes in the exception of CsHDA10 (Figure 7A). In this study, two kinds of synthetic precursors were unambiguously identified as OA and GPP, respectively, by comparing the retention times, adductions, and productions with those of authentic standards. We also tested the expression of representative genes in the cannabinoid metabolic pathway and the changes in the synthetic precursors OA and GPP. The expression of genes in the fatty acid synthesis pathway, OAC, OLS, and AAE, was generally downregulated (Figure 7C). However, the expression of PP1 (HMB-PP reductase1), PP2 (HMB-PP reductase2), and geranyl diphosphate synthase (GPPS) in the MEP pathway increased (Figure 7B). Correspondingly, the content of synthetic precursor OA was reduced, while GPP accumulated (Figure 7D). The above results indicated that TSA inhibited the expression of CsHDAC genes and altered the contents of representative genes and synthetic precursors in the cannabinoid synthesis pathway.