Peptide sequences of all Arabidopsis HD2 genes were obtained from the Phytozome (https://phytozome-next.jgi.doe.gov/) and a phylogenetic tree was generated using MEGA X (Kumar et al., 2018). Amino acid sequence identities and divergences were determined using the National Center for Biotechnology Information (NCBI) tool BLAST for nucleotide sequences (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Arabidopsis thaliana wild-type (WT) and T-DNA insertional mutants GK355-H03 (hd2a: At3g44750), SAIL1247-A02 (hd2b: At5g22650), SALK129799 (hd2c: At5g03740) and GK279-D04 (hd2d: At2g27840) were used in this study. All these genotypes were in accession Columbia (Col-0) background and have been used previously (Song et al., 2010; Buszewicz et al., 2016; Li et al., 2017). The T-DNA mutant lines were obtained from the Arabidopsis Biological Resource Center (ABRC). The double mutant lines used in this study were generated by crossing hd2 single mutants with each other in different combinations. All the T-DNA insertion mutant lines were verified by genotyping by PCR using primers listed in Supplementary Table 1 .

For plant growth experiments, all seeds were surface-sterilized in ethanol (70%) followed by bleach (25%) and stratified in the dark for 3 days at 4°C before sowing on ProMix-BX soil (Premier Horticulture, Québec) or on Murashige and Skoog (MS) growth medium (2.15 g/L) with or without additives such as antibiotics, plant growth regulators, or chemicals. Plants in the growth room were grown under constant 16 hours of daylight (long-day conditions) at 22°C with a relative humidity of 60%. Plants in the growth chambers were also grown under the same conditions, except dark cycles were at 18°C.

To generate HD2 overexpression (OE) vectors (CaMV35S:AtHD2), the full-length cDNA excluding the stop codon was amplified and cloned into the entry vector, pDONR221 (Invitrogen) using BP Clonase reaction mix (Thermo Fisher Scientific cat. 11789020) according to manufacturer’s protocol. The cloned vector was transferred into Escherichia coli strain DH5α. Subsequently, the cloned DNA was transferred from the entry vector to the destination vector, pEarleyGate101 (Hartley et al., 2000; Earley et al., 2006) using LR Clonase reaction mix (Thermo Fisher Scientific cat. 11791100) according to manufacturer’s protocol. The HD2-OE constructs were transferred to Agrobacterium tumefaciens strain GV3101 via electroporation (Wise et al., 2006). Primers used for gene cloning are listed in Supplementary Table 2 . Plant transformation of WT Arabidopsis was carried out using the floral dip method (Zhang et al., 2006). The putative transgenic seeds were grown on selection medium (MS agar medium containing BASTA (glufosinate ammonium 25 mg/L) and two highly expressed transgenic lines for HD2A and HD2C were selected.

To study mRNA expression of different genes in hd2 mutants and HD2 overexpression lines, quantitative real time PCR (RT-qPCR) was performed. RNA was isolated from leaves or seedlings using the Plant/Fungi Total RNA Purification Kit (Norgen) and cDNA was synthesized using iScript Reverse Transcription Supermix (Bio-Rad). The RT-qPCR was performed using a SsoFast EvaGreen Supermix kit (Bio-Rad) on CFX96 Real-time PCR detection system (Bio-Rad) following the manufacturer’s instructions. Data were analyzed using Bio-Rad CFX Manager 3.1 software. The expression levels were normalized to the ACTIN2 gene. The ΔΔCT method was applied to calculate fold change in the expression level (Livak and Schmittgen, 2001). Three independent RT-qPCR analysis were performed with three technical replicates for each biological replicate. All primers used for RT-qPCR analysis are listed in Supplementary Table 3 .

Nuclei protein was extracted following the histone extraction protocol (abcam). Briefly, 10-day old seedlings were ground and fixed in Triton Extraction Buffer [TEB: PBS containing 0.5% Triton X-100, 2 mM phenylmethylsulfonyl fluoride (PMSF)]. After incubating the mixture on ice for 10 minutes, samples were centrifuged at 2000 rpm for 10 minutes at 4°C and the supernatant was discarded. TEB buffer (100 µl) was added to the samples and the centrifugation step was repeated. Cell pellets were resuspended in 0.4 µl of 0.2 N HCl and samples were kept at 4°C for 4 hours. Samples were then centrifuged at 2000 rpm for 10 minutes at 4°C. The supernatant was collected and the protein concentration was estimated using Bradford assay reagents (Thermo Fisher Scientific) following the standard protocol.

To separate proteins by SDS PAGE, samples were denatured by adding 20 mM dithiothreitol (DTT) and heating at 95°C for about 10 minutes. Samples were loaded onto SDS-PAGE gels (Bio Rad) to separate the proteins. Proteins from gels were transferred to PVDF membranes (Bio-Rad) using a Trans-Blot Semi-Dry Electrophoretic Transfer Cell (Bio-Rad). The PVDF membranes were then incubated overnight at 4°C in anti-H3K9ac or Anti-H3 antibody solution (Cell Signaling Technology and Millipore), followed by incubation in the secondary antibody (Milipore Sigma) solution for 1 hour. Proteins blots were detected using the EZ-ECL Chemiluminescence Detection Kit (Biological Industries) as per the manufacturer’s instructions with a MicroChemi imager (DNR Bio-Imaging Systems).

Soil drying or Polyethylene glycol (PEG-6000) methods were employed for plant drought treatment (Jiang et al., 1995; Landjeva et al., 2008). For the soil drought method, four-week old plants in the drought group were withheld from watering for 7 days to establish a drought stress. For PEG-6000, 10-day old seedlings were exposed to PEG (10%) in MS media for 24-72 hours before sample collection for gene expression analysis.

For the drought stress tolerance assay, two-week old hd2 mutant or HD2-OE plants in drought-treated groups were withheld from watering for 11 and 13 days, respectively (Ge et al., 2018). After the drought stress period, plants were rewatered for two days to allow them to recover. Plants with at least two green turgid rosette leaves were considered to have survived. The survival rate for each genotype was determined based on the number of plants survived out of total number of plants grown.

Leaf relative water content (RWC) was measured from plant leaves as previously described (Barrs and Weatherley, 1962). A total of three healthy rosette leaves were excised from the plants and fresh weight (FW) of all the collected leaves was determined immediately. All the leaves were rehydrated by floating them on water surface for 4 hours and leaf turgid weight (TW) was obtained. All the turgid leaves were oven-dried at 55°C for 24 hours. Oven-dried leaves in the tubes were weighed again to determine the leaf dry weight (DW). Leaf RWC (%) was determined as [(FW – DW)/(TW – FW)] x 100.

To measure fresh leaf water loss in different genotypes, healthy leaves were detached from plants grown under control conditions and immediately weighed to measure their fresh weights. Leaves were weighed after every 30-minute interval for a total of 180 minutes. During each 30-minute interval, leaves were kept open in petri dishes at room temperature. Results were calculated as the percentage water loss compared to the leaf fresh weights.

To examine stomatal opening and closing in hd2 mutants and HD2-OE lines, stomatal aperture measurement assays were performed following the established protocol (Eisele et al., 2016; Scarpeci et al., 2017) with minor modifications. Briefly, healthy leaves were excised from plants and floated (with abaxial side down) in petri dishes containing MES/KOH stomatal opening buffer (10 mM MES, 5 mM KCl, 50 μM CaCl2, pH 6.15) for two hours and was used as control. For ABA treatment, the same steps were followed, except that, after two hours, 10 μM of ABA was added to the petri dishes containing MES/KOH buffer. After two hours of ABA treatment, leaves were taken out from the solution and quickly pat dried with napkins and prepared for microscopy. Preparation of leaf imprints to visualize stomata using the microscopy was performed by following the protocol described by Scarpeci et al. (2017) with minor modifications. Briefly, leaf imprints were obtained by applying a thin layer of transparent nail varnish on the abaxial leaf surface. The applied nail varnish layer was dried for about 30 minutes. The thin layer of nail varnish was transferred to a glass slide by gently pressing the leaf onto the slide. Stomata imprints were observed under a 40x objective lens using the EVOS™ XL Core Imaging System (AMEX1000). Stomatal pore widths and lengths were measured using the image processing software ImageJ (http://rsb.info.nih.gov/ij/). Results for stomatal aperture were shown as the ratio of width to length.

Primary root lengths and secondary root number were examined in different genotypes under control and different treatments. For osmotic stress in seedlings in MS media, PEG-6000 was used (Jiang et al., 1995; Landjeva et al., 2008). Seeds of hd2 mutants and HD2-OE lines were grown on agar plates containing only MS media or MS media supplemented with PEG-6000 (10%). For gibberellic acid treatments, GA3 and GA4 were added to the MS media at 1 μM and 10 μM concentrations. All the plates were placed vertically in the growth chamber to allow longitudinal root growth. Primary root length and number of secondary roots were examined on day 10 of germination for all the experiments.

To study the protein-protein interactions in vivo, yeast two-hybrid (Y2H) assays were performed following the instructions provided in the Matchmaker Gold Yeast Two-Hybrid System manual (Clontech) (Brückner et al., 2009). Using LR Clonase reaction mix, the cloned DNA from pDONR221 was transferred to destination vectors, pGBKT7 (bait vector), and the pGADT7 (prey vector). The Y2HGold strain was co-transformed with 5 µg of each of bait and prey vectors containing HD2 genes. For negative controls, the yeast strain was transformed with bait vector and empty prey vector (containing no gene of interest). The pGBKT7-53 and pGADT7-T plasmids provided with the kit (encode murine p53 and SV40 large T-antigen) were used as positive controls. The transformed Y2HGold cells were plated on minimal selective medium SD-Leu/Trp double dropout (DDO) and sub-cultured to selective medium SD-Leu/Trp/Ade/His quadruple dropout (QDO) and QDO supplemented with 40 mg/L X-Alpha Gal (5-Bromo-4-Chloro-3-indolyl a-D-galactopyranoside) and 200 µg/L Aureobasidin-A (QDO/X/A).

Bimolecular fluorescence complementation (BiFC) assays were performed to confirm the protein-protein interactions in a plant host (Tian et al., 2011). Using LR Clonase reaction mix, the cloned DNA from pDONR221 was transferred to destination vectors pEarleyGate201-YN and pEarleyGate201-YC to generate plasmid constructs for the BiFC assay. Both vectors were transferred to Agrobacterium tumefaciens strain GV3101 via electroporation. Agrobacterium cultures containing pEarleyGate201-HD2-YN and pEarleyGate201-HD2-YC constructs were co-infiltrated into Nicotiana benthamiana leaf epidermis. (Sparkes et al., 2006; Tian et al., 2011). For negative controls, pEarleyGate201-HD2-YN construct was co-infiltrated with an empty pEarleyGate202-YC vector. The YFP signal was observed after 48-72 hours using the Olympus Confocal Laser Scanning Microscope FV3000. The argon excitation laser wavelength was set at 514 nm to visualize the YFP signal.

For statistical analysis, Microsoft Excel 2016 (Microsoft Corp., Washington, USA) was used to calculate simple univariate statistics, such as means, standard deviations, and standard errors. The statistical analyses were performed using IBM SPSS Statistics version 25.0 (IBM Corp. Armonk, New York). A p-value < 0.05 (*) or 0.01 (**) was used to show a statistically significant difference.

First, we performed phylogenetic and sequence alignment analysis of Arabidopsis HD2 proteins. Analysis showed that HD2A and HD2B are the closest relatives and share the highest sequence similarity (50.5%) and the least divergence (0.493) with each other at the protein level ( Supplementary Figure S1 ). HD2A and HD2C shares the second highest similarity (40.5%) and the second least divergence (0.802), whereas HD2D is a distantly related member of the HD2 family and shares the least similarity and the maximum divergence with HD2A, HD2B and HD2C.

To examine how HD2 family genes respond to drought stress, we measured the expression levels of the HD2 family genes, HD2A, HD2B, HD2C, and HD2D under drought conditions. For soil drought, four-week old Arabidopsis wild-type (WT) plants were withheld from watering and leaf relative water content was measured to determine the water status of the plants ( Supplementary Figure S2A, B ). Leaf samples were collected daily from both the control group and the drought-treated group to examine gene expression by quantitative PCR (qPCR). The Responsive to desiccation 29A (RD29A) gene, which is induced under different abiotic stresses and frequently used as a drought responsive marker (Msanne et al., 2011; Bihmidine et al., 2013), was used to affirm the effectiveness of the drought treatment at the molecular level ( Supplementary Figure S2C ). RD29A showed significant upregulation starting on day 3 of drought treatment. When an effective drought stress initiating from day 3 was confirmed, we then examined the expression of HD2A, HD2B, HD2C and HD2D using samples from day 3 of the drought treatment. Analysis showed that all four HD2 genes were upregulated in response to drought ( Figure 1A ). Polyethylene glycol (PEG) is frequently used to mimic drought stress in plants and to study related responses (Jiang et al., 1995; Landjeva et al., 2008). To establish osmotic stress in MS media, seedlings were treated with PEG-6000 (10%) for 24 hours, and expression of HD2s was measured. Analysis showed that all HD2 genes were upregulated in response to PEG-induced osmotic stress ( Figure 1B ).

To investigate the relationship of different HD2s in drought stress response, we used hd2a, hd2b, hd2c, and hd2d T-DNA insertional mutants in our study. Genotyping by PCR was done to identify homozygous mutant lines for each HD2 gene ( Supplementary Figure S3 ). qPCR analysis was performed to confirm the knockout of HD2 gene expression at the mRNA level in hd2 single mutants. Analysis showed that the expression of the HD2A and HD2C was completely knocked out in the hd2a and hd2c mutant lines, respectively, whereas HD2B and HD2D genes retained their expressions at 21% and 10% of the WT levels in the hd2b and hd2d mutant lines, respectively.

HD2C has been shown to be involved in ABA and abiotic stress response and can mediate root growth under stress (Sridha and Wu, 2006; Luo et al., 2012c; Buszewicz et al., 2016; Chen et al., 2018). To study the relationship of HD2C with other HD2s, we crossed the hd2c mutant line with other hd2 single mutants and successfully generated double mutant lines hd2a.hd2c (Mac16) and hd2b.hd2c (Mbc68). We could not obtain seeds for hd2c.hd2d double mutant line in several attempts of crossing hd2c with hd2d. The Mbc68 double mutant line was morphologically similar to the single mutants hd2b and hd2c and did not show any significant change in primary root length as compared to the respective single mutants hd2b and hd2c ( Supplementary Figure S4 ). The Mac16 line was morphologically distinct from the respective single mutants hd2a and hd2c ( Supplementary Figure S5 ). Due to the relative paucity of information on hd2a.hd2c double mutants, we focused our attention on the Mac16 (hd2a.hd2c) double mutant line along with hd2a and hd2c single mutants for further study. The homozygosity of T-DNA insertions in the double mutant line Mac16 was confirmed by PCR genotyping ( Figure 2A ). The qPCR analysis of hd2 single and double mutants revealed that when the HD2A is knocked out, HD2C is upregulated in the hd2a mutant background. Similarly, when the HD2C is knocked out, HD2A is upregulated in the hd2c mutant line ( Figure 2B ). Next, we performed drought treatment assays to examine the survival of the single mutants hd2a and hd2c as well as the double mutant Mac16 under drought ( Figure 2C, D ). The single mutants hd2a and hd2c showed 79% and 67% survival, respectively, relative to WT. In contrast, the double mutant Mac16, showed only 47% survival relative to WT, which is also significantly lower than that observed in the single mutants. We also measured the fresh leaf water loss in these genotypes. The hd2a did not show significant change in the water loss, whereas the hd2c had increased water loss. The knockout of both HD2A and HD2C resulted in a significantly increased water loss in the Mac16 as compared to both single mutants and WT ( Figure 2E , Supplementary Figure S7 ).

Along with the double mutant Mac16, we also generated HD2A (OEA1 and OEA2) and HD2C (OEC1 and OEC2) overexpression (OE) lines ( Supplementary Figure S6 ) and studied whether HD2A or HD2C overexpression affects the survivability under drought stress ( Figures 3A, B ). While the WT plants showed 34% survival, the HD2A-OE lines, OEA1 and OEA2 showed increased survival of 52% and 56%, respectively, whereas HD2C-OE lines OEC1 and OEC2 showed 84% and 77% survival, respectively. Fresh leaf water loss analysis showed that both HD2A and HD2C overexpression lines had reduced water loss from the leaves as compared to WT ( Figures 3C, D , Supplementary Figure S7 ).

Accumulation of ABA in leaves under drought stress initiates a signalling cascade to induce stomatal closure to minimize transpirational water loss (Hopkins and Hüner, 2008). Stomata are closed to limit water loss and establish an equilibrium between water absorbed by roots and water loss by transpiration (Wilkinson and Davies, 2002). We examined stomatal closure in HD2-OE and hd2 mutant lines. The HD2A-OE lines did not show any differences in stomatal aperture under control conditions, whereas stomata closed completely under ABA-induced stress as compared to WT ( Figures 4A, B ). Stomata opened fully in HD2C-OE lines under control conditions and closed completely under ABA-induced stress as compared to WT ( Figures 4C, D ). Among the hd2 mutants, hd2a did not show any change in stomatal aperture under control conditions, whereas the hd2c and Mac16 mutants showed reduced stomatal aperture with respect to WT ( Figures 5A, B ). Under ABA-induced stress, the hd2c mutant showed larger stomatal aperture as compared to the hd2a and WT. Knockout of both genes in the Mac16 resulted in a significantly larger stomatal aperture as compared to both single mutants and WT. This suggests that HD2A and HD2C corelate in regulating stomatal closure to play a role in the drought stress response.

Accumulation of ABA in leaves under drought stress stimulates a core signalling pathway to induce stomatal closure by upregulating the expression of SLAC1 gene (Vahisalu et al., 2008; Munemasa et al., 2015; Zhang et al., 2016), whereas ABI1 and ABI2, the two protein phosphatases, are considered negative regulators of ABA-mediated stomatal closure. We examined the expressions of ABI1, ABI2, and SLAC1 genes in hd2 mutants under drought stress. Analysis of ABI1 and ABI2 gene expression ( Figure 5C ) showed that their expression was not changed in the hd2a mutant, while the hd2c mutant exhibited an upregulation of both genes. Knockout of both HD2A and HD2C in Mac16 resulted in a significant upregulation of ABI1 and ABI2 under drought stress as compared to the single mutants and WT. Gene expression analysis of SLAC1 ( Figure 5D ) showed that its expression was not affected in hd2a but decreased in the hd2c. Knockout of both HD2A and HD2C in Mac16 resulted in a drastic decrease of SLAC1 gene expression. Taken together, the HD2A and HD2C play a collective role in regulating stomatal opening and closure. Knockout of both genes led to a stronger phenotype in terms of drought tolerance and stomatal closure as compared to the single mutants, indicating an additive effect of HD2A and HD2C in drought stress response.

Roots also play a vital role in the strategy of plant response to drought. Our next objective was to study if HD2A and HD2C play a role in regulating root growth under drought stress. First, we measured the primary and secondary root growth of HD2-OE and hd2 mutant lines under control and drought stress conditions mimicked by PEG-6000. Overexpression of both HD2A and HD2C resulted in increased primary root length under control and osmotic stress conditions ( Figure 6 ). These overexpression lines also showed a higher number of secondary roots under control conditions. Conversely, primary root length was significantly decreased in the Mac16 (71%) and hd2c (87%) mutants under control conditions as compared to WT ( Figures 7A, B ), whereas root length did not change in the hd2a mutant. Application of osmotic stress (PEG) caused a significant decrease in the primary root lengths of both hd2a (86%) and hd2c (79%) and a further decrease in the Mac16 (59%) as compared to WT ( Figures 7C, D ). Regarding the secondary root number in hd2 mutants under control conditions, the hd2c did not show difference in the number of secondary roots, whereas the hd2a showed fewer secondary roots in comparison to WT. Knockout of both HD2A and HD2C resulted in a significant decrease in the number of secondary roots in the Mac16 line ( Figure 7E ).

The GA2ox enzymes play an important role in maintaining a dynamic homeostasis of bioactive GAs for normal plant root growth. We thus investigated whether the HD2A and HD2C can regulate the expression of the GA2ox genes. First, we performed western blot to examine global H3K9 acetylation levels in the hd2 mutants and HD2-OE lines. The H3K9ac/H3 levels were found to be significantly higher in the Mac16 (1.41) as compared to the hd2a (1.14), hd2c (1.18), and WT ( Figure 8A ). Conversely, lower levels of H3K9ac/H3 were found in both HD2A and HD2C overexpression lines ( Figures 8B, C ). Next, we examined the expression of GA2ox family genes in hd2 mutants and HD2-OE lines under control conditions. Expression of the GA2ox1 and GA2ox2 was upregulated many fold in the Mac16 as compared to the single mutants and WT ( Figure 8D ). The hd2a mutant did not show change in the expression of both genes. On the other hand, HD2A and HD2C overexpression caused a significant downregulation of both GA2ox1 and GA2ox2 ( Figures 8E, F ). Other GA2ox genes (GA2ox3, GA2ox4, GA2ox6) did not show significant changes in their expression levels. Previously, it has been shown that the transcriptional activation of GA2ox genes is associated with histone acetylation status triggered by ABA in response to internal developmental signals and abiotic stress conditions (Colebrook et al., 2014; Li et al., 2017; Chen et al., 2019). We analysed the expression of the GA2ox1 and GA2ox2 in WT under ABA and PEG treatments ( Figure 8G ). Analysis showed that the expression of both GA2ox1 and GA2ox2 was significantly upregulated in response to ABA and PEG. GA2ox2 showed many fold higher expression as compared to the GA2ox1. These results indicate that GA2ox1 and GA2ox2 are upregulated in drought stress response via ABA pathway and their expression is co-regulated by HD2A and HD2C, possibly through histone acetylation modification in response to ABA and osmotic stress.

Bioactive gibberellic acid GA4 serves as a substrate for GA2ox-mediated GA catabolism and GA3 is not considered its substrate. Upregulation of GA2ox genes leads to the lower levels of GA4 and GA1 in Arabidopsis (Yamauchi et al., 2007; Rieu et al., 2008; Li et al., 2017). We thus examined the primary root lengths in hd2 mutants under GA4 and GA3 treatments. The hd2 mutant lines did not show any difference in root phenotype under 1 μM of GA4 treatment, with a slight increase in root lengths under 10 μM of GA4 ( Figure 8H ). However, GA3 treatments (1 μM and 10 μM) in hd2 mutants led to a significant increase in primary root lengths of Mac16 as compared to the hd2a, hd2c, and WT ( Figure 8I ).

The additive effect of HD2A and HD2C in drought stress response and root growth led us to investigate if HD2A and HD2C can physically interact with each other. To study this protein-protein interaction, we performed yeast two-hybrid (Y2H) and bimolecular fluorescent complementation (BiFC) assays. In Y2H assays, the yeast strains containing constructs HD2A-AD or HD2C-AD along with an empty vector (pGBKT7) did not grow on selective medium QDO and QDO/X/A. However, when co-transformed with constructs expressing HD2A-AD and HD2C-BK (or HD2A-BK and HD2C-AD), yeast cells developed colonies on QDO and QDO/X/A, indicating a physical interaction of the HD2A and HD2C ( Figure 9A ). To further confirm this protein interaction in the plant host, BiFC assays were performed. The HD2A-YN or HD2C-YN construct did not show any YFP signal (yellow fluorescence) when co-transformed with an empty pEarleyGate202-YC vector. YFP signal in leaves was observed when both HD2A-YN and HD2C-YC constructs (or HD2A-YC and HD2C-YN) were infiltrated together, again indicating a physical interaction of HD2A and HD2C ( Figure 9B ).