Key Cannabis Salt-Responsive Genes and Pathways Revealed by Comparative Transcriptome and Physiological Analyses of Contrasting Varieties

: For the dissection and identiﬁcation of the molecular response mechanisms to salt stress in cannabis, an experiment was conducted surveying the diversity of physiological characteristics. RNA-seq proﬁling was carried out to identify differential expression genes and pathway which respond to salt stress in different cannabis materials. The result of physiological diversity analyses showed that it is more sensitive to proline contents in K94 than in W20; 6 h was needed to reach the maximum in K94, compared to 12 h in W20. For proﬁling 0–72 h after treatment, a total of 10,149 differentially expressed genes were identiﬁed, and 249 genes exhibited signiﬁcantly diverse expression levels in K94, which were clustered in plant hormone signal transduction and the MAPK signaling pathway. A total of 371 genes showed signiﬁcant diversity expression variations in W20, which were clustered in the phenylpropanoid biosynthesis and plant hormone signal transduction pathway. The pathway enrichment by genes which were identiﬁed in K94 and W20 showed a similar trend to those clustered in plant hormone signal transduction pathways and MAPK signaling. Otherwise, there were 85 genes which identiﬁed overlaps between the two materials, indicating that these may be underlying genes related to salt stress in cannabis. The 86.67% agreement of the RNA-seq and qRT-PCR indicated the accuracy and reliability of the RNA-seq technique. Additionally, the result of physiological diversity was consistent with the predicted RNA-seq-based ﬁndings. This research may offer new insights into the molecular networks mediating cannabis to respond to salt stress.


Introduction
Cannabis is one of the oldest cultivated crop, and has been grown for over ten thousand years; it is considered to have originated from the central Asia area, and then transmitted and planted around the world [1,2]. However, the cultivated area gradually decreased or disappeared in the mid-20th century [3], especially due to restrictions by the International Narcotics Convention [4,5]. Cannabis is famous as an source of tetrahydrocannabinol (∆9-THC), which is a primary cause of hallucinogenic effects [6]. Cannabis may be the most recognizable and controversial species around the world, because of its characterization as an illicit drug [4]. Another reason for the decrease in cultivation area is that the laws and rules related with cannabis do not distinguish the different forms of cannabis, resulting in restrictions to cultivating and scientific research, and not considering industrial cannabis which, without hallucinogenic material, could as an source of fiber, fuel and nutrition [7]. controls. Then, the samples were mixed at all stages of treatment, and the samples without salt stress were all subjected to transcriptome sequencing.

Physiological Analysis
The samples were collected at times of 6 h, 18 h, 24 h and 2 days after salt treatment, and the samples which had undergone treatment were the controls. The samples were quickly stored in liquid nitrogen. Proline contents and chlorophyll contents were determined with an ultraviolet spectrophotometer. Samples were ground in liquid nitrogen; then, 0.5 g of material was heated in a boiling water bath for 10 min after adding 5 mL sulfosalicylic acid solution. Subsequently, the 2 mL supernatant was extracted with the chromogenic reagent, which comprised ninhydrin and glacial acetic acid and bathe in boiling water for 40 min. The absorbance was determined at a wavelength of 520 nm, and finally, the proline content was determined from the absorbance through the standard curve [37].
The materials were mixed with an 80% acetone solution after grinding the sample in liquid nitrogen, and then shaking a few times over 48 h of lightless cultivation. The chlorophyll content was determined from the absorbance at a wavelength of 652 nm [38].
The relative electrical conductivity was measured with a conductometer. The leaves were cut into suitable strips, 0.1 g material was placed in 10 mL ddH 2 O for 12 h and the electrical conductivity was measured and classified as R1; next, the samples were heated in boiling water for 30 min, and the electrical conductivity measured and labelled as R2. The relative electrical conductivity was determined as R1/R2 × 100% [34].

RNA-Seq Library Construction and Sequencing
Total RNA was extracted using an RNA extraction kit (Sangon, Shanghai, China, SK1321) and genomic DNA was removed, following the manufacturer's procedure. The total RNA quantity and purity were evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Then, the cleaved RNA fragments were reversetranscribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). The cDNA libraries were sequenced using the Illumina PE150 platform (Illumina, San Diego, CA, USA).

Data Analysis
The quality of the raw data was validated based on the Phred quality score (Q). The raw data were filtered and trimmed for low-quality score reads, adaptor and primer sequences were removed, and raw data were filtered with Trimmomatic [39]. The reads containing ploy-N and the low-quality reads were removed; then, the clean reads were mapped to the reference genomes using HISAT2 [40]. To evaluate the assembly quality, all the paired-end reads were aligned back to these contigs using the Bowtie2 [41]. In-house Perl scripts were used to identify the length distribution of the transcripts, the number and the length of N50, and the average length, maximum length and number of contigs in different intervals.

Alignment and Annotations of New Transcripts
HISAT2 was utilized to map the clean reads to the reference genome. ANNOVAR was used to annotate the reads for alignment based on the reference genome, and the priority annotation standards were the exon region, the shear zone, intron region and the gene spacing [42]. At the same time, the distribution of the region was calculated. All the mapped reads were assembled with StringTie to construct the transcripts [43]; then, the transcripts were combined into a comprehensive transcript set of every sample with the StringTie Merge tool. The new genes and transcripts were identified by comparing the integrated transcript sets with the reference transcript, for which the genome was annotated to obtain information. The encoding capacity of the new transcription proteins was evaluated using CPAT [44].

Expression Analysis
The transcripts per million (TPM) model was used in this study to evaluate the expression of each gene based on transcript abundance [45], and the transcript accounts of each gene were identified through htseq-count [46]. The gene differential expression analysis was selected by Deseq2 v1.10.1, which estimates size factors, and the nbinom test was conducted with a p-value threshold of less than 0.05 and fold changes greater than 2 or less than 0.5 for significantly different expressions [47]. The genes detected in three replicates were considered to be the most effective expressed gene. Gene expression was identified through the average of three replicates for DEGs between treatment and control samples. Fold change refers to the expression ratio between the treatment group and the control group. By screening the log2 fold change and false discovery rate (FDR), genes with significant differences in expression between different treatments could be obtained. Generally, at a standard fold change of less than or equal to four and FDR less than 0.001, both upregulated and downregulated genes were expressed when plants were exposed to salt treatment. The pathway analyses of GO [48] and KEGG [49] enrichment of DEGs were carried out and illustrated in R. GO enrichment analyses and pathway annotation were based on the Gene Ontology Database and KEGG pathway, respectively.

qRT-PCR Analysis
A total of 15 DEGs were selected for evaluating the sequencing accuracy, which represented significant variation among the treatment and control samples. The 15 DEGrelated primers are detailed in Table 1, and two housekeeping actin genes, CsACT and CsEF2, were selected as the reference genes for normalization [50]. The qRT-PCR reactions were carried out in a total system volume of 20 µL (2 µL diluted cDNA, 1 µL forward primer, 1 µL reverse primer, 10 µL SYBR premix ex taq and 6 µL ddH 2 O). The parameters of the qRT-PCR program were set as follows: 95 • C for 3 min, followed by 30 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 20 s. Three technical replicates were calculated under the 2 −∆∆c method, and the average of three biological replicates represents their relative expression levels.

Physiological Analysis of Two Cannabis Seedling Varieties' Responses to Salt Stress
First, at the seedlings stage, of cannabis varieties K94 and W20 were subjected to salt stress treatment of 100 mmol·L −1 , which were sampled at 6 h, 18 h, 24 h, 2 days and 3 days after treatment, and sampling with no treatment was the control check (CK) in the greenhouse. The diverse variations in physiological and phenotypic responses to salt stress were investigated. In both cannabis varieties, the chlorophyll content decreased significantly (p ≤ 0.05) in the first 6 h after salt treatment; the K94 variety reached an extremely remarkable level (p ≤ 0.01), then gradually recovered from 6 h to 24 h, and then decreased with prolonged exposure to salt stress ( Figure 1a). Furthermore, the chlorophyll content reached its maximum levels 24 h after treatment in K94, whereas the maximum was reached after 18 h in W20. At 72 h after treatment, the chlorophyll content was decreased significantly in W20, but there was no significant change in K94. The results also indicated that photosynthesis was inhibited during the first 6 h, then recovered during the 6-24 h time period. Additionally, the degree of photosynthesis decreased significantly from 48 h to 72 h after treatment in W20, but not in K94. (c) the diversity in proline contents. Green represents the K94 variety; orange represents the W20 variety. The letters at the top of the histogram stand for the variation comparisons with the data before and/or after treatment: the same letters stand for data with no significant differences (p ≤ 0.05); two different letters represent differences in the data reaching extremely remarkable level (p ≤ 0.01).
RNA-seq reads generated data for cannabis varieties K94 and W20. To ensure the accuracy of the results, three sequenced samples were randomly selected at the same stage of sampling as a biological replicate. Four categories were severed from two different varieties; in total, 12 independent RNA libraries were prepared and sequenced. Overall, 32.76 Gb of clean data were generated, and the volume of effective data for each sample ranged from 4.06 Gb to 4.22 Gb after removing the adaptor and filtering the low-quality reads with Trimmomatic software. After trimming, the GC contents in samples ranged Green represents the K94 variety; orange represents the W20 variety. The letters at the top of the histogram stand for the variation comparisons with the data before and/or after treatment: the same letters stand for data with no significant differences (p ≤ 0.05); two different letters represent differences in the data reaching extremely remarkable level (p ≤ 0.01).
Although exhibiting diverse physiologies, the results indicated that the two cannabis varieties responded in a similar trend to salt tolerance, with no significant differences in chlorophyll contents at different stages of treatments between K94 and W20. Electrical conductivity increased significantly with the duration of stress exposure during the first day after treatment, and reached maximum values after 24 h: 92.26% in K94 and 96.70% in W20. This trend in increasing electrical conductivity may be caused by the cell membrane being severely injured with increasing the time of salt exposure from treatment to 24 h afterwards, resulting in the outflowing of osmotic substances. Subsequently, the value gradually decreased during the next two days, potentially because the plant recovered with the salt treatment ( Figure 1b). This trend indicated that the electrical conductivity of the two varieties was increasing after salt stress at first, reached the maximum after 24 h, and then decreased.
The proline contents exhibited different trends in the two varieties: it increased between 0 h (CK) and 6 h, reaching a peak at 6 h, decreased until 24 h, then rose and reached the maximum value at 48 h before decreasing in K94; however, the content of proline increased (from 0 to 24 h) and then decreased (from 24 to 72 h) in W20 (Figure 1c). Proline is one substance which can regulate abiotic stress in plants [51]. The diversities in proline content indicated that at the times of 48 h after treatment in variety K94 and 24 h after treatment in variety W20, there were sharp increases in the amount of proline content accumulated. Both reached extremely remarkable levels (p ≤ 0.01); then, the proline contents both decreased.
RNA-seq reads generated data for cannabis varieties K94 and W20. To ensure the accuracy of the results, three sequenced samples were randomly selected at the same stage of sampling as a biological replicate. Four categories were severed from two different varieties; in total, 12 independent RNA libraries were prepared and sequenced. Overall, 32.76 Gb of clean data were generated, and the volume of effective data for each sample ranged from 4.06 Gb to 4.22 Gb after removing the adaptor and filtering the low-quality reads with Trimmomatic software. After trimming, the GC contents in samples ranged from 43.52% to 44.14%, with an average of 44.87%; the Q30 index ranged from 91.96% to 92.35% ( Table 2). All clean reads were mapped to the genome of the assembly number was GCA_900626175.2 with TopHat2 software [52]; the results showed that the average coverages of the four categories were 84.45%, 83.61%, 84.12% and 85.20%, and the trimmed high-quality reads were mapped at unique positions in the gene set of the reference genome. Sample, sample name; Clean reads, the number of clean reads obtained after filtering; Total read pairs, total number of read pairs obtained; GC, the total number of G and C in clean bases constituting the percentage of the total number of bases; Q30, the percentage of bases whose phred QC value is greater than 30 in the RAW bases; Total mapped reads, statistics on the number of sequenced sequences that are located on a genome; in general, if there is no contamination and the reference genome selection is appropriate, the percentages of these data are greater than 70%; Unique mapped reads, statistical count of the number of sequencing sequences with unique alignment positions on the reference sequence; non-spliced reads, whole-segment alignment to exon sequencing sequence statistics; Multiple mapped reads, statistical enumeration of sequencing sequences with multiple alignment locations on the reference sequence.
A principal component analysis (PCA) was performed to verify the repeatability and reliability of the results; this indicated that the samples had extensive and diverse variations, and performing repetitions with the samples indicated relatively good repeatability of the analysis (Figure 2a). The results of expression pattern relationships between repetitions and treatments indicated that there was significant separation between the two cannabis varieties; additionally, there was significant separation between the samples before and after treatment. The repetitions between each treatment were clustered together, which indicated that the experimental procedure had relatively good repeatability and reliability ( Figure 2b). spliced reads, whole-segment alignment to exon sequencing sequence statistics; Multiple mapped reads, statistical enumeration of sequencing sequences with multiple alignment locations on the reference sequence.
A principal component analysis (PCA) was performed to verify the repeatability and reliability of the results; this indicated that the samples had extensive and diverse variations, and performing repetitions with the samples indicated relatively good repeatability of the analysis (Figure 2a). The results of expression pattern relationships between repetitions and treatments indicated that there was significant separation between the two cannabis varieties; additionally, there was significant separation between the samples before and after treatment. The repetitions between each treatment were clustered together, which indicated that the experimental procedure had relatively good repeatability and reliability (Figure 2b).  The transcript abundance of each expression genes was evaluated with the transcripts per million (TPM) model, and a total of 31,735 expression genes were identified. The expression levels among different samples were analyzed, with the results showing that there were rarely differences in the expression levels in all samples. The samples with the maximum expression values were W20 varieties before treatment (Figure 3a). The majority of transcripts in all samples were located in the exon region, as per the gene annotation ( Figure 3b). There was also a situation where some transcripts were located in the intergenic region, with the second highest percentage of genes located in these intergenic regions, which indicated that the cannabis genome was insufficient. Otherwise, it was also identified that some transcript annotations were in the intronic region, although the percentage was relatively low, which may have been the result of alternative splicing, because the main type of alternative splicing was intron retention in plants [53].
Agronomy 2021, 11, x FOR PEER REVIEW 8 the maximum expression values were W20 varieties before treatment (Figure 3a). The jority of transcripts in all samples were located in the exon region, as per the gene a tation (Figure 3b). There was also a situation where some transcripts were located in intergenic region, with the second highest percentage of genes located in these interg regions, which indicated that the cannabis genome was insufficient. Otherwise, it was identified that some transcript annotations were in the intronic region, although the centage was relatively low, which may have been the result of alternative splicing, bec the main type of alternative splicing was intron retention in plants [53]. There were also numerous reads which failed to map the cannabis genome; n transcripts were assembled with Stringtie v1.3.6 software, and a total of 16,149 novel scripts were identified, along with 11,304 annotations of 5392 known genes, accoun for 70.00% of new transcripts. This study also predicted the protein coding ability of transcripts by software CPAT v1.2.4, a total of 6032, accounting for 37.35% of novel scripts, were predicted to be lncRNAs. The novel transcripts were annotated by NR, KEGG and Swiss-Prot database, with the maximum proportion of annotation perfor with NR blast, which could annotate 89.50% of the new transcripts with the NR data There were also numerous reads which failed to map the cannabis genome; novel transcripts were assembled with Stringtie v1.3.6 software, and a total of 16,149 novel transcripts were identified, along with 11,304 annotations of 5392 known genes, accounting for 70.00% of new transcripts. This study also predicted the protein coding ability of new transcripts by software CPAT v1.2.4, a total of 6032, accounting for 37.35% of novel transcripts, were predicted to be lncRNAs. The novel transcripts were annotated by NR, GO, KEGG and Swiss-Prot database, with the maximum proportion of annotation performed with NR blast, which could annotate 89.50% of the new transcripts with the NR database, followed by Swiss-Prot (67.00%), GO (28.27%) and KEGG (27.44%) (Figure 4). of the annotation region. The expression density of genes in twelve samples was mapped using transcripts per million (TPM) reads. The abscissa was the sample name, and the ordinate was log2 (TPM + 1). The box chart for each region is divided into five statistics (from top to bottom are the maximum, upper quartile, median, lower quartile and minimum).
There were also numerous reads which failed to map the cannabis genome; novel transcripts were assembled with Stringtie v1.3.6 software, and a total of 16,149 novel transcripts were identified, along with 11,304 annotations of 5392 known genes, accounting for 70.00% of new transcripts. This study also predicted the protein coding ability of new transcripts by software CPAT v1.2.4, a total of 6032, accounting for 37.35% of novel transcripts, were predicted to be lncRNAs. The novel transcripts were annotated by NR, GO, KEGG and Swiss-Prot database, with the maximum proportion of annotation performed with NR blast, which could annotate 89.50% of the new transcripts with the NR database, followed by Swiss-Prot (67.00%), GO (28.27%) and KEGG (27.44%) (Figure 4).

The Expression Pattern of the Differential Expression Genes Analyzed in K94
The differential expression genes (DEGs) were identified with Deseq2 v1.10.1 software. A total of 249 differentially expressed genes with significant variation were identified before and after treatment with K94, in which 143 genes were upregulated and 106 genes were downregulated (Figure 5a). The number of upregulated DEGs was greater than the number of downregulated DEGs in the clustering analysis (Figure 5b). The clusters analysis of the samples based on the expression was divided into two categories, which divided the samples into controls and treatment of variety K94. All the DEGs' broad biological functions were determined and unique DEGs were identified by Gene Ontology (GO) enrichment analysis. The results based on the DEGs in the K94 variety showed that the vacuole, cytosol, mitochondrion and lysosome were enriched with the top cellular component, and DNA binding and RNA binding were enriched with the top molecular function (Figure 5c). Further KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were conducted; the results indicated that 18 pathways were significantly enriched, based on the 0.05 threshold Q-value, and the top 10 pathways are illustrated in Figure 5d. The pathways of peroxisome, plant hormone signal transduction, flavonoid biosynthesis and plant circadian rhythm were the most enriched cellular processes; environment information processing, metabolism and organismal systems were also enriched. In our study, environmental information processing was the focus; the plant hormone signal transduction pathway and MAPK signaling pathway were the most enriched pathways.
There was diverse variation in the expression profiles of DEGs which were related to the plant hormone signal transduction pathway through KEGG enrichment (Figure 6a). The results showed that nine DEGs exhibited significant variations in expression, where five of them were downregulated and four were upregulated. One of them was the novel gene, in which the expression was also significantly decreased after salt treatment. We also analyzed the expression variation of DEGs which were associated with the MAPK signaling pathway, and the results showed that six DEGs were significantly enriched in this pathway, where two of them were downregulated and four were upregulated. One of the DEGs was the novel gene. There were more upregulated genes than downregulated genes with the MAPK signaling pathway (Figure 6b); there were more downregulated genes in the pathway of plant hormone signal transduction. It was interesting that there were five overlapping DEGs, which was related to both pathways; these may be the pathways with the closest relationship to cannabis responses to salt stress.
showed that the vacuole, cytosol, mitochondrion and lysosome were enriched with the top cellular component, and DNA binding and RNA binding were enriched with the top molecular function (Figure 5c). Further KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were conducted; the results indicated that 18 pathways were significantly enriched, based on the 0.05 threshold Q-value, and the top 10 pathways are illustrated in Figure 5d. The pathways of peroxisome, plant hormone signal transduction, flavonoid biosynthesis and plant circadian rhythm were the most enriched cellular processes; environment information processing, metabolism and organismal systems were also enriched. In our study, environmental information processing was the focus; the plant hormone signal transduction pathway and MAPK signaling pathway were the most enriched pathways.  There was diverse variation in the expression profiles of DEGs which were related to the plant hormone signal transduction pathway through KEGG enrichment (Figure 6a). The results showed that nine DEGs exhibited significant variations in expression, where five of them were downregulated and four were upregulated. One of them was the novel gene, in which the expression was also significantly decreased after salt treatment. We also analyzed the expression variation of DEGs which were associated with the MAPK signaling pathway, and the results showed that six DEGs were significantly enriched in this pathway, where two of them were downregulated and four were upregulated. One of the DEGs was the novel gene. There were more upregulated genes than downregulated genes with the MAPK signaling pathway ( Figure 6b); there were more downregulated genes in the pathway of plant hormone signal transduction. It was interesting that there were five overlapping DEGs, which was related to both pathways; these may be the pathways with the closest relationship to cannabis responses to salt stress.

The Expression Pattern of DEGs Analyzed in W20
There were 456 DEGs identified, with significant variation in DEG expression after salt treatment in W20: 122 upregulated genes and 234 downregulated genes (Figure 7a). The upregulated and downregulated DEGs were divided into two categories based on the expression profile. The heatmaps of the DEGs' expression levels indicated that the sam-

The Expression Pattern of DEGs Analyzed in W20
There were 456 DEGs identified, with significant variation in DEG expression after salt treatment in W20: 122 upregulated genes and 234 downregulated genes (Figure 7a). The upregulated and downregulated DEGs were divided into two categories based on the expression profile. The heatmaps of the DEGs' expression levels indicated that the samples were clustered in two categories, dividing control and treatment samples. The expression profiles of downregulated DEGs in W20 were at extremely high levels (Figure 7b).
Analysis of GO enrichment (Figure 7c) showed that the intracellular cytoplasm, related with the cellular component; the catalytic activity, binding of which is related to molecular functions; and the metabolic process, i.e., the biological process, were significantly enriched. Combining this with the results of the KEGG enrichment (Figure 7d) indicated that the phenylpropanoid biosynthesis pathway, plant hormone signal transduction, pyrimidine metabolism and linoleic acid metabolism were significantly enriched.
DEGs related to the phenylpropanoid biosynthesis pathway (Figure 8a) and plant hormone signal transduction pathway (Figure 8b) exhibited expression patterns. There were nine DEGs related to the phenylpropanoid biosynthesis, three of which were upregulated genes, and six were downregulated genes in comparison to the expression levels between control and treatment samples in W20. There were seven DEGs related to the plant hormone signal transduction pathway, all of which were downregulated. The DEGs did not overlap between two pathways, which indicated the high diversity variations for responding to salt treatment in W20. were nine DEGs related to the phenylpropanoid biosynthesis, three of which were upregulated genes, and six were downregulated genes in comparison to the expression levels between control and treatment samples in W20. There were seven DEGs related to the plant hormone signal transduction pathway, all of which were downregulated. The DEGs did not overlap between two pathways, which indicated the high diversity variations for responding to salt treatment in W20.

Identification of DEGs between K94 and W20
There were 85 DEGs overlapping K94 and W20 among the control and treatment samples for responding to salt treatment (Figure 9a). The heatmap shows that some DEGs significantly reduced the expression levels among the control samples in comparison to treatment samples. There were more downregulated DEGs than the upregulated DEGs (Figure 9b). The downregulated DEGs in the W20-1 sample had the highest expression level in comparison with the other 11 samples. The oxidoreductase activity, heme binding, etc., related to the molecular function ( Figure 9c); and morphogenesis, cell growth, etc., related to the biological process, were identified through GO enrichment for the responses to salt treatment in K94 and W20 cannabis varieties. The signal transduction pathway, MAPK signaling pathway and plant hormone signal transduction were related to salt tolerance, based on the KEGG enrichment analysis (Figure 9d). The MAPK signaling pathway and plant hormone signal transduction were identified to overlap in K94 and W20.
The expression patterns of the DEGs which were related to the signal transduction pathway, MAPK signaling pathway, environmental information processing and plant hormone signal transduction are presented in Figure 10a-d. Interestingly, the DEGs which were related to the signal transduction pathway and MAPK signaling pathway were the same: three were upregulated and two were downregulated. The DEGs related to other pathway expression profiles were the same 5 genes, which showed that these DEGs with significantly diverse variations in expression levels which could enriched two different pathways. This result may indicate that these five DEGs which regulated the mechanism presented different capacities for responding to the salt treatment in cannabis.

Comparison of DEG Tag Data with qRT-PCR
To validate the expression profiles obtained by RNA-seq, real-time PCR was carried out to identify the 15 DEGs that exhibited significant variation in expression between control and treatment samples ( Figure 11). The DEGs which were selected had great variation in expression before and after the salt treatment. Overall, thirteen DEGs had the same variation trends in the TPM, with two other DEGs which had different patterns were identified in the RNA-seq profile analysis among the control and treatment samples. The 86.67% agreement of the RNA-seq with the qRT-PCR illustrated the high accuracy of RNA-seq, which indicated that the pathway and candidate DEGs which were identified were reliable for the salt stress response in cannabis. There were no expression levels detected for genes in some samples, indicating that the DEGs were not expressed at that stage. The expression levels of DEGs LOC115697737, LOC115698674, LOC115710226, LOC115716727 and LOC115700775 were accompanied by decreasing variation with the time of salt treatment; the expression levels were relatively low between 1 day and 2 days after treatment. Eight DEGs exhibited a gradual decline in expression level after the salt stress; the same was observed with the TPM level. There was no detected expression level of LOC115704163 at the time of 0 h or 6 h in W20, and 6 h in K94, indicating that the DEG was not expressed until responding to the salt stress in W20; the expression level of this DEG decreased first and then increased in K94, with the level especially low between 0 h and 6 h. With the same trend in other DEGs, the expression level in LOC115709853 similarly increased after salt treatment. The expression levels of LOC115715759 and LOC115707929 were not consistent with the TPM values. Interestingly, there were variations in the expression levels in LOC115715759; the expression was gradually increased in W20 after salt treatment, but the expression levels were not detected in K94.
son to treatment samples. There were more downregulated DEGs than the upregulated DEGs (Figure 9b). The downregulated DEGs in the W20-1 sample had the highest expression level in comparison with the other 11 samples. The oxidoreductase activity, heme binding, etc., related to the molecular function ( Figure 9c); and morphogenesis, cell growth, etc., related to the biological process, were identified through GO enrichment for the responses to salt treatment in K94 and W20 cannabis varieties. The signal transduction pathway, MAPK signaling pathway and plant hormone signal transduction were related to salt tolerance, based on the KEGG enrichment analysis (Figure 9d). The MAPK signaling pathway and plant hormone signal transduction were identified to overlap in K94 and W20.
The expression patterns of the DEGs which were related to the signal transduction pathway, MAPK signaling pathway, environmental information processing and plant hormone signal transduction are presented in Figure 10a-d. Interestingly, the DEGs which were related to the signal transduction pathway and MAPK signaling pathway were the same: three were upregulated and two were downregulated. The DEGs related to other pathway expression profiles were the same 5 genes, which showed that these DEGs with significantly diverse variations in expression levels which could enriched two different pathways. This result may indicate that these five DEGs which regulated the mechanism presented different capacities for responding to the salt treatment in cannabis.

Comparison of DEG Tag Data with qRT-PCR
To validate the expression profiles obtained by RNA-seq, real-time PCR was carried out to identify the 15 DEGs that exhibited significant variation in expression between control and treatment samples ( Figure 11). The DEGs which were selected had great variation in expression before and after the salt treatment. Overall, thirteen DEGs had the same variation trends in the TPM, with two other DEGs which had different patterns were identified in the RNA-seq profile analysis among the control and treatment samples. The 86.67% agreement of the RNA-seq with the qRT-PCR illustrated the high accuracy of RNA-seq, which indicated that the pathway and candidate DEGs which were identified were reliable for the salt stress response in cannabis. There were no expression levels detected for genes in some samples, indicating that the DEGs were not expressed at that stage. The expression levels of DEGs LOC115697737, LOC115698674, LOC115710226, LOC115716727 and LOC115700775 were accompanied by decreasing variation with the time of salt treatment; the expression levels were relatively low between 1 day and 2 days after treatment. Eight DEGs exhibited a gradual decline in expression level after the salt stress; the same was observed with the TPM level. There was no detected expression level of LOC115704163 at the time of 0 h or 6 h in W20, and 6 h in K94, indicating that the DEG was not expressed until responding to the salt stress in W20; the expression level of this DEG decreased first and then increased in K94, with the level especially low between 0 h

Discussion
Over 20% of the agricultural land around the world is affected by soil salinization [54], which seriously affects crop yields and security. Some researchers have divided plants into two categories based on their capacity of responding to salt stress: halophytic and glycophytic [55]. Those in the halophyte category could adapt to the salinized environment and have a high capacity to cope with soil salinization; in contrast, those in the glycophyte category are severely damaged by a salinized environment.
Under salt stress, plants are affected at every stage of growth. Salt toxicity to plants mainly includes osmotic and ionic stress, reactive oxygen species (ROS), and nutrient de- Figure 11. Expression profiles of DEGs which were selected based on qRT-PCR and the TPM. Green represents the expression of the K94 variety at different stages after treatment based on the qRT-PCR; orange represents the TPM expression of treatment and control samples based on the RNA-seq; blue-violet represents the expression of the W20 variety at different stages after treatment based on qRT-PCR.

Discussion
Over 20% of the agricultural land around the world is affected by soil salinization [54], which seriously affects crop yields and security. Some researchers have divided plants into two categories based on their capacity of responding to salt stress: halophytic and glycophytic [55]. Those in the halophyte category could adapt to the salinized environment and have a high capacity to cope with soil salinization; in contrast, those in the glycophyte category are severely damaged by a salinized environment.
Under salt stress, plants are affected at every stage of growth. Salt toxicity to plants mainly includes osmotic and ionic stress, reactive oxygen species (ROS), and nutrient depletion. When the salt content of soil is increased, the water potential of the soil is lower than that of the plant root tissue, resulting in plant dehydration [56][57][58]. In addition, osmotic stress could lead to stomatal closure; this phenomenon inhibits the absorption and transition of water [59], and also leading to reduced photosynthesis and the induced expression of many stress-responsive genes [60][61][62][63][64]. Ionic stress is mainly caused by the accumulation of sodium (Na + ) and chlorine (Cl − ) in cells, with variations in ionic stress mainly affecting Ca + signaling; some protein kinases are Ca + -dependent [65][66][67]. The toxicity of sodium is mainly due to the inhibitory effects on enzyme activity, which could negatively affect metabolism. The pathways affected by Na + accumulated under salt tolerance mainly include photosynthesis, the Calvin cycle, etc. Additionally, plants can regulate the harm caused by salt stress, such as by restricting Na + accumulation and stabilizing K + concentrations [26,62,68]. In addition, excess Na + accumulated in the cytoplasm can also affect the absorption and transportation of K + , Ca + and other mineral elements such as nitrogen, phosphorus and zinc. If there is excess Cl − in the cytoplasm, the absorption and transition capacities of key macronutrients, such as nitrogen and sulphone, are affected, which leads to deficiencies. Non-selective anion transporters such as Cl − mediate the uptake channels of NO 3− and SO 4 2− in plants [69][70][71]. Furthermore, salt stress induces variations in osmotic and ionic stress, and could also induce the accumulation of ROS, which severely damages the cell structure and macromolecules, such as DNA, lipids, enzymes, etc.
Soil salinity is a severe abiotic stress which severely restricts crop productivity [72,73]; for coping with the effects of ionic toxicity [74] and osmotic shock [75,76], plants have evolved a series of response mechanisms. With the development of salt stress research, many genes which are associated with salt tolerance have been identified or cloned in different plants. The HKT gene family [77], SNRK gene family [78], SPL gene family [79] and PGDH gene family [80], among others, have been identified to regulate salt stress. Correlational research has presented considerable progress in exploring the salt responses of major plants, such as Arabidopsis [81][82][83], soybean [84][85][86], rice [87][88][89] and maize [90][91][92]. Plant breeders and biotechnologists have tried to develop crop varieties with a strong capacity to respond to salt stress, mainly including improving the plant's salt tolerant and salt resistance. However, satisfactory results have not yet been achieved in agroecologically important crops due to many factors. Melatonin was discovered to accumulate in response to salt stress; this strategy could enhance the capacity of plant resistance to salt stress, predominantly through eliminating the reactive oxygen species and enhancing antioxidant enzyme activity [93].
Similarly to common crops, cannabis was able to accumulate high levels of NaCl in its roots and shoots in a short-term hydroponic experiment. Through the analysis of Illumina high-throughput mRNA sequencing, over 31,000 high-quality transcript contigs have been identified in cannabis for responding to salt stress. Furthermore, a total of 16,149 novel transcripts have been identified, including some candidate DEGs and pathways which respond to salt stress in cannabis. One key observation in this study was the response of different strains of cannabis to salt stress.
Transcriptome analysis showed that a large number of genes were upregulated or downregulated under salt stress: there were more upregulated DEGs than downregulated DEGs in cannabis varieties K94 and W20, with some overlapping DEGs also identified. In addition, there are significant differences in the pathways which control the efficiency of different cannabis varieties to respond salt stress. These conclusions were consistent with previous studies in rice [45,[83][84][85][86] and maize [94][95][96] responses to salt stress. The "plant hormone signal transduction" pathway, "MAPK signal pathway" and "phenylpropanoid biosynthesis" were detected to be significantly associated with the salt stress through the analysis of GO enrichment and KEGG enrichment based on the DGEs which were identified in K94 and W20 varieties. The diversity in molecular variation was strongly related with the phenotypical and physicochemical properties, namely, salt stress can induce increases in lipid peroxidation and ROS in cannabis, and upregulate the expression of antioxidant genes in plants to cope with ionic toxicity [82,[97][98][99][100].
The "MAPK signal pathway" is strongly involved in the response to salt stress in plants; the Ca + -dependent protein kinase CPK3 is required for the MAPK pathway in Arabidopsis [101]. Several modules of the MAPK pathway are active in receiving and transporting signals under salt stress, and these MAPK modules regulate the production, transportation and distribution of auxin, which enhances the capacity to tolerate salt stress in wild barley [102]. Overexpressing the GhCIPK6a gene could enhance the capacity to respond to salt stress, involving ROS scavenging in the MAPK pathway in cotton [103]. The "plant hormone signal transduction" pathway has also previously been reported to play an essential role. In pumpkins, enrichment analyses of the DEGs which are related to "plant hormone signal transduction" indicate that the ethylene signal transduction genes may be related to this pathway, exhibiting important effects in responding to salt stress [104]. In sophora, auxin, cytokinin and brassinosteroids have been identified to influence the recovery after salt stress through plant hormone signal transduction [105]. In cotton, the ABA-and GA-related DEGs of plant hormone signal transduction are regulated by mediating the expression of hormone-related DEGs to respond to salt stress [106].
Among the enrichment analyses of GO and KEGG in the K94 and W20 cannabis varieties, the results indicated that the "plant hormone signal transduction" pathway, "MAPK signal pathway" were both detected, and the DEGs related to these pathways were the main candidate genes for cannabis to respond to salt stress. However, other different pathways were enriched between the two varieties, showing different strategies in the capacity for resisting salt stress among different of cannabis varieties.
It is necessary to develop a comprehensive understanding of the physiology, biochemistry, metabolism and gene expression of cannabis under salt stress in order to improve yield. In this study, the capacity and expression variation of DEGs in salt stress responses were identified for two cannabis varieties with different genetic backgrounds. Moreover, more gene expression profiles were identified in cannabis, enabling further understanding of cannabis genomic genetics, and which will be beneficial for improving agricultural traits in cannabis.

Conflicts of Interest:
The authors declare no conflict of interest.