#### Sub-Cellular Localization of OoNAC72

The sequence analysis showed that OoNAC72 possessed a conserved nuclear localization signal (NLS, positions 74–88, 117–129 aa) (**Supplementary Figure S1**). Meanwhile, analysis by ProtComp v.9.0 indicated a high likelihood of nuclear localization for OoNAC72 protein. To confirm this prediction, an expression cassette fusing OoNAC72 with the eGFP protein was constructed. Then the fused protein was expressed transiently in Nicotiana benthamiana while the pCAMBI1302 eGFP functioned as a control. Fluorescence microscopy revealed that the OoNAC72-eGFP fusion protein was exclusively localized in the nucleus in the transformed cells, whereas the control eGFP was uniformly distributed throughout the cell (**Figure 1**). These results further confirmed that the OoNAC72 protein was a nuclear-localized protein.

### Transcriptional Activity Assay in Yeast Cells

The result of transcriptional activity analysis of OoNAC72 was illustrated in **Figure 2**. The vector pGBKT7 fused with OoNAC72 was used as the experimental set, and the empty pGBKT7 vector and pGAL4 were hired as the negative and positive controls, respectively (**Figure 2A**). All of the transformants grew well on selective SD/-Trp medium (**Figure 2B**), indicating that the three vectors were successfully transformed into the yeast cells. The GAL4-binding domain-OoNAC72 construct and pGAL4 grew well on SD/-Trp-His-Ade medium, while the transformants containing the pGBKT7 vector did not grow on the same medium (**Figure 2C**). The results of α-galactosidase activity assays showed that transformants containing GAL4 binding domain-OoNAC72 construct and pGAL4 appeared blue in color on SD/-Trp-His-Ade medium containing 5-bromo-4 chloro-3-indoxyl-α-D-galactopyranoside (X-α-Gal) (**Figure 2D**). These results indicated that OoNAC72 had transcriptional activity in yeast.

# Expression Patterns of OoNAC72 Under Various Treatments

independent replications. Asterisk indicates significant difference (∗P < 0.05; ∗∗P < 0.01) between transgenic lines and WT.

Quantitative real-time PCR was used to evaluate the expression patterns of OoNAC72 during dehydration (20% PEG-6000), salt (150 mM NaCl) and three hormone treatments (100 mM ABA, ETH, and GA) (**Figure 3**). When the seedlings were treated with dehydration, OoNAC72 mRNA abundance was slightly induced at 1 h, followed by progressive elevation until reaching the peak value at 3 h, which showed an approximately 11.9-fold increase relative to the initial level (**Figure 3A**). Treatment with 150 mM NaCl led to a quick accumulation of OoNAC72 mRNA level, which progressed until a maximum level was reached at 6 h, approximately 12.3-fold increase relative to the initial level (**Figure 3A**). Moreover, for the ABA treatment, the expression levels of OoNAC72 increased rapidly and reached the maximum level at 3 h, being 41.8-fold greater than the control (0 h), while there were no significant changes when treated with ETH and GA (**Figure 3B**). All the results suggested that OoNAC72 may respond to stress in O. ochrocephala by participating in ABA-dependent signal transduction pathways.

#### Overexpression of OoNAC72 Changed the Phenotype of Arabidopsis

To explore the function of OoNAC72 during the tolerance to abiotic stress, transgenic Arabidopsis over-expressing plants driven by the CaMV35S promoter were generated. We analyzed the growth status of two OoNAC72-OX transgenic lines (OX1, OX2) from the screening transgenic positive seedlings. Compared with wild-type plants, OoNAC72- OX transgenic lines (OX1, OX2) showed rosette leaves during the vegetative growth stage, and the number of leaves was significantly increased (**Figures 4A,B,D**); At the reproductive stage, OoNAC72-OX transgenic lines (OX1, OX2) showed significant delayed bolting and flowering (**Figures 4C,E**). OoNAC72-OX transgenic lines (OX1, OX2) showed no significant difference in leaf number, flowering plant height, dry and wet weight (**Figures 4D–G**). These

findings suggested that OoNAC72 was an important regulator of plant development.

### Overexpression of OoNAC72 Increases Tolerance to Salt and Osmotic Stresses Under Sterile Condition in Arabidopsis

To further investigate mechanisms of hypersensitivity to abiotic stress in OoNAC72-OX plants, we evaluated the stress tolerance of transgenic (OX1 and OX2 lines) and WT Arabidopsis in germination and root growth in seedlings (**Figure 5**). WT and OoNAC72-OX lines were no significant difference on 1/2 MS medium, whereas OX1 and OX2 lines better growth and longer root length under 150 mM NaCl, 200 mM NaCl, 200 mM mannitol and 250 mM mannitol stresses. These results suggested that OoNAC72-OX plants had improved tolerance to salt and drought stresses during seed germination and post-germinative growth periods.

#### Drought Tolerance and Salt Tolerance Phenotypes in Transgenic Plants Under Non-sterile Conditions and Analysis of Physiological Indices

To further elucidate the possible involvement of OoNAC72 in the response to drought and salt, we compared drought and salt tolerance of the OoNAC72-OX lines and WT plants at the vegetative growth stage (**Figure 6**). After 20 days of 150 mM NaCl stress treatment, the phenotypes of the WT and OX1 and OX2 lines began to show significant differences, in which 87.5% of the WT died, while the majority of OX1 and OX2 lines remained green,whereas OX1 and OX2 showed significantly higher survival rate 75 and 62.5% (**Figures 6A,C**). In order to

stresses. Seedlings were grown without water for 35 days and re-watered for 5 days. (C) Survival rate after recovery from drought and salt stress. Each experiment comprises 20 plants. (D–G) Proline and MDA content, SOD and POD activities in WT and OoNAC72-OX lines after drought and salt stress. Each data point is the mean of three biological replicates. Error bars indicate SD, and asterisks indicate a significant difference (∗P < 0.05; ∗∗P < 0.01) compared with WT.

detect the drought resistance of OoNAC72-OX lines, WT and OoNAC72-OX lines were simultaneously withholding water for 35 and 5 days after rehydration, the phenotypes of the WT and OoNAC72-OX lines began to show significant differences, the leaves of WT lines completely wilted and the plants survival rate was 3.33%. However, OX1 and OX2 lines survival rate were 87.5% (**Figures 6B,C**). Furthermore, there were no significant differences in proline and MDA content, POD and SOD activities between WT and OX1 and OX2 lines under normal condition (**Figures 6D–G**). Whereas, under salt and drought stresses, OX1 and OX2 lines exhibitedhigher levels of proline content, POD and SOD activities, and lower MDA level compared with those of WT (**Figures 6D–G**). Phenotypic analysis showed that the OoNAC72-OX plants had the highest drought and salt tolerance. These data demonstrated that the OoNAC72-OX plants exhibited increased tolerance to salt and drought stresses, thus we speculated OoNAC72 plays a critical role in O. ochrocephala response to salt and drought stresses.

#### Drought Resistance of WT and OoNAC72-OX Plants

To further investigate drought sensitivity of the OoNAC72-OX plants, the rate of nature water loss and stomatal apertures of leaves from the 25-day-old soil-grown WT and OoNAC72-OX plants were detected (**Figure 7**). After 3 h of air drying, the

leaves of OX1 and OX2 were slightly curled, while the WT plants were severely curled (**Figure 7A**). Additionally, OoNAC72-OX lines showed lower water loss rate compared with WT plants (**Figure 7B**), indicating that over expression of OoNAC72 had increased water retention capacity in Arabidopsis. After 10 days of drought control, stomatal apertures index of the transgenic and WT lines all appeared significant changes. OX1 and OX2 plants decreased from 0.40 and 0.49 to 0.09 and 0.08, which was significantly smaller than WT (from 0.46 to 0.27) (**Figures 7C,D**). WT and OoNAC72-OX plants were treated with 20% PEG-6000 simulated drought stress. ABA was detected at different time points. The ABA content of these OoNAC72-OX lines were significantly higher than that of the WT (**Supplementary Figure S3**). Those results indicated that OoNAC72-OX plants may reduce the loss of water by regulating the stomata closure to improve drought resistance.

# Altered Expression of Stress-Responsive Genes in Transgenic OoNAC72 Plants

To explore the underlying basis of this phenotype caused by drought or salt stress in the transgenic plants, eight stressresponsive genes: RD29A (Nakashima et al., 2006), RD29B (Nakashima et al., 2006), RD26 (Fujita et al., 2004), LEA14 (Jia et al., 2014), ANACOR19 (Tran et al., 2004), ZAT10 (Mittler et al., 2002), PP2CA (Yoshida and Hirayama, 2006) and NCED3 (Iuchi, 2001) were selected for expression pattern assays (**Figure 8**). Under normal condition, the expression levels of eight genes showed no significant difference between the OoNAC72- OX lines and WT. However, the expression levels of the eight genes in the OX1 and OX2 lines were significantly higher than those in WT plants under salt and drought conditions.

#### Increased ABA Sensitivity in OoNAC72-OX Plants

The ABA sensitivity of the OoNAC72-OX lines was assessed by analyses of seed germination and seedling growth (**Figure 9**). In the absence of ABA, no obvious difference was observed between transgenic and WT plants under normal growing conditions. However, when supplied with 1 or 3 µM ABA, seeds germination and seedling root length of the OX1 and OX2 lines were significantly inhibited compared with the WT lines (**Figures 9A– D**). Moreover, when 40-day-old plants were treated with 10 µM ABA, the stomatal apertures of both wild-type and OoNAC72- OX lines all happened to change (**Figure 9E**). The stomatal apertures index of WT plants decreased from 0.47 to 0.25, while the OX1 and OX2 lines dropped from 0.48 and 0.52 to 0.12 and 0.15 (**Figure 9F**). It indicated that the OoNAC72-OX lines were more rapid and variable in stomatal conductance. These

results indicated that overexpression of OoNAC72 gene led to increased ABA sensitivity, which resulted in retarded growth of transgenic plants.

# DISCUSSION

#### O. ochrocephala Can Be Exploited as a Pasture and Ecological Resource

O. ochrocephala is a perennial poisonous plant widely distributed in harsh environments, and livestock generally do not eat it when edible pasture is relatively abundant. In recent years, due to the emergence of extreme climatic and human overgrazing, O. ochrocephala has spread rapidly and become one of the main poisonous weeds in Western China. Studies have shown that O. ochrocephala is rich in protein and mineral elements and can be used as a good forage forb for livestock after detoxification (Shen and Mo, 2017). At the same time, as a highly resistant plant, O. ochrocephala has important ecological value in wind break and sand fixation, as well as biodiversity maintenance. In view of these characteristics, we focus on O. ochrocephala to elucidate the survival mechanism in adversity, and ultimately to provide a theoretical basis for its development and utilization.

# OoNAC72 Is a Stress-Related NAC Transcription Factor

In the present study, we characterized the OoNAC72 gene, a novel stress-related member of the NAC gene family in O. ochrocephala. OoNAC72 protein contains a typical NAC conserved domain located in the N-terminal region, which can be divided into five subdomains. As the transcriptional regulatory domain, the C-terminus is extremely variable. This result is in agreement with previous reports on other NAC TFs domains, such as PbeNAC1, ZmNAC55, CmNAC1 (Mao et al., 2016; Cao et al., 2017; Jin et al., 2017). Phylogenetic analysis confirmed the OoNAC72 had the highest homology with the two legumes, Medicago truncatula and Cicer arietinum. Subcellular localization and transactivation analysis jointly revealed that OoNAC72 was a transcription factor, which was consistent with previous research (Mao et al., 2016). Therefore, OoNAC72 act as a transcriptional activator, depending on interactions with other transcription factors or conformational changes. These results speculate that OoNAC72 is a novel stress-related member of the NAC gene family in O. ochrocephala.

# OoNAC72 Is Involved in Response to Diverse Stresses

In previous researches, many NAC genes had been identified as an important part of the progress of complex signal transduction when plants were subjected to various stresses, especially abiotic stress such as drought and salt stress (Puranik et al., 2012). In our study, we found that OoNAC72 can be induced by high salt, drought and ABA treatment, whereas not induced by GA and ETH. we speculated that OoNAC72 may respond to stress in O. ochrocephala by ABA-dependent signal transduction pathways.

# OoNAC72 Is Involved in Developmental Processes

Numerous reports have demonstrated that NAC TFs are involved in a number of biological processes, such as regulating the growth of plant cells (Hiroaki et al., 2010), seed development (Sperotto et al., 2009), embryonic development, fiber formation and development, cell differentiation and leaf senescence (Youn-Sung et al., 2007). In our study, we found that OoNAC72-OX (OX1, OX2, OX3) lines showed significant rosette leaves, delayed twitches and flowering compared with WT plants. However, the mechanism resulting from overexpression of OoNAC72 in regulating this phenotype is still unknown and waiting for further research.

#### Overexpression of OoNAC72 Improves Stress Tolerance in Arabidopsis

Overexpression or deletion mutation is the primary way to study the function of NAC transcription factors in adversity stress. In this study, several physiological changes of OoNAC72-OX lines seem to be involved in abiotic stress-resistant and molecular mechanisms. Firstly, in terms of phenotype, OoNAC72-OX lines appeared improved tolerance to salt and drought stresses during seed germination, post-germinative growth periods and vegetative growth stage. Consistent with the results of OoNAC72, CarNAC4-transgenic plants exhibited enhanced drought and salt tolerance than the WT plant, which were strongly demonstrated by both morphological and physiological changes (Yu et al.,

∗P < 0.05, ∗∗P < 0.01.

2016). Secondly, at the physiological and biochemical levels, proline and MDA content, POD and SOD activities in plants overexpressing OoNAC72 have been significantly induced under different stress conditions. Thirdly, at the molecular level, our data also demonstrated that over-expression of OoNAC72 enhanced the expression of stress-responsive genes such as RD29A, RD29B, RD26, LEA14, ANACOR19, ZAT10, PP2CA, and NCED8. Previous studies have confirmed that NAC proteins could bind to NACRS containing the core sequence "CGT[A/G]" to regulate gene expression through ABA dependent pathway for various stress response (Duval et al., 2002; Olsen et al., 2005a; Takasaki et al., 2010; Jensen et al., 2013). However, the genome of O. ochrocephala has not been sequenced and the OoNAC72 promoter sequence could not be cloned, and so it is impossible to predict these cis-acting elements associated with stress in the promoter, and resistance mechanisms need to be further investigated.

#### OoNAC72 Functions in ABA Signaling and Confers Drought Resistance to Transgenic Plants

Abscisic acid is a very important signaling molecule that respond to many adverse environmental stresses such as high salt, drought and extreme temperatures (Agarwal and Jha, 2010; Fujita et al., 2011). In plants, higher ABA sensitivity may irritate the stomata to maintain moisture and enhance drought resistance. In our study, ABA content of WT and OoNAC72-OX plants treated with 20% PEG-6000 simulating drought stress showed that ABA content of the OoNAC72-OX lines was significantly higher than that of the WT (**Supplementary Figure S3**). Furthermore, we found that the stomata in Arabidopsis overexpressing OoNAC72 had different degrees of closure after drought stress, indirectly mapping the relationship between drought stress and the ABA pathway. Similar results were found in Arabidopsis and maize (Ping-Li et al., 2007; Mao et al., 2016). Overexpression of KUP6 (K+ uptake transporter 6) in Arabidopsis increased the sensitivity of ABA to drought stress through faster stomatal closure, thereby improving the tolerance of transgenic plants to drought stress (Yuriko et al., 2013). Genes such as PP2Cs, SnRK2s, NAC and WRKY played a vital role in ABA signaling pathways induced by drought stress (Kim, 2014). Therefore, we supposed that OoNAC72 may enhance the drought resistance of transgenic Arabidopsis by participating in the ABA pathway.

# OoNAC72 Scavenges ROS Capability by Increasing ABA Content Under Stress

Under stress conditions, plants reduce the damage caused by stress-induced ROS by enhancing their antioxidant defense system (Oracz et al., 2009; Suzuki et al., 2012). MDA is one of the important parameters to measure the degree of oxidative damage in plant cells. While pro contributes to osmotic adjustment to effectively enhance the antioxidant system and reduce peroxidative damage (Zhuang et al., 2015a). In our study, transgenic plants produced less MDA and more pro under stress than wild-type plants, revealing that OoNAC72 overexpression enhanced the body's regulation of feedback regulatory substances.

Under drought stress, ABA acts as an upstream signaling of NO and they form cross-signal pathways co-regulate a balance between ROS (H2O2) and NO production (Zhuang et al., 2015b). Our resistance assays showed that the ABA contents, and SOD and POD activities were significantly induced by OoNAC72 under drought stresses condition. The increase in ABA content leads to an increase in the concentration of NO for activating NO signaling pathway. On the one hand, the low concentration of NO interacts with H2O2 to stimulate the antioxidant defense system. On the other hand, the increase of SOD and POD activity further enhances plant's ability to scavenge excessive H2O2. The synergy between NO signal and ROS signal promotes the dynamic balance between NO and H2O2, and ultimately improves plant resistance tolerance (Miller et al., 2010). These results indicate that showing the positive regulation of OoNAC72 on ABA signaling pathway under drought stress.

# CONCLUSION

Plant NAC transcription factors control diverse biological processes, such as differentiation, development and abiotic stress responses. In this study, we identified a gene encoding NAC72 type transcription factor from Oxytropis ochrocephala, and characterized its role. The nuclear localization and transcriptional activity of this protein indicate the possible function of OoNAC72 as a transcription factor. In addition, expression of OoNAC72 transcript was shown to be up-regulated in response to abiotic stresses and exogenous ABA. Furthermore, analyses of transgenic Arabidopsis expressing OoNAC72 also supported the involvement of OoNAC72 in drought and salinity responses as well as in the regulation of ABA-dependent processes.

# AUTHOR CONTRIBUTIONS

HG and YF conceived and designed the study. HG, XL, and FN performed the experiments. HG, XL, and QZ analyzed the data and wrote the manuscript. FN, QZ, NF, DC, and DM reviewed the manuscript. WH, BG, YW, and YF contributed reagents, materials, and fund support.

# FUNDING

This work was supported by Special Fund for Agro-scientific Research in The Public Interest (201203062), Key Research and Development Plan Project of Shaanxi Province (2018ZDXM-SF-016) and Opening Foundation of Key Laboratory of Resource Biology and Biotechnology in Western China (Northwest University), Ministry of Education (ZSK2016003) to YW, Bureau of Science and Technology in Shaanxi Province, Xi'an City, Beilin District (GX1702) to BG, Natural Science Foundation of Shanxi Province of China (2018JQ3029) and Natural Science Foundation of China (31702159) to YF, Key Research and Development Plan Project of Shaanxi Province (2018SF-034) and Scientific Research from Shaanxi Provincial Department of Education (17JK0243) to NF.

### ACKNOWLEDGMENTS

fpls-10-00890 July 10, 2019 Time: 17:16 # 13

We thank the reviewers for their constructive suggestions. We are also particularly grateful to other members of our research group for their helpful comments to improve the article.

#### REFERENCES

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00890/ full#supplementary-material

gene NCED3 in Arabidopsis thaliana. FEBS Open Bio. 3:321. doi: 10.1016/j.fob. 2013.07.006

transgenic Tamarix and Arabidopsis. Front. Plant Sci. 8:635. doi: 10.3389/fpls. 2017.00635

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Guan, Liu, Niu, Zhao, Fan, Cao, Meng, He, Guo, Wei and Fu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Roles of GmERF135 in Improving Salt Tolerance and Decreasing ABA Sensitivity in Soybean

Meng-Jie Zhao1† , Li-Juan Yin1†‡, Jian Ma2† , Jia-Cheng Zheng3 , Yan-Xia Wang4 , Jin-Hao Lan5 , Jin-Dong Fu1 , Ming Chen1 , Zhao-Shi Xu1,2,3 \* and You-Zhi Ma1 \*

#### Edited by:

1

Yan Guo, China Agricultural University, China

#### Reviewed by:

Xiuli Hu, Henan Agricultural University, China Yangrong Cao, Huazhong Agricultural University, China

#### \*Correspondence:

Zhao-Shi Xu<EMAIL_ADDRESS>You-Zhi Ma<EMAIL_ADDRESS>†These authors have contributed

equally to this work ‡Present address:

Li-Juan Yin, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China

#### Specialty section:

This article was submitted to Plant Abiotic Stress, a section of the journal Frontiers in Plant Science

Received: 22 November 2018 Accepted: 04 July 2019 Published: 23 July 2019

#### Citation:

Zhao M-J, Yin L-J, Ma J, Zheng J-C, Wang Y-X, Lan J-H, Fu J-D, Chen M, Xu Z-S and Ma Y-Z (2019) The Roles of GmERF135 in Improving Salt Tolerance and Decreasing ABA Sensitivity in Soybean. Front. Plant Sci. 10:940. doi: 10.3389/fpls.2019.00940

Institute of Crop Science, Chinese Academy of Agricultural Sciences (CAAS)/National Key Facility for Crop Gene Resources and Genetic Improvement, Key Laboratory of Biology and Genetic Improvement of Triticeae Crops, Ministry of Agriculture, Beijing, China, 2 Department of Agronomy, Jilin Agricultural University, Changchun, China, 3 College of Agriculture, Anhui University of Science and Technology, Fengyang County, China, 4 Hebei Academy of Agriculture and Forestry Sciences, Research Center of Wheat Engineering Technology of Hebei, Shijiazhuang, China, 5 College of Agronomy and Plant Protection, Qingdao Agricultural University, Qingdao, China

Abscisic acid (ABA) mediates various abiotic stress responses, and ethylene responsive factors (ERFs) play vital role in resisting stresses, but the interaction of these molecular mechanisms remains elusive. In this study, we identified an ABA-induced soybean ERF gene GmERF135 that was highly up-regulated by ethylene (ET), drought, salt, and low temperature treatments. Subcellular localization assay showed that the GmERF135 protein was targeted to the nucleus. Promoter cis-acting elements analysis suggested that numerous potential stress responsive cis-elements were distributed in the promoter region of GmERF135, including ABA-, light-, ET-, gibberellin (GA)-, and methyl jasmonate (MeJA)-responsive elements. Overexpression of GmERF135 in Arabidopsis enhanced tolerance to drought and salt conditions. In addition, GmERF135 promoted the growth of transgenic hairy roots under salt and exogenous ABA conditions. These results suggest that soybean GmERF135 may participate in both ABA and ET signaling pathways to regulate the responses to multiple stresses.

Keywords: ABA, ethylene-responsive factor, hypocotyl elongation, root growth, response mechanism, salt tolerance, soybean

# INTRODUCTION

As a vital hormone in plants, abscisic acid (ABA) is essential to a myriad of aspects of plant growth and developmental processes, including plant embryo development, seed maturation, fruit maturity, and stomatal movement (Finkelstein et al., 2002). ABA signal transduction has been studied for many years (Pandey et al., 2006, 2009; Shen et al., 2006; Ma et al., 2009; Park et al., 2009; Yang and Guo, 2018), and the widely accepted molecular mechanism is that the pyrabactin resistant/PYR-like/regulatory component of ABA receptor (PYR/PYL/RCAR) acts as an ABA receptor which can bind to ABA, and then binds to and inhibits protein phosphatases type 2C (PP2Cs) (Santiago et al., 2009). In addition, the activity of SNF1-related protein kinases 2 (SnRK2s) is enhanced and can phosphorylate ABRE binding factors (AREB/ABFs) to induce physiological and biochemical changes in the process of ABA response (Hubbard et al., 2010; Yoshida et al., 2015). A recent study showed that ABFs can directly bind to the promoters of group A PP2C genes, and

rapidly induce their expression on exogenous ABA treatments (Wang et al., 2018). However, the downstream molecular mechanism is yet not clearly understood.

Many transcription factors, such as MYC/MYB, bZIP/ABRE, AP2/ERF, and NAC, are regulated by ABA. The AP2/ERF family, the largest plant transcription factor family (Okamuro et al., 1997), can be divided into three subfamilies: the AP2, ERF, and RAV subfamilies. Among them, the ERF subfamily can specifically bind to GCC-box and/or the dehydrationresponsive element/C-repeat (DRE/CRT) cis-acting elements (Allen et al., 1998).

The AP2/ERF family is involved in responses to various abiotic stresses and exogenous hormones (Xu et al., 2008, 2011; Klay et al., 2018). For example, overexpressing TaERF1 enhances tolerance of physiological and environmental stresses, such as salt, drought, low temperature, exogenous ABA, ethylene (ET), salicylic acid (SA), and disease (Xu et al., 2007). Transgenic rice plants expressing OsERF922 displayed higher susceptibility to Magnaporthe oryzae and NaCl compared to the wild type (WT), while the knockout mutant and RNAi lines enhanced resistance to these stresses (Liu et al., 2012). A recent study described how ABA modulates the expression level of ERF family, showing that the expression levels of ABA-responsive genes such as RD22, LEA3, and PODs were up-regulated after overexpressing OsERF101 in rice, which enhanced its susceptibility to ABA (Jin et al., 2018).

Transcriptomic analysis of grapevine organs treated with or without ABA showed that ERF members were involved in differently expressed genes (DEGs), and the ERF subfamily had the most significant change compared to other transcription factors (Rattanakon et al., 2016). Although the ABA signal pathway has been extensively studied, it is yet unclear how regulation of the expression of downstream genes via ERF subfamily could enhance abiotic stress tolerances in soybean. To investigate whether the ERF subfamily is modulated by the ABA signal in soybean, we studied the response to exogenous ABA and identified the functions of GmERF135 in transgenic Arabidopsis and soybean hairy roots.

### MATERIALS AND METHODS

#### Plant Materials and Treatment

Soybean cultivar "Tiefeng 8" was sown in pots containing vermiculite and grown at 25◦C for 14 days and then treated with various abiotic stresses. For the various treatments, soybean plants were exposed in the air for rapid drought, dipped into 200 mM NaCl for salt stress, placed in 42◦C/4◦C chamber for high/low temperature stress, respectively. They were also placed in an airtight container filled with ET, or dipped in 100 µM ABA, 50 µM salicylic acid (SA), or 50 µM jasmonate (JA) for exogenous hormone stresses. Leaves of these plants were collected after 0, 0.5, 1, 2, 5, 12, and 24 h treatment and then immediately stored at −80◦C for RNA extraction.

# Gene Structure and Protein Domain of ERFs in Soybean

The whole genome data of the candidate ERF genes was obtained from JGI Glyma1.0 annotation1 (Goodstein et al., 2012). Gene structure was analyzed by submitting the CDS of the candidate genes and whole genomic DNA sequences to the Gene Structure Display Server (GSDS) website2 . Protein Fold Recognition Server (PHYRE2)3 was used for analysis of structural homology modeling of these genes, and DOG 2.0 was used to draw the protein domains.

# Quantitative Real-Time PCR (qRT-PCR)

Total RNA of the soybean plants was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The specific primer pairs for the 16 genes were designed by Primer Premier 5.0 according to the cDNAs. qRT-PCR was conducted using the SYBR Premix Ex TaqTM kit (TAKARA, Kusatsu, Japan) and the ABI Prism 7500 real-time PCR system (Thermo Fisher Scientific, Waltham, MA, United States). The 2−11CT method was used to conduct qRT-PCR analysis (Le et al., 2011). Soybean Actin (U60506) was used as an internal control and sequence of specific primers were shown in **Supplementary Table S1**.

#### Isolation and Promoter Analysis of GmERF135

The GmERF135 gene was amplified by PCR and the primers for cloning the GmERF135 were 50 - AATCATTATGTGTGGCGGTGCC-30 and 50 -TATTCCTCGC TAATCGAAACTCCAGAG-30 . The PCR product was then cloned into the pEASY-T1 vector (TransGen, Beijing, China). 1,886 bp promotor region of GmERF135 was cloned and submitted to PLACE4 and PlantCARE5 databases to analyze the putative cis-acting elements in the promoter region (Lescot et al., 2002).

# Subcellular Localization Analysis

The cDNA of GmERF135 were augmented with PCR, connected to the N-terminus of humanized green fluorescent protein (hGFP) reporter gene under the control of the double Cauliflower Mosaic Virus (2 × CaMV) 35S promoter, and a recombinant plasmid was obtained (Scott et al., 1999). The recombinant plasmid was introduced into onion epidermal cells, while the onion epidermal cells with hGFP vector acted as the control. Fluorescence microscopy was used to identify hGFP expression (Xu et al., 2007; He et al., 2016).

# Generation of Transgenic Arabidopsis and Stress Treatments

AT5G47230 is the orthologue of GmERF135 in Arabidopsis, which was named ERF2 and used to investigate function in

2http://gsds.cbi.pku.edu.cn/

1https://phytozome.jgi.doe.gov/pz/portal.html

3http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index

4http://www.dna.affrc.go.jp/PLACE/

5http://bioinformatics.psb.ugent.be/webtools/plantcare/html/

responses to various stresses. The seeds of two mutants, erf2- 1 and erf2-2 (SALK\_126889, SALK\_076967) were mutated via T-DNA insertion.

The cDNA of GmERF135 was obtained by using the specific primer pairs: 50 -TGATTACGCCAAGCTTATGTGTG GCGGTGCC-30 , 50 -CCGGGGATCCTCTAGAATCGAAACTC CAGAG-30 , and then were cloned into pBI121 under the control of the CaMV 35S promoter. The recombinant plasmid 35S::GmERF135 was sequenced and transformed into wild-type (WT) and the two mutants Arabidopsis lines using the vacuum infiltration method (Bechtold and Pelletier, 1998; Liu et al., 2013). T3 seeds of transgenic lines were selected for further analysis.

For root length assay, 30 seeds of each line were sown on 1/2 strength MS growth medium with or without 6% PEG, 75 mM NaCl, 1-aminocyclopropane-1-carboxylic acid (ACC), or in dark treatment for growth, respectively. At least 30 seedlings per line were randomly selected to measure root length. The cotyledon pieces of each line were recorded every 12 h. Each treatment contained three independent replicates.

### Soybean Hairy Root Induction and Stress Treatments

The Superroot of Lotus corniculatus and Cucumopine-type Agrobacterium rhizogene strain K599 with pGFPGUSPlus were provided by Professor Tian-Fu Han (CAAS, China). Seedling growth, rooting, hairy root induction, and hairy root transformation were performed as described by Chen et al. (2004, 2014)

GFP positive (GFP+) hairy roots were induced by K599 carrying the pGFPGUSPlus-GmERF135 binary vector. These hairy roots were cultured on 1/2 MS medium supplemented with 0, 50, 85, 120, or 150 mM NaCl for salt treatment, and 50, 100, or 150 µM ABA for hormone treatment, respectively. They were then incubated at 24◦C under a 16/8 h light/dark cycle condition for 2 weeks. After 24 h incubation at 105◦C, the increase in dry weight (30 roots per unit) was measured and recorded. Each treatment contained three independent replicates.

# RESULTS

### Molecular Characterization of Soybean Targeted ERF Genes

In a previous study, 160 non-redundant soybean ERFs were identified using the Pfam database and SMART program, and these soybean ERFs were clustered into eight groups (**Supplementary Table S2**; Zhao et al., unpublished). To comprehensively understand the responses of soybean ERFs to ABA, 16 ERF genes were selected for further investigation according to the phylogenetic tree. Gene structure analysis showed that six genes had no introns, including GmERF106, GmERF132, GmERF135, GmERF41, GmERF49, and GmERF84 (**Figure 1**). The remaining 10 ERF genes contained one intron, which was distributed in each cDNA region except for GmERF103, GmERF111, and GmERF15. The protein domains of each soybean ERF gene were drawn by DOG2.0 (**Figure 2**). All the proteins contained a conservative AP2/ERF domain, which was distributed in different positions of each protein sequence. The AP2/ERF domain of the proteins in each group displayed similar locations such as Group I, III, IV, V, VII, and VIII.

An expression pattern map of soybean ERFs was drawn based on the gene-chip data downloaded from the soybean genome database6 . As shown in **Supplementary Figure S1A**, the soybean ERFs were expressed at different levels in various tissues and organs. Among them, GmERF111 and GmERF135 showed high expression levels in almost all tissues and organs. GmERF75 showed high expression levels in root and nodule, and GmERF49 was highly expressed in nodule and 10 DAF pod shell. GmERF15 mainly expressed in young leaf. All the candidate genes had different expression patterns.

#### Isolation and Subcellular Localization Characterization of GmERF135

Abscisic acid plays an important role in plant growth and development which is closely related to quality and yield in plants. To investigate how ABA affected the expression pattern of the 16 soybean ERF genes, qRT-PCR was conducted. Almost all the ERFs were up-regulated by ABA, except for GmERF49 (**Figure 3**). The ERF gene GmERF135 was the highest expressed gene after exogenous ABA treatment except for GmERF75 which has been studied (Zhao et al., unpublished). The transcript level of GmERF135 was activated and showed a 17-fold increase within 2 h after treatment and then declined to normal. Considering the high expression level and ABA-responsive increase in plants, GmERF135 was selected for further study.

The subcellular localization assay was conducted to provide clues to understand intrinsic characteristics of cell activities. hGFP reporter gene fused to C-terminus of GmERF135, which could fluoresce when laser irradiated. The GmERF135::hGFP fused protein only fluoresced in nucleus, while the control fluoresced throughout the entire cell (**Figure 4**). This result showed that GmERF135 functions in nucleus.

# GmERF135 Promoter Region Comprises Diverse Stress-Responsive Elements

The promoter region is an important part of the gene which could regulate gene expression and control gene action. To investigate the potential regulation mechanism, the promoter region of GmERF135, which has 1886 bp length upstream of the start codon, was isolated. The PLACE and PlantCARE databases were used to analyze the putative regulatory elements in the promoter region. Several regulatory elements were identified to be involved in responses to abiotic stresses and plant hormones (**Table 1**).

The GmERF135 promoter region contains many ABA and stress responsive elements, including ABA-responsive elements (ABREs, 3 hits), MYBST1 core sequences (4 hits), MYB binding sites (MBS, 1 hit), and DPBF binding sites (2 hits) **(Table 1**). Except for ABREs, three hormone-responsive elements were predicted, including an ethylene- responsive element (ERE, 1 hit), MeJA-responsive element TGACGmotif (1 hit) and CGTCA-motif (1 hit), and GA-responsive

6https://soybase.org/soyseq

blue boxes represent introns, and black lines represent untranslated regions (UTRs).

element (GARE, 2 hits). In addition, TCT-motif (1 hit), G-box (2 hits), and Box 4 (1 hit) were found, which could act as light-responsive elements in the promoter region of GmERF135. These elements suggested that GmERF135 may be involved in responses to multiple abiotic stresses and exogenous hormones.

#### The Roles of GmERF135 in Responding to Multiple Stimuli

Quantitative Real-Time PCR was conducted to investigate the expression level of GmERF135 under abiotic stresses, such as drought, salt, high/low temperature, and exogenous hormones including ET, SA, and JA. GmERF135 was induced by almost all the abiotic stresses and exogenous hormones. Remarkably, the expression level of the gene was extremely up-regulated by drought, low temperature, and ET treatment. The transcription of GmERF135 peaked at 5 h under drought or low-temperature treatment which had 18-fold / 22-fold increases, respectively (**Figure 5**). The peak of GmERF135 transcription appeared at 12 h after treatment with ET, which showed a 28-fold increase. In contrast, GmERF135 rapidly reached the maximum transcript level within 0.5 h under the NaCl and JA treatments, and then declined to baseline after 1 h. Under the SA treatment, GmERF135 was

induced and reached the peak within 2 h (about 7-fold) and then declined to normal level in 24 h. These results showed that GmERF135 may be involved in responses to multiple stimuli.

#### GmERF135 Rescued Two erf2 Mutants and Affected Growth of the Root and Cotyledon

Total RNA was extracted for semi-quantitative PCR from hypocotyls, root, stem, and leaf tissues of normal growth soybean seedlings. Actin primers were used as a parallel reaction to normalize the added template amounts. GmERF135 was predominantly expressed in the leaf, less in the hypocotyl and stem, and very little in the root (**Supplementary Figure S1B**).

To assess whether GmERF135 could rescue the two mutants, the seeds of transgenic GmERF135::erf2 lines, WT, and two erf2 mutants were sown on the 1/2 MS medium for growth and the plant growth rate was recorded each 12 h. Six day later, the phenotypes were imaged (**Figures 6A,B**).

Growth rate assay showed GmERF135::erf2 lines displayed faster growth rate compared to the mutants **(Figure 6C**). WT plants required 2 d of growth to sprout the two pieces of the cotyledon while the two mutants required 4 days. After introduction of GmERF135 in erf2 lines, 3 days was required for the same process. The two mutants needed 6 days to reach 4-leaf stage, while the transgenic GmERF135::erf2 lines had similar growth rate with WT, which only needed 5 d. Root length assay showed there is an approximately 25% decrease in erf2 mutants compared to WT lines **(Figure 6D**). After overexpressing GmERF135 in the erf2 mutants, the lines had similar root lengths to WT and the phenotype of which was rescued. These results showed that overexpression of GmERF135 in the two erf2 mutants partly rescued two erf2 mutants and affected growth of root and cotyledon in Arabidopsis.

#### GmERF135 Promotes Plant Growth Under Drought and Salt Stresses in Arabidopsis

To determine whether GmERF135 confers abiotic stress tolerance to Arabidopsis plants, 3-day-old WT, erf2 mutants, and GmERF135 overexpression seedlings were transferred to 1/2 MS medium with or without 6% PEG, 75 mM NaCl, ACC, and dark environment grown for 3 days. The growth rate and the root length of 30 seedlings in each line were recorded (**Figures 7B,C**). Statistics showed that the average root length of erf2 was shorter than the WT and transgenic Arabidopsis plants under all conditions except for ACC treatment (**Figures 7A,B**). The transgenic plants displayed larger cotyledon than the WT and mutants under NaCl and PEG treatments, which suggested that GmERF135 enhanced the growth rate of Arabidopsis plants. Overexpression of GmERF135 in Arabidopsis resulted in a longer root compared to the WT and two erf2 mutants under 6% PEG, 75 mM NaCl, and dark conditions. No detectable difference was observed between the transgenic lines and WT under ACC conditions.

To further understood the response to various stresses in Arabidopsis, some marker genes, such as NCED3 (Ruggiero et al., 2004), RD29A (Yamaguchi-Shinozaki and Shinozaki, 1994), COR15A (Xu et al., 2011), DREB2A (Liu et al., 1998), ABA1, ABA2 (He et al., 2016), ACO4, and ACS2 (Kim et al., 2013) were selected for qRT-PCR (**Supplementary Figure S2**).

FIGURE 5 | Expression patterns of GmERF135 under multiple stimuli. Two-week soybean seedlings were given various abiotic stresses and exogenous hormones for 0, 0.5, 1, 2, 5, 12, and 24 h, which were used to extract RNA and obtained cDNA. The transcript levels of GmERF135 after multiple stresses were quantified by qRT-PCR, including drought, NaCl, high/low temperature, ET (ethylene), JA (jasmonate), and SA (salicylic acid). Data were shown referring to three biological replicates.

#### GmERF135 Improves Tolerance of Salt and ABA in Soybean

To investigate the roles of GmERF135 in soybean, the pGUS-GmERF135 expression vector was constructed and transformed into soybean hairy roots, which were studied on MS basal medium supplemented with different concentrations of ABA and NaCl. GmERF135 transgenic soybean hairy roots experienced greater growth under the treatments compared to the vector control, especially under the 100 µM ABA and 85 mM NaCl conditions (**Figure 8A**).

Dry weight measurement results showed that GmERF135 hairy roots increased by about 10% compared to the vector under the MS0 condition, while transgenic hairy roots displayed a 6.24–75% improvement compared to WT under NaCl and exogenous ABA at the different concentrations (**Figure 8B**). For NaCl treatment, there has a 6.24–75.03% increase in the average dry weight of hairy roots overexpressing GmERF135. For exogenous ABA, the transgenic hairy roots displayed a 36– 65% improvement compared to vector control under exogenous ABA treatment at different concentrations (**Figure 8B**). It is noteworthy that a total of 64.46% improvement in GmERF135 hairy roots was observed under the 100 µM ABA condition, and an improvement of 75.03% with 85 mM NaCl (**Figure 8B**). Results of dry weight measurements confirmed the conclusion that GmERF135 could promote plant growth under ABA and NaCl conditions.

# DISCUSSION

Abscisic acid is well-known to be involved in responses to multiple stresses, such as drought, salt, and cold, and to induce expression of stress-related genes in plants (Leung et al., 1997;

FIGURE 7 | GmERF135 transgenic Arabidopsis lines improved cotyledon growth under the NaCl and PEG conditions. (A) Phenotype of GmERF135 transgenic Arabidopsis lines under various stresses. Three-day-old GmERF135 transgenic Arabidopsis lines were transferred to 1/2 MS medium with 6% PEG, 75 mM NaCl, ACC, or grown in darkness. (B) Average root lengths of WT, erf2, and 35S::GmERF135 lines. Three-day-old seedlings were treated with various stresses for 3 days and the data were shown as the means ± SD of three biological replicates (at least 30 individuals per line). (C) Days needed for growing cotyledon. Statistic data were recorded every 12 h. Three biological replicates were performed.

Rohde et al., 2000; Rook et al., 2001; Xu et al., 2007). The interaction of the ABA and ET signal pathway is extremely complex and intricate. AREB/ABFs, members of the bZIP family of transcription factors, participate in the ABA signal pathway which can specifically bind to ABRE cis-element to modulate the expression of downstream target genes (Bonetta and McCourt, 1998; Leung and Giraudat, 1998; Li et al., 2012). Fujita et al. (2005) observed that AREB1 can act as a transcription activator in ABRE-dependent ABA signaling, which enhances drought tolerance in Arabidopsis. Overexpression of the AREB1 gene in Arachis hypogaea could enhance drought tolerance by modulating ROS scavenging and maintaining content (Li et al., 2013). Phosphorylated AREB/ABFs may bind to the ABRE ciselement of GmERF135 promotor region to activate its transcription level (Fujita et al., 2013), which could activate or repress the transcription of targeted genes. In our study, analysis of the GmERF135 promoter region showed that three ABREs were located at the promoter region, which may be binding sites of AREB/ABFs to affect the transcript of GmERF135. GmERF135 has a conversed AP2/ERF domain, which could specifically bind to GCC-box and/or dehydrationresponsive element/C-repeat (DRE/CRT) cis-acting elements to modulate the expression of pathogenesis- and abiotic stressesrelated genes (Hao et al., 1998; Chakravarthy et al., 2003; Zarei et al., 2011).

Previous studies showed there is an antagonistic interaction between ABA and JA-ET signaling pathways ABA for disease response (Beaudoin et al., 2000; Anderson et al., 2004; Martin-Rodriguez et al., 2011). A recent study showed that SlAREB1 could enhance the expression of ethylene biosynthetic genes ACS and ACO in tomato fruits, which are the two key genes of ET synthesis (Mou et al., 2018). In our study, qRT-PCR analysis showed that the transcription of both AtACO4 and AtACS2 increased (**Supplementary Figure S2**) in transgenic Arabidopsis lines, which suggested acceleration of ET production. ET acceleration could trigger a series of reactions of ET and activate the transcript of ERF responsive genes (Guo and Ecker, 2004). At the same time, the transcript levels of ABA1 and ABA2, key factors of ABA synthesis, were also upregulated after overexpressing GmERF135 in Arabidopsis lines. These results suggested that GmERF135 may participate in both ET and ABA signaling pathways, and the regulation between the two signaling pathways needs further research.

Except for ABREs, other stress-related elements were also distributed in the promotor region of GmERF135, such as the MYBST1 core sequence, MBS, and ERE (**Table 1**). The MYBST1 core sequence and MBS elements have been demonstrated to be involved in drought, low temperature, salt, and ABA stress responses (Abe et al., 2003). qRT-PCR showed that GmERF135 is induced by multiple stresses, including drought, salt, low temperature, exogenous ET, SA, and JA. These changes may be caused by some corresponding cis-elements in promoter regions which can be bound by specific transcription factors (Huang et al., 2012; Buscaill and Rivas, 2014). These results suggest that soybean ERF gene GmERF135 is a key factor which participates in multiple signaling pathways to regulate expression levels of stress-related genes.

# AUTHOR CONTRIBUTIONS

Z-SX and Y-ZM coordinated the project, conceived and designed experiments, and edited the manuscript. M-JZ and L-JY conducted the bioinformatic work, performed experiments and wrote the first draft. JM revised the manuscript and figures. J-CZ conducted the bioinformatic work. Y-XW, J-HL, and J-DF contributed with valuable discussions. MC provided analytical tools and managed reagents. All authors have read and approved the final manuscript.

# FUNDING

This research was supported by the National Natural Science Foundation of China (31871624), the National Transgenic Key Project of Ministry of Agriculture (2018ZX0800909B), and the Technological Innovation Projects of Modern Agriculture of Hebei Province (494-0402-YBN-RDC4).

# ACKNOWLEDGMENTS

We thank Professor Tian-Fu Han for providing the Superroot culture of L. corniculatus, A. rhizogenes strain K599, and the binary vector pGFPGUSPlus. We also thank Dr. Wen-Sheng Hou for plant material preparation, technical assistance, and soybean seeds.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00940/ full#supplementary-material

FIGURE S1 | Expression patterns of candidate soybean ERFs in different organs. (A) Expression patterns of candidate soybean ERFs in different organs. Normalized expression data for the soybean ERF genes were collected from the SoyBase (http://www.soybase.org/). Expression (vertical coordinates) is in transcripts per million (TPM). (B) Semi-quantitative PCR of GmERF135 in different organs.

FIGURE S2 | Express levels of various stresses-related genes after treatment. Two-week-old WT and transgenic Arabidopsis lines were used to extract total RNA. GmActin/AtActin was used for normalization. Data were shown as the means ± SDs of three experiments.

TABLE S1 | Specific primers of each gene for qRT-PCR. All the primers were designed via Primer Primer 5.0.

TABLE S2 | Classification of 160 soybean ERF genes.

#### REFERENCES

fpls-10-00940 July 19, 2019 Time: 15:31 # 10

response in reproductive tissues. Plant Mol. Biol. 98, 51–65. doi: 10.1007/ s11103-018-0762-5

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Zhao, Yin, Ma, Zheng, Wang, Lan, Fu, Chen, Xu and Ma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Transcriptome Analysis of the Hierarchical Response of Histone Deacetylase Proteins That Respond in an Antagonistic Manner to Salinity Stress

*Minoru Ueda1,2,3, Akihiro Matsui1,3, Shunsuke Watanabe4, Makoto Kobayashi5, Kazuki Saito5,6, Maho Tanaka1,3, Junko Ishida1,3, Miyako Kusano5,7, Mitsunori Seo4 and Motoaki Seki1,2,3,8\**

*1 Plant Genomic Network Research Team, RIKEN Center for Sustainable Resource Science, Yokohama, Japan, 2 Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST), Kawaguchi, Japan, 3 Plant Epigenome Regulation Laboratory, RIKEN Cluster for Pioneering Research, Wako, Japan, 4 Dormancy and Adaptation Research Unit, RIKEN Center for Sustainable Resource Science, Yokohama, Japan, 5 Metabolomics Research Group, RIKEN Center for Sustainable Resource Science, Yokohama, Japan, 6 Graduate School of Pharmaceutical Sciences, Chiba University, Chiba, Japan, 7 Graduate School of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Japan, 8 Kihara Institute for Biological Research, Yokohama City University, Yokohama, Japan*

#### *Edited by:*

*Vicent Arbona, University of Jaume I, Spain*

#### *Reviewed by:*

*Keqiang Wu, National Taiwan University, Taiwan Giorgio Perrella, ENEA - Centro Ricerche Trisaia, Italy*

*\*Correspondence:*

*Motoaki Seki<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Plant Abiotic Stress, a section of the journal Frontiers in Plant Science*

*Received: 18 January 2019 Accepted: 23 September 2019 Published: 18 October 2019*

#### *Citation:*

*Ueda M, Matsui A, Watanabe S, Kobayashi M, Saito K, Tanaka M, Ishida J, Kusano M, Seo M and Seki M (2019) Transcriptome Analysis of the Hierarchical Response of Histone Deacetylase Proteins That Respond in an Antagonistic Manner to Salinity Stress. Front. Plant Sci. 10:1323. doi: 10.3389/fpls.2019.01323*

Acetylation in histone and non-histone proteins is balanced by histone acetyltransferase and histone deacetylase (HDAC) enzymatic activity, an essential aspect of fine-tuning plant response to environmental stresses. HDACs in *Arabidopsis* are composed of three families (RPD3-like, SIRT, and HD-tuins). A previous study indicated that class I (HDA19) and class II (HDA5/14/15/18) RPD3-like family HDACs control positive and negative responses to salinity stress, respectively. Furthermore, quintuple *hda5/14/15/18/19* mutants (*quint*) exhibit salinity stress tolerance, suggesting that *hda19* suppresses the sensitivity to salinity stress present in quadruple *hda5/14/15/18* mutants (*quad*). In the present study, transcriptome analysis of the *quint* mutant was conducted to elucidate the hierarchical control of salinity stress response operated by RPD3-like family HDACs (HDA5/14/15/18/19). The analysis identified 4,832 salt-responsive genes in wild-type (Col-0), *hda19-3*, *quad*, and *quint* plants and revealed that 56.7% of the salt-responsive genes exhibited a similar expression pattern in both the *hda19-3* and *quint* plants. These results indicate that deficiency in HDA19 has a bigger impact on salinity stress response

than in class II HDACs. Furthermore, the expression pattern of genes encoding enzymes that metabolize phytohormones raises the possibility that a drastic change in the homeostasis of phytohormones, such as abscisic acid, brassinosteroid, and gibberellin, may contribute to increasing stress tolerance in *hda19-3* and *quint* plants. Among these phytohormones, abscisic acid accumulation actually increased in *hda19-3* and *quint* plants, and decreased in *quad*, compared with wild-type plants. Importantly, 7.8% of the salt-responsive genes in *quint* plants exhibited a similar expression pattern in *quad* plants, suggesting that some gene sets are regulated in an HDA5/14/15/18-dependent manner.

**324**

The transcriptome analysis conducted in the present study revealed the hierarchical and independent regulation of salt stress response that is mediated through HDA19 and class II HDACs.

Keywords: high salinity stress, epigenetics, histone acetylation, histone deacetylases, stress response

#### INTRODUCTION

Recent studies have demonstrated that a diverse type of histone modifications, such as acetylation, methylation, and phosphorylation, play a pivotal role in orchestrating plant response to environmental stresses (Kim et al., 2015). Acetylation levels are balanced by histone acetyltransferases (HATs) and histone deacetylases (HDACs). These enzymes have the ability to write or erase an acetylation mark, respectively. In general, the mark is positively correlated with mRNA expression. Notably, the mechanism by which histone acetylation governs plant response to environmental stresses through HDACs has been elucidated (Asensi-Fabado et al., 2016; Luo et al., 2017; Ueda et al., 2018b).

Three families of HDAC proteins (Reduced potassium dependency 3 (RPD3)-like, Silent Information Regulator 2 (SIRT), and HD-tuins) have been recognized in plants. The *Arabidopsis thaliana* genome encodes 18 genes that are members of the three HDAC families: 12 RPD3-like family proteins; 2 SIRT family proteins; and 4 HD-tuin family proteins. The RPD3-like family is further divided into three classes (I, II, and IV) based on their homology to yeast HDACs (Hollender and Liu, 2008). The RPD3-like family HDACs are pharmacological targets because their inhibition has potential application as a cancer treatment (Bolden et al., 2006; Seto and Yoshida, 2014). The potential value of their inhibition not only applies to cancer therapy but may also be applicable for increasing environmental stress tolerance in plants. *Arabidopsis* and cassava plants treated with HDAC inhibitors exhibit an increase in tolerance to salinity stress (Sako et al., 2016; Patanun et al., 2017; Nguyen et al., 2018). Pharmacological analysis indicated that inhibition of class I RPD3-like HDACs is essential for increasing salinity tolerance based on the evaluation of the survival of plants treated with class-selective inhibitors and subjected to salinity stress conditions (Ueda et al., 2017).

Consistent with the pharmacological studies, *Arabidopsis* plants deficient in class I RPD3-like HDACs (HDA6, HDA9, and HDA19) exhibit increased tolerance to environmental stresses, including different accessions of Col-0. A previous study reported that *hda19* plants (*hda19-3* and *hda19-5*) in the Col-0 background exhibit tolerance to salinity, drought, and heat stress (Ueda et al., 2017; Ueda et al., 2018a). In *hda19-3* plants, positive regulators, such as *NAC019* and *ABI5* in the ABA signaling pathway, that contribute to enhanced environmental stress tolerance are strongly upregulated (Ueda et al., 2017). HDA6 and HDA9 negatively regulate the expression of drought stress tolerance-related genes (class I RPD3-like HDAC) (Zheng et al., 2016; Kim et al., 2017). HDA9 is also involved in improving salt stress tolerance (Zheng et al., 2016). Collectively, the data indicate that inhibition of class I RPD3-like HDACs plays an essential role in improving environmental stress tolerance and that they function as a negative regulator of increased tolerance to environmental stresses, with the exception of the hypersensitive phenotype of the *hda6* mutant in response to salt stress (Chen et al., 2010). In contrast to the inhibition of class I RPD3-like HDACs, the multiple inhibition of class II RPD3-like HDACs (*quadruple hda5/14/15/18* mutant: *quad*) decreases salinity stress tolerance (Ueda et al., 2017). Thus, these previous studies imply that HDACs participate in an antagonistic manner in response to salinity stress.

Although the different classes of HDACs regulate salinity stress responses in a contrasting manner, the epistasis of *hda19-6* in *quad* has been observed. When a mutation of HDA19 is introduced in *quad*, the mutation converts the saltsensitive phenotype to a salt-tolerant phenotype (Ueda et al., 2017). The details of which quintuple *hda5/14/15/18/19* mutant (*quint*) acquired the tolerance to salinity stress, however, remain unknown. In the present study, a genomewide transcriptional profiling of *quint* was conducted to reveal the molecular mechanism underlying the epistasis of HDA19 against HDA5/14/15/18. The transcriptional profiling identified 4,823 salt-responsive genes in samples of whole seedlings that exhibited significant changes in expression in Col-0, *hda19-3*, *quad*, and *quint* plants. Further analysis indicated that 56.7% of the salt-responsive genes exhibited a similar expression pattern in both *hda19-3* and *quint* plants. Notably, genes for abscisic acid (ABA) biosynthesis were significantly upregulated in *quint* and *hda19-3* plants growing under non-stress and/or salinity stress conditions. In accordance with the transcriptional profiling, their ABA levels increased under salinity stress conditions. This was in addition to the upregulation of genes for transcriptional regulators, such as ABI5, reported in a previous study, suggesting that earlier expression of ABA signaling pathway genes contributes to increasing salinity tolerance in *hda19* and *quint*  plants. The results of the transcriptional profiling also indicate that an alteration of phytohormone homeostasis, in particular ABA, brassinosteroid (BR), and gibberellin (GA), occurs in *hda19* and *quint* plants.

The effects of salinity stress on plant growth could be divided into two types. Initially, osmotic stress has an impact on plant water uptake after onset of the salt stress, and, subsequently, specific ion toxicities cause to reduce cellular metabolism (Munns et al., 2006). Here, we mainly evaluate the effect of osmotic stress on plant growth because 5-day-old plants are subjected to higher salinity stress conditions (125 mM NaCl). The results indicate that 7.8% of the salt-responsive genes exhibit a similar expression pattern in *quint* and *quad* plants. Among them, the expression of *NAC016*, which is involved in leaf senescence (Kim et al., 2013), is strongly induced in *quad* plants and is considered to play a major role in producing the salt stress-sensitive phenotype. The upregulation of *NAC016*, however, was detected in *quint* mutant plants that exhibit similar tolerance to salinity stress to *hda19* plants, suggesting the presence of an HDA5/14/15/18 regulatory network independent from HDA19. The transcriptional profile conducted in the present study identified a set of genes that were regulated counteractively in *quad* and *quint* plants, which resulted in the enrichment of genes encoding enzymes involved in stress sensitivity (e.g., *IPT7* and *ISOPENTENYL TRANSFERASE 7*). The mode of action of HDA19 and HDA5/14/15/18 in salinity stress response and where crosstalk occurs between them during salinity stress response are discussed.

#### MATERIALS AND METHODS

#### Plant Material and Growth Conditions

*Arabidopsis thaliana* (L.) Heynh. (Columbia: Col-0), *hda19-3*, *hda5*/*14*/*15*/*18* (*quad*) mutant, and *hda5*/*14*/*15*/*18*/*19* (*quint*) mutant plants were used in this study. *quint* has the *hda19-6* allele. Details pertaining to the wild-type and mutant plants have been previously reported (Ueda et al., 2017). The T-DNA insertional mutants for the generation of the listed mutants were obtained from ABRC (Samson et al., 2002; Alonso et al., 2003). After surface sterilization with sodium hypochlorite, followed by two rinses with distilled water, seeds were floated on 1 ml liquid media (half-strength MS medium with 0.5% MES and 0.1% agar, pH 5.7) for 48 h at 4°C in 24-well tissue culture plates (TPP; Trasadingen, Switzerland). After germination, the seedlings in the 24-well tissue culture plates were placed on a cultivation rack at 22°C under a long-day photoperiod (16h:8h light/dark cycle) at 70–90 μE m− 2 s− 1. When plants were 5 days old (counted after seeds had germinated), the treated plants in each well were administered 25 μl 5 M NaCl (125 mM NaCl per well), while the controls received no salt treatment. Extraction of ABA, proline, and RNA from salt-treated and non-treated (control) plants was conducted for the following analysis at 2 h after the salt treatment was administered. In the evaluation of salinity tolerance in the *hda19-3*, *quad*, and *quint* mutants, 100 mM NaCl treatments were applied to 5-day-old plants (counted after seeds had germinated) in liquid culture and the percentage survival determined 4 days later (three biological replicates consisted of 10 plants; means ± SD). Plants with green true leaves were counted as a survival. Significant differences between the survival values of the experimental plants, relative to untreated or NaCl-treated wild-type plants, were determined using Student's *t* test (\*P < 0.05 and \*\*P < 0.01).

#### Microarray Analysis

Total RNAs were extracted from 5-day-old plants using an RNeasy Plant Mini Kit (Qiagen, Valencia, CA). All RNAs were further purified by incubation with RNase-free DNase I (Qiagen), according to the manufacturer's instructions. RNAs were reverse transcribed into cDNAs using 400 ng of total RNA. cDNA was labeled with a single color (Cy3) using a Quick Amp labeling kit (Agilent Technologies, Palo Alto, CA, USA) and hybridized to an *Arabidopsis* custom microarray (GEO array platform: GPL19830, Agilent Technologies) (Nguyen et al., 2015). Arrays were scanned with a microarray scanner (G2505B, Agilent Technologies). The resulting microarray data were deposited in and are available on the GEO website (GEO ID: GSE121225). The R program ver. 3.2.3 was used for the analysis of the microarray data. The fluorescence intensities of the microarray probes were normalized by quantile normalization using the limma package (Smyth, 2004). Genes with a significant change in expression were selected using the following criteria: an expression log\_2 ratio >0.5 and a controlled *p* value (false discovery rate, FDR; Benjamini and Hochberg, 1995) from a *t* test analysis <0.05. For construction of the heat map, a *Z* score was computed for each of the selected genes using gplots. For constructing the hierarchical cluster, pairwise distances between the expression data of significant expression changed genes were calculated using the "Euclidean" method and hierarchical clustering on this distance matrix was constructed using the "Ward" method. The resultant was divided into 10 groups by using the "cutree" function in R.

#### RT-qPCR Analysis

First-strand cDNA was synthesized from 250 ng total RNA with random primers. ReverTra Ace (TOYOBO) was used for the reverse transcription reaction according to the manufacturer's instructions. Transcript levels were assayed using THUNDERBIRD SYBR qPCR Mix (TOYOBO) and a StepOnePlus Real-Time PCR System (Applied Biosystems) according to the manufacturer's protocols. Gene-specific primers were designed using the PrimerQuest tool (http://sg.idtdna.com/primerquest/ Home/Index). Melting curve analyses were conducted to validate the specificity of the PCR amplification. At least three biological replicates were used in each reverse transcription–quantitative PCR (RT-qPCR) assay. *MON1* was used as a reference gene to normalize data because this gene has been demonstrated to be one of the best reference genes for *Arabidopsis* (Remans et al., 2008). The microarray data also support the use of this reference gene under the experimental conditions used in this study. The RT-qPCR scores and relevant primers are listed in **Supplementary Tables S1** and **S2**, respectively.

#### ABA Measurement

Endogenous ABA was extracted with 80% (*v*/*v*) acetonitrile containing 1% (*v*/*v*) acetic acid from whole wild-type and mutant seedlings after freeze drying. Hormone contents were determined using a UPLC-MS/MS system consisting of a quadrupole/timeof-flight tandem mass spectrometer (Triple TOF 5600, SCIEX, Concord, Canada) and a Nexera UPLC system (Shimadzu Corp., Kyoto, Japan) as described previously (Kanno et al., 2016).

#### Proline Measurement

Quantification of proline content was conducted as described in Kusano et al. (2007). Briefly, each frozen sample was extracted with a concentration of 25 mg fresh weight (FW) of whole seedling per milliliter of extraction medium (methanol/

chloroform/water, 3:1:1 *v*/*v*/*v*) containing 10 stable isotope reference compounds. An equivalent of the 27.8 µg FW of the derivatized sample was injected to Leco Pegasus HT gas chromatography–time-of-flight mass spectrometry (LECO, St. Joseph, MI, USA). Data processing was performed using ChromaTOF version <IP_ADDRESS> (LECO, St. Joseph, MI, USA). A calibration curve was generated by analyzing the derivatized proline at concentrations of 0.1, 0.5, 1.0, 3.0, 5.0, and 10 ng/µl of injection buffer, respectively. We applied specific ions at *m/z* 142 and 216 to quantify proline content. For normalization, 15 ng/µl of [13C5]-proline (Cambridge Isotope Laboratories, Andover, MA, USA) was used. We chose *m/z* 146 and 220 for it.

#### Statistical and GO Enrichment Analysis

Changes in gene expression and ABA and proline accumulations were statistically analyzed with one-way ANOVA from data obtained from three or four biological replicates. P < 0.05 was considered as significant. Gene ontology (GO) analyses of upregulated genes in *hda19* (Group 1 in **Figure 1C**) and *quad* (Group 2 in **Figure 1C**) based on results obtained from the above microarray analysis were carried out using PANTHER (Mi et al., 2013).

### Accession Numbers

Sequence data from this article can be found in the *Arabidopsis* Genome Initiative or GenBank/EMBL databases under the following accession numbers: *ABA2*; AT1G52340, *ABI5*; AT2G36270, *GA 2ox7*; AT1G50960, *GA 20ox1*; AT4G25420, *IPT7*; AT3G23630, *MON1*; AT2G28390, *NAC016*; AT1G34180, *NCED4*; AT4G19170, *UGT73C5*; AT2G36800.

# RESULTS

#### Global Transcriptome Analysis in Response to Salinity Stress

Under salinity stress conditions, the *hda19-3* and *quad* mutant plants showed a higher and a lower survival rate than wildtype (Col-0) plants, respectively. Similar to the *hda19-3* mutant plants, the *quint* plants exhibited a comparable survival rate (**Figure 1A**). To understand responses to salt stress and elucidate the details of the mechanism underlying increased tolerance to salinity stress in *quint*, transcriptional profiling using a genome-wide microarray was conducted for *quint* (*hda19-6*, an allele of *hda19*), *hda19-3*, and *quad* mutant plants, as well as wild-type plants, growing under normal and salinity stress conditions. The analysis identified 4,823 saltresponsive genes in whole seedling samples that exhibited significant expression changes in at least one of the genotypes (one-way ANOVA with Benjamini–Hochberg correction [FDR] < 0.05; **Figures 1B**, **C**). A principal component analysis (PCA) of gene expression was conducted to understand the global expression patterns of the 4,823 salt-responsive genes in the mutant and wild-type plants under salt stress and non-stressed conditions. The results indicated that the *quad*  mutant does not cluster with the wild type anymore after salt stress (open circles in **Figure 1B**). Furthermore, the clustering between *hda19-3* and *quint* is very similar between the control (filled circles in **Figure 1B**) and salt (open circles in **Figure 1B**) conditions. These results indicate that the global expression pattern in the *quint* mutant in response to salinity stress was similar to the response of the *hda19-3* mutant.

The heat map illustrates that there are significant differences in the transcriptional accumulation of the 4,823 salt-responsive genes under non-stress and salinity stress conditions in wildtype plants and each mutant genotype. A hierarchical cluster analysis revealed that differentially expressed genes in each condition can be divided into ten groups (Group 1: 729, Group 2: 281, Group 3: 903, Group 4: 183, Group 5: 607, Group 6: 498, Group 7: 304, Group 8: 239, Group 9: 193, and Group 10: 881 genes; **Figure 1C** and **Supplementary Tables S3–S12**). The lack of HDAC (*HDA19* or *HDA5/14/15/18*) expression had no significant effect on the expression of the collective 1,401 genes in Groups 3 and 6. The remainder of the 4,823 salt-responsive genes was affected under non-stress and/or salt stress conditions by the genetic defects. As indicated by the PCA, *hda19* is clearly pleiotropic over the *quad*. A total of 56.7% of the salt-responsive genes (Groups 1, 2, 5, 8, and 10) exhibited a similar expression pattern in both *hda19-3* and *quint* plants. In contrast, only 7.8% of the salt-responsive genes (Groups 4 and 9) exhibited a similar expression pattern in *quad* and *quint* plants. Notably, Group 7 enriched genes were expressed additively in *quint* plants relative to *hda19-3* and *quad* plants, suggesting that *HDA19* and *HDA5/14/15/18* redundantly act on mRNA expression. Importantly, however, HDA19 and HDA5/14/15/18 mainly regulate a set of genes in response to salinity stress in an independent manner.

#### Upregulation of ABA Signaling Pathway-Related and ABA Metabolism-Related Genes That Are Independently Regulated by HDA19 in *hda19-3* and *quint* Mutants

In general, histone deacetylases catalyze the deacetylation of histone tails, leading to DNA compaction tightening, inducing a transcriptional inactivation; hence, one would expect that mutation of HDACs has the opposite effect: transcriptional activation. GO enrichment analysis revealed that the majority of the genes in Group 1 (**Figure 1C**) containing upregulated genes under salinity stress conditions in *hda19* and *quint* compared with wild-type plants belongs to GO:0009611 (response to wounding), GO:0009414 (response to water deprivation), GO:0009651 (response to salt stress), and GO:0009737 (response to abscisic acid; **Supplementary Table S13**). Consistent with our previous study, transcription factors such as ABI5, which are considered to be a major positive regulator in the ABA signaling pathway, were strongly upregulated in both *hda19-3* and *quint* plants subjected to salinity stress conditions (Group 1 in **Figure 1C** and **Figure 2A**; **Supplementary Table S1**). In addition to the transcriptional factors involved in ABA signaling, the transcriptional profiling analysis revealed that the ABA biosynthesis genes *ABA2* and

conditions. (A) The survival rate (in percent) of each plant was evaluated 4 days after treatment with 100 mM NaCl or without NaCl (means ± SD; *n* = 3, where each biological replicate was a collection of 10 plants). *P* values were calculated using Student's *t* test (\**P* < 0.05 and \*\**P* < 0.01). (B) Principal component analysis (PCA) based on whole-genome transcriptome analysis. *Filled* and *open circles* reflect the expression profiles under normal and salinity growth conditions, respectively. (C) Hierarchical cluster analysis of salt-responsive genes. The genome-wide mRNA profiles determined by microarray analysis were obtained from 5-day-old plants treated with or without 125 mM NaCl for 2 h. Transcript data were generated from three biological replicates. The heat map represents the *Z* score, with *bars* showing values from −2 to 2. *Red* represents upregulated genes, while *blue* represents downregulated genes. Genes with a significant change in expression were selected using the following criteria: an expression log\_2 ratio >0.5 and a controlled *p* value (FDR; Benjamini and Hochberg, 1995) from a *t* test analysis <0.05.

*NCED4* were also upregulated in *quint* and/or *hda19-3* plants under non-stress and/or stress conditions (**Figures 2A**, **B** and **Supplementary Table S14**). Their upregulation was also confirmed by RT-qPCR (**Figure 2A** and **Supplementary Table S1**). In contrast, the upregulation of *ABA2* and *NCED4* was not observed in the transcriptional profiling analysis of *quad* plants (**Figure 2A** and **Supplementary Table S1**). Consistent with the transcriptome analysis, relative to wildtype plants, the ABA levels increased in *hda19-3* and *quint*

plants. Decreased accumulation of ABA in *quad* under salinity stress conditions was detected, although our transcriptome analysis could not find any signs of the decline to the best of our knowledge (**Figure 2C** and **Supplementary Table S1**). Thus, it appears that HDA19 exerts its anti-stress effect through its regulation of ABA signaling, which operates in a class II HDAC-independent manner.

ABI5 functions at the core of ABA signaling, and the induction of ABI5 expression causes growth arrest and the induction of

some of the LATE EMBRYOGENESIS ABUNDANT (LEA) proteins, such as RAB18, COR6.6, COR15A, EM1, and EM6, that prevent protein aggregation resulting from water loss at all stages of plant growth (Skubacz et al., 2016). Only the induction of *COR6.6* was almost significant (*p* < 0.058) under normal growth conditions and was clearly significant under salinity stress conditions (*p* < 0.016). Nonetheless, *EM6* was one of the direct targets (Finkelstein and Lynch, 2000), while the expression of the remaining four *LEA* genes did not increase significantly (**Supplementary Table S15**). These data suggest that an atypical activation of ABA signaling, coordinated by an absence of *HDA19* expression, may occur in *hda19* and *quint* plants.

To characterize the phenotype of mutants, we measured the accumulation of proline, known to act as a compatible osmolyte to counteract salinity stress (Szabados and Savoure, 2010). Salinity stress induced proline accumulation in all lines, and, under salinity stress conditions, there were significant declines in proline accumulation in *hda19-3* and *quint* plants compared with wild-type and *quad* plants (**Figure 2D**).

#### Altered Expression of Genes Encoding Enzymes Related to Phytohormone Biosynthesis or Catabolism

Phytohormones play an important role in fine-tuning abiotic stress response in plants (Qin et al., 2011). GA is believed to counteract the effect of ABA. Thus, fine-tuning the ABA/ GA ratio is essential to regulate both plant development and plant response to environmental stress (Vanstraelen and Benkova, 2012). Microarray data indicated that the expression of genes involved in GA biosynthesis was altered in *hda19-3* and *quint* plants. *GA2ox7*, a GA catabolic gene that degrades bioactive GA under salinity stress conditions (Magome et al., 2008; Colebrook et al., 2014), was significantly suppressed in *hda19-3* and *quint* plants under salinity stress conditions (**Figures 3A**, **C** and **Supplementary Table S1**). The expression of *GA20ox1*, which encodes a key oxidase enzyme in the biosynthesis of GA (De Vleesschauwer et al., 2012), was significantly upregulated in *hda19-3* and *quint* plants relative to wild-type plants under non-stress conditions (**Figures 3A**, **C** and **Supplementary Table S1**).

BR/ABA antagonism is known to play a role in balancing growth and stress response in plants (Clouse, 2016). BR homeostasis appeared to be altered in *hda19-3* and *quint* plants. UGT73C5, which catalyzes BL-23-*O*-glucosylation and a subsequent reduction in BR activity (Poppenberger et al., 2005), was significantly upregulated under salinity stress conditions (**Figures 3B**, **C** and **Supplementary Table S1**). These data indicate that the homeostasis of phytohormones, in particular ABA and ABA-antagonistic phytohormones such as GA and BR, may be altered in *hda19-3* and *quint* plants.

#### Presence of an HDA19-Independent Pathway and a Pathway Counteractive to HDA19 Coordinated by Class II HDACs

We have not detected any explainable genes from Group 2 (**Figure 1C**) containing upregulated genes under salinity stress conditions in *quad* compared with wild-type plants by the GO enrichment test (**Supplementary Table S16**). However, *NAC016*, whose expression is involved in senescence (Kim et al., 2013), was strongly induced in *quad* plants (Ueda et al., 2017). *Arabidopsis* plants overexpressing *NAC016* (*ANAC016-OX*) rapidly turn white when subjected to either salt or oxidative stress (Kim et al., 2013). Consistent with a previous study (Ueda et al., 2017), the *quad* plants in the current study were sensitive to salinity stress (**Figure 1A**). Unexpectedly, the significant induction of *NAC016* was retained in *quint* plants, although *quint* plants exhibited an elevated level of salinity tolerance that is similar to *hda19-3* plants (**Figure 4A**, Group 9 in **Figure 1C**, **Supplementary Table S1**). This suggests that there is an independent mechanism regulating *NAC016* expression in HDA5/14/15/18 (*quad*) plants that does not involve HDA19.

Transcriptional analysis detected a set of genes that were regulated counteractively in *hda19* alleles (*hda19-3* and *quint*) and *quad* (Groups 2 and 8 in **Figure 1C**; **Supplementary Tables S4** and **S10**) plants, which could account for the discrepancies in the upregulation of *NAC016* and the salttolerant phenotype in *quint*. Deletions of *HDA5/14/15/18* activate the expression of the cytokinin (CK) biosynthetic gene, *ISOPENTENYL TRANSFERASE7* (*IPT7*), resulting in increased levels of CK (Yanai et al., 2005) (**Supplementary Table S4**). Overproduction of endogenous cytokinin and decreased levels of CK have been reported to decrease and increase salt tolerance, respectively (Nishiyama et al., 2012; Wang et al., 2015). Downregulation of *IPT7* was observed in the current study by the introduction of the *HDA19* mutation into *quad* plants under salinity stress conditions (**Figure 4A** and **Supplementary Table S1**). This counteractive expression of *IPT7* may reflect the opposite salt sensitivity and salttolerant phenotypes observed in *quad* and *quint* plants, respectively (**Figure 4B**).

# DISCUSSION

ABA signaling is one of the major phytohormonal pathways responsible for increasing tolerance to a variety of environmental stresses, including salinity and drought (Larosa et al., 1987; Mantyla et al., 1995; Bari and Jones, 2009; Vishwakarma et al., 2017). The current study indicated that ABA signaling is strongly induced in *hda19* mutant plants. There is no clear morphological difference between *quad* and wild-type plants. On the other hand, *hda19-3* and *quint* showed lower germination rates under non-stress conditions, although statistical significance was not detected (**Figure 1A**). A wide range of developmental abnormalities such as growth inhibition was observed in *hda19-3* (Hollender and Liu, 2008). Some of them might be explainable by the enhanced acclimation of ABA. The transcriptome analyses revealed that an HDA19 defect primarily activates the ABA signaling pathway, which is independently regulated from class II HDACs. ABA signaling plays a pivotal role in the growth reduction of plants during the early phase of plant response to salinity stress (Ismail et al., 2014). In contrast, deficiencies in class II HDAC genes appear to result in a malfunction in the growth reduction because expansin 12 (a cell wall-loosening protein; **Supplementary Table S4**) (Lee et al., 2001) and TRH1 (K+ transporter for root hair elongation; **Supplementary Table S4**) (Desbrosses et al., 2003) are activated in response to salinity stress conditions in *quad* plants as in the case of *IPT7* (**Figure 4**). In *hda19-3* and *quint* plants, the upregulation of these genes could not be detected (**Supplementary Tables S4** and **S10**). The HDA19 defect, in which HDA19 is not produced,

same as shown in Figure 2A. (B) A simplified model for the counteractive and hierarchical control of the expressions of *NAC016* and *IPT7* genes mediated through HDA19 and HDA5/14/15/18 deficiencies.

may induce an earlier expression of genes involved in the ABA signaling pathway, including ABA metabolism and other transcriptional factors.

BR homeostasis appeared to be altered in *hda19* and *quint* plants. UGT73C5, which catalyzes BL-23-*O*-glucosylation leading to reduced BR activity (Poppenberger et al., 2005), was significantly upregulated under salinity stress conditions (**Figure 3C**). Under the experimental conditions used in the current study, ABA signaling may be activated due to the reduced accumulation of active BR, which results in increased salinity stress tolerance. In support of this premise, *ABI5* expression has been reported to directly suppress BR signaling through the BRASSINAZOLE RESISTANT 1 (BZR1) transcription factor (Yang et al., 2016). The upregulation of *ABI5* expression under salinity stress conditions may indicate the reduced accumulation of active BRs in *haa19* and *quint* plants. A previous study also indicated that BR signaling is involved in the suppression of GA biosynthetic genes, such as *GA20ox* and *GA3ox3*, and the induction of a GA inactivation gene, such as *GA2ox*, which leads to the suppression of GA signaling (De Vleesschauwer et al., 2012). In the current study, the transcriptome analysis indicates that GA signaling may be activated in *hda19* and *quint* plants because suppression of *GA2ox7* and the induction of *GA20ox1* were detected (**Figure 3C**). Considering the above data, a reduced accumulation of active BRs may trigger a change in the expression pattern of genes encoding phytohormonecatabolizing enzymes.

In contrast to ABA metabolism, a defect in GA metabolism generally enhances tolerance to abiotic stresses, such as salinity and drought stress (Achard et al., 2006; Colebrook et al., 2014). The results of the present study indicate that GA metabolism is activated in *hda19* and *quint* plants based on the transcriptome analysis (**Figure 3C**). These plants, however, exhibit tolerance to salinity stress. Furthermore, *hda19* plants also exhibit tolerance to drought stress (Ueda et al., 2018a). If the hormone network is altered in *hda19-3* or *quint* plants, the altered hormone network (coincidental activation of antagonistic hormone signaling of ABA and GA) may allow plants to alleviate growth arrest induced by ABA signaling. Notably, *ga2ox7* mutants that accumulate active GA exhibit less growth retardation than Col wild-type plants of *Arabidopsis* subjected to high-salinity stress (Magome et al., 2008). ABA stimulates generally an increase in proline accumulation (Szabados and Savoure, 2010). ABA increased in *hda19* and *quint*; nonetheless, a lower proline accumulation was found in them than in wild-type plants and *quad* (**Figure 2D**). The unusual response in proline accumulation might be due to the alteration in phytohormone homeostasis under salinity stress conditions. Further studies are needed to reveal how epigenetic regulation mediated by HDA19 controls stress responses and crosstalk between phytohormones.

We confirmed the altered mRNA expressions of *ABI5*, *ABA2*, *NCED4*, *GA 2ox7*, *GA 20ox1*, *UGT73C5*, *NAC016*, and *IPT7* genes under normal and/or salinity growth conditions by microarray and RT-qPCR analysis. Among them, we found tissue- or organspecific expressions, except for *ABA2* and *GA 2ox7* genes, by *Arabidopsis* eFP Browser (Winter et al., 2007) (**Supplementary Figure S1**). In the study, we have not revealed whether their alterations are due to the derepression of their tissue-specific expression or activation through the inhibition of HDAC activity. A detailed analysis of their expression is needed to understand how HDAC inhibition alters salinity stress response.

HDA19 mainly appears to mediate histone H3 acetylation (Zhou et al., 2005; Krogan et al., 2012). In class II HDACs, HDA5 participates in histone H3 acetylation (Luo et al., 2015). HDA15 participates in histone H3 and H4 acetylation (Liu et al., 2013; Gu et al., 2017). HDA14 is located in chloroplasts, and its localization to chromatin has not been documented (Alinsug et al., 2012). Synergistic upregulation of HDA19 and HDA5/14/15/18 deficiencies observed in *quint* plants may be due to the redundancy in histone or non-histone acetylation coordinated by both HDA19 and class II HDACs. As shown in **Figure 4A**, the expression of *NAC016* is perhaps independent of HDA19, as the majority of the responses described in the manuscript seem to be otherwise. In the case of HDA19 and HDA5/14/15/18-dependent regulation (HDA19- and HDA5/14/15/18-dependent expressions, categories 1, 2, 5, 8, and 10 and categories 4 and 9, respectively, in **Figure 1C**), levels in the acetylation of histones around target sites or proteins, including non-histones, may be independently coordinated by these HDACs. In actuality, multiple acetylations of non-histone proteins have been documented using HDAC inhibitors or mutants (Hao et al., 2016; Hartl et al., 2017). For example, it remains unclear whether the altered expression of genes for BR, GA, and other hormone pathways is due to the HDA19 action or a reflection of the effect on ABA signaling. Further analysis to determine whether HDACs directly regulate the gene expression of candidate signaling pathways through histone acetylation and/or are involved in the acetylation of key enzymes involved in hormonal crosstalk, as highlighted in the current study, will elucidate the detailed processes in stress responses in plants that are coordinated by HDA19 and HDA5/14/15/18.

#### AUTHOR CONTRIBUTIONS

MU and MSek designed the experiments. MU, AM, SW, MKo, MT, and JI conducted the experiments. MU, AM, SW, MKo, MKu, KS, and MSeo analyzed the data. MU, AM, SW, MKo, MKu, MSeo and MSek wrote the manuscript.

# FUNDING

This work was supported by grants to MU and MSek from the RIKEN and grants to MSek from Japan Science and Technology Agency (JST) [Core Research for Evolutionary Science and Technology (CREST, grant number JPMJCR13B4)], and KAKENHI on Innovative Areas (grant no. 16H01476, 18H04791, and 18H04705) of the Ministry of Education Culture, Sports and Technology of Japan.

#### ACKNOWLEDGMENTS

The authors would like to show their appreciation to Ms. C. Torii, Ms. K. Mizunashi, and Ms. A. Sato for their technical support. They are also grateful to Ms. S. Nakae, Ms. K. Kaneko, and Ms. F.

#### REFERENCES

Sakai in the support unit for Bio-Material Analysis and RIKEN CBS Research Resources Division for technical help with the nucleotide sequencing analyses.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.01323/ full#supplementary-material

FIGURE S1 | eFP diagrams of *ABI5*, *ABA2*, *NCED4*, *GA 2ox7*, *GA 20ox1*, *UGT73C5*, *NAC016*, and *IPT7* genes

and RAB18 proteins in *Arabidopsis thaliana. Plant Physiol.* 107, 141–148. doi: 10.1104/pp.107.1.141

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2019 Ueda, Matsui, Watanabe, Kobayashi, Saito, Tanaka, Ishida, Kusano, Seo and Seki. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*


• First chromosome section

A precedence-matrix (P) is used to derive feasible disassembly sequences. It describes the order of precedence of the disassembly steps and can be obtained from carrying out manual disassembly experiments or by using the computer-aided design (CAD) models of the product to be disassembled [28]. Another approach to determining possible precedence relationships in an automatic way based on computer vision was presented in [13]. This can avoid the disadvantages of the two mentioned methods: while manual pre-processing is error-prone and time-consuming, pre-processing using CAD data is often inaccurate. This is because CAD models can rarely still reliably describe the product at their EoL

phase, for example, due to corrosion or changes that have been applied to the product during the usage phase. We used the disassembly precedence graph (DPG) to generate the P matrix for the use case in this paper. It allows a simple one-to-one comparison of all components. If a component i is the predecessor of the component j, the value Pij gets the value 1, otherwise 0. The numerical method to derive feasible disassembly sequences from the P matrix is illustrated in Figure 5.

**Figure 4.** Disassembly strategy optimizer using a genetic algorithm.

```
*$
$
*!
#%

,=
A%	$

;
A	6
M-
2
/O%1

>$
=

7-
A
,

M-
2
/
O1
&77

,=,

E
7-
+

=

7-

M
N-
"
MN-
"
```
**Figure 5.** Method to generate feasible disassembly sequences from a given precedence matrix.

• Second chromosome section

The transfer of battery components into a circular economy takes place via the implementation of circular economy strategies. Here, there are different strategies, which mainly differ in preparing the components and the application field. The waste hierarchy of the European Commission contains five priorities [29]. In this work, we focus on the priorities on the top, as they should be preferred in a circular economy. The highest one is prevention, for example, by extending the life of products through predictive maintenance or simple repair operations. The second priority is repair for reuse, followed by recycling. These priorities can be achieved by applying diverse circularity strategies. Potting et al. [30] identified ten strategies and divided them into three categories: (1) smarter product use and manufacture, (2) extend lifespan of product and its parts, and (3) useful application of materials. Here, we consider only the strategies, which can be applied at the EoL phase. Circularity solutions during the stages of design, manufacturing, and use are ignored. In this scope, the EoL strategies can be reduced to four strategies: reuse, remanufacturing, repurposing, and recycling. The integration of these strategies in the waste hierarchy and the material flows after their application are illustrated in Figure 6. In the literature, different definitions can be found for these strategies [24]. In this paper, we adopt the definitions of Potting et al. [30]. Reuse means the utilization of a discarded battery or a set of its components by another user in the automotive field. Remanufacturing is the treatment of battery parts so that they can at least meet the requirements of newly manufactured products and their utilization for the manufacturing of EVBs. If the battery or its parts are reconditioned and used in another application field, such as stationary energy storage, this is called repurposing. Recycling is the recovery of materials. In this case, pure and high-quality materials for use in the automotive sector should be targeted.

Figure 6 shows that disassembly is of central importance for the implementation of all these circularity strategies. The components of the EVBs can be allocated to different strategies depending on their condition. However, the selection of the optimal circular economy strategy does not depend exclusively on this. Other factors, such as the disassembly costs and the market constraints, such as the potential revenues, also play an important role. Therefore, we consider the selection of the EoL strategy as a part of the disassembly planning. This paper assumes that the feasibility of a circular economy strategy for a given component depends on its condition. That is why the second section of the chromosome must fulfill the condition constraints defined by a condition vector S containing the feasible circular economy strategies CESi for every component i—see Equation (1). Here, the priority of the circular economy strategies is taken into account. If a component i is in excellent condition, it can be allocated to all strategies (reuse, remanufacturing, repurposing, and recycling). In this case, CESi is assigned the value 1. If CESi equals 2, the reuse option will be excluded. Part i can be neither reused nor remanufactured if CESi equals 3. CESi is assigned the value 4 if recycling is the only possible recovery option.

$$\mathbf{S} = \begin{pmatrix} \text{CES}\_1 \\ \dots \\ \text{CES}\_n \end{pmatrix}; \text{CES}\_i \in \{1; 2; 3; 4\} \tag{1}$$

**Figure 6.** Overview about the considered end-of-life strategies.

The S vector cannot just be seen as a collection of testing results. Other factors can play a role in determining the potentially possible circular economy strategies for the different battery parts, such as the employees' experience in the disassembly factory. Testing results include the state of health (SoH) and state of charge (SoC) of battery cells and modules and additional parameters for the rest of the components.

• Third chromosome section

This section consists of a single gene and is used to define the disassembly depth. An EVB could be entirely disassembled by separating all its parts. However, this approach is neither economically nor environmentally practical in an industrial context [16]. Therefore, EVBs are more likely to be subject to incomplete disassembly. Here, there are two methods to perform incomplete disassembly: (1) the selective method and (2) the unrestricted method. The selective method means that specific components are selected to be disassembled. Subsequently, the disassembly planner needs to calculate a strategy for the optimal extraction of these parts. Here, the high-value strategy and the high-impact strategy can be distinguished [31]. For EVBs, the removal of the modules could present a high-value disassembly strategy. The high-impact strategy applies when, for example, a module with safety risks is identified and has to be replaced before reusing the battery. In contrast, no target components are selected in the unrestricted incomplete disassembly. The disassembly planner can freely calculate the optimal disassembly strategy based on an objective function. This method is considered in this paper. Thereby, the gene representing

the disassembly depth is randomly generated with values between 0 (no disassembly) and the maximum number of components n while generating the initial population.

In Fehler! Verweisquelle konnte nicht gefunden werden, a disassembly precedence graph of a theoretical product, the associated precedence matrix, and the entire structure of a possible chromosome depending on a given condition vector are presented. In this case, one feasible disassembly sequence is 1-2-3-5-6-4, possible circular economy strategies at the component level are 1-2-4-1-2-4, and the disassembly is complete, since the last gene matches the component number.