Investigation of a Novel Salt Stress-Responsive Pathway Mediated by Arabidopsis DEAD-Box RNA Helicase Gene AtRH17 Using RNA-Seq Analysis

Previously, we reported that overexpression of AtRH17, an Arabidopsis DEAD-box RNA helicase gene, confers salt stress-tolerance via a pathway other than the well-known salt stress-responsive pathways. To decipher the salt stress-responsive pathway in AtRH17-overexpressing transgenic plants (OXs), we performed RNA-Sequencing and identified 397 differentially expressed genes between wild type (WT) and AtRH17 OXs. Among them, 286 genes were upregulated and 111 genes were downregulated in AtRH17 OXs relative to WT. Gene ontology annotation enrichment and KEGG pathway analysis showed that the 397 upregulated and downregulated genes are involved in various biological functions including secretion, signaling, detoxification, metabolic pathways, catabolic pathways, and biosynthesis of secondary metabolites as well as in stress responses. Genevestigator analysis of the upregulated genes showed that nine genes, namely, LEA4-5, GSTF6, DIN2/BGLU30, TSPO, GSTF7, LEA18, HAI1, ABR, and LTI30, were upregulated in Arabidopsis under salt, osmotic, and drought stress conditions. In particular, the expression levels of LEA4-5, TSPO, and ABR were higher in AtRH17 OXs than in WT under salt stress condition. Taken together, our results suggest that a high AtRH17 expression confers salt stress-tolerance through a novel salt stress-responsive pathway involving nine genes, other than the well-known ABA-dependent and ABA-independent pathways.


Introduction
RNA helicases (RHs) are present in most prokaryotic and eukaryotic organisms, which catalyze the unwinding of DNA or the secondary structure of RNA, and thus, play essential roles in almost every aspect of genetic processes such as replication, transcription, translation, repair, and recombination [1,2]. RHs have been classified into six superfamilies (SF1-SF6) based on specific motif sequences and domain structures. Superfamily II (SF2), the largest helicase family, mainly consists of the DEAD-box RHs, which are named after the strictly conserved sequence "Asp-Glu-Ala-Asp" (D-E-A-D) [1,2]. Fifty-eight DEAD-box RHs have been identified in Arabidopsis [1,3].

Transcriptomic Profiling of AtRH17 OXs
In a previous study, we reported that overexpression of AtRH17 confers salt stress-tolerance at the seedling and mature plant stages. Interestingly, the well-known ABA-dependent and ABA-independent salt stress-responsive genes such as RD29A, RAB18, RD29B, RD22, COR47, DREB2A, and DREB2B showed similar or lower expression levels in AtRH17 OXs than in WT under salt stress conditions [17], implying that salt stress-tolerance of AtRH17 OXs is mediated by an uncharacterized pathway or mechanism other than the well-known stress-responsive pathways. In this study, to clarify the regulatory mechanism of salt stress-tolerance of AtRH17 OXs, RNA-Seq was performed using 10-day-old WT and AtRH17 OX whole seedlings. The mapping of RNA-Seq reads to the Col-0 genome was successful, with a mapping rate of 97.47% to 98.01% (Table 1). The number of mapped reads ranged from 26.4 to 30.3 million (Table 1). Genes having very low abundance were removed from the analysis, leaving 22,884 genes for further analysis.

Analysis of Up-and Downregulated Genes in AtRH17 OXs
A total of 397 differentially expressed genes (DEGs) between WT and AtRH17 OX plants were identified with >3-fold differences in expression and a p-value < 0.05. Among these, 286 genes were upregulated and 111 genes were downregulated in AtRH17 OXs (Figure 1a). The distribution of upand downregulated genes is shown using volcano plot ( Figure 1b). Differentially expressed genes in AtRH17-overexpressing transgenic plants (OXs) compared to wild type (WT) from the RNA-Sequencing (RNA-Seq) analysis. (a) Numbers of genes more than 3-fold up-and downregulated in AtRH17 OXs compared to WT. (b) Volcano plot of differentially expressed genes (DEGs) identified between WT and AtRH17 OXs (p-value < 0.05 and log2 ratio ≥ 1). The upregulated genes are represented by a red dot and downregulated genes by a green dot, grey dots indicate genes that are not differentially expressed.
The analysis of GO terms with three categories, namely, biological processes, molecular functions, and cellular components, using the selected genes is useful in predicting the altered biological and molecular processes. Therefore, 397 DEGs were subjected to GO enrichment analysis to determine their functional significance in AtRH17 OXs. The upregulated genes were enriched in response to toxic substance, toxin catabolic process, glutathione metabolic process, and oxidation-reduction process in the biological process categories of GO annotation; glutathione transferase activity, lipid binding, heme binding, flavin adenine dinucleotide binding, and electron carrier activity in the molecular function categories; and extracellular region, apoplast, extracellular space, RNA polymerase II transcription factor complex, and cell wall in the cellular component categories (Supplementary Figure S1). The downregulated genes were enriched for plant-type cell wall organization, hydrogen peroxide catabolic process, and response to brassinosteroid in the biological process categories of GO annotation; structural constituent of cell wall, transcription factor activity, peroxidase activity in the molecular function categories; and cell wall, extracellular region, and plant-type cell wall in the cellular component categories (Supplementary Figure S2). GO enrichment analysis showed that many genes, which are involved in various biological and molecular processes such as signaling, secretion, detoxification, transcription, and stress responses, were up-and/or downregulated in AtRH17 OXs.
We further analyzed the biological process categories of GO annotation to decipher the biological functions of the DEGs in AtRH17 OXs. The upregulated genes were highly enriched in the biological process of oxidation-reduction (34 genes), defense response to fungus (14 genes), defense response (14 genes), response to oxidative stress (12 genes), response to salt stress (11 genes), response to toxic substance (10 genes), and response to water deprivation (10 genes) ( Table 2), demonstrating that AtRH17 might be involved in abiotic and biotic stress responses. In addition, AtRH17 might be involved in the metabolic processes of sugar compounds and homeostasis ( Table 2). The GO analysis of the downregulated genes revealed that the biological processes in which the downregulated genes are involved were mostly associated with oxidation-reduction (11 genes), defense response (seven genes), plant-type cell wall organization (six genes), response to oxidative stress (five genes), and hydrogen peroxide catabolic process (four genes) ( Table 3), indicating that AtRH17 might negatively regulate oxidative stress-related genes. Moreover, AtRH17 can have a minor role in negative regulation of plant growth-related genes based on the enriched GO terms such as response to brassinosteroid, response to gibberellin, unidimensional cell growth, and response to ethylene (Table 3).
To identify the functional pathways in which AtRH17 is involved, KEGG pathway analysis was performed. Both the up-and downregulated genes were involved in metabolic pathways, biosynthesis of secondary metabolites, and phenylpropanoid biosynthesis (Figures 2 and 3). In contrast, only the upregulated genes were involved in glutathione metabolism and flavonoid biosynthesis, whereas the downregulated genes were involved in plant hormone signal transduction and MAPK signaling pathway (Figures 2 and 3), indicating that AtRH17 might be involved in several functional pathways through the up-and/or downregulation of various related genes. We also performed hierarchical clustering to classify and identify the relationship among DEGs. The up-and downregulated genes were classified into individual hierarchies. The upregulated genes were classified into six groups, whereas the downregulated genes were classified into four groups ( Figure 4).
To validate the RNA-Seq results by quantitative RT-PCR, we selected five upregulated genes, namely, TSPO, LEA4-5, ABR, LEA18, and DIN2/BGLU30 and four downregulated genes, namely, PIL1, FRF4, MYB108, and NAC019. Thereafter, we analyzed the expression of the selected genes in 10-day-old WT and AtRH17 OX seedlings. The expression of all five upregulated genes was significantly higher in AtRH17 OXs than in WT (Figure 5a-e). Especially, TSPO expression was more than 5-fold higher in AtRH17 OXs than in WT, whereas LEA4-5 expression was 2.3-fold higher in AtRH17 OXs than in WT (Figure 5a,b). The expression levels of ABR, LEA18, and DIN2/BGLU30 were more than 3-fold higher in AtRH17 OXs than in WT (Figure 5c-e). In contrast, the four downregulated genes showed lower expression in AtRH17 OXs than in WT (Figure 5f-i). PIL1 expression in AtRH17 OXs was only almost half of that in WT (Figure 5f). FRF4 and MYB108 expression in AtRH17 OXs were only one-fifth of that in WT (Figure 5g,h). Moreover, NAC019 expression in AtRH17 OXs was one-tenth of that in WT ( Figure 5i). These results are consistent with the RNA-Seq results.
(i), and AtRH17 (j) in 10-day-old WT and AtRH17 OX seedlings. GAPc was used as an internal control. Transcript levels in WT were set as 1. Three independent reactions were performed for each technical replicate. Two technical replicates were performed for each biological replicate. At least two biological replicates showed similar results, with one shown here. Error bars represent standard deviation (n = 6 reactions) and * indicates t-test p < 0.05.

Identification of the Salt Stress-Tolerance Pathway in AtRH17 OXs
To elucidate the salt stress-tolerance mechanism in AtRH17 OXs, we isolated 19 upregulated genes enriched in the salt stress-related GO terms such as "response to salt" (GO:000965) and "response to water deprivation" (GO:0009414) ( Table 4), indicating that these genes might function in salt stress-tolerance in AtRH17 OXs. To see the expression patterns of the 19 genes under stress conditions, we performed expression analysis using Genevestigator under drought, osmotic, and salt stress conditions. Five experiments were used for each stress condition. We divided the 19 genes into three groups depending on their expression patterns. In the first group, eight genes, namely, CHIA, CSD1, PER23, MYB12, ZFHD1, ANNAT7, GLP9, and MYB29, did not respond to any of the stresses ( Figure 6). In the second group, LOX2 and GSTU17 showed high expression levels under drought condition, whereas their levels did not increase under osmotic and salt stress conditions ( Figure 6). In the third group, nine genes, namely, DIN2/BGLU30, LEA18, GSTF6, GSTF7, LTI30, TSPO, HAI1, LEA4-5, and ABR, showed high expression levels under all three examined stress conditions (Figure 6), demonstrating that they might be involved in salt stress-tolerance in AtRH17 OXs (Table 5).  Next, we analyzed the expression of genes selected based on the Genevestigator analysis in WT and AtRH17 OXs under salt stress condition to verify whether these genes are involved in the salt stress response of AtRH17 OXs. LEA4-5 and ABR, late embryogenesis abundant (LEA) protein genes, and TSPO were selected for expression analysis among the above nine genes. The biological functions of LEA proteins in osmotic stress-tolerance are well known [18,19]. TSPOs have also been shown to be responsive to salt and osmotic stresses and ABA [20,21]. Quantitative RT-PCR analysis showed that the expression levels of LEA4-5, ABR, and TSPO were higher in AtRH17 OXs than in WT under salt stress conditions (Figure 7). Interestingly, the differences in expression levels of LEA4-5 and TSPO between WT and AtRH17 OXs were the highest under early salt stress conditions such as after 1 and 2 h of NaCl-treatment. Subsequently, the difference decreased gradually, and no difference was observed at 8 h after the NaCl-treatment (Figure 7a,b). However, the expression of ABR was almost similar in WT and AtRH17 OXs until 2 h after the NaCl-treatment. Thereafter, the expression in AtRH17 OXs significantly increased to more than that in WT until 8 h after the NaCl-treatment (Figure 7c), indicating that LEA4-5 and TSPO might function in the early salt stress response, whereas ABR might be involved in the late salt stress response. Protein network analysis of the nine salt stress-responsive genes using Cytoscape String Apps revealed that the nine genes were tightly connected and interacted with each other (Figure 8), indicating that the salt stress-tolerance in AtRH17 OXs might be the result of precise upregulation of these nine genes.

Discussion
In this study, we investigated the regulatory mechanism of the salt stress-tolerance of AtRH17 OXs using RNA-Seq analysis. AtRH17 is a member of the DEAD-box RHs, which function in RNA metabolism, ribosome biogenesis, and transcriptional and translational regulation [1,2,17].
RNA-Seq, a high-throughput and next-generation sequencing technology, has allowed for rapid analysis of large genomic datasets and quantification of transcriptomes [26]. In addition, RNA-Seq analysis can be used to identify and quantify transcripts without prior information of genes and can provide information regarding alternative splicing and sequence variations in identified genes [26,27]. Global gene expression patterns have been determined using RNA-Seq analysis of samples at different developmental stages, in response to different stimuli or in different genotypes [28,29]. RNA-Seq analysis of AtRH17 OXs revealed that 286 genes were upregulated and 111 genes were downregulated in AtRH17 OXs relative to WT (Figure 1). The upregulated and downregulated genes were classified into six and four individual groups by hierarchy analysis, respectively ( Figure 4). GO annotation enrichment analysis showed that AtRH17 might regulate genes that are involved in various functions such as signaling, secretion, detoxification, and transcription, depending on enriched GO terms such as response to toxic substances, toxin catabolic process, glutathione metabolic process, oxidation-reduction process, transcription factor activity, RNA polymerase II complex, and extracellular region ( Supplementary Figures S1 and S2). KEGG pathway analysis showed that AtRH17 might be involved in metabolic pathways, biosynthesis of secondary metabolites, plant hormone signal transduction, and MAPK signaling pathway via the regulation of downstream genes functioning in these pathways (Figures 2 and 3).
Previously, we reported that the expression levels of the well-known ABA-dependent and ABA-independent salt stress-responsive genes such as RD29A, RAB18, RD29B, RD22, COR47, DREB2A, and DREB2B were not higher in AtRH17 OXs than in WT under salt stress condition [17]. Using RNA-Seq analysis, we confirmed our previous results that the expression levels of these seven genes were not significantly increased in AtRH17 OXs (data not shown), indicating that the salt-tolerance in AtRH17 OXs is mediated through pathway(s) other than the well-known ABA-dependent and/or ABA-independent stress-responsive pathways.
Analysis of the upregulated genes in AtRH17 OXs using GO annotation and Genevestigator revealed that nine genes, namely, LEA4-5, GSTF6, DIN2/BGLU30, TSPO, GSTF7, LEA18, HAI1, ABR, and LTI30, are salt stress-responsive genes and might be involved in the salt-tolerance in AtRH17 OXs (Table 5). These nine genes are tightly inter-connected ( Figure 8). Among these, LEA4-5, TSPO, and ABR were selected for quantitative RT-PCR analysis because their expression increased under all abiotic stress conditions examined in Genevestigator analysis and they showed high expression levels in AtRH17 OXs ( Figure 6 and Table 5). Three genes showed higher expression levels in AtRH17 OXs than in WT under salt stress conditions (Figure 7). LEA4-5 is a member of LEA proteins. The C-terminal region of LEA4-5 is responsible for its antioxidant activity and in the scavenging of metal ions under stress conditions, whereas the N-terminal can function as a chaperone in the folding of enzyme proteins and in preventing their unfolding that results in protection of the enzyme activity [19,30]. TSPO, belonging to the Trp-rich sensory protein/peripheral-type benzodiazepine receptor group protein, interacts with and regulates PIP2;7, a plasma membrane aquaporin, in the endoplasmic reticulum and Golgi membrane during abiotic stress conditions [20,21]. The interaction between TSPO and PIP2;7 triggers the reduction of PIP2;7 channels present in the plasma membrane and reduces the water transport activity [20,21]. ABR is a Ser/Thr kinase and interacts with phospholipase D-derived phosphatidic acid (PA) [31,32]. ABR increases the ABR-dependent phosphorylation of PIN2, which activates auxin efflux, alters auxin accumulation, and promotes root growth under salt stress conditions [31,32]. T-DNA insertional mutants of ABR are hypersensitive to salt stress in primary root elongation and do not respond to PA [31,32]. These results suggest that the upregulated LEA4-5, TSPO, and ABR confer salt stress-tolerance in AtRH17 OXs.

Plant Materials and Growth Conditions
All plant materials used in this study were of the Arabidopsis thaliana accession Col-0 background. Sterilized seeds were incubated in the dark for 2-3 days at 4 • C and then germinated, grown on half-strength Murashige and Skoog (MS) agar plates, supplemented with 1.5% (w/v) sucrose and B5 vitamin. Seedlings were grown under short-day (SD) conditions (8 h light:16 h dark photoperiod) at 22 • C.

Plant Stress Treatment for Quantitative RT-PCR
For the analysis of salt stress-responsive gene expression under salt stress conditions, 10-day-old WT and AtRH17 OX seedlings grown under SD conditions were placed on filter papers soaked in MS solution containing 150 mM NaCl. After 0, 1, 2, 4, and 8 h, the seedlings were harvested; seedlings harvested at 0 h were used as the control.

Quantitative RT-PCR
Quantitative RT-PCR analysis was performed using 2× POWER SYBR Green PCR Master mix (Applied Biosystems, Foster, CA, USA) with diluted cDNA as the template. The C t (cycle at the threshold) value was set to be constant throughout the study and corresponded to the log-linear range of PCR amplification. The normalized amount of target reflected the relative amount of target transcripts with respect to that of the endogenous reference gene, GAPc. The amplification conditions included an initial denaturation at 95 • C for 10 min, followed by repeated cycles at 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 30 s. The primers used are listed in Supplementary Table S1.

Library Preparation and RNA-Sequencing
The RNA quality was assessed with Agilent 2100 bioanalyzer using the RNA 6000 Nano Chip (Agilent Technologies, Santa Clara, CA, USA), and RNA quantification was performed using a ND-2000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA). For the control and test RNAs, a library was constructed using QuantSeq 3 mRNA-Seq Library Prep Kit (Lexogen, Inc., Vienna, Austria), according to the manufacturer's instructions. In brief, 500 ng total RNA was hybridized to an oligo dT primer containing an Illumina-compatible sequence at its 5 -end and reverse transcription was performed. After degradation of the RNA template, second strand synthesis was initiated using a random primer containing an Illumina-compatible linker sequence at its 5 -end. The double-stranded library was purified of the reaction components by magnetic separation. The library was amplified to add the complete adapter sequences required for cluster generation, and the finished library was purified from the PCR components. High-throughput single-end 75 sequencing was performed using NextSeq 500.

RNA-Seq Data Analysis
QuantSeq 3 mRNA-Seq reads were aligned using Bowtie2 [33]. Bowtie2 indices were either generated from the genome assembly sequence or representative transcript sequences for aligning to the genome and transcriptome. The alignment file was used for assembling the transcripts, estimating their abundance, and detecting the differential expression of genes. DEGs were determined based on the counts from unique and multiple alignments using coverage in Bedtools [34]. The Read Count data were processed based on the quantile normalization method employing EdgeR within R [35] using Bioconductor [36]. GO annotation enrichment was performed using DAVID [37] with default parameters. KEGG pathway analysis was conducted using KEGG Mapper [38] Gene clustering was performed using MeV ver. 4.9.0 [39]. Protein network analysis was performed using String Apps of Cytoscape ver. 3.7.2 [40].

Statistical Analysis
Statistical analysis was performed using one-way analysis of variance and Dunnett's post-hoc test, as implemented in IBM SPSS v.23 (IBM Corp., Armonk, NY, USA).

Conflicts of Interest:
The authors declare no conflicts of interest.