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

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).

Subcellular localization assays

The CDSs of DoHDA7, DoHDA9, DoHDA10, DoHDT3, DoHDT4, DoSRT1 and DoSRT2 without a termination codon were cloned into binary vector pCAMBIA1302 with the green fluorescent protein (GFP) (GenBank Accession No. AF234298) with the Nco I restriction enzyme. The resulting binary vectors, pCAMBIA1302-DoHDA7, pCAMBIA1302-DoHDA9, pCAMBIA1302-DoHDA10, pCAMBIA1302-DoHDT3, pCAMBIA1302-DoHDT4, pCAMBIA1302-DoSRT1 and pCAMBIA1302-DoSRT2, were transferred into Agrobacterium tumefaciens strain EHA105 (Shanghai Weidi Biotechnology Co. Ltd., Shanghai, China), according to the manufacturer’s protocol. Then, using Agrobacterium-mediated transformation, these recombinant binary vectors were transformed into epidermal cells of onion (Allium cepa L. “Red Sun”) (Huang et al., 2011). A Zeiss Model Axio Imager A2 upright fluorescence microscope (Carl Zeiss, Oberkochen, Germany) was used to observe gene expression in transformed onion epidermal cells. In order to show the nucleus, according to the manufacturer’s protocol, the epidermal cell layers were stained with 2-(4-amiylphenyl)-6-indolecarbamoyl dihydrochloride (DAPI) staining solution (Beyotime Biotech Inc., Shanghai, China). The positive control was onion epidermal cells with an empty pCAMBIA1302-GFP vector. ZEN 2011 software (Carl Zeiss) was used to analyze cellular localization.

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₃), 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™ 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™ qPCR SYBR^® 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₂H₂ 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 Agrobacterium-mediated 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₃, 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₃ treatment (Fig. 9). The results showed that after 4 h of GA₃ 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₃ treatment, DoHDA10 and DoHDT4 were highly up-regulated in leaves. At the same time, cis-acting elements involved in the GA₃-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₃-, SA-, ABA- and MeJA-mediated signal transduction pathways, DoHDA10 might be involved in GA₃-, 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 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 in the regulation of salt and drought stress responses during the growth of D. officinale plants.