Meta-Analysis of Salt Stress Transcriptome Responses in Different Rice Genotypes at the Seedling Stage

Rice (Oryza sativa L.) is one of the most important staple food crops worldwide, while its growth and productivity are threatened by various abiotic stresses, especially salt stress. Unraveling how rice adapts to salt stress at the transcription level is vital. It can provide valuable information on enhancing the salt stress tolerance performance of rice via genetic engineering technologies. Here, we conducted a meta-analysis of different rice genotypes at the seedling stage based on 96 public microarray datasets, aiming to identify the key salt-responsive genes and understand the molecular response mechanism of rice under salt stress. In total, 5559 genes were identified to be differentially expressed genes (DEGs) under salt stress, and 3210 DEGs were identified during the recovery process. The Gene Ontology (GO) enrichment results revealed that the salt-response mechanisms of shoots and roots were different. A close-knit signaling network, consisting of the Ca2+ signal transduction pathway, the mitogen-activated protein kinase (MAPK) cascade, multiple hormone signals, transcription factors (TFs), transcriptional regulators (TRs), protein kinases (PKs), and other crucial functional proteins, plays an essential role in rice salt stress response. In this study, many unreported salt-responsive genes were found. Besides this, MapMan results suggested that TNG67 can shift to the fermentation pathway to produce energy under salt stress and may enhance the Calvin cycle to repair a damaged photosystem during the recovery stage. Taken together, these findings provide novel insights into the salt stress molecular response and introduce numerous candidate genes for rice salt stress tolerance breeding.


Introduction
Since plants are sessile, they are forced to continuously face a multitude of biotic and abiotic stresses during their lifespan. Abiotic stresses, such as salt, drought, and low temperature, seriously threaten the growth and agricultural productivity of plants [1]. Of the abiotic stresses, salt stress is a serious threat to crops' yield worldwide. According to the FAO Land and Plant Nutrition Management Service, at least 6% of the world's land is affected by salt stress to different degrees [2]. In recent years, environmental pollution and climate change have intensified the adverse effects of salt stress through raising soil salinity [3]. Previous investigations have clarified that salt stress impairs plants in the form of osmotic stress and oxidative stress [4]. Osmotic stress breaks the selective permeability of the cell membrane and ion homeostasis due to a high-salinity external condition and an excessive accumulation of Na + and Cl − in plant cells, which seriously impede the absorption of water and

Plant Materials
The 'Nipponbare' rice (O. sativa ssp. japonica) was chosen for the quantitative real-time RT-PCR (qRT-PCR) verification of randomly selected DEGs. After 2 days of germination in water at 37°C, seeds were grown in containers with sponges as supporting materials in Yoshida solution with 60% relative humidity and with a light and temperature regime of 14 h/10 h, light/dark, 30°C/22°C. Three-leaf stage seedlings were transferred to 200 mM NaCl Yoshida solution for salt treatment. Then, the roots of treatment/control seedlings were collected at 0, 3, 6, 12, and 24 h for RNA extraction. In this study, 0 h was the control group, and 3 h, 6 h, 12 h, and 24 h were treatment groups. Three biological replicates were adopted for each group. For each biological replicate, 15 seedlings were collected and mixed. Totally, 15 RNAs (three control groups and 12 treatment groups) were extracted using the TRIzol method and all RNAs were reverse-transcribed into cDNAs using the PrimeScript RT reagent Kit (TakaRa, Dalian, China).

Data Collection and Meta-Analysis
The series matrix file of the TNG67 genotype (rice subspecies indica, salt-tolerant, 48 datasets of the roots and shoots) and the TCN1 genotype (rice subspecies japonica, salt-sensitive, 48 datasets of the roots and shoots) were obtained from the Gene Expression Omnibus (GEO) repository (accession number GSE76613; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76613). Within the two genotypes, 0-h, 3-h, 24-h, and recovery 24-h roots and shoots from three-leaf stage seedlings after 250 mM NaCl treatment were used in this study. Detailed information (experiment design, transcriptome analysis, array information, data processing, and platform ID) of GSE76613 can be obtained from the GEO repository, and this information is partly summarized in Table S1. Raw data analysis and ID conversion were performed using the R language package [25,26]. A difference analysis was carried out with the limma package in the Bioconductor package (http://www.bioconductor. org/) [26]. DEGs were considered those with a |log 2 fold change| >1.5 after normalization to the control and with significant results for the t-test (p value <0.05) based on six replicates (3 biological repeats x 2 technical repeats) for each treatment compared with the control [18].

Gene Ontology (GO) Enrichment Analysis
To understand DEGs' functions, a GO enrichment analysis of core salt-responsive DEGs (common DEGs at 3 h and 24 h of both genotypes in shoots or roots) was implemented by the GOseq R package (https://www.bioconductor.org/packages/release/bioc/html/goseq.html) based on the Wallenius noncentral hyper-geometric distribution [27]. GO terms with a corrected p value of less than 0.001 were considered significantly enriched.

TF, TR, and PK Identification and MapMan Analysis
For the identification of TFs, the DEG sequences were searched against the Plant transcription factor database (PlantTFDB 4.0, http://planttfdb.cbi.pku.edu.cn/) with an E-value cut off of ≤10 −5 [28]. For the identification of transcriptional regulators (TRs) and protein kinases (PKs), all of the core salt-responsive sequences in roots and shoots were analyzed by iTAK software (http://itak.feilab.net/ cgi-bin/itak/index.cgi) [29]. The average log 2 fold change values of 3 biological repeats x 2 technical repeats for each treatment compared with the control at different points were displayed using MapMan 3.6.0 (https://mapman.gabipd.org/) [18,30].

Sequence Alignment and Gene Comparison
Genes that have previously been identified as salt-responsive genes were obtained based on previous papers [8,23,31]. Gene ID conversion between different genome versions (such as 9311 and Nipponbare) was conducted by diamond software (https://ab.inf.uni-tuebingen.de/software/) using the blastp method with the following parameters: max-target-seqs 1, evalue 1e-10 [32]. Common DEGs are shown by a Venn diagram.

DEG Mapping on the Previously Identified Salt-Stress-Related QTL Intervals
A total of 17 salt-stress-related QTLs of rice were downloaded from the Gramene QTL database (http://archive.gramene.org/db/) [23]. All DEGs were mapped on these QTLs according to sequence and QTL location information.

Quantitative Real-Time PCR (qRT-PCR) Validation of DEGs
In this study, six DEGs were randomly selected for the verification of the DEG results. Primers of these genes were designed by Primer 5.0 in specific regions or 3'-UTR regions (the primers in Table S2). The qRT-PCR reaction ( ubiquitin fusion protein) was used as an internal control. The average threshold cycle (Ct) from three biological replicates was used to determine the fold change of gene expression by the 2 −∆∆CT method [33].

Identification of DEGs Involved in Salt Stress
Comparisons of gene expression in the shoots from the two genotypes (TNG67 and TCN1) (Table S3) revealed a total of 1241 DEGs in TCN1 and 904 DEGs in TGN67 at 3 h of salt treatment (S3) and two genotypes shared 586 DEGs ( Figure 1A). Following 24 h after salt stress (S24), 1519 and 1319 DEGs were identified in TCN1 and TGN67, respectively, and 738 DEGs were shared by two genotypes ( Figure 1A). At S3 and S24, two genotypes shared 377 DEGs ( Figure 1A). These genes showed similar expression patterns at S3 and S24 in these two genotypes ( Figures S2 and S3, Tables S4 and S5), suggesting that these genes were core salt-responsive genes in shoots between the two genotypes. After 24 h recovery after salt stress (ReS24), 395 DEGs were identified in TCN1, while 436 DEGs were identified in TGN67, and the two genotypes had 180 DEGs in common ( Figure 1A). In roots, 993 and 1448 DEGs were identified at S3 (714 DEGs in common), while 2493 and 2360 DEGs were identified at S24 (1605 DEGs in common) ( Figure 1B). At S3 and S24, a total of 488 DEGs were shared by the two genotypes ( Figure 1B). These genes displayed similar expression patterns at S3 and S4 in both genotypes ( Figures S2 and S4, Tables S4 and S6), indicating that these genes were core salt-responsive genes in the roots of these two genotypes. Interestingly, a total of 123 DEGs were found between core salt-responsive genes in shoots and core salt-responsive genes in roots ( Figure S1 and Table S4). As these genes also showed similar expression patterns in the shoots and roots of the two genotypes at S3 and S24, these genes play important roles during the salt stress process. At ReS24, in roots, 1037 DGEs were identified in TCN1 as well as 1057 DEGs in TGN67, and the two genotypes shared 353 DEGs ( Figure 1B). Under salt stress, the number of DEGs identified in roots was greater than that in shoots within these two genotypes (Figure 1), suggesting that the salt stress has a more broad effect on gene expression in roots than in shoots. We also observed that the number of DEGs at S24 was greater than that at S3 in shoots and roots in both genotypes. biological replicates was used to determine the fold change of gene expression by the 2 −ΔΔCT method [33].

Identification of DEGs Involved in Salt Stress
Comparisons of gene expression in the shoots from the two genotypes (TNG67 and TCN1) (Table S3) revealed a total of 1241 DEGs in TCN1 and 904 DEGs in TGN67 at 3 h of salt treatment (S3) and two genotypes shared 586 DEGs ( Figure 1A). Following 24 h after salt stress (S24), 1519 and 1319 DEGs were identified in TCN1 and TGN67, respectively, and 738 DEGs were shared by two genotypes (Figure 1A). At S3 and S24, two genotypes shared 377 DEGs ( Figure 1A). These genes showed similar expression patterns at S3 and S24 in these two genotypes ( Figures S2 and S3, Tables S4 and S5), suggesting that these genes were core salt-responsive genes in shoots between the two genotypes. After 24 h recovery after salt stress (ReS24), 395 DEGs were identified in TCN1, while 436 DEGs were identified in TGN67, and the two genotypes had 180 DEGs in common ( Figure 1A). In roots, 993 and 1448 DEGs were identified at S3 (714 DEGs in common), while 2493 and 2360 DEGs were identified at S24 (1605 DEGs in common) ( Figure 1B). At S3 and S24, a total of 488 DEGs were shared by the two genotypes ( Figure 1B). These genes displayed similar expression patterns at S3 and S4 in both genotypes ( Figures S2 and S4, Tables S4 and S6), indicating that these genes were core salt-responsive genes in the roots of these two genotypes. Interestingly, a total of 123 DEGs were found between core salt-responsive genes in shoots and core salt-responsive genes in roots ( Figure S1 and Table S4). As these genes also showed similar expression patterns in the shoots and roots of the two genotypes at S3 and S24, these genes play important roles during the salt stress process. At ReS24, in roots, 1037 DGEs were identified in TCN1 as well as 1057 DEGs in TGN67, and the two genotypes shared 353 DEGs ( Figure 1B). Under salt stress, the number of DEGs identified in roots was greater than that in shoots within these two genotypes (Figure 1), suggesting that the salt stress has a more broad effect on gene expression in roots than in shoots. We also observed that the number of DEGs at S24 was greater than that at S3 in shoots and roots in both genotypes. Figure 1. A Venn diagram of differentially expressed genes (DEGs) in shoots of TNG67 and TCN1 seedlings at different points (exposed to salt stress for 3 h or 24 h and allowed to recover for 24 h) (A); a Venn diagram of DEGs in roots of TNG67 and TCN1 seedlings at different points (B).
The heat map of all DEGs revealed that all DEGs showed similar expression patterns in these two genotypes (Figure 2), indicating that the main response mechanism of salt stress was the same in both genotypes. In addition, all DEGs can be grouped into two groups (I and II). Genes in Group I showed higher expression levels in roots than in shoots, while genes in Group II had higher expression levels in shoots than in roots ( Figure 2). The core salt-responsive genes in shoots and roots also showed different expression patterns ( Figures S3 and S4). These results suggest that the The heat map of all DEGs revealed that all DEGs showed similar expression patterns in these two genotypes (Figure 2), indicating that the main response mechanism of salt stress was the same in both genotypes. In addition, all DEGs can be grouped into two groups (I and II). Genes in Group I showed higher expression levels in roots than in shoots, while genes in Group II had higher expression levels in shoots than in roots ( Figure 2). The core salt-responsive genes in shoots and roots also showed different  Figures S3 and S4). These results suggest that the salt-response mechanisms were different between shoots and roots. Thus, the core salt-responsive genes in shoots and roots were annotated by a GO enrichment analysis, respectively. salt-response mechanisms were different between shoots and roots. Thus, the core salt-responsive genes in shoots and roots were annotated by a GO enrichment analysis, respectively. To verify the reliability of DEGs in this study, six salt-responsive DEGs with significantly upregulated expression levels under salt stress were randomly chosen for qRT-PCR ( Figure S5). The qRT-PCR result revealed that these six genes were significantly upregulated after salt stress.

GO Enrichment Analysis of Core Salt-Responsive Genes
To understand the core salt-responsive genes' functions in shoots or roots, a GO enrichment analysis of these genes was performed. A total of 30 and 25 terms were enriched in shoots and roots, respectively (Tables 1 and 2). Shoots and roots shared seven terms, including 'response to cadmium ion (GO:0046686)', 'response to water deprivation (GO:0009414)', 'response to cold (GO:0009409)', 'response to salt stress (GO:0009651)', 'hyperosmotic salinity response (GO:0042538)', 'response to abscisic acid (GO:0009737)', and 'response to wounding (GO:0009611)'. These terms are all related to abiotic stress, especially salt stress. As shown in Table 1 and Table 2, the remaining terms involve multiple biological processes. These results revealed that the salt-response mechanisms were different in shoots and roots. Shoots involve an abscisic acid response, the hydrogen peroxide catabolic process, the glucose catabolic process, the ethylene biosynthetic process, and the others, To verify the reliability of DEGs in this study, six salt-responsive DEGs with significantly upregulated expression levels under salt stress were randomly chosen for qRT-PCR ( Figure S5). The qRT-PCR result revealed that these six genes were significantly upregulated after salt stress.

GO Enrichment Analysis of Core Salt-Responsive Genes
To understand the core salt-responsive genes' functions in shoots or roots, a GO enrichment analysis of these genes was performed. A total of 30 and 25 terms were enriched in shoots and roots, respectively (Tables 1 and 2). Shoots and roots shared seven terms, including 'response to cadmium ion (GO:0046686)', 'response to water deprivation (GO:0009414)', 'response to cold (GO:0009409)', 'response to salt stress (GO:0009651)', 'hyperosmotic salinity response (GO:0042538)', 'response to abscisic acid (GO:0009737)', and 'response to wounding (GO:0009611)'. These terms are all related to abiotic stress, especially salt stress. As shown in Tables 1 and 2, the remaining terms involve multiple biological processes. These results revealed that the salt-response mechanisms were different in shoots and roots. Shoots involve an abscisic acid response, the hydrogen peroxide catabolic process, the glucose catabolic process, the ethylene biosynthetic process, and the others, while roots mainly involve nitrate transport, lipid transport, selenate transport, the salicylic-acid-mediated signaling pathway, and the others. Table 1. The biological process categories of Gene Ontology (GO) annotation of core salt-responsive genes in shoots under salt stress. The GO enrichment analysis of core salt-responsive genes was implemented by the GOseq R packages, and GO terms with p <0.001 were adopted in this study. Red terms represent common terms of shoots and roots.

MapMan Analysis of All Core Salt-Reponsive Genes of Shoots and Roots
For a better understanding of core salt-stress-response mechanisms, all core salt-responsive genes of shoots and roots (742 genes) were visualized by a biotic stress overview in MapMan 3.6.0. This result showed that the salt-stress-response mechanisms in rice were very complicated and involved multiple signal transductions, MAPK, TFs, and defense genes ( Figure 3). For example, lots of signaling receptor kinase genes, signaling calcium genes, signaling phosphoinositide genes, and signaling G-protein genes were significantly upregulated. Similarly, Auxins, ABA, Ethylene, and jasmonic acid (JA) signaling genes also showed significant upregulations. In addition, we found that ERF-, bZIP-, WRKY-, and MYB-TFs were significantly upregulated. Defense genes, including heat shock protein genes, secondary metabolite genes, and the others, also showed upregulation.
lignin biosynthetic process (GO:0009809); 3.23% 0.000737 Note ： A total of 465 core salt-responsive genes in roots were enriched on GO terms. Cluter_frequency = enriched gene numbers of each term/465.

MapMan Analysis of All Core Salt-Reponsive Genes of Shoots and Roots
For a better understanding of core salt-stress-response mechanisms, all core salt-responsive genes of shoots and roots (742 genes) were visualized by a biotic stress overview in MapMan 3.6.0. This result showed that the salt-stress-response mechanisms in rice were very complicated and involved multiple signal transductions, MAPK, TFs, and defense genes (Figure 3). For example, lots of signaling receptor kinase genes, signaling calcium genes, signaling phosphoinositide genes, and signaling G-protein genes were significantly upregulated. Similarly, Auxins, ABA, Ethylene, and jasmonic acid (JA) signaling genes also showed significant upregulations. In addition, we found that ERF-, bZIP-, WRKY-, and MYB-TFs were significantly upregulated. Defense genes, including heat shock protein genes, secondary metabolite genes, and the others, also showed upregulation.  Figure S6. ABA, abscisic acid; MAPK, mitogen-activated protein kinase; SA, salicylic acid; JA, jasmonic acid; HSPs, heat shock proteins.

Identification of TFs, TRs, and PKs Involved in Salt Stress
In order to identify more salt-responsive TFs, TRs, and PKs, we conducted an identification analysis of all core salt-responsive genes of shoots and roots (742 genes) using iTAK software. In this study, a total of 555 DEGs-TFs were identified (74.80%). These contained 47 Table S7). In addition, 23  Figure S6. ABA, abscisic acid; MAPK, mitogen-activated protein kinase; SA, salicylic acid; JA, jasmonic acid; HSPs, heat shock proteins.

Identification of TFs, TRs, and PKs Involved in Salt Stress
In order to identify more salt-responsive TFs, TRs, and PKs, we conducted an identification analysis of all core salt-responsive genes of shoots and roots (742 genes) using iTAK software. In this study, a total of 555 DEGs-TFs were identified (74.80%). These contained 47 families, of which the top seven families are NAC (66, 8.9%), MYB (60, 8.1%), bHLH (49, 6.6%), ERF (47, 6.3%), C2H2 (38, 5.1%), FHR1 (22, 3.0%), and HSF (22, 3.0%), respectively ( Figure 4A and Table S7). In addition, 23 PK subfamilies and four TR families were identified ( Figure 4B). Among them, the majority of PK subfamilies belonged to the receptor-like kinases (RLK/Pelle) family. In TRs, three members of HMG, two members of GNAT, two members of AUX/IAA, and one member of PHD were identified. PK subfamilies and four TR families were identified ( Figure 4B). Among them, the majority of PK subfamilies belonged to the receptor-like kinases (RLK/Pelle) family. In TRs, three members of HMG, two members of GNAT, two members of AUX/IAA, and one member of PHD were identified.

DEGs Involved in the Ca 2+ Signal Transduction Pathway
The Ca 2+ signal transduction pathway widely exists in eukaryotes [34]. Ca 2+ acts as a second messenger through binding to Ca 2+ sensors, causing a series of downstream reactions. At present, Ca 2+ sensors can be divided into three types in plants: calmodulins (CAM/CML), calcium-dependent protein kinases (CDPKs), and calmodulins B-like proteins (CBLs) [35]. In this study, 63 DEGs related to the Ca 2+ signaling pathway were identified, including genes encoding CMLs, CBLs, CDPKs, Ca 2+ -transporting ATPase, H + -ATPase, cation/Ca 2+ exchangers (CCXs), SOS1, HKT, AKT1, calcineurin B-like-interacting protein kinases (CIPKs), and ABI ( Figure 5A). Expression patterns revealed that three genes (LOC_Os10g25010, LOC_Os06g14030, and LOC_Os01g41510) were upregulated in shoots and roots of both genotypes at all time points and that three genes (LOC_Os12g12730, LOC_Os01g43410, and LOC_Os01g45990) were downregulated in shoots and roots of both genotypes at all time points. The remaining genes showed variable expression profiles under salt stress.

DEGs Involved in the Ca 2+ Signal Transduction Pathway
The Ca 2+ signal transduction pathway widely exists in eukaryotes [34]. Ca 2+ acts as a second messenger through binding to Ca 2+ sensors, causing a series of downstream reactions. At present, Ca 2+ sensors can be divided into three types in plants: calmodulins (CAM/CML), calcium-dependent protein kinases (CDPKs), and calmodulins B-like proteins (CBLs) [35]. In this study, 63 DEGs related to the Ca 2+ signaling pathway were identified, including genes encoding CMLs, CBLs, CDPKs, Ca 2+ -transporting ATPase, H + -ATPase, cation/Ca 2+ exchangers (CCXs), SOS1, HKT, AKT1, calcineurin B-like-interacting protein kinases (CIPKs), and ABI ( Figure 5A). Expression patterns revealed that three genes (LOC_Os10g25010, LOC_Os06g14030, and LOC_Os01g41510) were upregulated in shoots and roots of both genotypes at all time points and that three genes (LOC_Os12g12730, LOC_Os01g43410, and LOC_Os01g45990) were downregulated in shoots and roots of both genotypes at all time points. The remaining genes showed variable expression profiles under salt stress.

DEGs Involved in the ABA Signal Transduction Pathway
In total, 17 DEGs associated with the ABA signaling transduction pathway were identified: two PYL genes, six PP2C genes, five SnRK2 genes, and four ABF genes ( Figure 5B). The expression pattern results revealed that two PYL genes showed different expression patterns. LOC_Os03g18600 was downregulated under salt stress, while LOC_Os05g39580 was upregulated in TCN1 shoots under salt stress. Interestingly, these two genes both showed upregulation in TGN67 roots during the recovery process. Compared with PYL genes, the majority of PP2C genes showed upregulation in shoots and roots under salt stress. Similarly, the majority of SnRK2 and ABF genes were also upregulated in shoots and roots under salt stress.

DEGs Involved in the MAPK Cascade Pathway
In this study, 24 DEGs related to the MAPK cascade pathway were identified: three MPK genes, one MEK gene, and 20 MAPKKK genes ( Figure 5B). These genes showed unequal expression patterns. For example, LOC_Os01g50400 (MAPKKK) was upregulated in shoots and roots under salt stress. However, LOC_Os07g43900, LOC_Os01g54480, LOC_Os01g66860, and LOC_Os02g39560 showed downregulation at all tested points.

DEGs Involved in the ABA Signal Transduction Pathway
In total, 17 DEGs associated with the ABA signaling transduction pathway were identified: two PYL genes, six PP2C genes, five SnRK2 genes, and four ABF genes ( Figure 5B). The expression pattern results revealed that two PYL genes showed different expression patterns. LOC_Os03g18600 was downregulated under salt stress, while LOC_Os05g39580 was upregulated in TCN1 shoots under salt stress. Interestingly, these two genes both showed upregulation in TGN67 roots during the recovery process. Compared with PYL genes, the majority of PP2C genes showed upregulation in shoots and roots under salt stress. Similarly, the majority of SnRK2 and ABF genes were also upregulated in shoots and roots under salt stress.

DEGs Involved in the MAPK Cascade Pathway
In this study, 24 DEGs related to the MAPK cascade pathway were identified: three MPK genes, one MEK gene, and 20 MAPKKK genes ( Figure 5B). These genes showed unequal expression patterns. For example, LOC_Os01g50400 (MAPKKK) was upregulated in shoots and roots under salt stress. However, LOC_Os07g43900, LOC_Os01g54480, LOC_Os01g66860, and LOC_Os02g39560 showed downregulation at all tested points.

Key TFs and Functional Proteins Related to Salt Stress
We found that NAC (NAM, ATAF, and CUC) family genes formed the largest salt-responsive TF family of this study, and the majority of NAC family genes showed upregulation under salt stress (Figures 3 and 6). Of these NAC genes, 14 have been reported in previous studies ( Figure 6A) and their functions were associated with plant senescence delay [36], abiotic stress responses [37], cellulose synthesis [38], and others [39,40]. These results imply that the NAC family is a multifunctional family and also plays vital roles in regulating plant salt stress tolerance. Thus, NAC family genes may be good candidate genes in salt stress tolerance for genetic engineering breeding.
Previous reports have indicated that many functional proteins, including AQPs (aquaporins), HSPs (heat shock proteins), LEA (late embryogenesis abundant) proteins, F-box proteins, transporter proteins, and other functional proteins are differentially expressed under salt stress and play key roles in regulating plant salt tolerance [24,41]. In this study, seven HSPs and eight AQPs were identified in DEGs ( Figure 6B). In addition, 22 of the previously identified functional proteins also were found in this study's DEGs.

Key TFs and Functional Proteins Related to Salt Stress
We found that NAC (NAM, ATAF, and CUC) family genes formed the largest salt-responsive TF family of this study, and the majority of NAC family genes showed upregulation under salt stress (Figures 3 and 6). Of these NAC genes, 14 have been reported in previous studies ( Figure 6A) and their functions were associated with plant senescence delay [36], abiotic stress responses [37], cellulose synthesis [38], and others [39,40]. These results imply that the NAC family is a multifunctional family and also plays vital roles in regulating plant salt stress tolerance. Thus, NAC family genes may be good candidate genes in salt stress tolerance for genetic engineering breeding.
Previous reports have indicated that many functional proteins, including AQPs (aquaporins), HSPs (heat shock proteins), LEA (late embryogenesis abundant) proteins, F-box proteins, transporter proteins, and other functional proteins are differentially expressed under salt stress and play key roles in regulating plant salt tolerance [24,41]. In this study, seven HSPs and eight AQPs were identified in DEGs ( Figure 6B). In addition, 22 of the previously identified functional proteins also were found in this study's DEGs.

TRs and PKs Involved in Salt Stress
TRs and PKs play vital roles in rice tolerance to abiotic stresses [24]. However, they have been rarely mentioned in previous salt stress studies. In this study, many TRs and PKs were differentially expressed under salt stress ( Figure 6C). For instance, one-third of the PK genes showed upregulation under salt stress. In TRs, AUX/IAA and GNAT genes showed upregulation under salt stress, while HMG and PHD genes were downregulated under salt stress.

Key Genes Associated with Higher Salt Stress Tolerance in TNG67
To investigate the possible reason for different salt stress tolerances between TCN1 and TNG67, metabolic pathways were visualized by a metabolism overview in MapMan 3.6.0. For shoots, these results revealed a greater enhancement of gene expression associated with light reactions in TNG67 than in TCN1 under salt stress (Figure 7). For roots, there was a greater enhancement of gene expression associated with starch, sucrose, and fermentation in TNG67 than in TCN1 under salt stress ( Figure 8). Interestingly, we observed that Calvin cycle genes had higher expression levels in TNG67 than in TCN1 during the recovery process ( Figure 8). Thus, we speculated that TNG67 shifted to a fermentation pathway to produce energy for growth under salt stress and produced lots of starch and sucrose for coping with osmotic stress. In addition, TNG67 may enhance the Calvin cycle to repair a damaged photosystem during the recovery stage.

TRs and PKs Involved in Salt Stress
TRs and PKs play vital roles in rice tolerance to abiotic stresses [24]. However, they have been rarely mentioned in previous salt stress studies. In this study, many TRs and PKs were differentially expressed under salt stress ( Figure 6C). For instance, one-third of the PK genes showed upregulation under salt stress. In TRs, AUX/IAA and GNAT genes showed upregulation under salt stress, while HMG and PHD genes were downregulated under salt stress.

Key Genes Associated with Higher Salt Stress Tolerance in TNG67
To investigate the possible reason for different salt stress tolerances between TCN1 and TNG67, metabolic pathways were visualized by a metabolism overview in MapMan 3.6.0. For shoots, these results revealed a greater enhancement of gene expression associated with light reactions in TNG67 than in TCN1 under salt stress (Figure 7). For roots, there was a greater enhancement of gene expression associated with starch, sucrose, and fermentation in TNG67 than in TCN1 under salt stress ( Figure 8). Interestingly, we observed that Calvin cycle genes had higher expression levels in TNG67 than in TCN1 during the recovery process ( Figure 8). Thus, we speculated that TNG67 shifted to a fermentation pathway to produce energy for growth under salt stress and produced lots of starch and sucrose for coping with osmotic stress. In addition, TNG67 may enhance the Calvin cycle to repair a damaged photosystem during the recovery stage.  Figure 7 is provided in Figure S7.

The Ca2+ Signaling Pathway, ABA Signaling Transduction Pathway, and MAPK Cascade Pathway Play Important Roles under Salt Stress Conditions
Previous genetic and biochemical studies have revealed that many genes participate in salt resistance via hormone and Ca 2+ signal pathways, TFs, ion metabolism, ion transfer, nitrogen metabolism, and secondary metabolism [3,8,24]. In this study, these genes were also identified. A total of 63, 17, and 24 DEGs were identified in Ca 2+ signaling pathway, the ABA signaling transduction pathway, and the MAPK cascade pathway, respectively. These results support previous findings that these pathways can form a close-knit signaling network and play vital roles in  Previous genetic and biochemical studies have revealed that many genes participate in salt resistance via hormone and Ca 2+ signal pathways, TFs, ion metabolism, ion transfer, nitrogen metabolism, and secondary metabolism [3,8,24]. In this study, these genes were also identified. A total of 63, 17, and 24 DEGs were identified in Ca 2+ signaling pathway, the ABA signaling transduction pathway, and the MAPK cascade pathway, respectively. These results support previous findings that these pathways can form a close-knit signaling network and play vital roles in plant salt stress tolerance by interacting with each other or starting up downstream factors, such as TFs, TRs, and functional protein genes (Figure 3) [41][42][43].

Many TFs and PKs May Be Good Candidate Genes for Plant Salt Stress Tolerance Breeding
Previous studies have reported that NAC, MYB, bHLH, and AP2/ERF are associated with salt stress tolerance [44]. Of these TFs, some TFs are associated with stress-signal pathways. For example, Zhang et al. (2017) reported that several bHLH-TFs are associated with the ABA signal pathway [44]. Jakoby et al. (2002) reported that bZIP-TFs are downstream factors of the ABA signal pathway and regulate the pivotal cell process in the response to salt stress [45]. Zhang et al. (2012) reported that HD-Zip-TFs affect ABA biosynthesis and regulate rice salt tolerance via the ABA signal pathway [46]. Similarly, we also found that TFs participate in the plant stress response via the ABA signal pathway. In this study, we found that LOC_Os02g17500 (encoding a GRAS TF) was associated with SnRK (LOC_Os03g27280, key genes for the ABA signal pathway). However, beyond that, TFs are also regulators of salt stress. For example, Jiang et al. (2017) reported that WRKY-TFs are important regulators for salt stress tolerance [47]. OsMYB2 is an important regulator for salt stress in rice [48]. Some previously reported TFs (associated with abiotic stress) were found to be DEG-TFs in this study. These genes may also play vital roles in the salt stress tolerance of rice, such as LOC_Os03g48780, LOC_Os05g49730, LOC_Os04g56430, LOC_Os03g48750, and LOC_Os09g31031 (Table S7). For instance, several reported TFs in salt stress were present in the DEG-TFs identified in this study, such as Oshox22 (LOC_Os04g45810, HD-ZIP) [46] and OsHsfB2b (LOC_Os08g43334, HSP) [49]. In this study, many TFs were identified and can be important candidate genes for salt stress breeding. Among them, two HD-ZIP genes (LOC_Os06g46740 and LOC_Os02g43330) showed high upregulation (LogFC >2.0, Table S7). Similarly, two MYB genes (LOC_Os01g54030 and LOC_Os08g39730) also were highly upregulated (LogFC >3.0, Table S7). These genes can be important for rice salt stress tolerance.
RLKs (receptor-like protein kinases), one of the largest gene families in plants, play vital roles in the regulation of plant developmental processes, signaling networks, and disease resistance [50]. Previous studies have reported that many RLKs have been proved to be involved in abiotic stress responses, including the ABA response, Ca 2+ signaling, and antioxidant defense [50]. Vaid et al. (2015) reported that Pisum sativum LecRLKs (PsLecRLKs) were upregulated under salt stress and PsLecRLK-overexpressing plants had a greater tolerance to salt stress than wild-type plants due to ROS-scavenging enzymes, reducing ROS accumulation and leading to lower membrane damage [51]. Li and Sun (2014) found that SIT1, a LecRLK gene mainly expressed in root epidermal cells, mediated rice salt sensitivity [52]. In this study, many TR genes were upregulated under salt stress. RLKs may be good candidate genes for rice salt stress tolerance breeding.

Common DEGs between This Study and Previous Studies Contain Important Salt-Responsive Genes
Common DEGs from multiple salt stress studies can help us to identify core salt-responsive genes. In this study, we conducted a comparative analysis of this study with two previous salt stress studies. This result showed that a total of 1738 and 1080 common DEGs were found with the Wang et al. study [8] and the Zhou et al. study [23], respectively ( Figure 9A,B, Tables S8 and S9). Next, 450 common DEGs were found in these three studies ( Figure 9C). Based on MSU7.0 description and expression data, we noticed that several genes may be good candidate genes for future salt stress tolerance breeding. For example, LOC_Os01g11730 (a GDSL-like lipase/acylhydrolase gene), LOC_Os05g31670 (encoding AWPM-19-like membrane family protein), LOC_Os01g21420 (encoding pre-mRNA-splicing factor, SF2), and LOC_Os05g31020 (encoding eukaryotic peptide chain release factor subunit 1-1) showed significant upregulation under salt stress (Table S10). Moreover, we found that 100 salt-responsive genes were reported in previous studies (Table S11). Of these genes, 34 were identified as DEGs in this study (Table S12). This result revealed that many unreported salt-responsive genes were found in this study.

Multiple DEGs Mapped to the Previously Identified Salt-Stress-Related QTL Intervals and Could Be Good Candidate Genes for Rice Salt Stress Tolerance Breeding
The previously identified salt-stress-related QTLs can help us to target candidate genes. A total of 17 QTLs were obtained from the Gramene QTL database [23]. The DEGs of this study were mapped to 12 QTLs (Table S13). Among them, AQGR001, AQEM001, AQEM008, and AQEM007 had the greatest number of co-localized DEGs: 122, 93, 93, and 64, respectively ( Figure 10). A previous study reported that AQEM002 and AQGR001 are the most important QTLs related to rice salt stress tolerance at the seedling stage [12,53]. In this study, 35 and 122 DEGs, respectively, were co-localized on these two QTLs. These DEGs could be good candidate genes for rice salt stress tolerance breeding, such as LOC_Os01g64360 (one MYB gene, in AQGR001), LOC_Os01g64790 (one AP2/ERF gene, in AQGR001), LOC_Os01g62760 (one PP2C gene, in AQGR001), LOC_Os06g39040 (one dehydrogenase/reductase SDR family gene). Interestingly, only one gene (LOC_Os03g41064) was mapped on AQGR002, which encodes one natural-resistance-associated macrophage protein. Thus, LOC_Os03g41064 may play a vital role under salt stress.

Multiple DEGs Mapped to the Previously Identified Salt-Stress-Related QTL Intervals and Could Be Good Candidate Genes for Rice Salt Stress Tolerance Breeding
The previously identified salt-stress-related QTLs can help us to target candidate genes. A total of 17 QTLs were obtained from the Gramene QTL database [23]. The DEGs of this study were mapped to 12 QTLs (Table S13). Among them, AQGR001, AQEM001, AQEM008, and AQEM007 had the greatest number of co-localized DEGs: 122, 93, 93, and 64, respectively ( Figure 10). A previous study reported that AQEM002 and AQGR001 are the most important QTLs related to rice salt stress tolerance at the seedling stage [12,53]. In this study, 35 and 122 DEGs, respectively, were co-localized on these two QTLs. These DEGs could be good candidate genes for rice salt stress tolerance breeding, such as LOC_Os01g64360 (one MYB gene, in AQGR001), LOC_Os01g64790 (one AP2/ERF gene, in AQGR001), LOC_Os01g62760 (one PP2C gene, in AQGR001), LOC_Os06g39040 (one dehydrogenase/reductase SDR family gene). Interestingly, only one gene (LOC_Os03g41064) was mapped on AQGR002, which encodes one natural-resistance-associated macrophage protein. Thus, LOC_Os03g41064 may play a vital role under salt stress.

Multiple DEGs Mapped to the Previously Identified Salt-Stress-Related QTL Intervals and Could Be Good Candidate Genes for Rice Salt Stress Tolerance Breeding
The previously identified salt-stress-related QTLs can help us to target candidate genes. A total of 17 QTLs were obtained from the Gramene QTL database [23]. The DEGs of this study were mapped to 12 QTLs (Table S13). Among them, AQGR001, AQEM001, AQEM008, and AQEM007 had the greatest number of co-localized DEGs: 122, 93, 93, and 64, respectively ( Figure 10). A previous study reported that AQEM002 and AQGR001 are the most important QTLs related to rice salt stress tolerance at the seedling stage [12,53]. In this study, 35 and 122 DEGs, respectively, were co-localized on these two QTLs. These DEGs could be good candidate genes for rice salt stress tolerance breeding, such as LOC_Os01g64360 (one MYB gene, in AQGR001), LOC_Os01g64790 (one AP2/ERF gene, in AQGR001), LOC_Os01g62760 (one PP2C gene, in AQGR001), LOC_Os06g39040 (one dehydrogenase/reductase SDR family gene). Interestingly, only one gene (LOC_Os03g41064) was mapped on AQGR002, which encodes one natural-resistance-associated macrophage protein. Thus, LOC_Os03g41064 may play a vital role under salt stress. Figure 10. Co-localization of DEGs onto the previously detected quantitative trait loci (QTLs) responsible for salt treatment in rice. Red words denote previously detected QTLs, green numbers denote the DEGs of this study on previously detected QTLs.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Venn diagram of salt-responsive DEGs in shoots and in roots. Figure S2: Heat map of salt-responsive genes (123 common DEGs in shoots and roots in both genotypes). Figure S3: Heat map of salt-responsive genes in shoots (254 common DEGs in shoots in both genotypes). Figure S4: Heat map of salt-responsive genes in roots Figure 10. Co-localization of DEGs onto the previously detected quantitative trait loci (QTLs) responsible for salt treatment in rice. Red words denote previously detected QTLs, green numbers denote the DEGs of this study on previously detected QTLs.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2223-7747/8/3/64/s1, Figure S1: Venn diagram of salt-responsive DEGs in shoots and in roots. Figure S2: Heat map of salt-responsive genes (123 common DEGs in shoots and roots in both genotypes). Figure S3: Heat map of salt-responsive genes in shoots (254 common DEGs in shoots in both genotypes). Figure S4: Heat map of salt-responsive genes in roots (365 common DEGs in shoots in both genotypes). Figure S5: qRT-PCR results of DEGs in 'Nipponbare'. Figure S6: MapMan biotic stress overview of core salt-responsive genes (742 DEGs) in shoots after 3 h, 24 h, and recovery 24 h in shoots, as well as 3 h, 24 h, and recovery 24 h in roots. Table S1: Detailed information of GSE76613. Table S2: Primers of qRT-PCR for DEG validation. Table S3: All expression datasets of DEGs under salt and recovery conditions in shoots and roots of the two genotypes. Table S4: Expression datasets of Figure  S2 DEGs. Table S5: Expression datasets of Figure S3 DEGs. Table S6: Expression datasets of Figure S4 DEGs. Table S7: Expression datasets of TF genes (these genes are identified in 742 core salt-responsive genes).  Table S11: Details of previously reported salt responsive genes (100 genes). Table S12: The common genes between this study and previously reported salt responsive genes. Table S13. Co-localization of DEGs onto the previously detected QTLs responsible for salt treatment in rice.
Author Contributions: W.K. performed the experiments, analyzed the data, prepared figures and tables, and wrote the manuscript. Y.L. conceived and designed the experiments as well as modified the paper. H.Z., Z.G., X.F., X.D., and T.S. analyzed parts of the data and prepared parts of figures and tables. All authors read and approved the final manuscript.