Comparative Transcriptomics Reveals Osmotic and Ionic Stress Key Genes Contributing To Salinity Tolerance Differential in Two Pak Choi Cultivars

Background: Pak choi is an important leafy vegetable crop. Salinity is among the most harmful agents that negatively inuence pak choi yield. However, the mechanism of salinity tolerance in pak choi has not been well understood. In this study, the root transcriptomics of two cultivars differing in salinity tolerant, Shanghaijimaocai (S, salinity tolerant) and Te’aiqing (T, salinity sensitive), were investigated under 0 and 100 mM NaCl treatments. Results: Using de novo assembly, 214,952 assembled unigenes were generated. Totals of 6765, 2454, 2451 and 5798 differentially expressed unigenes (DEUs) were identied in comparison of S 100 /S 0 , T 100 /T 0 , S 0 /T 0 and S 100 /T 100 , respectively. Shanghaijimaocai is more sensitive to NaCl stress than Te’aiqing in terms of root transcriptomics. Based on GO and KEGG pathway analysis, several osmotic and ionic stress-related genes including MP3K18, PYL8, PP2C15/16/49, ARF2, bHLH112, bZIP43, COL5, CDF1/3, ERF25/60, HSFA6, MYBS3/59/92/CCA1/PHL5, POD21, GOLS7, CIPK4/7/12, NHX7, SLAH1 and ALMT13, showed higher expression in Shanghaijimaocai than in Te’aiqing. These genes, therefore, might be contributed to the difference in salinity tolerance. Moreover, the physiological shift of peroxidase activity was in accordance with dynamic transcript prole of the relevant unigenes. Conclusions: We determined digital expression prole and discovered a broad survey of unigenes associated with the difference of salinity tolerance between Shanghaijimaocai and Te’aiqing. These ndings would be useful for further functional analysis as potential targets to improve resistance to salinity stress via genetic engineering.

from the raw reads. Then, the clean reads from all the eight samples were assembled using Trinity pairedend assembly method [58] with default settings for all parameters, and mapped to transcripts. The transcripts which would be further clustered and eliminated into non-redundant unigenes using TIGR Gene Indices clustering tools (TGICL) [59]. For the study of multi-samples, the previous step using the TGICL for down-stream analysis of the unigenes was repeated. The processed unigenes were referred to as "All-unigenes." The All-unigenes were divided into two classes. One is clusters, several unigenes with over 70% similarity among them (pre x is ''CL'', followed by the gene family ID), and the other unigenes were singletons (pre x is ''Unigene'').

Identi cation and functional enrichment analysis of DEGs
All clean reads were mapped to the unigene sequences using Bowtie2 program (v2.2.5, http://bowtiebio.sourceforge.net/bowtie2/index.shtml). Expected count tables for each gene were obtained by RSEM (v1.2.8, http://deweylab.biostat.wisc.edu/rsem/rsem-calculate-expression.html). Gene expression levels of each library were normalized to FPKM (fragments per kilobase of transcript per million mapped reads). Based on the negative binomial distribution, read counts were used to determine the differential expression analysis among four groups (S 0 , S 100 , T 0 and T 100 , two biological replicates per group) by the DEseq 2 package as described by Love et al. [60]. The gene with |log 2 fold-change|≥2 and P adj -value ≤ 0.05 was considered as differentially expressed between two groups. The identi ed DEGs were subsequently carried into GO and KEGG pathway enrichment analysis with Phyper in R package by a corrected P-value ≤ 0.05.

RT-qPCR analysis
Total RNA was extracted from three biological replicates of roots samples of both pak choi cultivars exposed to 0 and 100 mM NaCl treatments, respectively. Then the rst strand cDNA synthesis was performed using Prime Script ® RT reagent Kit (Takara, Dalian, China). RT-qPCR was conducted using a SYBR Premix EX TaqKit (Takara) in a 20μl reaction mixture on an ABI7300 (Applied Biosystems, Foster City, CA, USA) according to the method described by Yu et al. [61]. The primers of all selected genes were designed using Beacon Designer 7.0 software (Premier Biosoft International, USA) (Additional le 2: Table   S9). The equation ratio=2 -ΔΔCT was applied to calculate the relative expression level of the selected genes using Actin gene as an internal control [62]. Data statistical analysis with IBM SPSS Statistics Version 25 software (Armonk, NY, USA) was performed using Duncan's multiple range test at the P <0.05 level of signi cance.

Results
Physiology assay Salinity tolerance coe cient worked well in screening the cultivars for salinity-tolerance based on salinity tolerant indices. It well showed the results of the comparison between the control and salinity-treated cultivars for each index. When pak choi seedlings were cultivated under NaCl conditions, the LFW, LDW, RFW and RDW in two cultivars showed considerable variability in salinity tolerance coe cient (Additional le 1: Figure S1). Compared with the values under normal conditions, the values of LFW, LDW, RFW and RDW in Shanghaijimaocai were higher under 100 mM NaCl treatment (STC>1), while in Te'aiqing, they were lower (STC<1). Furthermore, the values of LFW, LDW, RFW and RDW in two cultivars were lower under 200 mM NaCl treatment (STC<1). Therefore, 100 mM NaCl was used for physiological and transcriptome analyses in the present study.
Since K + and Na + uptake and sequestration are considered as among the key components differentiating between sensitive and tolerant genotypes, the K + and Na + concentrations in leaves of two pak choi cultivars were measured (Fig. 1). The Na + concentrations in two cultivars were no difference in control.
The control leaves K + concentration showed higher quantity in Te'aiqing than in Shanghaijimaocai, whereas, salinity-stressed leaves K + and Na + concentrations showed signi cantly higher quantity in Shanghaijimaocai than in Te'aiqing. Moreover, the cultivar Shanghaijimaocai has higher K + /Na + value than the cultivar Te'aiqing under NaCl treatment. These ndings indicate that Shanghaijimaocai differed in the properties of Na + (K + ) translocation and accumulation from Te'aiqing.
To investigate whether the cellular antioxidant defence system was activated, the activities of key antioxidant enzymes including superoxidase dismutase (SOD), peroxidase (POD), catalase (CAT) and ascorbate peroxidase (APX) were tested. In Shanghaijimaocai, under NaCl stress, the SOD (Fig. 2a) and APX ( Fig. 2d) activity obviously decreased compared with the control, whereas POD (Fig. 2b) and CAT (Fig. 2c) activities signi cantly increased. In Te'aiqing, NaCl stress induced a signi cant increase in POD activity compared with the control, while the SOD, CAT and APX activities were unchanged (Fig. 2).
Sequencing, de novo transcriptome assembly and annotation To obtain an overview of pak choi transcriptome in roots, and identify candidate genes involved in response to NaCl exposure, eight cDNA libraries prepared from two replicates per treatment were sequenced using a BGISEQ-500 platform. Approximately 46.63−49.08 million raw reads were produced for the eight libraries through RNA sequencing (Table 1). After removing low quality reads and reads containing adapter sequences, a total of 87.64, 88.25, 87.13 and 85.44 million clean pair-end reads remained for S 0 , S 100 , T 0 and T 100 , respectively. The clean reads from each library were assembled and then merged into 214,952 all-unigenes as the reference transcripts of pak choi (Table 2). Of all-unigenes, the total length was 267,453,003 bp, mean length was 1244 bp, the N50 was 1888 bp, which were longer than that obtained in pak choi [28,29], implying a good assembled quality of the transcriptome in the present study. The length distribution of all unigenes were as shown in Additional le 1: Figure S2, the length of assembled unigenes were mostly ranged from 200 to 500 nt accounted for 32.23%, 500 to 2000 nt accounted for 48.43% and 41,582 unigenes (19.34%) with length > 2000 nt. These results suggested that the quality of this assembled unigenes was credible.

Identi cation of differentially-expressed unigenes (DEUs)
Pairwise comparison analysis for each gene was performed between Shanghaijimaocai and Te'aiqing (S 0 /T 0 and S 100 /T 100 ) or between control and NaCl-treated samples in each cultivar (T 100 /T 0 and S 100 /S 0 ). DEUs were identi ed by the threshold of |log 2 fold-change| ≥ 2 and P adj -value ≤ 0.05. A total of 12,836 unigenes were differentially regulated in the four comparisons (Additional le 2: Table S1). Under control condition, 5798 DEUs were identi ed between two cultivars, while this value reached 2451 DEUs under NaCl-treated condition (Fig. 3a). Of which, only 354 unigenes were common ones (Fig. 3b). Compared with the control, 6765 DEUs including 2298 up-regulated and 4467 down-regulated unigenes were differentially regulated in Shanghaijimaocai in NaCl exposure. Whereas, only 2454 NaCl-treated DEUs were identi ed in Te'aiqing, including 1195 up-regulated and 1259 down-regulated unigenes. Among them, only 242 unigenes were common salinity-responsive genes in two cultivars (Fig. 3b).
Furthermore, the results of signi cantly enriched GO terms (P adj -value ≤ 0.05) are listed in Additional le 2: Table S2.  Table S2).

KEGG metabolic pathway analysis of DEUs
Based on Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology-Based Annotation System, a total of 2700 (S 100 /S 0 ), 1086 (T 100 /T 0 ), 2428 (S 0 /T 0 ) and 1118 (S 100 /T 100 ) DEUs were respectively assigned to 127, 125, 129 and 126 pathways (Additional le 2: Table S3). The four comparisons differed from each other in metabolic pathways of DEUs (Fig. 4). For the salinity-responsive DEUs, the top ve pathways in Shanghaijimaocai were ribosome, nitrogen metabolism, linoleic acid metabolism, alpha-Linolenic acid metabolism and carbon xation in photosynthetic organisms; (Fig. 4a); glucosinolate biosynthesis, Ubiquinone and other terpenoid-quinone biosynthesis, Isoquinoline alkaloid biosynthesis, Phenylalanine, tyrosine and tryptophan biosynthesis and Carotenoid biosynthesis were the top ve pathways in Te'aiqing ( Fig. 4b). Among these pathways, they were all signi cantly enriched pathways (P adj -value ≤ 0.05) For DEUs between Shanghaijimaocai and Te'aiqing, ribosome, photosynthesis, oxidative phosphorylation, carbon xation in photosynthetic organisms and phagosome were the top ve pathways under control condition, and only the ribosome and photosynthesis was the predominantly enriched pathway (Fig. 4c). Under salinity-treated condition, the top ve pathways were circadian rhythm-plant, glucosinolate biosynthesis, 2-Oxocarboxylic acid metabolism, amino sugar and nucleotide sugar metabolism and C5-Branched dibasic acid metabolism; among these pathways, they were all signi cantly enriched pathways ( Fig. 4d).

Reverse transcription-quantitative PCR (RT-qPCR) validation
To verify the differential expression patterns of DEUs identi ed in our transcriptome data, a total of ten candidate unigenes were selected for RT-qPCR analysis. The results indicated that the expression patterns of nine out of ten DEUs showed good agreement with the RNA-seq-based gene expression patterns, and only one gene (MYBS3) was not well consistent with the results of sequencing (Fig. 5). For example, the NaCl-responsive up-regulated DEUs also indicated relatively high expression level in Shanghaijimaocai under NaCl treatment in RT-qPCR analysis (Fig. 5). For DEUs between Shanghaijimaocai and Te'aiqing, the expression levels of NHX7, SLAH1, ALMT13, MYB59, CDF3, ERF60, PYL8, CIPK7 and POD21 genes were higher in Shanghaijimaocai than that in Te'aiqing under NaCl treatment (Fig. 5). These results suggested that the data detected by RNA-Seq analysis of this study provide reliable inference.

Discussion
Comparison of gene expression in roots between Shanghaijimaocai and Te'aiqing RNA sequencing (RNA-seq) analysis has been widely utilized to reveal the mechanism of salinity tolerance in some plant species [30][31][32]. To elucidate the mechanisms involved in cultivar differences in salinity tolerance in pak choi, root gene expression pro les of Shanghaijimaocai and Te'aiqing were rst analyzed under control and NaCl-treated conditions by RNA-Seq. Totally, 6765, 2454, 5798 and 2451 unigenes exhibited signi cant differential expression in comparison of S 100 /S 0 , T 100 /T 0 , S 0 /T 0 and S 100 /T 100 , respectively (Fig. 3). Among them, NaCl-induced 2454 unigenes in the roots of Te'aiqing, whereas 6765 unigenes were identi ed between 0 and 100 mM NaCl treatments in Shanghaijimaocai, which was 2.8-fold higher than that of Te'aiqing ( Fig. 3a), thereby suggesting that Shanghaijimaocai is more sensitive to salinity stress than Te'aiqing. Furthermore, between the two cultivars, 67.6% (NaCl-free) and 57.0% (NaCl-treated) DEUs showed higher expression in Shanghaijimaocai (Fig. 3a). Moreover, only 242 unigenes were common salinity-responsive genes in two cultivars (Fig. 3b). These results suggested that the two cultivars differed in the molecular mechanisms of roots adaption to salinity stress tolerance.
Potential DEUs may determine differences between Shanghaijimaocai and Te'aiqing in salinity-tolerance

Osmotic stress-related DEUs
High osmotic pressure on plant roots is one of the problems plants face when growing in salinized soils.
The rst response of plants to salinity is signal perception, and relay of this information via signal transduction pathways that eventually re-establish cellular osmotic and ionic homeostasis [1,4,7].
Ca 2+ /CDPK signaling pathways have been shown to be involved in maintaining osmotic homeostasis [4,7]. In the current study, 52 DEUs were identi ed to be involved in CDPK signaling pathways (Additional le 2: Table S4). Among of these unigenes, the majority of unigenes encoding CDPKs (e.g. CPK14, CPK28 and CPK32) were down-regulated by NaCl treatment in Shanghaijimaocai (85.2%) and Te'aiqing (88.2%), respectively (Table 4). Moreover, most unigenes (71.4%) encoding CDPKs showed lower expressions in Shanghaijimaocai than in Te'aiqing under NaCl treatment (Table 4). Thus, the CDPK may not be the key pathway to help ght against salinity stress in the two pak choi cultivars.
MAPK cascades and phosphatidic acid are also involved in the regulation of osmotic, ionic, and ROS homeostasis under salinity stress [4,7,14]. In the present study, 27 MAPK unigenes were found to expression altered in four comparisons (Additional le 2: Table S4). Of these, more importantly, one transcript encoding MAP3K18 was up-regulated by NaCl treatment in Shanghaijimaocai, but in Te'aiqing, it was unchanged, and it showed higher in Shanghaijimaocai than in Te'aiqing under NaCl treatment ( Table 5). The overexpression of the AtMAP3K18 enhanced the resistance to drought in transgenic Arabidopsis plants [33]. Thus, NaCl induced up-regulation of MAP3K18 may contribute to salinity tolerance in Shanghaijimaocai.
Overexpression of AtPYL8 enhances the resistance to drought and osmotic stress by positively regulating ABA signaling in transgenic Arabidopsis [37]. The PP2C15 (CL18346.Contig10), PP2C16 (Unigene22896) and PP2C49 (CL1555.Contig17) were found to be down-regulated by NaCl treatment in Shanghaijimaocai, while in Te'aiqing, they were unchanged, as well as they showed lower expressions in Shanghaijimaocai than in Te'aiqing under NaCl treatment (Table 5). These results indicated that the PYL8, PP2C15, PP2C16 and PP2C49 genes may be the critical factors in determining the contrasting salinity tolerance in two cultivars.
TFs, which are crucial components in osmotic stress mediated signaling pathways, are generally phosphorylated by protein kinases and then activate the expression of numerous downstream salinity stress-responsive genes [1,4,13,16,17]. In this study, we identi ed 547 DEUs mapped to 21 TF families, among which 284 were salinity-responsive DEUs, 154 were DEUs between Shanghaijimaocai and Te'aiqing, while the remaining 103 were not only salinity-responsive but also differentially expressed between the two cultivars (Additional le 2: Table S5). These results indicated that the two pak choi cultivars differed in the expression of TFs that may contribute to the difference in salinity tolerance.
Under osmotic stress, ROS as toxic products, which result in oxidative damage and cell death [1].
Antioxidant defence system can alleviate oxidative damage in response to salinity/oxidative stress [1,4,21]. Our physiological work showed that the antioxidative enzymes, POD (Fig. 2b) and CAT (Fig. 2c) activities signi cantly increased by NaCl in Shanghaijimaocai, while in Te'aiqing, only POD activity was increased by NaCl (Fig. 2). Moreover, the POD, CAT and APX activities were obviously higher in Shanghaijimaocai than in Te'aiqing under 100 mM NaCl treatment. Meanwhile, our transcriptome analysis also showed that 70 antioxidative enzymes-related DEUs including APXs, CATs, PODs, SODs, GPXs and GRs, non-enzymatic antioxidants-related DEUs including GSTs and GSHs (Additional le 2: Table S7) were identi ed in pak choi. Of these, the majority of PODs were up regulated by NaCl treatment in both Shanghaijimaocai and Te'aiqing. It was well consistent with our physiological data (Fig. 2b). Importantly, we also found that the all POD unigenes showed higher in Shanghaijimaocai than in Te'aiqing (Additional le 2: Table S7). Of these, POD21 (Unigene3772) gene was up regulated by NaCl in Shanghaijimaocai, while in Te'aiqing, it was unchanged ( Table 5), suggesting that POD21 may contribute to the cultivar difference in pak choi salinity tolerance.

Ionic stress-related DEUs
Ionic stress (also called "ionic imbalance" or "ion toxicity") is another of the problems plants face when growing in salinized soils. Ionic stress (mainly Na + ) triggers an increase in cytosolic Ca 2+ , and thereafter, Ca 2+ -binding proteins further activate downstream signaling pathways [4,7]. The Ca 2+ -SOS signaling is well known as the central salinity excretion system that helps plants against Na + accumulation [48]. The genes for SOS signaling pathway (CBLs-CIPKs as SOS3-SOS2) have shown to provide enhanced resistance to salinity via positive regulation of the plasma membrane-localized Na + /H + antiporter SOS1/NHX7 [7,49,50]. We found that the majority of CBL unigenes were down-regulated by NaCl treatment in both Shanghaijimaocai (71.4%) and Te'aiqing (100%) ( Table 4; Additional le 2: Table S4). The CBLs were not DEUs as judged by our criteria between the two cultivars under NaCl treatment (Table  4). This result suggested that CBLs may not be involved in the difference of salinity tolerance between the two cultivars. Most unigenes (62.5%) belonging to CIPK genes (CIPK3, CIPK4, CIPK7, CIPK10, CIPK11, CIPK12, CIPK 21 and CIPK25) were found to be up-regulated by NaCl treatment in Shanghaijimaocai, while in Te'aiqing only one CIPK1 transcript was up-regulated (Table 4). Moreover, CIPK4 (Unigene6888), CIPK7 (CL13846.Contig1, CL13846.Contig6) and CIPK12 (CL16987.Contig7) showed higher expressions in Shanghaijimaocai than in Te'aiqing (Table 5). Thus, we speculate that three CIPKs (CIPK4, CIPK7 and CIPK12) might contribute to the difference of salinity tolerance between two pak choi cultivars.
The CBL-CIPK modules can also interact with and regulate the activity of many ion transporters [9,10].
Previous studies have reported that a number of transporters can control Na + , Cland K + transport, which are essential for ion homeostasis to improve salinity tolerance in plants [50,51]. Based on the physiological data (Fig. 1), we reasoned that differential root uptake and translocation capacity of K + and Na + may be main factors determining the contrasting of salinity tolerance in Shanghaijimaocai and Te'aiqing. Within the pak choi transcriptome, we identi ed 128 transporter-encoding DEUs, among which 65 (including 24 up and 41 down-regulated) were salinity-responsive DEUs in Shanghaijimaocai, and only 31 (including 16 up and 15 down-regulated) were salinity-responsive DEUs in Te'aiqing (Additional le 2: Table S6). Some transporters including NHXs, ALMTs, MSLs and SLAHs have been veri ed to involve in the cellular Na + and Clexclusion or sequestration in plants, and thus enhanced the tolerance to salinity [52][53][54]. In this study, 33 DEUs that showed similarity to NHXs (6 unigenes), ALMTs (4 unigenes), MSLs (11 unigenes) and SLAHs (12 unigenes) (Additional le 2: Table S6). Of them, most ALMT13 (100%), MSLs (75%), SLAHs (100%) and NHXs (60%) unigenes showed higher expressions in Shanghaijimaocai than in Te'aiqing. We speculate that the salinity tolerance mechanism of Shanghaijimaocai may be the inherent characteristics. Moreover, NHX7 (CL850.Contig14), ALMT13 (CL1489.Contig2) and SLAH1 (CL10603.Contig2), whose expressions were not affected by NaCl in Te'aiqing, were signi cantly upregulated by NaCl treatment in Shanghaijimaocai, and their expressions were higher in Shanghaijimaocai than in Te'aiqing under NaCl treatment (Table 5). NHX7, ALMT13 and SLAH1 of pak choi were homologous with AtNHX7, AtALMT13 and AtSLAH1 in Arabidopsis, respectively. Overexpression of the plasma membrane-localized transporter genes AtSOS1/AtNHX7 [10,52], (mediating Na + e ux), and AtSLAH1 [54], (mediating Clextrusion), resulted in increased salinity tolerance in transgenic plants. These results indicated that NHX7 and SLAH1 could contribute to elevate salinity tolerance in Shanghaijimaocai. Besides SLAC1, ALMT transporters localize to the plasma and/or vacuolar membrane and might also contribute to Clexclusion or sequestration [53], indicating that ALMT13 may be responsible for the Cldetoxi cation, and consequently, increasing salinity tolerance in Shanghaijimaocai.
The cellular balance between K + and Na + is among the key salinity-tolerance mechanisms during salinity stress [55,56]. In this study, 11 DEUs (homologous with KEA6, AKT1 and POT9) were identi ed as K + transporter or K + channel protein (Additional le 2: Table S6). Of them, two AKT1 were down-regulated by NaCl treatment in Shanghaijimaocai (Unigene19806) and Te'aiqing (CL4.Contig8), respectively. Three POT9 were not affected by NaCl treatment in Shanghaijimaocai roots, as well as they were lower expression in Shanghaijimaocai than in Te'aiqing roots under NaCl treatment. These results indicated that AKT1 and POT9 may not involve in the difference of salinity-tolerance between the two cultivars. Two KEA6 (CL6150.Contig22, CL6150.Contig12) and one KAT1 (CL2130.Contig1) were speci cally upregulated by NaCl treatment in Shanghaijimaocai, but unchanged between the two cultivars exposed to NaCl. This result indicated that KEA6 and KAT1 may be related to NaCl-response in roots of Shanghaijimaocai, but not be involved in the difference of salinity-tolerance between the two cultivars.
Proton pumps including plasma membrane H + -ATPase, vacuolar membrane H + -ATPase and H +pyrophosphatase, play an essential role in cellular Na + extrusion or sequestration into vacuoles [10,23].

Conclusions
Our study is the rst systematic report of root transcriptome changes between two contrasting pak choi cultivars in response to NaCl stress. The sequence was assembled into 214,952 unigenes. Next, a total of 6765, 2454, 2451 and 5798 DEUs were identi ed during S 100 /S 0, T 100 /T 0 , S 0 /T 0 and S 100 /T 100 comparison, respectively. The two cultivars differed in the molecular mechanisms in response to NaCl stress. Shanghaijimaocai is more sensitive to NaCl stress than Te'aiqing in terms of root transcriptomics. Moreover, many candidate functional genes including four protein phosphorylation cascades (CIPK4/7/12, MP3K18), four ABA signaling genes (PYL8, PP2C15/16/49), 14 TF genes (ARF2, bHLH112, bZIP43, COL5, CDF1/3, ERF25/60, HSFA6 and MYBS3/59/92/CCA1/PHL5), three ion transporter genes (NHX7, SLAH1 and ALMT13), one antioxidant enzyme gene (POD21) and one osmoprotectant-related gene (GOLS7), which are involved in salinity tolerance mechanism, were found to be up-regulated by NaCl treatment in Shanghaijimaocai, as well as they showed higher expressions in Shanghaijimaocai than in Te'aiqing under NaCl stress. Thereby, these genes might be involved in cultivar difference in salinity tolerance of pak choi (Fig. 6), and these should be studied further. Our ndings could provide information for further investigations of key genes and molecular mechanisms of salinity tolerance in pak choi.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
All the authors have signed the consent form.

Competing interests
The authors declare that they have no competing interests for this research.  Tables   Table 1 Overview of raw and clean reads in two pak choi cultivars exposed to 0 or 100 mM NaCl for two Note: GC, is short for guanine-cytosine