Identification of histone deacetylase genes in Dendrobium officinale and their expression profiles under phytohormone and abiotic stress treatments

The deacetylation of core histones controlled by the action of histone deacetylases (HDACs) plays an important role in the epigenetic regulation of plant gene transcription. However, no systematic analysis of HDAC genes in Dendrobium officinale, a medicinal orchid, has been performed. In the current study, a total of 14 histone deacetylases in D. officinale were identified and characterized using bioinformatics-based methods. These genes were classified into RPD3/HDA1, SIR2, and HD2 subfamilies. Most DoHDAC genes in the same subfamily shared similar structures, and their encoded proteins contained similar motifs, suggesting that the HDAC family members are highly conserved and might have similar functions. Different cis-acting elements in promoters were related to abiotic stresses and exogenous plant hormones. A transient expression assay in onion epidermal cells by Agrobacterium-mediated transformation indicated that all of the detected histone deacetylases such as DoHDA7, DoHDA9, DoHDA10, DoHDT3, DoHDT4, DoSRT1 and DoSRT2, were localized in the nucleus. A tissue-specific analysis based on RNA-seq suggested that DoHDAC genes play a role in growth and development in D. officinale. The expression profiles of selected DoHDAC genes under abiotic stresses and plant hormone treatments were analyzed by qRT-PCR. DoHDA3, DoHDA8, DoHDA10 and DoHDT4 were modulated by multiple abiotic stresses and phytohormones, indicating that these genes were involved in abiotic stress response and phytohormone signaling pathways. These results provide valuable information for molecular studies to further elucidate the function of DoHDAC genes.


INTRODUCTION
Higher plants have a sessile lifestyle and often suffer from exposure to abiotic stresses such as drought, cold or high salinity (Shao, Wang & Tang, 2015). Various stress-inducing genes and signaling factors with different functions participate in stress responses, and the expression of these genes usually depends on chromatin remodeling (Kim et al., 2010).
The eukaryotic genome is stored in a nuclear protein complex composed of nucleosomes (Teif & Clarkson, 2019). In the nucleosome, two copies of histone H2A, H2B, H3, and H4 form a histone octamer, which is surrounded by DNA 147 base pairs long (Luger et al., 1997). The basic residues such as lysine and arginine in the N-terminal tail of each histone are covalently modified by acetylation to regulate the transcription of genes wrapped around the core histone (Srivastava, Singh & Dubey, 2016). Since a report that indicated that gene activity is related to histone acetylation, the modification of histone acetylation and deacetylation have received increasing attention (Grunstein, 1997). Generally, the acetylation of lysine residues in H3 and H4 N-tails neutralizes the positive charge of the histone tails, thereby reducing their affinity for DNA strands, hence histone acetylation mediated by histone acetyltransferases (HATs), which connect acetyl groups to histones, promotes loosening of the chromatin structure and activates gene transcription (Berger, 2007;Kim et al., 2015;Shahbazian & Grunstein, 2007). Conversely, histone deacetylation mediated by histone deacetylases (HDACs), which remove an acetyl group from histones, facilitates the formation of compact chromatin leading to the repression or silencing of genes (Hollender & Liu, 2008).
Dendrobium Swartz is the second largest genus of the Orchidaceae, which consists of nearly 1,600 accepted species (Juswara, Schuiteman & Champion, 2019). Among them, Dendrobium officinale Kimura et Migo is a popular tonic and traditional Chinese medicine (TCM) with high commercial value, containing compounds with antioxidant and antitumor activity (Li et al., 2011;Teixeira da Silva & Ng, 2017;Wu et al., 2014). In addition to polysaccharides thought to be among the major active ingredients of D. officinale, a wide variety of low molecular weight compounds including bibenzyls, phenanthrenes, and flavanones have been detected by phytochemical analyses (Chen et al., 2012). D. officinale grows well under a suitable environment (Ding et al., 2018;Yang et al., 2014). In a natural environment, D. officinale is in an epiphytic state, often growing on humid rocks in mountain climates at an altitude of 500-1600 m, or on the tree trunks in virgin forests with warm and moist environments (Zeng et al., 2019). Based on its growth habits, D. officinale is susceptible to abiotic stresses such as high and low temperatures, drought and salinization, resulting in low natural reproduction rate and slow growth (Zeng et al., 2019).
Phytohormones not only mediate developmental processes in plants but can also sustainably alleviate adverse effects of abiotic and biotic stresses in plants (Wolters & Jurgens, 2009). For instance, SA improved plant abiotic stress tolerance via SA-mediated control of major plant metabolic processes by strengthening salinity and drought stress tolerance (Khan et al., 2015). The JA signaling pathway plays a key role in coordinating the activation of biosynthetic routes involved in defense-related metabolites (Ballare & Austin, 2019). Secondary metabolites, including terpenes, phenolics, compounds with nitrogen, and others, which are biosynthetically derived from primary metabolites, have also been established as beneficial for plants' tolerance to biotic and abiotic stresses (Ramakrishna & Ravishankar, 2011). Environmental stresses trigger the biosynthesis of certain secondary metabolites in plants which are important candidates for human nutrition (Ashraf et al., 2018). Furthermore, epigenetic regulation of gene expression is important for a plant's adaptation to environmental changes (Zheng et al., 2020). Therefore, understanding how histone acetylation modification is involved in the responses of D. officinale to environmental stresses will contribute significantly to our understanding of the molecular mechanisms underlying epigenetic regulation in this orchid. However, studies on the evolutionary relationships and functional characteristics of HDAC genes in D. officinale (DoHDAC genes) are scarce.
The genome and transcriptome sequences of D. officinale are now available, allowing the DoHDAC genes to be isolated and identified (Shen et al., 2017). In this study, fourteen DoHDAC genes were identified from the D. officinale genome and their structure, phylogeny, conserved motifs, and putative promoter were analyzed. Subsequently, the subcellular localization and expression patterns of selected DoHDAC genes under hormone, salt, drought and cold treatments were analyzed. These results provide a foundation for further clarifying the functions of DoHDAC proteins.

Identification and bioinformatics analysis of DoHDAC genes
To identify potential DoHDAC genes, the coding sequences (CDSs) of AtHDAC genes were used for a BLASTP algorithm-based query against the D. officinale genome database (https://www.ncbi.nlm.nih.gov/genome/?term=Dendrobium, NCBI Biosample ID: SAMN03083363). Genes were identified by a hidden Markov model search based on the HDAC domain using the Pfam protein domain database. To determine the sequence identity among HDACs, full-length amino acid sequences of HDAC proteins were aligned and compared using DNAMAN 8 software (Lynnon Biosoft, San Ramon, CA, USA). For phylogenetic analysis, the amino acid sequences of HDAC proteins from D. officinale, A. thaliana and O. sativa in FASTA format were aligned with Clustal X 2.0 (Larkin et al., 2007), and sequence identity was calculated by UniProt BLAST (http://www.uniprot.org/blast/), and then a neighbor-joining phylogenetic tree was constructed using the MEGA 7 program (Kumar, Stecher & Tamura, 2016) with a bootstrap analysis of 1000 replicates. Nucleus localization signals (NLSs) were found in amino acid sequences of DoHDAC proteins at http://localizer.csiro.au/ (Sperschneider et al., 2017). Recognizable conserved domains of DoHDAC proteins were identified with EBL-EBI (http://www.ebi.ac.uk/interpro/search/sequence-search) and also verified by the CDD database (https://www.ncbi.nlm.nih.gov/cdd). The domain architecture was drawn using DOG2.0 software (Ren et al., 2009). The protein sequence motifs of DoHDAC proteins were identified using MEME (http://memesuite.org/tools/meme) (Bailey et al., 2009). Gene structure analysis of the DoHDAC genes was conducted with GSDS (Gene Structure Display Server, http://gsds.gao-lab.org/) (Hu et al., 2014). The functional interacting networks of functional proteins were integrated in STRING (https://string-db.org/) based on an A. thaliana association model with the confidence parameter set at a threshold of 0.700 and no more than 20 interactors (Szklarczyk et al., 2019). DNA sequences about 2000 bp long upstream of the initiation codon ATG of DoHDAC genes were regarded as putative promoters. Cis-elements in these promoters were analyzed using the online program PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The predicted cis-acting elements of promoters of DoHDAC genes were illustrated by TBtools software (Chen et al., 2020).

Plant materials and treatments
Based on previously reported concentrations (He et al., 2017), D. officinale plantlets with a stem length of about three cm growing on half-strength Murashige and Skoog (MS; (Murashige & Skoog, 1962) (1/2 MS) agar medium containing 2% sucrose were separately exposed to one of several abiotic stresses: 100 µM gibberellic acid (GA 3 ), 100 µM ABA, 100 µM SA, 200 µM methyl jasmonate (MeJA), 15% polyethylene glycol (PEG)-6000 (all compounds: Sigma-Aldrich, Shanghai, China), 250 mM NaCl (Guangzhou Chemical Reagent Factory, Guangzhou, China) and 4 • C , separately. Five D. officinale plantlets were used for each treatment. Plant tissues, including roots, stems and leaves, were harvested at different selected time points after applying the above treatments, then immediately frozen in liquid nitrogen for further gene expression experiments. Based on the transcriptome data, the expression patterns of DoHDAC genes in different tissues were analyzed, and a heatmap was generated with the heatmap software package on the BMKCloud platform (http://www.biocloud.net).
About 100 mg of samples, including leaves, roots and stems from previously collected D. officinale plants, were placed in a pre-chilled mortar and ground with liquid nitrogen to a fine powder. Each finely powdered sample was transferred to a 2 mL microcentrifuge tube containing 1.5 mL of pre-warmed extraction buffer (50 mM Tris-HCl, 20 mM EDTA, 1 M NaCl, 2% SDS and 4% β-mercaptoethanol; pH 8.0). Samples were incubated in a water bath at 65 • C for 15-20 min. The tube was centrifuged at 13,000 g for 10 min at 4 • C and the supernatant was transferred to a sterile 2 mL centrifuge tube. An equal volume of water-saturated phenol (pH 4.5) was added, vortexed for 2 min, then kept at room temperature for 10 min. The tube was centrifuged at 13,000 g for 10 min at 4 • C, the supernatant was transferred to a sterile 2 mL centrifuge tube, then equal volumes of chloroform and isoamyl alcohol (24/1, v/v) were added and mixed thoroughly. The tube was again centrifuged at 13,000 g for 10 min at 4 • C, the supernatant was transferred to a sterile two mL centrifuge tube, and an equal volume of isopropanol was added. The solution in the centrifuge tube was mixed, placed in a refrigerator at −20 • C for 10 min, and then centrifuged at 13,000 g for 10 min at 4 • C. The supernatant was discarded, the pellet was washed with 75% chilled ethanol, and the tube was centrifuged at 13,000 g for 15 min at 4 • C. The RNA pellet was dried at room temperature for 5 min. The dried RNA pellet was resuspended in 50 µL of diethyl phosphorocyanidate (DEPC) in distilled water. According to the manufacturer's instructions, samples were treated with Recombinant DNase I (RNase-free) (Takara Bio Inc., Kusatsu, Shiga, Japan) to remove genomic DNA, and stored at −80 • C. The NanoDrop One/OneC micro nucleic acid protein concentration analyzer (Thermo Fisher Scientific, Waltham, MA, USA) was used to assess the purity and concentration of RNA samples. Agarose gel electrophoresis was used to survey the RNA integrity by Clinx GenoSens gel documentation system (Clinx Science Instruments, Shanghai, China).

Quantitative real-time reverse transcriptase-PCR analysis
Following the manufacturer's instructions, the first-strand cDNA was synthesized using the GoScript TM Reverse Transcription System (Promega, Madison, WI, USA). Quantitative real-time reverse transcriptase-PCR (qRT-PCR) assays of three biological replicate samples were performed on a LightCycler 480 system (Roche, Basel, Switzerland) using the Aptamer TM qPCR SYBR R Green Master Mix (Tianjin Novogene Bioinformatics Technology Co. Ltd., Tianjin, China). The reaction conditions were as follows: 95 • C for 5 min, followed by 40 cycles of 95 • C for 10 s, and 60 • C for 30 s. D. officinale actin (NCBI accession number JX294908) was used as an internal control gene based on the advice of He et al. (2015). The relative expression levels of genes were calculated using the 2 − CT method (Livak & Schmittgen, 2001). The expression of genes was significantly altered upon different treatments when the fold change was greater than 2 (Li et al., 2016). The gene-specific primers for qRT-PCR were designed by Primerquest (https://sg.idtdna.com/Primerquest/Home/Index) and are listed in Table S1.

Identification and phylogenetic analysis of the DoHDAC genes
In this study, 14 DoHDAC genes were identified in the D. officinale genome. The CDS length of DoHDAC genes was diverse, ranging from 267 bp to 2085 bp ( Table 1). Analysis of gene structure helps to understand gene function and evolution (Feng et al., 2016). Hence, the structure of the DoHDAC genes was analyzed. Our results showed that their conserved coding regions consisted of different numbers of exons, while DoHDA10 and DoHDT4 have only one exon (Fig. 1). Amino acid sequence alignment of the 14 DoHDAC proteins showed that, except for DoHDA4 and DoHDA10, the eight histone deacetylases with amino acid sequences between 200 aa and 312 aa all contained a typical histone deacetylase catalytic domain. DoHDA4 and DoHDA10 showed significantly low sequence identity with other RPD3/HDA1 family proteins in D. officinale (Fig. 2), suggesting that they might play some unique roles in certain cellular events. Functional domain analysis revealed that DoHDA6 possessed a RanBP2-type zinc finger domain (Fig. 2). DoHDT3 contained a conserved domain with a pentapeptide (MEFWG) on the N-terminus and a C 2 H 2 zinc finger domain on the C-terminus of the protein (Fig. S1, Fig. 2). In addition, DoSRT1 and DoSRT2 possessed a single copy of the sirtuin-type HDAC domain (Fig. 2). To investigate the conserved nature of motifs in DoHDAC proteins, we investigated 10 motifs among the 14 DoHDAC proteins. The width of these 10 motifs ranged from 21 (motif 6, 8 and 9) to 50 (motif 3, 4 and 5) amino acids (Fig. 3). Closely related proteins, such as DoHDA3, DoHDA7, DoHDA9, DoHDA1, DoHDA6, DoHDA2, DoHDA5, DoHDT3 and DoHDT4, had similar conserved motifs. For instance, most of the D. officinale RPD3/HDA1 superfamily proteins contained HDACs domain motifs 1 and 2. However, the HDAC domain motif was not found in DoSRT1 and DoSRT2. We constructed a protein-protein interaction network of the 14 DoHDAC proteins using STRING (http://string-db.org/) based on an A. thaliana association model (Fig. 4). DoHDA9 showed a high degree of amino acid sequence similarity to AtHDA6, thus it might interact with the flavin-containing amine oxidoreductase family protein (FLD), which might be histone demethylase that promotes flowering by repressing FLOWERING LOCUS C (FLC) ( Table S3, Fig. 4, Table S4). Moreover, DoHDA9 might interact with the WD40 repeat-containing protein HOS15 in response to abiotic stress (Table S4). To study the evolutionary relationships among DoHDAC proteins, a phylogenetic tree was constructed by aligning the amino acid sequences of histone deacetylases of four, ten and one HDAC proteins from A. thaliana, O. sativa, and Z. mays using MEGA 7.0. In the phylogenetic tree, 14 DoHDAC proteins were divided into three families: RPD3/HDA1, SIR2, and HD2 (Fig. 5). DoHDA9, AtHDA6 homologue, shared approximately 71.1% amino acid sequence similarity with AtHDA6. In addition, DoHDA7, AtHDA19 homologue, had approximately 80% amino acid sequence similarity with AtHDA19.

Subcellular localization analysis of selected DoHDAC genes
To verify the subcellular localization of certain DoHDAC proteins, an Agrobacteriummediated transient transformation system was used in onion epidermal cells. The expression of GFP fused to DoHDA7, DoHDA9, DoHDA10, DoHDT3, DoHDT4, DoSRT1 and DoSRT2 was tracked by the GFP marker signal. The blue fluorescence pattern of DAPI-stained cells completely merged with the green fluorescence pattern, showing the nuclear localization of DoHDAC proteins in onion epidermal cells (Fig. 6). As shown in Fig. 6, the positive control was distributed throughout the onion cells. In contrast, DoHDA7, DoHDA9, DoHDA10, DoHDT3, DoHDT4, DoSRT1 and DoSRT2 were localized in the nucleus.

Analysis of cis-acting elements in the putative promoters of DoHDAC genes
A cis-acting element is a non-coding part of DNA, is usually limited to the 5 upstream region of a gene, and is responsible for transcriptional regulation. To better predict gene functions, we identified the cis-acting elements within the putative ∼2 kb promoter region of DoHDAC genes. In addition to the core cis-acting elements, such as the TATA-box and the CAAT-box, several other cis-elements related to plant growth or development, stress responses, and hormone responses were detected (Table S2, Fig. 7). Cis-acting elements related to growth or development included light-responsive elements such as Box4, G-box, GT1-motif, Sp1, and TCT-motif, elements expressed only in the endosperm such as the AACA_motif and GCN4_motif, and elements under circadian control or showing meristem-specific activation such as CAT-box and NON-box. Cis-acting elements related to phytohormone responses included GA-responsive elements such as GARE-motif, P-box and TATC-box, MeJA-responsive elements such as CGTCA-motif and TGACG-motif, ABA-responsive element (ABRE), and SA-responsive element such as TCA-element. Cis-acting elements related to abiotic stresses included drought-responsive elements such as MBS and a low-temperature-responsive (LTR) element. For example, both DoHDA3 and DoHDT3 had cis-acting elements involved in MeJA-responsiveness, low-temperature responsiveness and drought inducibility. Therefore, based on the large number of cis-acting elements related to these stresses, the D. officinale plantlets were treated with SA, MeJA, ABA, GA 3 , cold stress and drought stress.

Expression analysis of DoHDAC genes in different tissues
Since high-throughput sequencing has been performed on D. officinale tissues at different developmental stages, the publicly available RNA-seq data obtained is a useful resource for studying gene expression profiles. We used previously reported transcriptome data (Zhang et al., 2017) to examine the expression pattern of DoHDAC genes in different tissues. The heat map (Fig. 8) revealed that three DoHDAC genes (DoHDA3, DoHDA7, and DoHDT3) were highly expressed with FPKM values > 25 in all detected tissues, while DoHDA10 showed a low expression with FPKM values < 5 in all the detected tissues. In addition, some genes displayed a tissue-specific expression pattern. For example, DoHDA5 and DoSRT1 were most highly expressed in the white parts of root, while DoHDA6 and DoHDT4 were expressed predominantly in flower buds (Fig. 8). These results suggested that DoHDAC genes might play an important role in the development of different organs.

Expression analysis of DoHDAC genes in response to exogenous phytohormones
Hormones play a vital role in the vegetative and reproductive growth of plants. Analysis of the cis-acting elements in promoters indicates that the DoHDAC genes might respond to a variety of stresses and signaling molecules. Therefore, based on potential cis-acting  elements, we selected some genes and used qRT-PCR to detect their expression patterns under various phytohormone treatments. We analyzed the expression patterns of 12 DoHDAC genes in different subfamilies after GA 3 treatment (Fig. 9). The results showed that after 4 h of GA 3 treatment, 11 genes in roots were highly up-regulated, DoHDA8 and DoHDA10 were the highest (> 20-fold), and DoSRT2 was highly up-regulated in stems. In addition, after 24 h of GA 3 treatment, DoHDA10 and DoHDT4 were highly up-regulated in leaves. At the same time, cis-acting elements involved in the GA 3 -response were found in  these four genes (Table S2). Most selected DoHDAC genes were induced by exogenous SA treatment (Fig. 10). All genes in roots were highly up-regulated after 2 h under SA treatment, while the expression of DoHDA10 and DoHDT4 in leaves was highly up-regulated after 24 h of SA treatment. Moreover, cis-acting elements in response to SA were found in these two genes (Table S2). ABA is a widely distributed plant hormone in land plants and plays an important role in plant development (Finkelstein, Gampala & Rock, 2002). After 4 h of ABA treatment, 11 DoHDAC genes, especially DoHDA8, DoHDA10 and DoHDT4, were up-regulated in roots (Fig. 11), and cis-acting elements in response to ABA were found in the promoters of DoHDA8 and DoHDT4 (Table S2). After 24 h of treatment with MeJA, 12 DoHDAC genes in roots were up-regulated and two DoHDAC genes were down-regulated (Fig. 12). Interestingly, we found that after 2 h of MeJA treatment, DoHDA10 and DoHDT4 were significantly up-regulated in stems and leaves (> 9-fold), and cis-acting elements in response to MeJA were found in the promoters of these two genes (Table S2). These results indicate that DoHDT4 might be involved in GA 3 -, SA-, ABA-and MeJA-mediated signal transduction pathways, DoHDA10 might be involved in GA 3 -, SA-and MeJA-mediated signal transduction pathways, and DoHDA8 might be involved in ABA-mediated signal transduction pathways.

Expression analysis of DoHDAC genes in response to abiotic stresses
HDAC proteins appear to activate positive and negative responses in salinity stress, with some discrepancies. For example, AtHDA9 and AtHD2D negatively regulate the salt response, whereas AtHDA6 and AtHD2C positively regulate it (Ueda et al., 2017). After salt (NaCl) treatment, most of the selected DoHDAC genes were highly up-regulated in leaves, especially DoHDA10 and DoHDT4, which showed similar expression patterns (Fig. 13). After PEG treatment, most genes were highly up-regulated in roots, stems and leaves, especially DoHDA10 and DoHDT4, which showed similar expression patterns (Fig. 14). In addition, cis-acting elements involved in responsiveness to drought were found in the promoters of DoHDA10 and DoHDT4 (Table S2). The drought stress response in plants is related to the status of histone acetylation. In response to drought stress, the level of acetylation of histone H3K9 increased in drought-responsive genes (Kim  DoHDA2,DoHDA5,DoHDA6,DoHDA7,DoHDA8,DoHDA9,DoHDA10,DoHDT4,DoSRT1,DoSRT2, respectively. Error bars represent the standard error of the means of three independent replicates. Expression levels were calculated relative to 0 h of each organ. Full-size DOI: 10.7717/peerj.10482/ fig-9 et al., 2008). Drought stress significantly induced four HAT genes in rice (OsHAC703, OsHAG703, OsHAF701, and OsHAM701) (Fang et al., 2014). However, in this study, the up-regulated expression of DoHDAC genes under drought stress may be related to the inhibited expression of drought-insensitive genes. AtHD2C-overexpressing A. thaliana plants showed ABA insensitivity, reduced transpiration, and enhanced tolerance to drought stress (Sridha & Wu, 2006). After low-temperature treatment, the selected DoHDAC genes were weakly up-regulated in roots, stems and leaves (Fig. 15). This is consistent with the expression of most HDAC genes in rice being regulated by salt and drought, but less affected by cold (Hu et al., 2009). Therefore, DoHDA10 and DoHDT4 might play an important role  DoHDA4,DoHDA5,DoHDA6,DoHDA7,DoHDA8,DoHDA9,DoHDA10,DoHDT4,DoSRT1,DoSRT2,respectively. Error bars represent the standard error of the means of three independent replicates. Expression levels were calculated relative to 0 h of each organ. Full-size DOI: 10.7717/peerj.10482/ fig-11 and plant responses to abiotic stress (Liu et al., 2014). Plant HDAC proteins have been identified and characterized in some plant species such as A. thaliana, rice, soybean, tomato and others. For instance, a total of 18 AtHDAC proteins have been identified, 12 of which belong to the RPD3/HDA1 family, four belong to the HD2 family, and two belong to the SIR2 family (Pandey et al., 2002;Wang et al., 2014). Furthermore, at least 18 HDAC proteins have been identified in the rice genome, 14 of which belong to the RPD3/HDA1 family, two belong to the HD2 family, and two belong to the SIR2 family (Hu et al., 2009). In this study, we identified 14 HDAC proteins in D. officinale by using bioinformatics analysis. However, the number of members of the HDAC proteins in D. officinale identified was less than that in A. thaliana and rice. The other DoHDAC proteins need to be identified and characterized in a future study. Histone post-translational modifications (PTMs) are a major regulator of gene transcription. A number of PTMs are located within interfaces between H3/H4 tetramers and/or H2A/H2B dimers, such as the acetylation of H4(K91) (Zhang et al., 2003). Shortly after synthesis, histones H3 and H4 are acetylated by histone acetyltransferase 1/2 in the cytoplasm and then imported into the nucleus together with the histone acetyltransferase 1 complex. In the nucleus, H3-H4 dimers or tetramers are recognized by chaperones such as CAF-1, Asf1 and Hif1, perhaps through their acetylation markers and are deposited together with histone H2A-H2B dimers onto newly replicated DNA to form nucleosomes, and then level of DoHDA7 was up-regulated in stems and leaves (Fig. 13). The expression of HDAC genes occurs in response to stresses and is regulated by stress-related hormones like SA, JA or ABA (Fu, Wu & Duan, 2007;Hu et al., 2009). ABA is the plant hormone most directly involved in stress signal transduction (Xiong, Schumaker & Zhu, 2002). JA generally inhibits plant growth and promotes defense responses against insect pests and pathogens while ABA is involved in the water stress response, regulating plant water balance and osmotic stress tolerance (Ma et al., 2013). After ABA treatment, the expression levels of AtHD2C, which is a synonym of AtHDT3, OsHDT701, OsHTD702, OsSRT701 and OsSRT702 were down-regulated (Fu, Wu & Duan, 2007;Sridha & Wu, 2006). In contrast, the expression levels of DoHDA10, DoHDT4, and DoSRT1 were significantly up-regulated in roots after ABA treatment and in leaves after NaCl treatment (Figs. 11 & 13), so they may affect the expression level of some stress-responsive genes. In addition to their functions in plant development and stress responses, HDAC proteins are also involved in cellular processes such as cell death. OsSRT1 in rice and NtHD2a and NtHD2b in tobacco act as negative regulators of cell death (Bourque et al., 2011;Huang et al., 2007). Programmed cell death (PCD) in the maize aleurone layer is induced by GA, and HDAC activity is required for this process. HDAC activity gradually increased relative to histone acetyltransferase (HAT) activity, leading to a global decrease in histone H3 and H4 acetylation levels during PCD after treatment with GA (Hou et al., 2017). These data further prove that the DoHDAC genes play an important role in the growth and development of D. officinale and in response to environmental stresses. Further research is required to identify proteins that interact with DoHDAC proteins and their global targets in D. officinale.

CONCLUSIONS
We comprehensively identified and analyzed the HDAC genes from D. officinale. We performed a phylogenetic analysis, as well as analyses of conserved motifs, cis-acting elements, gene structure, protein interactions, and subcellular localization of some DoHDAC genes. These results help to elucidate the classification and function of the DoHDAC genes. Expression pattern analysis in different tissues showed that the DoHDAC genes were widely expressed in roots, stems, leaves and flower buds. The expression pattern of DoHDAC genes under abiotic stresses and various phytohormone treatments show that some DoHDAC genes were modulated by abiotic stresses such as salt, cold and drought, as well as plant hormones such as GA 3 , SA, ABA and MeJA. In summary, this information about the HDAC genes in D. officinale helps to reveal their role in epigenetic regulation of plant growth, development, and response to stresses.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by grants from the Natural National Science Foundation of China (No. 31871547) and the Science and Technology Project Foundation of Guangzhou