Comprehensive Analysis of the Histone Deacetylase Gene Family in Chinese Cabbage ( Brassica rapa ): From Evolution and Expression Pattern to Functional Analysis of BraHDA3

: Histone deacetylases (HDACs) are known as erasers that remove acetyl groups from lysine residues in histones. Although plant HDACs play essential roles in physiological processes, including various stress responses, our knowledge concerning HDAC gene families and their evolutionary relationship remains limited. In Brassica rapa genome, we identiﬁed 20 HDAC genes, which are divided into three major groups: RPD3/HDA1, HD2, and SIR2 families. In addition, seven pairs of segmental duplicated paralogs and one pair of tandem duplicated paralogs were identiﬁed in the B . rapa HDAC (BraHDAC) family, indicating that segmental duplication is predominant for the expansion of the BraHDAC genes. The expression patterns of paralogous gene pairs suggest a divergence in the function of BraHDACs under various stress conditions. Furthermore, we suggested that BraHDA3 (homologous of Arabidopsis HDA14 ) encodes the functional HDAC enzyme, which can be inhibited by Class I/II HDAC inhibitor SAHA. As a ﬁrst step toward understanding the epigenetic responses to environmental stresses in Chinese cabbage, our results provide a solid foundation for functional analysis of the BraHDAC family. rapa remains relatively unknown. In this study, we identiﬁed 20 putative BraHDAC genes divided into three major classes: RPD3/HDA1, HD2, and SIR2 families. Most HDAC gene duplications in B. rapa appeared to have been caused by segmental duplication, which occurred between 12.73 MYA and 42.04 MYA. Several BraHDACs were up- or down-regulated by abiotic stress treatments, indicating that the variation of BraHDAC transcription levels may contribute to alteration in histone H3 acetylation levels in response to various stresses. Our results provide a solid foundation for further understanding of the underlying evolutionary mechanisms in the HDAC family and


Introduction
Epigenetic regulation, including DNA and histone modifications, plays an essential role in controlling chromatin structure, with the nucleosome as the basic organization unit [1]. In eukaryotes, post-transcriptional histone modifications such as acetylation, methylation, ubiquitination, and sumoylation enable the regulation of mRNA expression by recruiting histone writers, erasers, and readers [2]. Among these modifications, histone acetylation in different lysine residues of histones by histone acetyltransferases (HATs) is required for opening the chromatin structure to allow binding of RNA polymerase II and transcription factors; in contrast, histone deacetylases (HDACs) maintain the homeotic balance of histone acetylation by removing the acetyl groups from hyperacetylated histones [1,3]. The highly dynamic balance of the histone acetylation and deacetylation is mostly maintained by the physical and functional interplay between HAT and HDAC activities [3], and serves as a major epigenetic regulatory mechanism for gene transcription and in turn governs numerous physiological processes [4].
In higher plants, HDACs are generally classified into three distinct groups: the reduced potassium dependence 3/histone deacetylase 1 (RPD3/HDA1) family, silent information regulator 2 (SIR2) family, and histone deacetylase 2 (HD2) family [5]. Since 18 HDACs were identified in the Arabidopsis thaliana genome [6], several HDACs have been identified and isolated from rice [7], soybean [8], grape [9], litchi [10], maize [11], and cotton [12]. These indicate that genome-wide mining and comparative analysis are powerful approach to identify genes and understand the evolutionary relationships of genes. Genetic and physiological studies on Arabidopsis have shown that HDA6, HDA9, HDA15, HDA19, HD2C, and AtSRT1 negatively regulate stress tolerance by regulating gene expression through histone acetylation [13]. The Arabidopsis HDA6 mutants showed enhanced induction of genes involved in the acetate biosynthesis pathway, which resulted in enhanced drought tolerance [14]. Similarly, the loss-of-function of HDA9 showed increased tolerance against salt and drought stresses [15]. In addition, HDA19-deficiency enhanced tolerance to highsalinity, drought, and heat stresses in Arabidopsis [16]. In contrast, HDA19 overexpression in Arabidopsis enhanced disease resistance against Alternaria brassicicola via increased transcription of jasmonic acid and ethylene responsive genes [17]. In rice plants, overexpression of OsHDA705, a member of the RPD3/HDA1 family, resulted in decreased salt and abscisic acid (ABA) stress tolerance during seed germination, whereas its overexpression resulted in enhanced osmotic stress resistance during the seedling stage [18]. The decrease in levels of histone H4 acetylation because of the overexpression of OsHDT701 (rice HD2 family) caused reduced resistance against rice pathogens Magnaporthe oryzae and Xanthomonas oryzae pv oryzae [7], whereas the overexpression of HDT701 in transgenic rice leads to increased tolerance to NaCl and polyethylene glycol stresses [19]. Furthermore, mitogenactivated protein kinase 3 phosphorylates the histone deacetylase HD2B to control the transcription of biotic stress response genes [20]. Collectively, these indicate the importance, diversity, and specificity of HDACs in response to different stimuli.
Despite the knowledge concerning HDACs, the evolutionary relationships and functional information regarding the HDAC family in Brassica rapa, a crop species of economic importance, remain largely unknown. In this study, we aimed at genome-wide identification of the HDAC gene family in B. rapa genome. Expansion and evolutionary mechanisms of this gene family were studied to explore the evolutionary relationships among B. rapa HDAC family (BraHDAC). Furthermore, we analyzed the expression pattern of BraHDACs under different abiotic stress condition, providing valuable information for further study of BraHDAC gene functions.

Identification and Characterization of HDAC Family in B. rapa
Phytozome v12.1 (https://phytozome.jgi.doe.gov; Brassica rapa FPsc v1.3) was used to identify HDACs in B. rapa genome using BLASTp with Arabidopsis HDACs as queries. The conserved domains, molecular weight, phylogenetic, and isoelectric point (pI) analyses and subcellular localization of putative BraHDACs has been conducted by Eom and Hyun [21]. A nomenclature system for BraHDACs (generic name: BraHDA1 to 16 for the RPD3/HDA1 family, BraHDT1 and 2 for the HD2 family, and BraSRT1 and 2 for the SIR2 family) was used as described by Hollender and Liu [6].

Plant Growth and Treatment Conditions
Chinese cabbage (cultivar Chunkwang) seeds were grown in a growth chamber under long-day conditions at 22 • C and 60% relative humidity. After 4 weeks, Chinese cabbage plants were treated at 40 • C, and samples were harvested at the times indicated in Figure 2.
Heat treatment was carried out three times, resulting in three biological replicates.

The Analysis of BraHDAC Expression
Total RNA was extracted, quantified, and qualified using nanodrop. cDNA was synthesized using ReverTra Ace qPCR RT Master Mix with gDNA Remover. Quantitative RT-PCR (qRT-PCR) was performed using the SYBR Green Real-time PCR Master Mix. The expression levels of each gene were normalized to the internal reference gene actin, and then was expressed relative to its value at 0 h. The specific primer pairs were designed using GenScript Real-time PCR Primer Design tool, and are listed in Table S1.

Production of Recombinant BraHDA3 in E. coli, and the Analysis of HDAC Activity
The PCR product of BraHDA3 was ligated into the pCR-Blunt, and subcloned into the Escherichia coli expression vector pPROEX-Htb vector that produces a histidine-fused translation product. To produce his tagged-BraHDA3 fusion proteins, E. coli strain BL21 harboring pPROEX-BraHDA3 plasmid was incubated at 37 • C. Once the optical density at 600 nm reached 0.5, fusion proteins were produced by 0.1 mM isopropyl βd-1thiogalactopyranoside (IPTG) treatment for 16 h at 30 • C. The fusion protein was purified using Ni-NTA affinity purification, and the HDAC activity was determined using HDAC Activity Colorimetric Assay Kit (BioVision, Milpitas, CA, USA). HDAC activity was calculated based on a standard curve generated with Ac-Lys-pNA (deacetylated standard).

Statistical Analyses
The experimental results are expressed as means ± SE. SPSS Statistics version 25 (IBM Microsoft, New York, NY, USA) was used to carry out statistical analyses (Duncan multiple range tests). P values < 0.05 were regarded as statistically significant difference.

Identification of Histone Deacetylases in Chinese Cabbage
Since HDACs have been introduced as erasers of epigenetic acetylation marks on lysine residues in histones, a number of plant HDAC genes have been identified from a variety of plant species [3]. From the B. rapa genome, we identified 20 HDAC members (Table 1). Among them, BraHDA6 was the largest protein, with 651 amino acids (AA), whereas BraHDA1 was identified as the smallest protein with 254 AA. The molecular weight of BraHDACs varied according to protein size, ranging from 28.8 KDa to 72.1 KDa, and their pI varied from 4.52 (BraHDT2) to 8.96 (BraSRT1). Subcellular locations as predicted by WoLF PSORT indicated that BraHDACs were potentially localized in the cytosol, nucleus, chloroplast, mitochondria, and peroxisome (Table 1). Similar to the Arabidopsis and soybean HD2 proteins [23,24], the B. rapa HD2 family were predicted as nuclear proteins (Table 1). In addition, BraSRT1 and BraSRT2 appeared to be localized in the mitochondria and nucleus, respectively, similar to the rice SRT (silent information regulator) family [25]. Taken together, the various localization patterns of BraHDACs indicated that they may have distinct roles.  [11,12]. As shown in Figure 1A, BraHDACs were divided into three major classes: RPD3/HDA1, HD2, and SIR2 families. The B. rapa genome encodes 16 proteins (Figure 1 and Figure S1) clustered in the RPD3/HDA1 family contained the conserved Hist_deacetyl domain. In addition, two SIR2 family HDACs with highly conserved SIR2 domain and two HD2 family HDACs with conserved zinc finger C2H2 type domain (ZnF_C2H2) were identified. Another type of zinc finger domain, ZnF_RTZ (SM000547), has been found in BraHDA07, indicating that BraHDT1, BraHDT1, and BraHDA07 have high binding affinity to DNA or modulate protein-protein interactions mediated by zinc finger domains [26]. The presence of domains similar to those in other plant HDACs suggested that all putative BraHDACs belong to the histone deacetylase family.

Gene Duplication Analysis of BraHDACs
Gene duplication has long been regarded as a key evolutionary process for acquiring new genes and genetic novelty [27]. It could occur on any of these scales: whole genome duplication, segmental duplication, and single gene duplication including tandem duplica-tion [28]. It has been shown that some HDACs expanded in the Gossypium and soybean genomes partly because of segmental duplication events [12,24]. In the B. rapa genome, 20 BraHDACs were asymmetrically distributed on 7 chromosomes (Table 1). Chromosome A06 contains the highest density of HDACs with six members (BraHDA8 to 13), followed by chromosome A01 with four genes (BraHDA1 to 3, BraSRT1), whereas only one BraHDAC was present on chromosome A05 (Table 1). A tandem duplication is characterized as two or more adjacent homologous genes located physically "close" to one another; tandem duplicated pairs are paralogous genes on the same chromosome within a 50 kb physical distance [29], indicating that only one set (BraHDA8 and 9) was tandemly duplicated on chromosome A06 separated by 14,097 bp (Table 1) with 85.7% and 75.4% of sequence identity at the nucleotide and protein levels. To further determine whether segmental duplication events also contributed to the expansion of BraHDACs, we analyzed the PGDD database combined with phylogenetic analysis. As shown in Figure 1B, seven duplicated BraHDAC gene pairs were identified, indicating that segmental duplication is a major contributory event leading to the expansion of BraHDAC genes in B. rapa genome. The highest frequency of segmental duplication events occurred between chromosomes A06 and 09, which contained two segmentally duplicated gene pairs.
To examine the evolutionary selection process of eight duplicated gene pairs in B. rapa genome, we analyzed the evolutionary constraint (Ka/Ks). Ka/Ks ratios ranged from 0.0895 to 0.8713, and the average Ka/Ks value of the duplicated gene pairs was 0.404 (Table 2). This indicates that these duplicated gene pairs were subjected to negative selection (purifying selection), similar to duplicated HDAC gene pairs in Gossypium and soybean genomes [12,24]. The segmental duplications of HDAC genes in B. rapa genome occurred between 12.73 MYA and 42.04 MYA, with the average value being 29.19 MYA. The Ks of tandem duplications of BraHDAC genes was 0.1640, dating the duplication event at 5.47 MYA (Table 2). Brassicaceae, including Arabidopsis and Brassica, have been shown to have undergone three rounds of whole genome duplication, including γ triplication (135 MYA) and β (90-100 MYA) and α (24-40 MYA) duplications [30][31][32]. In addition, B. rapa exhibits this complex history, with the addition of a whole-genome triplication occurring between 13 and 17 MYA [33], indicating that the segmental duplications of BraHDAC genes, except BraHDA 5/7 and BraHDA 12/15 gene pairs, occurred before wholegenome triplication of the Brassica genome.

Effects of Heat Stress on Acetylation Status of Histone H3 in Chinese Cabbage
Heat stress induces excessive generation of reactive oxygen species, resulting in alterations in plant growth, development, photosynthesis, physiological processes, and production yield [34]. To evaluate the efficacy of heat treatment (40 • C), the physiological response of Chinese cabbage including lipid peroxidation and proline content was analyzed. As shown in Figure 2A, heat stress induced the accumulation of malondialdehyde (MDA), which is an indicator of lipid peroxidation [35]. In higher plants, the accumulation of proline, an amino acid and a compatible solute, is a common phenomenon in response to various environmental stresses [36]. We observed an increase in proline levels in response to heat stress ( Figure 2B). The accumulated proline purportedly acts as a compatible osmolyte, molecular chaperone, cytosolic pH buffer, cell redox balancer, free radical scavenger, and stabilizer for subcellular structures during various stresses, especially osmotic and salt stresses [36,37], and improves resistance to various stresses [38]. However, proline accumulation in response to heat stress leads to increased production of reactive oxygen species in the mitochondria and inhibits ABA and ethylene biosynthesis, resulting in increased sensitivity to heat stress [36]. Taken together, these physiological changes suggested that the heat stress treatment was effective, and the leaves of heat-treated plants were harvested for further analyses.
Agriculture 2021, 11, x FOR PEER REVIEW 7 of 10 patterns were altered in response to different stresses including drought, salt, and wounding ( Figure S2), indicating that BraHDACs had diversified substantially, which might be due to the different divergent fates after duplications.

Histone Deacetylase Activity of BarHDA3
Although many fundamental questions remain unanswered regarding the physiological function of plant HDACs in response to environmental stresses, it has been shown that RPD3/HDA1 family including Arabidopsis HDA5/14/15/18/19 positively or negatively regulated to the environmental-stress tolerances [41]. Among them, AtHDA14 controlled the activation state of RuBisCO under low-light condition, and is involved in biosynthesis of melatonin, which plays an important role in plant response to environmental stresses [42,43]. Under heat stress condition, exogenous melatonin improved the antioxidant defense systems, resulting in the suppression of heat stress-induced damage [44]. In Chinese cabbage, BraHDA3, homologous of AtHDA14 ( Figure S1), dramatically downregulated in response to heat stress ( Figure 2C). This indicates that down-regulation of BraHDA3 might be required for the expression of heat response genes, although the physiological function of BraHDA3 needs to be characterized. To analyze the function of BraHDA3, recombinant BraHDA3 (BraHDA3-His) was prepared using E. coli expression system. As shown in Figure 3, purified recombinant BraHDA3 protein gave an activity of 18.1 µM/100 µg of protein. When HDAC inhibitor SAHA (5 µM) was mixed with purified recombinant BraHDA3 protein, the HDAC activity decreased to 3.8 µM/100 µg of protein. This indicates that BraHDA3 encodes the functional HDAC enzyme, which can be inhibited by Class I/II HDAC inhibitor SAHA. To gain insight into the possible involvement of BraHDACs during their response to heat stress, the expression profiles of BraHDACs were analyzed. As shown in Figure 2C, exposure of Chinese cabbage plants to heat stress induced or repressed the expression of BraHDACs. After 48-h exposure to heat stress, expression levels of BraHDA2, BraHDA4, BraHDA5, BraHDA6, BraHDA7, BraHDA11, BraHDA12, BraHDA15, BraHDA16, BraHDT2, and BraSRT2 had increased, whereas those of other BraHDACs decreased, remained unchanged, or were undetected. In addition, the transcript level of BraHDA13, BraHDA15, BraHDA16, BraHDT2, and BraSRT2 was up-regulated in response to drought, salt, and wounding, whereas the decreasing level of BraHDA2, BraHDA3, BraHDA9, and BraHDA14 transcripts was observed ( Figure S2). Under heat-stress condition, BraHDA4 and BraHDA14 exhibited a different expression pattern, although BraHDA4 and BraHDA14 were clustered into same sub-group in RPD3/HDA1 family, indicating that the duplication of those genes resulted in divergence of their expression pattern and function. Similar with BraHDACs, several Arabidopsis HDACs and maize HDACs in the class II of RPD3/HDA1 family were induced by heat treatment [11,39]. In addition, Arabidopsis, maize and cotton HDACs were differentially expressed in response to hormones and stresses, suggesting the functional divergence and specific regulation of HDACs in response to a particular stress [11,12,39]. According to the abiotic stress induced expression pattern of BraHDACs, we found that BraHDA8 was not expressed in leaves and on abiotic stresses. Based on transcriptome analysis, transcript level of BraHDA8 has been detected only in the silique [40], indicat-ing that BraHDA8 should be expressed in specific tissue. In addition, the expression pattern of paralogous gene pair (BraHDA4/14) diverged, whereas three paralogous gene pairs (BraHDA4/6, BraHDA5/7, and BraHDA12/15) exhibited similar expression patterns ( Figure 2C) under heat stress condition. However, these expression patterns were altered in response to different stresses including drought, salt, and wounding ( Figure S2), indicating that BraHDACs had diversified substantially, which might be due to the different divergent fates after duplications.

Histone Deacetylase Activity of BarHDA3
Although many fundamental questions remain unanswered regarding the physiological function of plant HDACs in response to environmental stresses, it has been shown that RPD3/HDA1 family including Arabidopsis HDA5/14/15/18/19 positively or negatively regulated to the environmental-stress tolerances [41]. Among them, AtHDA14 controlled the activation state of RuBisCO under low-light condition, and is involved in biosynthesis of melatonin, which plays an important role in plant response to environmental stresses [42,43]. Under heat stress condition, exogenous melatonin improved the antioxidant defense systems, resulting in the suppression of heat stress-induced damage [44]. In Chinese cabbage, BraHDA3, homologous of AtHDA14 ( Figure S1), dramatically down-regulated in response to heat stress ( Figure 2C). This indicates that down-regulation of BraHDA3 might be required for the expression of heat response genes, although the physiological function of BraHDA3 needs to be characterized. To analyze the function of BraHDA3, recombinant BraHDA3 (BraHDA3-His) was prepared using E. coli expression system. As shown in Figure 3, purified recombinant BraHDA3 protein gave an activity of 18

Conclusions
Despite the knowledge concerning HDACs, evolutionary and functional information regarding HDACs in B. rapa remains relatively unknown. In this study, we identified 20 putative BraHDAC genes divided into three major classes: RPD3/HDA1, HD2, and SIR2 families. Most HDAC gene duplications in B. rapa appeared to have been caused by segmental duplication, which occurred between 12.73 MYA and 42.04 MYA. Several BraHDACs were up-or down-regulated by abiotic stress treatments, indicating that the variation of BraHDAC transcription levels may contribute to alteration in histone H3 acetylation levels in response to various stresses. Our results provide a solid foundation for further understanding of the underlying evolutionary mechanisms in the HDAC family and will provide the basis for future research on functional characterization of the BraHDAC family in response to environmental stresses.

Conclusions
Despite the knowledge concerning HDACs, evolutionary and functional information regarding HDACs in B. rapa remains relatively unknown. In this study, we identified 20 putative BraHDAC genes divided into three major classes: RPD3/HDA1, HD2, and SIR2 families. Most HDAC gene duplications in B. rapa appeared to have been caused by segmental duplication, which occurred between 12.73 MYA and 42.04 MYA. Several BraHDACs were up-or down-regulated by abiotic stress treatments, indicating that the variation of BraHDAC transcription levels may contribute to alteration in histone H3 acetylation levels in response to various stresses. Our results provide a solid foundation for further understanding of the underlying evolutionary mechanisms in the HDAC family and will provide the basis for future research on functional characterization of the BraHDAC family in response to environmental stresses.