Comparative Transcriptome Analysis of NaCl and KCl Stress Response in Malus Hupehensis Rehd. Provide Insight into the Regulation Involved in Na+ and K+ Homeostasis

Background: Apple is among the most widely cultivated perennial fruit crops worldwide. It is sensitive to salt stress, which seriously affects the growth and productivity of apple trees by destroying the homeostasis of Na + and K + . Previous studies focused on the molecular mechanism underlying NaCl stress. However, signaling transduction under KCl stress has not been thoroughly studied. Results: We comprehensively analyzed the salt tolerance of Malus hupehensis Rehd., which is a widely used rootstock in apple orchards, by using RNA-Seq. Roots and leaves were treated with NaCl and KCl. Based on mapping analyses, a total of 762 differentially expressed genes (DEGs) related to NaCl and KCl stress in the roots and leaves were identied. The Gene Ontology (GO) terms were enriched in ion transmembrane transporter and oxidoreductase activity under NaCl and KCl stress. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways focused on the plant hormone signal transduction and mitogen-activated protein kinase signaling pathway. We also screened out 28 candidate genes from 762 DEGs and veried their expression by quantitative reverse transcription polymerase chain reaction (qRT-PCR). All of these enriched genes were closely related to NaCl and KCl stress and take part in mediating the Na + and K + homeostasis in M. hupehensis. Conclusions: This transcriptome analysis provides a valuable resource for elucidating the signaling pathway of NaCl and KCl stress and is a substantial genetic resource for discovering genes related to the NaCl and KCl stress response of M. hupehensis.


Background
Soil salinization is one of the most serious abiotic stress worldwide [1]. To date, more than 20% of the world's land was affected by salinization, and the trend is constantly expanding [1,2]. Apple (Malus domestica Borkh.) is one of the most valuable horticultural fruit crops cultivated worldwide. It is sensitive to salt stress, which seriously affects the growth and productivity of apple trees [3]. Therefore, strategies to improve the salt tolerance of M. hupehensis, a widely used apple rootstock in China, should be explored. However, the studies on the physiological and molecular mechanisms of plants' response to salt stress focus on model plants, such as Arabidopsis, rice, and tobacco [4][5][6]. The defense response mechanism of M. hupehensis to salt stress remains unclear.
The in uence of salt stress to the plants has been explored. Salt stress is usually caused by high concentrations of soluble sodium (Na + ) and chloride (Cl − ) ions in soil [7]. Salt stress causes primary damage, which are osmotic and ionic stress [8]. Salt stress can reduce soil water potential and change the cell osmotic potential, thereby resulting in physiological drought [9]. Excessive external and internal Na + accumulation in the cytoplasm disrupts the Na + /K + homeostasis and evokes ion toxic effects [10].
Oxidative damage is usually triggered by osmotic and ionic stress. Salt stress can also produce a large amount of reactive oxygen species (ROS), including hydroxyl radicals, hydrogen peroxide, and superoxide anion accumulation in cells; this condition can severely damage plants [11,12].
differentially expressed genes and the signaling pathways responding to NaCl and KCl stress were identi ed. In addition, we compared the genes and pathways between the NaCl and KCl stress and found 28 key regulation and function genes. This nding would be crucial in K + and Na + homeostasis of the M. hupehensis seedlings under salt stress. This study contributes to the discovery of genes related to NaCl and KCl stress and further explores the molecular mechanisms regulating salt tolerance in M. hupehensis.

Results
Analysis of RNA-Seq datasets As shown in Figure S1, the M. hupehensis seedlings were wilting seriously after salt treatment for one and 6 h compared with the control. The 30 libraries of RNA-Seq were shown in Table 1. After removing the low-quality reads and adapter sequences, 1,702,856,596 clean reads and 255.42G clean bases were obtained. In addition, the average of sample GC content was 48.61%. The sequencing error rate was only 0.03% of every sample, and the average Q20 and Q30 was 96.43% and 90.78%, respectively (Table S1). The clean reads mapped in M. hupehensis genome's ratios ranged from 32.33-84.71% among the 30 RNA-Seq libraries, and 31.54-82.68% unique reads were mapped to the reference genome (Table S2). The read1 mapped ratios were similar with read2 mapped ratios, and the positive mapped ratios were also similar with the negative mapped ratios. These results indicated that the root and leaf sample libraries data could be used for further analysis. at 6 h after KCl and NaCl treatments were more than one hour in the roots and leaves. These results indicated that the roots may be more responsive than the leaves under salt stress. DEGs quantity also increased with the extension of treatment time.
The Venn diagrams showed the distribution of similarly regulated genes in the treatment and comparison groups in the roots and leaves to further survey DEGs, which respond to different salt stress (Fig. 3). A total of 11,455, 16,516, 11,917, and 15,234 DEGs were observed in K1R, K6R, Na1R, and Na6R treatment groups, compared with the control group, respectively (Fig. 3). In the roots, 7144 DEGs were noted coresponding to NaCl stress at one and 6 h with 7528 DEGs to KCl stress. We also found 4278 DEGs collectively responded to KCl and NaCl stress at one and 6 h in the roots (Fig. 3a). In the leaves, 3009 DEGs were found responding to NaCl stress at one and 6 h, and 3830 DEGs were changed to responding to KCl stress at one and 6 h. Among these DEGs, 2295 DEGs were collectively responded to KCl and NaCl stress at one and 6 h (Fig. 3b). Moreover, 762 genes responded together to NaCl and KCl stress in the roots and leaves (Fig. 3c).

GO enrichment analysis of DEGs related to NaCl and KCl stress
To identify the DEGs functions responding to NaCl and KCl stress, we successfully analyzed the DEGs enrichment on the basis of the GO classi cations. The DEGs were separated into three categories: molecular function (MF), cellular component (CC), and biological process (BP) ( Tables 2 and 3). The functions of DEGs were similar between NaCl and KCl stress in the roots and leaves, indicating the similar signaling pathways responding to NaCl and KCl stress. 31 GO terms enriched in the three categories in the roots and leaves, and they may play major roles in plant response to NaCl and KCl stress. In the category of molecular function, the GO terms were enriched in ion (such as potassium ion, metal ion, and cation) transmembrane transporter activity, ion (such as iron and calcium) binding, oxidoreductase activity, ion channel activity, hydrolase activity, and antiporter activity. In cellular component, membrane protein complex, endomembrane system, and organelle membrane were signi cantly enriched. The enriched biological process includes signal transduction, cation and metal ion transport, response to stress, chemical and oxidative stress, cellular component organization or biogenesis. In addition, two new GO terms (ADP binding and hydrolase activity, hydrolyzing O-glycosyl compounds) were enriched in the leaves but not in roots, indicating the unique mechanism responding to salt stress in the leaves.    In the roots, plant hormone signal transduction pathway was the largest group with 138 and 179 DEGs enriched after KCl treatment for one and 6 h, whereas 141 and 151 DEGs were enriched after NaCl treatment for one and 6 h, respectively. This pathway was also enriched in the leaves for multiple DEGs.
These results suggested that plant hormone may play an important role in response to salt stress.

Identi cation of novel transcripts related to NaCl and KCl stress
To identify the key genes related to the NaCl and KCl stress, we screened them out from the 762 candidate genes, which responded together to NaCl and KCl stress in the roots and leaves by RNA-Seq ( Fig. 3). 28 candidate genes were identi ed through GO enrichment analysis and KEGG pathway enrichment analysis (Fig. 4). As shown in Fig. 4 In the second category ( Fig. 6), the kinase MdMAPKKKa, transcription factors of MdWRKY28, MdNAC56, MdMYB108, MdERF019 and MdERF109 were induced by NaCl and KCl stress in three rootstocks, resembling with RNA-SEq. In the hormone signaling pathway (Fig. 7a), the MdPP2C24 and MdPP2C51 were continue signi cantly upregulated under NaCl and KCl treatments in M26, M. hupehensis and M. zumi. This outcome was consistent with the RNA-Seq results. The expression of MdCYP82 was only induced by KCl stress in M. hupehensis. In three rootstock, MdPYL4 was upregulated in roots but downregulated in leaves. In contrast, the antioxidant enzyme genes MdPER19 and MdAPXT were induced by NaCl and KCl stress in the roots but reduced in the leaves (Fig. 7b). In accordance with RNA-Seq, the MdHIPP26, and MdPFK2 expression was all signi cantly induced by NaCl and KCl treatments in M. hupehensis and M. zumi, whereas no change in M26 (Fig. 7b).

Discussion
Increased salinity in soil has severely hampered the development of apple plantations and resulted in economic losses for orchardists [3,38,39]. High concentrations of NaCl in the soil could reduce soil water potential, change the cell osmotic potential [9], produce a mass of ROS [11,12], and disrupt the uptake of K + into plant cells, adversely affecting Na + /K + homeostasis and many metabolic pathways [8]. Excess KCl stress has a strong impact on plant growth, development, and decrease of the antioxidant enzyme activity [40]. The KCl stress could affect the Na + /K + homeostasis and cause the ion toxicity [41]. When the stress of NaCl and KCl were at the same concentration, the KCl stress caused more serious damage to plants than the NaCl stress [42]. To date, the molecular mechanism underlying plants' response to NaCl and KCl stress has only been reported in model plants. However, such mechanism was still unclear in woody plants. In our study, we performed transcriptome sequencing under NaCl and KCl stress in M. hupehensis to explore the key genes, functional families, and signaling pathways related to NaCl and KCl stress.
At present, the molecular mechanism of plants responding to NaCl stress had been reported. Three primary signaling pathways (CDPK, SOS, and MAPK pathways) were considered the main ways responding to NaCl stress in plants [20].  [25,27,28]. In our study, most genes responding to NaCl stress in M. hupehensis belong to the ion transporter families, including potassium ion, metal ion, and cation transmembrane transporters. This result indicated that the ion transport would be the most important mechanism responding to NaCl stress in M. hupehensis (Figs. 4 and 5). In addition, the MAPK pathway regulates the production of compatible osmolytes and antioxidants and may also participate in cell cycle regulation under NaCl stress [43]. The expressions of MdMAPKKKa and the antioxidant enzyme genes MdPER19 and MdAPXT in M. hupehensis were induced by NaCl stress in the roots (Fig. 5). In the woody plant pear, a RNA-Seq analysis had been conducted to explore the KEGG pathways related to NaCl stress [44]. The pathway of plant-pathogen interaction, circadian rhythm, diterpenoid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, steroid biosynthesis, and RNA degradation were enriched under NaCl stress in pear [44]. However, our study found that the KEGG pathways focused on plant hormone signal transduction, plant-pathogen interaction, MAPK signaling pathway, RNA transport, and protein processing in the endoplasmic reticulum in M. hupehensis (Tables 4 and 5). This result determines the difference between the pear and apple responding to NaCl stress.
To date, the mechanism underlying the response of woody plants to KCl stress has not been reported. Our RNA-Seq research found that the GO terms focused on the ion transporter activity (ion transmembrane transporter activity, calcium ion binding, ion binding, ion channel activity, and cation transport) and response to oxidative stress (oxidoreductase activity, hydrolase activity, response to stress, and signal transduction), as shown in Tables 2 and 3. The KEGG pathway of plant hormone signal transduction, carbon metabolism, plant-pathogen interaction, endocytosis, and MAPK signaling pathway were signi cantly enriched under KCl stress in the roots and leaves at one and 6 h (Tables 4 and 5). These results indicated that the ion transport, response to oxidative stress, plant hormone signal transduction, and MAPK signaling pathway are the most important mechanisms related to KCl stress in M. hupehensis.
The interconnected mechanism of plants in response to NaCl and KCl stress was the balance of Na + /K + in the cytoplasm [16,17,24]. Cytosolic Na + /K + ratio is one of the most important features conferring plant salt tolerance, and this trait is a potential screening tool for plant breeders [14]. We identi ed 28 key genes from the 762 candidate genes simultaneously related to NaCl and KCl stress to explore the crossing key genes regulating the Na + and K + homeostasis under NaCl and KCl stress (Fig. 4). The candidate genes were classi ed into six categories, namely, ion transporters, transcription factors, hormone signal pathway, antioxidant enzymes, MAPK signal pathway, and others.
The ion transporters were the most important families regulating the Na + and K + homeostasis. These ion transporters include Ca 2+ , K + , and H + transporters. AtCCX2 (cation/calcium exchanger) is directly involved in the control of Ca 2+ uxes between the endoplasmic reticulum and the cytosol, thereby promoting the tolerance of NaCl stress [45]. AtACA2 (calcium-transporting ATPase) overexpression could enhance the tolerance of NaCl stress by depletion of cytosolic Ca 2+ [46]. In our study, calcium transport genes, including MdACA2 and MdACA13, were signi cantly upregulated under NaCl stress and changed their expression under KCl stress (Fig. 5). MdCCX2 expression was only induced at 1 h in the leaves under NaCl stress. However, its expression was signi cantly changed in the roots and upregulated in the leaves under KCl stress. This result indicates the more important role of the gene under KCl stress than under NaCl stress. The MdCCX2 function under KCl stress still needs further research. Under excess KCl stress, K + transporter and K + channel are of great importance of regulating the K + concentration [47]. The GO term "potassium ion transmembrane transporter activity" was enriched, and the expression of potassium transporter gene MdPOT3, MdPOT6, potassium channel gene MdSKOR, K + /H + antiporter gene MdNHX2, and two-pore potassium channel gene MdTPK1 were signi cantly changed under NaCl and KCl stress.
SKOR was the K + release channel functioned to transport K + outside of the plasma membrane [48]. TPK1 was reported to be responsible for the K + compartmentation to different organelles [49]. NHX2 was K + /H + exchangers essential for active K + uptake at the tonoplast [34]. We inferred that the enrichment of potassium ion transmembrane transporters and the change of the potassium transporter genes could regulate the cytosolic Na + /K + ratio through expelling K + out of the cells or compartmentalizing K + in different organelles under NaCl and KCl stress. In addition, H + was an important cation involved in the resistance to salt stress. The expression of cation/proton exchanger gene MdCAX5 and cation/H + antiporter gene MdCHX15 was signi cantly induced by NaCl and KCl stress in the roots and leaves.
Overexpressing GmCAX1 in Arabidopsis could enhance the NaCl tolerance by accumulating less Na + in cells [50]. Knockout chx17 mutants accumulated less K + in the roots in response to salt stress and potassium starvation compared with the wild type, indicating the role of this gene in K + acquisition and homeostasis [51]. Our results indicated that the Ca 2+ , K + , and H + transporters are the crucial protein families mediating the cytosolic Na + and K + homeostasis under NaCl and KCl stress.
In the long-term process of evolution, plants have derived complex signaling pathways related to NaCl and KCl stress [43,52]. In our study, 73-85 genes were enriched in MAPK signaling pathway under NaCl and KCl stress. Our qRT-PCR results also indicated that the MdMAPKKKa expression was signi cantly induced by NaCl and KCl stress. We inferred that the MAPK signaling pathway would be the most important pathway related to NaCl and KCl stress in apple plants. In the signal transduction pathways, the transcription factors could not only receive the signal upstream but also transmit the signal to the downstream genes and played a key role in signal transduction [53]. In previous reports, overexpression of SpERFs [8], MdNAC047 [3], MdMYB44 [54], and MdSIMYB1 [55] could signi cantly improve NaCl tolerance of transgenic plants. However, the key TFs responding to KCl stress had not been reported. Our results found out the key TFs related to NaCl and KCl stress, including MdWRKY28, MdNAC56, MdMYB108, and MdERF109, with the signi cantly increased expression under NaCl and KCl stress (Fig. 5). We hypothesized that these key TFs could receive the MAPK signals and regulate the function genes for mediating the Na + and K + homeostasis under NaCl and KCl stress. The relationships between the TFs, MAPKs, and the function genes under salt stress still need further research.
Different from the RNA-Seq results related to NaCl stress in pear [44], our study found that the KEGG pathways focused on plant hormone signal transduction in M. hupehensis. Among the plant hormones related to salt stress, ABA was the most crucial one [8]. The ABA-dependent signaling pathways were involved in mediating salt stress in various plants [8]. Previous studies reported that the ABA-dependent signal transduction was regulated by protein phosphatases [56]. Singh et al (2015) reported that PP2C, as the protein phosphatases, negatively regulated the ABA signal transduction and was induced by salt stress [56]. In our study, MdPP2C51 and MdPP2C24 were signi cantly upregulated by NaCl and KCl stress, indicating their potential role in regulating the ABA signal transduction in M. hupehensis. Moreover, the MdPYL4 expression was signi cantly suppressed by NaCl and KCl stress in the leaves (Fig. 5). PYL, as the ABA receptor, could regulate ABA signaling by inhibiting PP2C under salt stress [57]. The NaCl and KCl stress also cause oxidative damage to plants [12], and antioxidant enzymes play an important role in alleviating ROS toxicity [58]. Our GO terms of oxidoreductase activity, hydrolase activity, response to stress, and signal transduction were signi cantly enriched under salt stress (Tables 2 and 3). We identi ed two key antioxidant enzyme genes, namely, peroxidase gene MdPER19 and L-ascorbate peroxidase gene MdAPXT, which were downregulated by NaCl and KCl stress in the leaves and upregulated in the roots. Overexpression of ascorbate peroxidase AtAPX in tobacco enhances NaCl stress tolerance [59]. The function of MdPER19 and MdAPXT under KCl stress requires further research.

Conclusion
In our study, we comprehensively analyzed the NaCl and KCl tolerance of Malus hupehensis Rehd. by using RNA-SEq. 28

RNA isolation
Total RNA was isolated using Trizol Reagent (Invitrogen, Carlsbad, USA) and evaluated by 1% agarose gel electrophoresis. The RNA purity and integrity were checked using NanoPhotometer spectrophotometer (Implen, München, Germany) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, USA). Qubit® RNA Assay Kit in Qubi 2.0 Flurometer (Thermo Fisher, Waltham, USA) was used to measure the RNA concentration.
Library construction and transcriptome sequencing RNA per sample (3 μg) was used as input material for RNA sample preparation. The sequencing library was generated using the NEBNext® UltraTM RNA Library Preparation Kit (NEB, Ipswich, USA). In addition, an index code was added to the attribute sequence for each sample [60]. Poly-T oligo-attached magnetic beads were utilized to purify mRNA from total RNA. Fragmentation was performed using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X) (NEB, Ipswich, USA).
Random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-) (Sangon, Shanghai, China) was used to synthesize rst strand cDNA. DNA Polymerase I and RNase H were used to degrade RNA and second strand cDNA. NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization after adenylation of 3' ends of DNA fragments. The library fragments were puri ed with the AMPure XP system (Beckman Coulter, Beverly, USA) to select cDNA fragments of preferentially 250~300 bp. Then 3 μl USER Enzyme (NEB, Ipswich, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase (NEB, Ipswich, USA), Universal PCR primers and Index (X) Primer. Finally, PCR products were puri ed by using AMPure XP system and library quality was assessed on the Agilent

Quality control
Raw data in the fastq format (original reads) were rst processed by an internal perl script. In this step, clean data (clean reads) were obtained by deleting the reads containing the adapter, the reads containing the ploy-N, and the low quality reads from the original data. Q20, Q30, and GC content of the raw data were calculated. All downstream analyses were based on high-quality cleaning data.

RNA-Seq reads mapping
Reads were aligned to the apple genome (https://www.rosaceae.org/species/malus/all) by using hierarchical indexing for spliced alignment of transcripts (Hisat2 v2.0.5), which initially removed a portion of the reads according to the quality information accompanying of each read. Then, the reads were mapped to the reference genome. Hisat2 can generate a database of splice points based on it, and the gene model annotation le was better mapped than the other non spliced mapping tools.
Quanti cation of transcript abundance and differential expression analysis The feature count v1.5.0-p3 was used to count the readings mapped to each gene. The fragments per kilobase of transcript per million mapped reads (FPKM) of each gene was calculated based on the length of the gene and the reading count mapped to the gene, which was used to determine the expression levels for mRNAs by calculating FPKM.
Differential expression analysis was performed using the DESeq2 R package (1.16.1) which provides statistical routines for determining differential expression in digital gene expression data by using a model based on the negative binomial distribution [61]. Differential expression analysis of two conditions was carried out using the edgeR R package (3.18.1). The resulting P-values were tted using Benjamini and Hochberg's approach for controlling the false discovery rate [62]. Genes with a |log 2 (fold change)|>0 and padj<0.05 in sequence counts between libraries were considered signi cantly differentially expressed using DESeq2 compared with control.
GO and KEGG enrichment analysis of DEGs GO enrichment analysis of DEGs was implemented by the cluster Pro ler R package, in which the gene length bias was corrected [63]. GO terms with corrected P-value <0.05 were considered signi cantly enriched by DEGs. KEGG is a database resource which for understanding high-level functions and utilities of the biological system (cell, organism and ecosystem) from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies (http://www.genome.jp/kegg/). ClusterPro ler R package was used to test the statistical enrichment of differential expression genes in KEGG pathways [64].
Validation of RNA-Seq gene expression by using qRT-PCR qRT-PCR analysis was performed to validate the expression levels of candidate DEGs revealed by the transcriptome data. The samples were obtained and conserved in liquid nitrogen after 100 mM NaCl and 100 mM KCl treatments for 0, 1h, 6h, 5d and 10day respectively. Total RNA was extracted with RNAprep Pure Plant Plus Kit (Tiangen, Beijing, China), and the rst-strand cDNA was synthesized from 2μg of total RNA using a 5X All-In-One RT MasterMix (abm, Sydney, Australia). All primers were synthesized at Sangon Biotech (Sangon, Shanghai, China) Co., Ltd. SYBR (Roche, Basel, Swiss) was used as the intercalating dye, and the β-actin gene of Malus was used as the inner reference gene. Three replicates were analyzed using a LightCycler® 480 (Roche, Basel, Swiss) instrument. The 10 μl reactions was 95°C for 10 min, followed by 45 cycles of 95°C for 10 s, 58°C for 10 s and 72°C for 10 s, nally 72°C for 10 min.
The relative expression of candidate DEGs was calculated using the 2 −ΔΔCt method. The experiment was repeated three times.

Experimental design and statistical analysis
All experiments were repeated thrice according to a completely randomized design. The data were analyzed by ANOVA followed by Fisher's least signi cance difference or Student's t-test analysis. Statistically signi cant differences were indicated by P-value <0.05 and extremely signi cant differences were indicated by P-value <0.01. Statistical computations were conducted using SPSS (IBM, Armonk, USA).

Availability of data and materials
The data sets supporting the results of this article were included within the article and its additional les.
The reads produced in this study had been deposited in the National Center for Biotechnology Information (NCBI) SRA database with accession number PRJNA588566. Access to the data was available upon publication at http://www.ncbi.nlm.nih.gov/sra/.

Competing Interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest. Authors' Contributions C. W. and Y. T. planned and designed the research. Y. L. and X. Z. performed experiments, conducted eldwork and analyzed data etc. Y. L., X. Z., and C. W. wrote the manuscript. S. Y. and C. M. contributed to edit the manuscript. All authors read and approved the nal manuscript.       DEGs analyses between treatments and control in apple. The x-axis represents the log2fold change, and the y-axis represents the −log10 (padj). Red dots represent the upregulated, and green dots represent the