Transcriptome Analysis Reveals Regulatory Framework for Salt and Drought Tolerance in Hibiscus hamabo Siebold & Zuccarini

: Hibiscus hamabo Siebold & Zuccarini ( H. hamabo ) is tolerant to salt and drought conditions, but the molecular mechanisms that underlie this stress tolerance remain unclear. In this study, the transcriptome of H. hamabo roots was investigated under NaCl or PEG treatment. A total of 20,513 and 27,516 signiﬁcantly changed known genes at 6 h and 24 h, respectively, were detected between the salt or drought treatments and the control libraries. Among these, there were 3845 and 7430 overlapping genes under the two stresses at 6 h and 24 h, respectively. Based on the analysis of enriched KEGG pathways and clustering of expression patterns, the DEGs that were continuously up- or down-regulated under both salt and drought treatments were mainly enriched in MAPK signaling pathway, transcription factors, transporters and other pathways. The transcriptome expression proﬁles of H. hamabo provide a genetic resource for identifying common regulatory factors involved in responses to different abiotic stresses. In addition, the identiﬁed factors may be useful to developing genetic breeding strategies for the Malvaceae.


Introduction
Water-deficit and salinity are abiotic stresses that severely impact various cellular and whole plant processes, influencing plant yield and quality [1]. These stress conditions can reduce the availability of water for basic cellular functions and affect the maintenance of turgor pressure [2], while salinity can cause an ionic imbalance [3]. Longer periods and a higher degree of water-deficit and salinity eventually generate oxidative stress in plants and damage cell components, ultimately affecting growth [2]. In response to drought and salt stresses, therefore, plants exhibit complex mechanisms [4]. Abscisic acid (ABA), Ca 2+ and reactive oxygen species (ROS), key signal transduction components, accumulate in plants to regulate downstream stress-related transcription factors through signaling cascades that regulate functional genes and improve stress resistance [5,6]. For example, mitogen activated protein kinase (MAPK) module, which was mainly composed with MAPK (MAPK/MPK), MAPK kinase (MAPKK/MKK), and MAPK kinase kinase (MAP-KKK/MAP3K), played diverse roles in eukaryotic cell signaling and was activated by stresses like drought, salt and then triggered stress responses [7,8]. It has been found that drought and salt could activate MKK1 which could induce H 2 O 2 via CAT pathway to enhance drought and salt tolerance in Arabidopsis (Arabidopsis thaliana L.) [7]. Under salt stress, MEKK1 may activate MKK2, which would further phosphorylate MPK6. Then MPK6 could activate MYB41 to response to salt stress in Arabidopsis [7]. The MYB transcription factor family is widely distributed in plants and has been reported playing roles in regulating plant-specific biological processes, as responsing and creasing tolerance to biotic stresses [9]. For example, AtMYB15 was found to increase salt and drought tolerance in Arabidopsis [10]. However, there may be differences in the adaptive mechanisms used by different plants under salt and drought stresses, and these require further exploration, especially in non-model plants.
RNA sequencing (RNA-Seq) technology is an efficient and economical platform that uses high-throughput sequencing techniques and is applied to analyse gene expression at the transcriptome level under various stresses [11]. For example, stress-regulated genes and their corresponding mechanisms have been investigated in Taxodium hybrid "Zhongshanshan" [12], Iris lactea var. chinensis [13] and Paspalum vaginatum [14].
Hibiscus hamabo Siebold & Zuccarini (Malvaceae) is naturally distributed in coastal sands near sea level in China, Japan and Korea and is cultivated in India and the Pacific islands [15]. As a semi-mangrove plant, H. hamabo has high tolerance to tidal saline soil, immersion and arid environments and is considered to be one of the best choices in increasing the vegetation in higher salinity landscapes [16,17]. Although there are several studies on the physiological responses of H. hamabo to salinity and drought [15,16], the molecular mechanisms underlying its stress tolerance remain unknown.
In the current study, a comprehensive transcriptome analysis was performed on H. hamabo roots under salt and drought stress to investigate its tolerance mechanisms. This work may provide useful information and potential genetic resources for the improvement of the tolerance characteristics in the Malvaceae.

Sampling
H. hamabo seeds were collected from Nanjing's Sun Yat-Sen Memorial Botanical Garden, treated with concentrated sulfuric acid for 15 min and rinsed thoroughly with sterile water. The treated seeds were sown into trays containing a mixture of peat and perlite (1:1, v/v) in an illuminated incubator with day/night temperatures of 26 • C/20 • C and a 14 h/10 h photoperiod. When the seedlings had grown 6-8 true leaves, the roots were washed carefully with sterilized water and the seedlings were transferred into triangle bottles (100 mL) filled with 1/4 MS nutrient solution. These were kept under identical conditions for 10 days in the same illuminated incubator. Then, the seedlings were transferred into new triangle bottles with 1/4 MS nutrient solution containing 200 mM NaCl or 15% PEG 6000. The roots of the seedlings exposed to the PEG and NaCl treatments for 0, 6, and 24 h were sampled. For each treatment, nine samples collected from three triangle bottles were used as three biological replicates, respectively. All of the collected samples were frozen in liquid nitrogen immediately after harvest and stored at −80 • C until RNA extraction.

Illumina Sequencing for Transcriptome Analysis
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and the purity and integrity of the RNA were monitored using agarose gels. A NanoPhotometer ® spectrophotometer (IMPLEN, Palo Alto, CA, USA) and RNA Nano 6000 Assay Kit for the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA) were also used. Sequencing libraries were generated using a NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Lpswich, MA, USA), following the manufacturer's recommendations, with a total of 2 µg RNA per sample. Index codes were added to attribute sequences to each sample. CK corresponded to treatment for 0 h, P6 and P24 corresponded to treatments with the drought stress solution (15% PEG 6000) for 6 h and 24 h, and CL6 and CL24 corresponded to treatments with the salt solution (200 mM NaCl) for 6 h and 24 h. After cluster generation, the library preparations were sequenced on an Illumina HiSeq 2500 system. Size of reads of paired sequencing was 150 bp. And the sequencing indicator was set as 8 Gb clean data for each sample.

Bioinformatic Analysis of RNA-Seq Data
Clean reads were obtained by removing adapter sequences, low-quality reads and uncertain bases (N). Then, the reads were aligned to the H. hamabo reference genome sequence (not published, sequenced with PacBio sequencing by Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, China) using Hisat2 and Bowtie2 (v 2.2.2) with default parameters. Gene and isoform expression values were normalized as fragments per kb per million reads (FPKM) referenced to the distribution and degree of coverage that passed quality controls for alignment. Principal component analysis (PCA) was performed using the "principal" function in R software (v. 3.0.1). The differentially expressed genes (DEGs) were screened with the threshold of a false discovery rate < 0.05 and fold change ≥ 2 using the edgeR package. Gene expression pattern analysis was performed using Short Time-series Expression Miner software (STEM) on the OmicShare tools platform [18]. The functional analysis of the DEGs was performed using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Heat maps representing the gene expression levels were created using TBtools [19].

Quantitative Real-Time PCR (RT-qPCR)
To verify the sequencing results in this study, nine functional genes were selected randomly to detect their expression levels using qRT-PCR (Table S1). The first-strand cDNA was synthesized using a SuperScript II reverse transcriptase kit (TaKaRa, Dalian, China) according to the manufacturer's instructions with the same RNA sample as used for RNA-seq. The qRT-PCR was performed with a StepOnePlus RealTime PCR system (Applied Biosystems). The real-time PCR program in the study was as follows: 10 min at 95 • C, 40 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 30 s; a process at 60-95 • C was carried out to generate the melting curve. The reaction mixture in a 20 µL volume contained the following: 10 µL of 2 × SYBR Green Master Mix (Bimake, TX, USA), 1 µL of diluted cDNA and 0.4 µL of each primer (10 µM). Three repetitions were performed for each qRT-PCR analysis. The relative gene expression level was calculated using the 2 −∆∆Ct method [20]. The transcript level of 18S rRNA served as the internal standard (Table S1).

Summary of H. hamabo Sequences
To obtain an overview of the H. hamabo gene expression profiles under salt and drought conditions, cDNAs from salt-or drought-treated H. hamabo seedlings were sequenced. After assembly, low complexity, low quality and repeat reads were deleted ( Table 1). The clean reads were aligned to the reference genome sequence, and the results showed that the total mapped ratio of all samples was higher than 64.00%, indicating that although there was a little contamination of phycophyta from hydroponic environment, the mapped ratio was relatively high, which can be used for subsequent analysis ( Table 2).

Global Comparisons of Salt-Treated and Drought-Treated Transcriptomes
To understand how H. hamabo genes respond to salt or drought stress, we first compared the mRNA populations with PCA ( Figure 1a). Transcriptomes within the P group and CL group were likely to show similarity in total gene expression, respectively, so these formed two separate groups that were separated from each other and the CK.   Figure S1 and Table S2). Meanwhile, 3356 and 5139 genes were specifically regulated under drought stress at 6 h and 24 h. Based on KEGG analysis, these genes were significantly enriched in 17 pathways, including signal transduction, environmental information processing, phenylpropanoid biosynthesis, MAPK signaling pathway, organismal systems, environmental adaptation, chaperones and folding catalysts, biosynthesis of other secondary metabolites and others (Figure 1b, Supplement Figure S2 and Table S3). It's worth noting that there were 3845 and 7430 overlapping genes under the two stresses at 6 h and 24 h, respectively (Figure 1b and Supplement Table S4). The results indicated that the responsive transcripts under salt differed from those during drought, while they may share some mechanisms dependent on key genes.

Validation of Gene Expression Levels with RT-qPCR Analysis
To validate the reliability of the sequencing, nine functional genes were selected randomly, and the expression patterns were tested by RT-qPCR (Figure 2). The expression levels of these genes (Hha049880, Hha006182, Hha009560, Hha063473, Hha043982,  Hha0411252, Hha026816, Hha016168 and Hha044473) showed general agreement between The differences in gene expression at different time points under salt treatment were relatively larger than those seen during drought treatment. The comparison of DEGs responding to salt and drought was further demonstrated using a Venn analysis (Figure 1b). A total of 20,513 and 27,516 significantly changed known genes at 6 h and 24 h, respectively, were detected between the salt or drought treatment and the control libraries. Among these, 13,302 and 14,947 genes, which were significantly enriched in the pathways as biosynthesis of other secondary metabolites, phenylpropanoid biosynthesis, lipid metabolism, carbohydrate metabolism, metabolism, transcription factors, transporters and plant hormone signal transduction, were specifically regulated under salt stress at 6 h and 24 h, respectively (Figure 1b, Supplement Figure S1 and Table S2). Meanwhile, 3356 and 5139 genes were specifically regulated under drought stress at 6 h and 24 h. Based on KEGG analysis, these genes were significantly enriched in 17 pathways, including signal transduction, environmental information processing, phenylpropanoid biosynthesis, MAPK signaling pathway, organismal systems, environmental adaptation, chaperones and folding catalysts, biosynthesis of other secondary metabolites and others (Figure 1b, Supplement Figure S2 and Table S3). It's worth noting that there were 3845 and 7430 overlapping genes under the two stresses at 6 h and 24 h, respectively (Figure 1b and Supplement Table S4). The results indicated that the responsive transcripts under salt differed from those during drought, while they may share some mechanisms dependent on key genes.

Validation of Gene Expression Levels with RT-qPCR Analysis
To validate the reliability of the sequencing, nine functional genes were selected randomly, and the expression patterns were tested by RT-qPCR ( Figure 2). The expression levels of these genes (Hha049880, Hha006182, Hha009560, Hha063473, Hha043982, Hha0411252, Hha026816, Hha016168 and Hha044473) showed general agreement between the results of the Illumina sequencing analysis and RT-qPCR ( Figure 2).

Time Course RNA-seq Analysis in Response to Salt and Drought
To excavate stress-related genes under both salt and drought treatments, we classified overlapping genes into 30 clusters according to their expression patterns with STEM software, and 13 clusters showed significant enrichment trends (Figure 3a; Table S5).
Clusters 1, 2 and 3, consisting of 2779 DEGs, were mainly down-regulated under salt or drought treatment (Figure 3a). Based on the KEGG enrichment analysis, the DEGs in Clusters 1, 2 and 3 were mainly enriched in cytoskeleton proteins, carbohydrate metabolism, chromosome and associated proteins, environmental information processing, MAPK signaling pathway and other pathways related to signal transduction, biosynthesis and metabolism (Figure 3b). Among them, six DEGs (Hha087483, Hha065651, Hha003634,  Hha019082, Hha19915 and Hha065146) belonged to the ABA receptor PYR/PYL family, which act upstream in the MAPK signaling pathway (Supplement Figure S3). There were 593 genes in Cluster 5 that showed a downward trend under drought but not under salt stress, and they were mainly enriched in certain metabolism-related pathways. In addition, a total of 324 DEGs, which were enriched in transcription factors, transporters and other pathways, including R2R3-MYB (Hha082080, Hha096127 and Hha072588), ATHB

Time Course RNA-seq Analysis in Response to Salt and Drought
To excavate stress-related genes under both salt and drought treatments, we classified overlapping genes into 30 clusters according to their expression patterns with STEM software, and 13 clusters showed significant enrichment trends (Figure 3a; Table S5).
Clusters 1, 2 and 3, consisting of 2779 DEGs, were mainly down-regulated under salt or drought treatment (Figure 3a). Based on the KEGG enrichment analysis, the DEGs in Clusters 1, 2 and 3 were mainly enriched in cytoskeleton proteins, carbohydrate metabolism, chromosome and associated proteins, environmental information processing, MAPK signaling pathway and other pathways related to signal transduction, biosynthesis and metabolism (Figure 3b). Among them, six DEGs (Hha087483, Hha065651, Hha003634, Hha019082, Hha19915 and Hha065146) belonged to the ABA receptor PYR/PYL family, which act upstream in the MAPK signaling pathway (Supplement Figure S3). There were 593 genes in Cluster 5 that showed a downward trend under drought but not under salt stress, and they were mainly enriched in certain metabolism-related pathways. In addition, a total of 324 DEGs, which were enriched in transcription factors, transporters and other pathways, including R2R3-MYB (Hha082080, Hha096127 and Hha072588), ATHB (Hha057357, Hha026637), OCT (organic cation/carnitine transporter, Hha088838), ABCB (ABC transporter B family member, Hha063491) and other genes, showed a trend of continuous up-regulation in response to salt and drought in group 13. Although they behaved differently in response to drought, under salt treatment, the DEGs in Clusters 8, 9 and 12 were only down-regulated at 6 h and Cluster 4 and 11 were only down-regulated at 24 h. These findings verified the results of the PCA.

Discussion
Studies on the mechanisms of plant responses to drought and salt stress have mainly focused on model plants, in which many regulatory networks are known to be involved in these processes [5]. Recently, several woody plants were used as materials to reveal important mechanisms and identify functional genes associated with stress tolerance.

Discussion
Studies on the mechanisms of plant responses to drought and salt stress have mainly focused on model plants, in which many regulatory networks are known to be involved in these processes [5]. Recently, several woody plants were used as materials to reveal important mechanisms and identify functional genes associated with stress tolerance. Comprehensive analysis of the regulatory genes in Zygophyllum xanthoxylum under salt and osmotic tolerance showed that hormone signaling pathways, ubiquitin-mediated proteolysis machinery, protein kinases and transcription factors such as NAC were significantly changed to participate in the stress response [21]. In Camellia sinensis, a total of 2131 overlapping DEGs, which were involved in galactosyl transferase activity, tetrapyrrole binding, hydrolase activity, and biosynthetic pathways of polyphenol and caffeine, were identified under drought and salt stress conditions [22]. To learn more about how H. hamabo adapts to severe environments, we analysed the transcriptome data of seedling roots exposed to drought stress and salt treatments for 6 h and 24 h, based on previous studies [21,23]. Many genes, which were involved in biosynthesis, metabolism, transporters, plant hormone signal transduction and other critical pathways, were found specifically regulated under salt stress in this study. And thousands of genes participating in regulation of signal transduction, environmental information processing, phenylpropanoid biosynthesis, MAPK signaling pathway, environmental adaptation and others, were specifically regulated under drought stress. We also found that thousands of up/down-regulated genes related to MAPK signal transduction, certain transcription factors, transport proteins and other pathways were found to be involved in both drought and salt-stress responses in H. hamabo. These overlapping genes will become the focus of this and future research. These results provided a molecular basis for understanding the specific drought-and salt-tolerance mechanisms of H. hamabo.
MAPK cascades are reported to be involved in signal transduction under salt or osmotic stress [2]. MKK4-MPK3 and MKKK20-MPK6 cascades were found in Arabidopsis to mediate responses to osmotic stress [24]. The MKK1-MPK4 cascade is involved in saltstress signal transduction in rice by regulating the expression of transcription factors [25]. In our studies, the MAPK signal pathway was enriched under both salt and drought stresses, indicating that it is important in the tolerance mechanisms of H. hamabo. The MAP3K17/18-MKK3-MPK1/2/7/14 cascade is activated by the PYR/PYL/RCAR-SnRK2-PP2C ABA core signaling module following drought stress [2]. The six DEGs in H. hamabo belonging to the PYR/PYL family in this pathway were continuously down-regulated in response to salt and drought, which was similar to previous studies [26], indicating that abiotic stress may increase the ratio of PP2Cs: PYR/PYLs and activate the ABA signaling pathway and affect the expression of PYL, PP2C and SnRK2, leading to the activation of the MAP3K17/18-MKK3-MPK1/2/7/14 cascade to improve the tolerance of seedlings.
Signaling cascades can target transcription factors such as MYB, WRKY and NAC, which regulate downstream stress-related genes [27,28]. The MYB family is known to be involved in plant responses to environmental stress [28,29]. For example, CaMYB85 can induce drought and salt tolerance by activating stress-responsive genes in transgenic Arabidopsis [30]. OsMYB6-overexpressing rice shows higher tolerance to drought and salt than controls by increasing proline content, activating CAT and SOD and reducing the REL and MDA content [28]. Many R2R3-MYB TFs in Arabidopsis have been demonstrated to participate in responses to various stresses [10]. Based on the results of this study, 3 R2R3-MYB TF genes in H. hamabo may function as a stress-responsive transcription factor that played a positive regulatory role in response to drought and salt stress resistance. As AtMYB108 was implicated in both biotic and abiotic stresses, Hha072588, the homologous genes of AtMYB108 in H. hamabo, may be used as a candidate gene for the molecular breeding of tolerant varieties [31].
Transporters embedded in the plasma membrane, cytoplasm or nucleus are associated with various signaling pathways [32]. In this study, 10 DEGs presenting upward trends consistently under salt and drought treatment were significantly enriched in transporters including OCTs and ABCB. Localisation of OCT proteins at the tonoplast and regulation of the gene expression under stress conditions suggested a specific role in plant adaptation to environmental stress [33]. As a result, OCTs may play an important role in protecting H. hamabo from abiotic stress damage. ABC transporters, localized at the plasma membrane, have been reported to regulate the influx and efflux of ABA, which is involved in drought stress tolerance, and export ABA from shoots and roots to the guard cells [34]. Under salt stress, ABC transporters take part in regulating the concentration of Na + in the cell by utilizing the power of ATP to increase the uptake of osmoprotectants [34]. Hence, the expression level of ABC transporter family members in H. hamabo was upregulated in the roots under both salt and drought treatments, while its role under different conditions may be different.

Conclusions
In summary, genes related to stress responses were identified in the roots of H. hamabo on a whole-transcriptome scale using an Illumina sequencing platform. A total of 20,513 and 27,516 significantly changed known genes at 6 h and 24 h, respectively, were detected under salt or drought treatment. Among these, there were 3845 and 7430 overlapping genes under the two stresses at 6 h and 24 h, respectively. According to the analyses of enriched KEGG pathways and clustering of expression patterns based on these genes, MAPK signal transduction, certain transcription factors, transport proteins and other pathways were found to be involved in both drought and salt-stress responses in H. hamabo. These related overlapping genes will also become the focus of future research. These results provided a reference for identifying common factors in plant responses to different abiotic stresses and may be useful to develope genetic breeding strategies for the Malvaceae.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12040454/s1, Table S1: Primers for RT-qPCR, Table S2: Specifically regulated genes under salt stress, Table S3: Specifically regulated genes under drought stress, Table S4: Regulated genes under both salt and drought stress, Table S5: The KEGG results of DEGs in the 13 clusters, Figure  S1: KEGG enrichment of specifically regulated genes under salt stress, Figure S2: KEGG enrichment of specifically regulated genes under drought stress, Figure S3: The expression levels of PYR/PYL genes in MAKP pathway.