Physiology assay
Salinity tolerance coefficient 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 coefficient (Additional file 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 significantly 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 findings 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 significantly increased. In Te’aiqing, NaCl stress induced a significant 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 S0, S100, T0 and T100, 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 file 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.
Due to the lack of a reference genome for B. rapa L. ssp. Chinensis, BLASTx (E-value ≤ 10-5)searches were carried out to perform functional annotations with assembled unigenes against seven public databases, including NCBI non redundant protein (NR), Gene Ontology (GO), Clusters of Orthologous Groups (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), Protein family (Pfam), nucleotide databases (NT) and SwissProt public protein databases (Additional file 1: Figure S3). Out of the 214,952 unigenes, 186,301 (86.67%) unigenes were matched with the public databases, suggesting that this project has captured a majority of the pak choi transcriptome. Among these annotated unigenes, 171,134 (79.61%) unigenes annotated in the NR database, 133,774 (62.23%) unigenes annotated in the KEGG database, 142,027 (66.07%) unigenes annotated in the COG database, 118,909 (55.32%) unigenes annotated in the GO database, and 129157 (60.09%) Unigenes annotated in the SwissProt database (Additional file 1: Figure S3). Furthermore, the CDSs of 141,274 (75.83%) unigenes were predicted based on the SwissProt and Pfam public protein databases (Table 3). Of these CDSs, 3189 (0.02%) were longer than 3000 bp, 100,248 (70.96%) ranged from 500bp to 3000 bp, and 37,837 (26.78%) were shorter than500 bp (Additional file 1: Figure S4).
Identification of differentially-expressed unigenes (DEUs)
Pairwise comparison analysis for each gene was performed between Shanghaijimaocai and Te’aiqing (S0/T0 and S100/T100) or between control and NaCl-treated samples in each cultivar (T100/T0 and S100/S0). DEUs were identified by the threshold of |log2 fold-change| ≥ 2 and Padj-value ≤ 0.05. A total of 12,836 unigenes were differentially regulated in the four comparisons (Additional file 2: Table S1). Under control condition, 5798 DEUs were identified 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 identified 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).
Gene ontology (GO) analysis of DEUs
Based on GO functional annotation, a total of 8457 (65.89%) DEUs including 4314(S100/S0), 1693 (T100/T0), 3698 (S0/ T0) and 1666 (S100/T100) DEUs, were assigned into 54 Go terms consisting of 25 biological process, 16 cellular component and 13 molecular function subcategories (Additional file 1: Figure S5a-d). Among these GO terms, the top five abundant categories of biological process (‘cellular process’, ‘metabolic process’, ‘biological regulation’, ‘response to stimulus’ and ‘regulation of biological process’), cellular component (‘cell’, ‘cell part’, ‘membrane’, ‘membrane part’ and ‘organelle’) and molecular function (‘binding’, ‘catalytic activity’, ‘transporter activity’, ‘structural molecule activity’ and ‘transcription regulator activity’), were similar in the four comparisons.
Furthermore, the results of significantly enriched GO terms (Padj-value ≤ 0.05) are listed in Additional file 2: Table S2. For salinity-responsive DEUs, the enriched GO terms for DEUs of Shanghaijimaocai were assigned into 65 enriched GO terms, including 27 biological process (‘translation’, ‘oxylipin biosynthetic process’, ‘response to oxidative stress’, etc.), 8 cellular component (‘ribosome’, ‘extracellular region’ and ‘mucus layer’, etc.), and 30 molecular function (‘structural constituent of ribosome’, ‘zinc ion transmembrane transporter activity’ and ‘chitin binding’, etc.); meanwhile, the DEUs in T100/T0 were assigned into 45 enriched GO terms, including 23 biological process (‘L-proline biosynthetic process’, ‘response to abscisic acid’ and ‘lipid transport’, etc.), 3 cellular component (‘chloroplast inner membrane’, ‘integral component of membrane’ and ‘vacuole’), and 19 molecular function (‘glutamate 5-kinase activity’, ‘glutamate-5-semialdehyde dehydrogenase activity’ and ‘transaminase activity’, etc) subcategories (Additional file 2: Table S2). For DEUs between Shanghaijimaocai and Te’aiqing, the DEUs in S0/T0 were assigned into 42 enriched GO terms, including 14 biological process (‘translation’, ‘microtubule-based process’ and ‘cyanide metabolic process’, etc.), 10 cellular component (‘ribosome’, ‘cytosolic small ribosomal subunit’ and ‘lysosome’, etc.), and 18 molecular function (‘structural constituent of ribosome’, chitin binding’ and ‘cellulose binding’, etc.); the DEUs in S100/T100 were assigned into 21 enriched GO terms, including 12 biological process (‘phosphorelay signal transduction system’, ‘circadian rhythm’, ‘cellular ion homeostasis’, etc.) and 9 molecular function (‘protein serine/threonine phosphatase activity’, ‘UDP-glucose 4-epimerase activity’ and ‘voltage-gated anion channel activity’, etc.) subcategories (Additional file 2: 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 (S100/S0), 1086 (T100/T0), 2428 (S0/T0) and 1118 (S100/T100) DEUs were respectively assigned to 127, 125, 129 and 126 pathways (Additional file 2: Table S3). The four comparisons differed from each other in metabolic pathways of DEUs (Fig. 4). For the salinity-responsive DEUs, the top five pathways in Shanghaijimaocai were ribosome, nitrogen metabolism, linoleic acid metabolism, alpha-Linolenic acid metabolism and carbon fixation 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 five pathways in Te’aiqing (Fig. 4b). Among these pathways, they were all significantly enriched pathways (Padj-value ≤ 0.05)
For DEUs between Shanghaijimaocai and Te’aiqing, ribosome, photosynthesis, oxidative phosphorylation, carbon fixation in photosynthetic organisms and phagosome were the top five pathways under control condition, and only the ribosome and photosynthesis was the predominantly enriched pathway (Fig. 4c). Under salinity-treated condition, the top five 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 significantly enriched pathways (Fig. 4d).
DEUs related to signal intermediates
Stress sensing and signal transduction are an important process for plants to combat both osmotic and ionic stress. In the present study, a total of 212 DEUs involved in SOS, ABA, MAPK and CDPK signaling pathway were identified, including calcineurin B-like proteins (CBLs), CBL-interacting serine/threonine-protein kinases (CIPKs), serine/threonine-protein kinases (SRK2s), abscisic acid receptor pyrabactin resistance like proteins (PYLs), protein phosphatase 2Cs (PP2Cs), MAPKs/MPKs, MAPK kinase kinases (MAP3Ks), and CDPKs/CPKs (Additional file 2: Table S4). Among them, 96 were salinity-responsive DEUs, 59 were DEUs between Shanghaijimaocai and Te’aiqing, while the remaining 57 were salinity-responsive DEUs, which also expressed differentially between two cultivars (Additional file 2: Table S4).
Furthermore, a significant distinction in the expression profile was also observed among different gene families. All unigenes belonging to PYL gene family was notably up-regulated by NaCl treatment in both Shanghaijimaocai and Te’aiqing (Table 4). In contrast, most unigenes belonging to CBL and CPK/CRK families were obviously down-regulated by NaCl treatment in both Shanghaijimaocai and Te’aiqing. MPK family in Shanghaijimaocai and CIPK in Te’aiqing showed mostly down-regulation after NaCl treatment (Table 4), while CIPK in Shanghaijimaocai showed mostly upregulation. Under NaCl exposure, most unigenes belonging to PYLs, CIPKs, PP2Cs and SRK/KINs showed higher expression in Shanghaijimaocai than that in Te’ aiqing, while the majority of DEUs, such as CPK/CRK, MAP3K and MPKs were lower in Shanghaijimaocai than in Te’aiqing (Table 4). Moreover, the up-regulated expression of signal intermediates in S100/S0 and S100/T100 were also analyzed, including CIPKs (Unigene6888, CL13846.Contig6, CL13846.Contig1), PYL8 (CL16764.Contig5, CL16764.Contig3), and MAP3K (CL9384.Contig1) (Table 5). The PP2Cs (CL18346.Contig10, CL1555.Contig17, Unigene22896) downregulation in S100/S0 and S100/T100 were also found (Table 5).
Identification of DEUs associated with TFs
The transcriptome data identified 547 DEUs mapped to 21 different TF families, such as APETALA2/ethylene responsive factors (AP2/ERFs), basic leucine zippers (bZIPs), basic helix–loop–helixs (bHLHs), zinc finger proteins CONSTANS-LIKE (C2C2-CO-likes), ethylene-responsive factors (ERFs), golden2-likes (G2-likes), myeloblastosis proteins (MYBs), no apical meristem domain-containing proteins (NACs), nuclear transcription factor Y subunits (NF-Ys), WRKYs and zinc finger proteins (C2H2s, C3Hs and Dofs) (Additional file 2: Table S5). Among them, 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 file 2: Table S5).
Meanwhile, most unigenes homologous to C3Hs, Dofs, G2-likes, MYBs, and NF-Ys were notably up-regulated by NaCl treatment in both Shanghaijimaocai and Te’aiqing (Table 4). Besides the aforementioned five TF families, five families (ARRs, B3s, bZIPs, C2C2-Cos and MADSs) in Shanghaijimaocai and two families (NACs and WRKYs) in Te’ aiqing showed mostly up-regulation after NaCl treatment (Table 4). In contrast, six TF families (AP2/ERFs, bHLHs, GRASs, HSFs, Tifys and WRKYs) in Shanghaijimaocai and AP2/ERFs in Te’ aiqing showed mostly down-regulation after NaCl treatment (Table 4). Under NaCl exposure, most unigenes homologous to AP2/ERFs, ARFs, bHLHs, bZIPs, C2C2-CO-likes, Dofs, HSFs and MYBs showed higher expression in Shanghaijimaocai than that in Te’aiqing, while the majority of DEUs, such as ARRs, B3s, C3Hs and NACs were lower in Shanghaijimaocai than in Te’aiqing (Table 4). Moreover, the expression abundance of up-regulated TFs in S100/S0 and S100/T100 were also analyzed, including eight MYBs, three Dofs, AP2/ERFs and C2C2-CO-likes which contained two unigenes each, and ARF, ARR, bHLH, bZIP and HSF contained one unigene each (Table 5).
DEUs involved in ions transport
According to GO functional annotation, 91 DEUs were identified to highly similar with transporters that may be involved in the translocation of Na+, Cl- and K+ in plants, including potassium transporter 9 (POT9), potassium channels (AKTs, KATs), K+ efflux antiporter 6 (KEA6), cyclic nucleotide-gated ion channels (CNGCs), sodium/hydrogen exchangers (NHXs), sodium/proton antiporter1 (NHD1), cation/H(+) antiporters (CHXs), plasma membrane (PM) H+-ATPases (AHAs, PMAs, HA1), vacuolar V-type proton H+-ATPases (VHAs), pyrophosphate-energized vacuolar membrane proton pumps (AVPs), aquaporins (NIPs, PIPs, TIPs, SIP1-2), Aluminum-activated malate transporters (ALMTs), Chloride channel protein CLC-b (CLC-B), mechanosensitive ion channel proteins (MSLs) and S-type anion channel SLAHs (SLAHs) (Additional file 2: Table S6). These transporters were divided into seven categories (I-VII) based on their expression patterns (Additional file 2: Table S6). The 12 unigenes in the first category and the five unigenes in the second category were specifically up-regulated in Shanghaijimaocai and Te’aiqing roots by NaCl treatment, respectively. 15 unigenes in the third category and four unigenes in the fourth category were specifically down-regulated by NaCl treatment in Shanghaijimaocai and Te’aiqing roots, respectively. The expression of 14 unigenes in the fifth category was higher in Shanghaijimaocai roots than that in Te’aiqing. The sixth category had 9 unigenes expressed lower in Shanghaijimaocai than that in Te’aiqing. Categories I to IV were salinity-responsive DEUs whereas categories V and VI were DEUs between two cultivars. The 32 unigenes in category VII were salinity-responsive DEUs between Shanghaijimaocai and Te’aiqing (Additional file2: Table S6). Among of them, the expression of PIP2-1 (Unigene18298) gene was higher in Shanghaijimaocai than in Te’aiqing, while it did not change under NaCl treatment in both cultivars (Additional file 2: Table S6). More importantly, three major transporter genes NHX7 (CL850.Contig14), ALMT13 (CL1489.Contig2) and SLAH1 (CL10603.Contig2) were significantly up-regulated by NaCl treatment in Shanghaijimaocai, and their expressions were higher in Shanghaijimaocai than in Te’aiqing under NaCl treatment (Table 5).
Reverse transcription-quantitative PCR (RT-qPCR) validation
To verify the differential expression patterns of DEUs identified 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.