Comparative transcriptome analysis reveals significant differences in the regulation of gene expression between hydrogen cyanide- and ethylene-treated Arabidopsis thaliana

Hydrogen cyanide (HCN) is a small gaseous molecule that is predominantly produced as an equimolar co-product of ethylene (ET) biosynthesis in plants. The function of ET is of great concern and is well studied; however, the function of HCN is largely unknown. Similar to ET, HCN is a simple and diffusible molecule that has been shown to play a regulatory role in the control of some metabolic processes in plants. Nevertheless, it is still controversial whether HCN should be regarded as a signalling molecule, and the cross-talk between HCN and ET in gene expression regulation remains unclear. In this study, RNA sequencing (RNA-seq) was performed to compare the differentially expressed genes (DEGs) between HCN and ET in Arabidopsis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were subsequently performed to investigate the function and pathway enrichment of DEGs. Parts of key genes were confirmed by quantitative real-time PCR. The results showed that at least 1305 genes and 918 genes were significantly induced by HCN and ET, respectively. Interestingly, a total of 474 genes (|log2 FC| ≥1) were co-regulated by HCN and ET. GO and KEGG analyses indicated that the co-regulated genes by HCN and ET were enriched in plant responses to stress and plant hormone signal transduction pathways, indicating that HCN may cooperate with ET and participate in plant growth and development and stress responses. However, a total of 831 genes were significantly induced by HCN but not by ET, indicating that in addition to ET, HCN is in essence a key signalling molecule in plants. Importantly, our data showed that the possible regulatory role of a relatively low concentration of HCN does not depend on ET feedback induction, although there are some common downstream components were observed. Our findings provide a valuable resource for further exploration and understanding of the molecular regulatory mechanisms of HCN in plants and provide novel insight into HCN cross-talk with ET and other hormones in the regulation of plant growth and plant responses to environmental stresses.


Background
Hydrogen cyanide (HCN) is a natural metabolite in bacteria, fungi and plants that has received special attention from scientists since the beginning of the nineteenth century, due to its toxic effect on living organisms [1,2]. Acute or chronic exposure to HCN can lead to intoxication, mild to severe illness, and in extreme cases even death in humans and animals because a certain dose of HCN inhibits the activity of metalloenzymes, principally cytochrome c oxidase, the final enzyme in the respiratory electron transport chain [3]. In plants, HCN is an equimolar by-product of ethylene (ET) biosynthesis via the ACC pathway [1]. In certain physiological states, such as fruit ripening and flower senescence, and in many environmental conditions, such as flooding and chilling, ET biosynthesis is greatly induced and at the same time, HCN accumulats rapidly [4]. In addition, HCN is released from cyanogenic lipids and cyanogenic glycosides during tissue disruption, infection, or cyanogenic plant food processing such as cassava roots (Manihot esculenta). There are more than 3000 species of higher plants including ferns, gymnosperms and angiosperms produce cyanogenic glucosides that are actively cleaved to produce HCN [5].
ET is certainly one of the most important plant hormones and is involved in numerous biological processes, such as root growth and stem elongation, fruit ripening, senescence, response to pathogens, response to gravity or submergence stress [6,7]. In contrast to the roles of ET in plants, HCN is generally regarded as a phytotoxic agent and plays a protective effect in plants against predators such as herbivores [8,9]. However, abundant information indicates that HCN, apart from being toxic, plays a regulatory (perhaps signalling) function in many physiological processes, e.g. seed germination, nitrate assimilation or in plant responses to some environmental stimuli [2,[10][11][12]. It was shown that exogenous 1 mM HCN treatment significantly promoted seed dormancy breakage and stimulated germination by inducing the production of reactive oxygen species (ROS) or feedback regulation of ET synthesis and signalling transduction [13][14][15]. Chivasa and Carr (1998) reported that potassium cyanide (KCN; 0.5 mM) pre-treatment could protect plants against tobacco mosaic virus (TMV) infection [16]. Additionally, the research of Seo et al. (2011) indicated that exogenous cyanide (0.5 or 1.0 mM) treatments contributed to the resistance of rice to blast fungus [11]. Our previous study demonstrated that a lower concentration of cyanide (20 μM KCN) could enhance cucumber seedlings against abiotic stress, such as drought stress and salinity stress [12].
Despite the well-described effects of HCN on germination and stress acclimation, its molecular regulatory mechanism remains largely unknown. It is possible that HCN plays a positive role in ET biosynthesis feedback regulation under conditions where rapid HCN accumulation occurs [17]. Some evidence proposed that the beneficial effect of HCN in seed germination might be due to feedback regulation of ET synthesis and signalling transduction [13], while other evidence hypothesized that HCN might be an important defence signal rather than ET in plant resistance to biotic stress [11].
HCN may play a dual role in plants, depending on its concentration [2]. HCN may be used in defense against herbivores at high toxic concentrations and may have a regulatory function at lower concentrations. Generally, most HCN produced in plants is detoxified quickly by the key enzyme β-cyanoalanine synthase (CAS) [18]. The remaining HCN at a lower level, probably at a non-toxic concentration, may act as a signalling molecule involved in the control of some metabolic processes in plants [2]. In Arabidopsis, it was shown that the CAS is encoded by a small gene family of three members, CYS-C1 (At3g61440), CYS-D1 (At3g04940), and CYS-D2 (At5g28020) [19]. The most abundant CAS enzyme is CYS-C1, which is localized in the mitochondria and contributes most of the CAS activity in root and leaf tissue [20]. Interestingly, Garcia et al. (2010) mentioned that the accumulation of HCN within the mitochondrion in the Arabidopsis cys-c1 mutant can act as an inhibitor of root hair development but is not toxic to the plant [21]. They also found that HCN accumulation in cys-c1 mutant plants presented an increased susceptibility to the necrotrophic fungus Botrytis cinerea and an increased tolerance to the biotrophic Pseudomonas syringae pv tomato DC3000 bacterium and Beet curly top virus, suggesting that HCN might act by stimulating the salicylic acid (SA)-dependent signalling pathway of the plant immune system [22,23]. However, how HCN participates in plant growth and stress responses is still largely unknown.
The aim of this study, therefore, is to uncover the regulatory role of HCN in Arabidopsis and compare the significant differences in gene expression regulation with ET. RNA sequencing (RNA-seq) methods were used to analyse the differentially expressed genes (DEGs) after treatment with HCN and ET, compared to the control. The data in this study provide a valuable resource for further exploration and understanding of the detailed molecular mechanisms of HCN in plant growth and response to environmental stress, and also provide novel insight into the cross-talk between HCN, ET and other hormones in Arabidopsis.

Illumina sequencing and gene function annotation
In this study, a total of nine cDNA libraries from HCN-treated, ET-treated and control (CK) seedlings were prepared and subjected to Illumina deep sequencing, with each group created in triplicate. After removing the adaptors and low-quality sequences, the Illumina sequencing generated 467,217,912 sequence reads and 58.4 Gb of sequence data ( Table 1). The good sequence and good ratio of each library was over 99%. In addition, the GC content of each library was approximately 45%, and CycleQ30% was greater than 90% for each library. Thus, the quality and accuracy of the sequencing data were sufficient for further analysis. All reads were aligned to the Arabidopsis genome and the results showed that most of the reads matched Arabidopsis genomic locations. All the unigenes matched previously described sequences with approximately 92% coverage.
To elucidate potential gene functions, the gene annotation was carried out against KOG, SwissProt, TrEMBL, GO, KEGG databases. As shown in Fig. 1 (Fig. 1).

Differentially expressed genes (DEGs) analysis
To compare the upregulated or downregulated genes under the conditions of HCN and ET treatment, the differentially expressed genes (DEGs) were analysed according to the results of RNA-seq. As shown in Fig. 2a, a total of 6719 and 5230 DEGs were detected in the samples of CK vs HCN and CK vs ET, respectively. By comparison, a total of 3512 DEGs were upregulated and 3207 DEGs were downregulated in the HCN-treated samples (Fig. 2a), while 2807 DEGs upregulated and 2423 DEGs downregulated in ET-treated samples when compared with the control samples (Fig. 2a). If the absolute log 2 fold change (log 2 FC) ≥1 were used to judge the significance of differences in gene expression, a total of 1305 DEGs and 918 DEGs were detected in HCN-treated and ET-treated samples, respectively (Fig.  2b). There were 567 DEGs up-regulated and 738 DEGs down-regulated by HCN treatment, and 488 DEGs up-regulated and 430 DEGs down-regulated by ET treatment (Fig. 2b). These data indicated that although HCN is a co-product of ET, as expected it plays important roles in gene regulation in Arabidopsis thaliana.
To confirm the reliability of the transcriptome sequencing, parts of DEGs regulated by HCN or ET were investigated by qRT-PCR. As shown in Additional file 1: Figure  S1, the results of qRT-PCR were generally consistent with the transcriptome data. In addition, linear regression analysis of the correlation between qRT-PCR and RNA-seq showed an R 2 (goodness-of-fit) value of 0.9352 and a corresponding slope of 0.9322, suggesting a strong positive correlation between the qRT-PCR and transcriptome data.

Genes commonly regulated by HCN and ET
As mentioned above, HCN and ET may play regulatory roles in gene expression in Arabidopsis. Therefore, it is interesting to know how many genes were commonly regulated by both HCN and ET. The results showed that a total of 3474 DEGs were commonly regulated by HCN and ET, of which 1643 DEGs and 1626 DEGs were commonly upregulated and downregulated, respectively (Fig. 3a). As shown in Fig. 3b, a total of 474 genes were  Figure S2. Interestingly, among all of the 474 co-regulated genes (|log 2 FC| ≥1), 184 genes upregulated and 267 genes downregulated by both HCN and ET, whereas only 23 genes were downregulated by HCN but upregulated by ET. The details of the top 10 commonly upregulated genes are listed in Table 2. Of these, the gene expression of AT2G15020 was upregulated by 23-fold and 6-fold with treatment of HCN and ET, respectively, compared to the CK.

Genes exclusively induced by HCN or ET
In addition to the co-regulated genes, the results showed that a total of 831 genes (|log 2 FC| ≥1) were exclusively induced by HCN, including 383 genes that were upregulated and 448 genes that were downregulated by HCN (Fig. 4a).
In comparison, there were approximately 444 genes (|log 2 FC| ≥1) specially induced by ET, including 281 genes that were upregulated and 163 genes that were downregulated by ET (Fig. 4b). The details of the top 10 significantly upregulated genes by HCN or ET are listed in Table 3. It was shown that the DEGs significantly induced by HCN but not by ET included the genes AT1G56600 (galactinol synthase 2, GolS2) and AT5G37990 (S-adenosyl-L-methionine-dependent methyltransferase superfamily protein), which were 15-fold (log 2 FC = 3.91) and 12.9-fold (log 2 FC = 3.69) higher than in the CK. The DEGs significantly induced by ET but not by HCN included AT3G11340 (UDP--Glycosyltransferase superfamily protein, UGT76B1) and AT4G25490 (C-repeat/DRE binding factor 1, CBF1), which were approximately 12.8-fold (log 2 FC = 3.7) and 9.6-fold (log 2 FC = 3.3) higher than in the CK, respectively.
GO and KEGG analysis show that the DEGs regulated by HCN are enriched in the plant hormone signal transduction pathway To further investigate the functions of DEGs regulated by HCN or ET, GO and KEGG analyses were carried out. According to the GO analysis, the DEGs enriched in HCN-treated samples were assigned 720 GO terms (P < 0.05) and 219 GO terms (FDR < 0.05), while a total of 637 GO terms (P < 0.05) and 149 GO terms (FDR < 0.05) were significantly observed in ET-treated samples (Fig. 5a). There were 327 GO terms (P < 0.05) and 105 GO terms (FDR < 0.05) assigned to both HCN and ET (Fig. 5b, c). It should be noted that 481 GO terms of HCN regulated DEGs and 411 GO terms of ET regulated DEGs were associated with biological process (BP) (Fig. 5 d, e). In addition, there were no significant differences in the top 30 GO terms between HCN and ET, but the number of DEGs in each terms was significantly different between them (Fig. 6). The top 30 common GO terms assigned to BP between CK vs HCN and CK vs ET are shown in Additional file 1: Table S2, which shows that HCN-and ET-regulated DEGs were mainly enriched in plant responses to stimuli, environmental stress and hormones. These data further indicated that HCN and ET may be cooperatively involved in regulating plant growth and development, and plant resistance to environmental stresses. KEGG analysis showed that the DEGs from CK vs HCN and CK vs ET were assigned to 24 and 21 KEGG pathways (P < 0.05), respectively (Fig. 7a). The common KEGG terms (P < 0.05) between CK vs HCN and CK vs ET were assigned to 10 pathways (Fig. 7a). Of these, the common pathway with the greatest number of DEGs was plant hormone signal transduction (ko04075) (Fig.  7c), consistent with the results of the GO analysis. Further analysis revealed that 103 DEGs (CK vs HCN) and 94 DEGs (CK vs ET) were assigned to the plant hormone signal transduction pathway (Fig. 7d). Among them, a total of 61 DEGs were co-regulated by both  Figure S3 and Table S3). Interestingly, the genes co-regulated by HCN and ET were mainly associated with the auxin signalling transduction pathway, including the family genes from SMALL AUXIN UP RNAs (SAURs), GRETCHEN HAGEN3s (GH3s) and AUXIN/INDOLE ACETIC ACID (AUX/IAA) (Additional file 1: Table S3). Moreover, the common DEGs regulated by HCN and ET were also involved in the abscisic acid (ABA) signalling pathway (e.g., HAB1; hypersensitive to ABA1), cytokinin (CTK) signalling pathway (e.g., ARR6; cytokinin response regulator 6), gibberellin (GA) signalling pathway (e.g., RGL2; RGA-like 2) and brassinosteroid (BR) signaling pathway (e.g., BSK2; brassinosteroid signaling kinase 2) (Additional file 1: Table S3). These data further indicated that in addition to ET, HCN also engages in cross-talk with other plant hormone signalling molecules.
The effect of HCN treatment on the gene expression of ethylene biosynthesis and signalling pathway As mentioned previously, HCN is a co-product of ET synthesis and probably plays a role in plant growth and stress tolerance via the feedback-inducing synthesis of ET [13,15,17]. Therefore, the DEGs from CK vs HCN were analysed based on the GO annotation and KEGG pathway. As shown in Fig. 8a, there were no significantly up-regulated genes observed in the ET biosynthetic process (GO:0009693) with HCN treatment. In contrast,  the gene expressions of 1-amino-cyclopropane-1-carboxylate synthase (e.g., ACS7) and 1-amino-cyclopropane-1-carboxylate oxidase (e.g., ACO1) were downregulated by HCN treatment, especially for the transcript of ACO1, a key gene for ET synthesis, which was downregulated more than 3.5-fold when compared to the CK. Further analysis of the ET-mediated signalling pathway (GO:0009873) showed that a total of 25 genes were significantly regulated by HCN, including 4 genes that were upregulated and 21 genes that were downregulated (Fig. 8b). For instance, the ET response factors including ERF1B, ERF5 and ERF6 were markedly down-regulated by HCN (Fig. 8b).
The effect of HCN treatment on the gene expression of plant mitochondrial respiration HCN poisoning is known to block the mitochondrial respiratory chain electron transport system by affecting complex IV [24]. Thus, the genes related to the mitochondrial respiratory chain (GO:0005746) were analysed according to the DEGs and GO enrichment analysis data. Interestingly, no significant DEGs (|log 2 FC| ≥1) were found in HCN-treated samples compared to the CK (Table 4). Further analysis of all DEG isoforms (|log 2 FC| < 1) showed that 3 genes were upregulated by HCN, namely, AT5G25450 (cytochrome bd ubiquinol oxidase; fold change = 1.72), AT3G10370 (FAD-dependent oxidoreductase family protein; fold change = 1.41) and AT3G27240 (Cytochrome C1 family; fold change = 1.24) ( Table 4). In contrast, the RNA-seq data showed that ET treatment significantly induced at least 11 genes (|log 2 FC| ≥1) related to the mitochondrial respiratory chain, including cytochrome c oxidase and NADH dehydrogenase.
In addition, it is noteworthy that HCN treatment induced, but not significantly, the gene expression of alternative oxidase 1A (AOX1A; fold change = 1.58) ( Table 5), which is the HCN-insensitive protein that mediates HCN-resistant respiration. In contrast, the RNA-seq data showed that AOX1D (fold change = 0.17) was significantly downregulated by HCN. Moreover, HCN induced, but not significantly, the gene expression of cysteine synthase D2 (CYSD2; fold change = 1.38), which is involved in HCN detoxification. Importantly, the expression of genes involved in ROS production was not significantly affected by HCN treatment (Additional file 1: Table S5). In comparison, ET treatment induced, but not significantly, the expressions of AOX1A (fold change = 1.97) and AOX1C (fold change = 1.73) but slightly reduced the expression of CYSD1 (fold change = 0.78).
Comparative analysis of stress-related genes between HCN and ET treatment As described above and previous studies, treatment with HCN at lower concentrations could enhance plant stress tolerance [12,25]; thus, the DEGs related to stress response were analysed. As shown in Fig. 9a, a total of 315 DEGs (|log 2 FC| ≥1) related to stress (GO:0006950) were found in HCN-treated samples, while a total of 208 DEGs (|log 2 FC| ≥1) were found in ET-treated samples. Further analysis showed that a total of 203 and 92 DEGs in HCN-treated samples were assigned to abiotic stress (GO:0009628) and biotic stress (GO:0009607) (Fig. 9b,  c), respectively, indicating that HCN plays a more important role in plant response to abiotic stress than biotic stress. In ET-treated samples, a total of 141 and 60 DEGs were assigned to abiotic stress and biotic stress, respectively (Fig. 9b, c). A total of 117 DEGs related to stress were commonly regulated by HCN and ET. Of these, a total of 80 DEGs related to abiotic stress and 27 DEGs related to biotic stress were regulated by both HCN and ET (Fig. 9b, c). Further analysis revealed that,  Table S6. Among them, it is interesting to note that the stress-related gene with the highest induced expression level regulating by both HCN and ET was AT3G09440 (heat shock protein 70-3, HSP70-3), which belongs to HSP70 family proteins and is essential for plant development and plant resistance to environmental stress [26]. In addition to HSP70-3, six more HSP protein transcripts were induced by both HCN and ET (Additional file 1: Table S6). Moreover, the gene expression of UDP-Glycosyltransferase superfamily protein (UDP72B1), which is involved in metabolizing xenobiotica (chloroaniline and chlorophenol) and salt stress responses [27,28], was significantly induced by HCN (fold  Table S6).
As shown in Fig. 9, there were 198 DEGs related to stress that were significantly induced by HCN but not by ET. Of these, the top 10 DEGs included the AT1G56600 (galactinol synthase 2, GOLS2), AT3G22840 (Chlorophyll A-B binding family protein, ELIP1), and AT5G36910 (Thionin 2.2, THI2.2), where the expressions were upregulated by 15-fold, 10-fold and 3.7-fold, respectively, when compared to the CK (Additional file 1: Table S7). Interestingly, there were 91 DEGs related to stress that were significantly induced by ET but not by HCN. Among them, the gene with the highest expression level induced by ET was AT3G11340 (UDP-Glycosyltransferase superfamily protein, UGT76B1), which showed approximately 13-fold higher than the CK (Additional file 1: Table S7). These data indicated that HCN and ET play different roles in regulating plant responses to environmental stress, apart from the possible synergistic relationship.

Discussion
Since the confirmation that HCN is formed as a co-product during ET biosynthesis, questions have been raised regarding the physiological significance of this metabolite in plants [9]. It has been proposed that HCN may be a phytotoxic agent at higher concentration and may have a regulatory function at lower concentration [2]. However, the toxic or regulatory function for HCN in plant metabolism remains controversial. In the present study, RNA-seq method was used to reveal the regulatory role of HCN in Arabidopsis thaliana and the significant differences in the regulation of gene expression between HCN-and ET-treated seedlings. The results showed that a total of 1305 DEGs and 918 DEGs were significantly (|log 2 FC| ≥1) regulated by HCN and ET, respectively. A total of 474 genes (|log 2 FC| ≥1) were commonly regulated by HCN and ET. These findings, in agreement with previous studies, indicate that HCN is possibly an important signalling molecule involved in the control of different metabolic, physiological and developmental processes in plants, rather than a waste co-product of ET biosynthesis [2,9]. Since a large number of genes were co-regulated by HCN and ET, suggesting that there may be a synergistic relationship in the regulation of cellular processes between them, as previous observations have shown that they elicit some similar physiological responses including alleviation of seed dormancy and enhancement of plant resistance to environmental stress [14,22].
According to the results of GO and KEGG analysis, we found that the DEGs regulated by HCN and ET were mainly enriched in plants response to stimuli and hormones (Figs. 6 and 7). Interestingly, HCN and ET co-regulated DEGs were mainly associated with auxin signalling transduction, including the family genes from SAURs, AUX/IAA and GH3s (Additional file 1: Table  S3). Auxin is identified as a plant growth hormone that plays a critical role in leaf growth and development and is also implicated in plant defence signalling pathways [29,30]. SAUR genes are rapidly upregulated in response to auxin; however, their functions remain elusive [29]. In comparison, more functions have been uncovered about the Aux/IAA family genes. It has been demonstrated that Aux/IAA genes are required for stress tolerance. The Aux/IAA proteins are auxin-sensitive repressors that mediate diverse physiological and developmental processes in plants [31,32]. Shani et al. (2017) reported that promoting the transcription of IAA5, IAA6 and IAA19 resulted in enhanced tolerance to abiotic stress, whereas recessive mutations in these IAA genes resulted in decreased tolerance to stress conditions [33]. Interestingly, our RNA-seq data showed that, in addition to IAA19 gene (fold change = 1.8), IAA6 and IAA17 were upregulated 3.2-fold and 2.1-fold by HCN, respectively (Additional file 1: Table S4). Thus, these findings encourage us to determine whether these genes contribute to HCN-mediated plant stress tolerance in the future.
In addition to auxin signaling pathway, our data showed that the genes regulated by HCN were also involved in ABA, GA, and BR signalling pathway (Additional file 1: Table S3). It has been demonstrated that the phytohormone ABA serves as an endogenous messenger in the biotic and abiotic stress responses of plants [34]. Under non-stress conditions, ABA is also required to fine-tune growth and development. It is evident that  It is noteworthy that the gene expression of RGL2, which encodes a DELLA protein, was significantly upregulated by HCN and ET 3.8-fold and 2.5-fold higher than the CK, respectively. DELLA proteins act as repressors of GA-dependent processes to retard plant growth [35]. However, under stress conditions, deceleration of the growth is regarded as one of the strategies that help to improve the survival of plants [36,37]. In addition, it is likely that the ability to reduce cell growth under unfavourable conditions may not only allow conservation of energy for defence purposes but also limit the risk of heritable damage [2]. ET signalling was found to provide salt stress and cold stress tolerance by enhancing the function of DELLAs [38,39]. Therefore, it seems that at least part of the growth regulatory action and stress acclimation of HCN is through cross-talk with DELLA proteins, such as RGL2. DELLA proteins have been demonstrated to interact with multiple hormone pathways and they have been regarded as key components of the plant growth regulation and stress response network [40,41]. Thus, more research is necessary to uncover whether DELLAs are involved in the interactions of HCN with other hormones in the future. In this study, parts of genes related to the BR signalling pathway were regulated by both HCN and ET (Additional file 1: Table S3), thus it is worth focusing on the cross-talk among them. BRs are an important group of plant steroid hormones involved in numerous aspects of plant life, including growth, development and response to various stresses [42]. It was shown that BR binds to the extracellular domain of the cell-surface receptor kinase BRASSI-NOSTEROID INSENSITIVE1 (BRI1) to activate BRI1 kinase activity [43]. BRI1 activation involves the recruitment of the co-receptor kinase BRI1-ASSOCIATED RE-CEPTOR KINASE1 (BAK1) and disassociation of the inhibitory protein BRI1 KINASE INHIBITOR1 (BKI1). In addition, it was suggested that BSKs are the substrates of BRI1 kinase that activate downstream BR signal transduction such as BRI1-SUPPRESSOR1 (BSU1). BSU1 can activate BRASSINAZOLE RESISTANT1 (BZR1) and BRI1-EMS-SUPPRESSOR1 (BES1) indirectly by inactivating the kinase BRASSINOSTEROID INSENSITIVE2 (BIN2), which is a negative regulator in the BR signalling pathway [44]. However, the interactions between BRs and ET are still largely unknown. Some studies suggest that BRs positively influence ET biosynthesis through the regulation of ACS and ACC oxidase activities [45], whereas how ET affects the BRs synthesis and signalling transduction remain unclear. In the present study, the RNA-seq data showed that gene expression of the BSK2 was upregulated while the BKI1, BIN2 and BZR1 were downregulated by HCN and ET when compared to the CK (Additional file 1: Table S3). Consequently, it is likely that both HCN and ET positively affect BR signal transduction to some extent.
Strikingly, the genes associated with ET biosynthetic pathway were downregulated by HCN based on the RNA-seq data in this study (Fig. 8), although several previous studies have demonstrated that seed dormancy removal by 1 mM HCN involves modifications in the ET biosynthetic pathway [13,15]. However, our findings are consistent with the report of Garcia et al. (2010) that HCN accumulation in cys-c1 mutant plants showed reduced ET production and repressed expression of several genes related to ET signalling and metabolism, such as ACS6, ERF6, and ERF105, when compared to wild-type plants [21]. Consequently, it seems that the role of HCN in plants does not rely on the ET feedback effect. In fact, our hypothesis is also in agreement with the results of Seo et al. (2011), showing that exogenous cyanide (KCN) but not ET complements blast fungus resistance in ACS/ACO knockdown rice plants [11]. In addition, it should be noted that Oracz et al. (2008) stated that the expression of the transcription factor ERF1 was markedly stimulated by 1 mM HCN during the release of sunflower seed dormancy but, it did not significantly affect ET production or the expression of genes involved in ET biosynthesis or the first steps of the ET signalling pathway [15]. Furthermore, there is no direct evidence that ERF1 is regulated by HCN or other molecular signals such as ROS during seed germination because a higher concentration of HCN (1 mM) has been shown to promote ROS accumulation apparently [10,14]. Importantly, our data showed that the transcription of NADPH oxidase (NADPHox) genes were not significantly affected by 20 μM HCN treatment. As shown in Additional file 1: Table S5, the fold changes of NADPH/respiratory burst oxidase protein D (RbohD) and RbohF in CK vs HCN were 0.77 and 1.38, respectively. Consistent with our findings, Arenas-Alfonseca et al. (2018) reported that the effect of HCN on root hair elongation is independent of H 2 O 2 production and direct NADPH oxidase inhibition [46]. Taken together, we speculate that HCN itself, especially at lower concentrations, should serve as a signalling molecule in plants but more studies are needed to reveal its target proteins (receptors) and its cross-talk with ET under certain conditions.
Our data suggest that a lower concentration of HCN (20 μM) treatment may not significantly affect mitochondrial respiration because there were no significant inhibition of gene expression (|log 2 FC| ≥1) associated with the mitochondrial respiration chain (Table 4). This result is in agreement with the previous study showed that transient accumulation of HCN in the cys-c1 mutant did not alter mitochondrial respiration rates in Arabidopsis seedlings [21]. Notably, it has been proposed that HCN is involved in the induction of the AOX gene [2,16,47], which mediates cyanide-resistant respiration and whose expression may play positive roles in plant resistance to biotic and abiotic stress [48][49][50]. However, it appears that the induction of AOX expression by HCN is an indirect effect rather than a direct effect, according to our current data from RNA-seq. Since a variety of studies have demonstrated that the main functions of AOX are to maintain the mitochondrial redox state and decrease the production of ROS [50][51][52], it is possible that the observed induction in the expression of AOX by HCN under stress conditions might be associated with some other signal molecules, such as ROS, because the transient accumulation of HCN, during an extreme ET burst in stressed tissue, can significantly block cellular respiration and induce large amounts of ROS production [2,12]. In other words, the regulatory actions of HCN on the plant mitochondrial respiration pathway, depending on its concentration, should be further confirmed by molecular and genetics methods.
In contrast to HCN treatment, ET treatment significantly induced several genes related to the mitochondrial respiratory chain, including cytochrome c oxidase and NADH dehydrogenase, indicating that ET has a positive role in the regulation of mitochondrial respiration. Additionally, it should be noted that ET treatment upregulated the gene expression of AOX (HCN resistance gene) rather than CAS (HCN detoxification gene), suggesting that AOX gene is probably regulated by ET but CAS gene is generally regulated by HCN (especially at a higher concentration). Considering that a large amount of HCN is generated along with ET biosynthesis in vivo but this is not the case when ET is applied externally, which may explain why the HCN detoxification genes (CAS family) were not significantly induced by ET treatment. Likewise, the induction of AOX by ET probably contributes to the alleviation of respiration inhibition by HCN burst during ET biosynthesis in plants.
In plants, ET has been regarded as a stress-hormone besides its roles in regulation of plant growth and development [53]. Notably, it has been speculated that sub lethal levels of HCN may trigger many events which lead to the acclimation of plants growing under adverse conditions [2], although the mechanism is still unclear. In this study, we found that a total of 117 DEGs related to stress were co-regulated by HCN and ET (Fig. 9), demonstrating that HCN and ET may synergistically regulate plant stress response. Moreover, the data showed that there were 198 stress-related DEGs were significantly induced by HCN but not by ET, thus it is necessary to further investigate who is the key component of HCN-induced plant stress resistance in the future. It is noteworthy to mention that GolS2, whose over-expression has been shown to increase plant tolerance to salt, chilling, and high-light stress [54][55][56], was upregulated by 15-fold by HCN when compared to the CK. Taji et al. (2002) stated that overexpression of GolS2 in transgenic Arabidopsis caused an increase in endogenous galactinol and raffinose, and showed reduced transpiration from leaves to improve drought tolerance [55]. The results from Sengupta et al. (2015) proposed that the stress-inducible GolS2 plays a key role in the accumulation of galactinol and raffinose under abiotic stress conditions, which may function as osmoprotectants in drought-stress tolerance of plants [57]. Similarly, Selvaraj et al. (2017) reported that over-expression of AtGolS2 was able to confer drought tolerance and increase grain yield in two different rice (Oryza sativa) genotypes under dry field conditions [58]. Given that GolS2 is involved in stress acclimation and markedly responds to HCN treatment, it appears that GolS2 is probably one of the key candidate genes contributing to HCN-mediated plant stress adaptation. Similarly, it is noteworthy that the gene of THI2.2, whose expression was markedly induced by HCN (3.7-fold) but not by ET, might be another key gene that is associated with HCN-induced plant pathogen resistance. It has been demonstrated that THI2.2 is expressed at a low basal level in seedlings and rosette leaves, encodes a PR (pathogenesis-related) protein and belongs to the PR-13 family [59,60]. However, the function of THI2.2 remains unclear although it was predicted to be involved in the defence response. Interestingly, it was shown that salicylate, ethephon, methyl jasmonate, and silver nitrate did not affect the transcript level of the THI2.2 gene [60], which is consistent with our findings that ET treatment did not induce its transcript at all. In this case, it is necessary to study whether THI2.2 is a key member of HCN-induced plant disease resistance in the future.
Recently, Garcia et al. (2019) reported that the altered immune response observed in the HCN accumulated Arabidopsis mutant (cas-c1) through posttranslational modification of proteins by S-cyanylation, which is involved in the regulation of primary metabolic pathways, such as glycolysis, and the Calvin and S-adenosylmethionine cycles [61]. Of these, a set of 163 proteins susceptible to S-cyanylation included the PEPTIDYL-PROLYL CIS-TRANS ISOMERASE 20-3 (CYP20-3) and ENOLASE 2 (ENO2). Here, we found that CYP20-3 and ENO2 were induced but not significantly by HCN, where the expression levels were approximately 1.47-fold and 1.43-fold higher than those in CK (Additional file 1: Table S8). In addition, HCN treatment regulated a large number of DEGs enriched in the cysteine and methionine metabolism pathways ( Fig. 7c; Additional file 1: Table S9). However, more experiments are required to uncover how HCN participates in the regulation of these primary carbon metabolisms, and the identification of S-cyanylated proteins should be further considered as it is beneficial for elucidating the HCN signalling mechanism [61].
In addition to the above-mentioned genes, another striking gene regulated by HCN that should be mentioned is AT2G15020, whose gene expression was upregulated by 23-fold. However, the function of the AT2G15020 is unknown currently. Therefore, further research is needed in the future to decipher the role of AT2G15020 in plants and its possible mechanism of responding to HCN.

Conclusion
In this study, we focused on the regulation of gene expression in Arabidopsis by HCN and ET treatment. The transcriptome sequencing data indicated that HCN should be recognized as an important signal molecule, rather than being simply considered a toxic by-product of ET biosynthesis. Here, we found that a large number of genes were regulated by HCN, and some of these genes were co-regulated with ET. The DEGs of CK vs HCN were associated with plant growth and development and plant response to stress. In addition, HCN-induced gene expression might be, at least partly, shared with other plant hormone signal transduction pathways. However, the cross-talk between HCN, ET and other hormones is required to perform further validation experiments with some key genes discovered in this study.
In addition, HCN and ET are small gaseous molecules with similar chemical properties and simple structures that are generated at the same time. However, there was no experimental data, including the data from the RNA-seq determination in this study, indicating that HCN shares the same receptor(s) as ET. Since HCN is a simple, small and diffusible molecule, it is highly improbable that its transduction involves specific receptor(s) [2]. Thus, the question arises as to what might be the receptor(s) of HCN and how HCN affects gene expression in plants. Genetic analysis and further physiological studies will improve our understanding of how HCN is perceived and transduced into specific downstream responses.

Plant materials and treatments
The Arabidopsis Columbia ecotype (Col-0) seeds used in this study were obtained from Arabidopsis Biological Resource Center. For growth under sterile conditions, seeds were surface sterilized with 25% (v/v) commercial bleach, followed by six washes with sterile distilled water. The seeds were sown onto half-strength Murashige and Skoog (MS)-containing 0.8% agar plates with 10 g/L sucrose. Seedlings were transferred to soil and grown in growth chambers with 16 h of light (approximately 120 μmol m − 2 s − 1 ) at 22°C, 8 h of dark at 18°C, and 70% relative humidity. For transcriptome analysis, approximately 4-weeks old seedlings were used for the following treatments.
HCN treatment was carried out according to the method of Bogatek and Lewak (1988) with some modification [62]. Arabidopsis seedlings were placed in a glassy closed container (20 L) and the HCN inside was released from another closed round-bottom flask where the HCN was produced by acidifying the 10 mL 1 mM KCN solution with 10 mL of lactic acid (10%, v/v). The HCN contents were detected with a cyanide gas detector (GT-901), which has been calibrated before use. The reaction was terminated while the final concentration of cyanide in glassy container reached 20 μM HCN, which is a relatively lower concentration and we found this concentration helps to improve the plants tolerance to environmental stress [12].
ET treatment was carried out according to the method of Khan et al. (2008) with some modification [63]. Arabidopsis seedlings were placed in a glassy closed container (20 L) and then the ET donor: 500 ppm ethephon (2-chloro-ethylphosphonic acid) was applied to generate ET gas. The control seedlings were placed in a similar size of glassy closed container (20 L). After all materials exposure to light (120 μmol m − 2 s − 1 ), 22°C for 2 h, the containers were opened and gaseous HCN and ET released. Then the HCN-treated, ET-treated and the control seedlings were collected for the total RNA extraction.

RNA extraction and Illumina sequencing
Total RNA was extracted using an EasyPure Plant RNA Kit (Transgene Co., China) and at least 5 individual seedlings of each treatment were homogenized with liquid nitrogen, and three biological replicates for each treatment were carried out in this study. The RNA quantities were analyzed by an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Equal quantities of total RNA from the three biological replicates were pooled before used for cDNA library construction and sequencing. cDNA libraries were constructed using the Illumina RNA-seq kit (Illumina, USA). Solexa adapters were then ligated to the ends of the cDNA fragments for Solexa sequencing. High-throughput sequencing was performed using Illumina Hiseq 2500 (Illumina, USA). Reads were 125 bases in length and generated from each end of the DNA fragments in paired-end sequencing.

Gene annotation and differential expression analysis
The adaptor sequences and low-quality sequence reads were removed from the data sets. Raw sequences were transformed into clean tags after data processing. Gene Ontology (GO; http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.kegg.jp) analyses were subsequently performed. GO is an internationally standardized gene function classification system for comprehensively describing the properties of genes and their products in any organism. The basic unit of GO is the GO term and each GO term belongs to a type of ontology. In gene expression profiling analysis, GO enrichment analyses of functional significance was performed using hypergeometric testing to map all differentially expressed genes to terms in the GO database by looking for GO terms that were significantly enriched in a given differentially expressed genes (DEGs) relative to the genome background. According to GO analysis, all unigenes were divided into three groups: molecular function (MF), cellular component (CC) and biological process (BP).
KEGG is the major public pathway-related database. Different genes usually cooperate with each other to exercise their biological functions. Pathway-based analysis helps the user to further understand the biological functions of specific genes. Pathway enrichment analysis identifies significantly enriched metabolic pathways and signal transduction pathways in DEGs by comparing them to the whole-genome background [64].
The expression level of each gene was measured as the normalized number of matched clean tags. The normalization method of reads per kilobase per million mapped (RPKM) was used in this study. The false discovery rate (FDR) method was used to tune the threshold P value.

Real-time quantitative PCR analysis
In order to validate the results from transcriptome sequencing analysis, part of genes were confirmed by quantitative real-time PCR (qRT-PCR). All the Primers are listed in Additional file 1: Table S1. qRT-PCR reactions were prepared with the SYBR Green Master Mix Reagent (Applied Biosystems), following the manufacturer's instruction. Reactions were carried out in Applied Real-Time System (ABI7500). The thermal cycling profile consisted of an initial denaturation at 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 30 s. All samples were performed in triplicate and relative expression levels were calculated using the delta-delta Ct method of the system. In this study, ACTIN2 (AT3G18780) gene was used as internal control.

Additional file
Additional file 1: Table S1. All primers for qRT-PCR. Figure S1. qRT-PCR analysis of DEGs from CK vs HCN and CK vs ET. (a) The top 5 genes co-regulated by HCN and ET were determined by qRT-PCR. The top 5 genes exclusively regulated by HCN (b) or ET (c) were determined by qRT-PCR. Expression ratios (FPKM fold change) obtained from transcriptome data (green) and qRT-PCR (red). (d) Lineage analysis between the transcriptome and qRT-PCR. Figure Table S3. The coregulated DEGs related to plant hormone signaling transduction pathway by HCN and ET. Table S4. DEGs related to Auxin/IAAs that were regulated by HCN and ET. Table S5. DEGs related to ROS production in CK vs HCN. Table S6. Top 10 co-upregulated DEGs related to stress by HCN and ET. Table S7. Top 10 DEGs related to stress that were exclusively induced by HCN or ET. Table S8. Parts of DEGs related to post-translation regulation by S-cyanylation in CK vs HCN. Table S9. DEGs enriched in cysteine and methionine metabolism that were regulated by HCN.