Sodium and potassium in roots and shoots
The sodium (Na+) and potassium (K+) contents and Na+/K+ ratio of roots and shoots were evaluated by ANOVA and comparison of means (Fig. S1, Table S2). Salinity treatment and increasing its duration from 6 to 54 hours, led to increased Na+ and decreased K+ in both organs of both genotypes. Further, their Na+/K+ ratio increased significantly under salinity treatment. While there were no significant differences between the genotypes under the control condition, we observed lower Na+, higher K+ and lower Na+/K+ ratio in both organs of CSR28 compared to those of IR28 under salinity stress. In particular, salinity stress in the roots resulted in a reduction of Na+ concentration in CSR28 by 20.6% and 40.4% at 6-hour and 54-hour timepoints, respectively, as compared to IR28, while the reductions in the shoots were 23.8% and 46%, respectively. Conversely, the absorption of K+ in the roots of CSR28 increased 10.4% and 34.3% at 6 hours and 54 hours after salt treatment as compared to IR28, while the increases in shoots were 17.2% and 41.4%, respectively. The Na+/K+ ratio in both organs showed significant difference between the two genotypes only under high salinity, as the ratio of salt-sensitive (IR28) was significantly higher than that of the tolerant one (CSR28).
Global view of gene expression
In this study, single-end RNA-Seq reads were generated from 48 samples, including two organs (root and shoot), two contrasting genotypes (salt-tolerant CSR28 and salt-sensitive IR28), two treatments (control and 150 mM salinity), two sampling times (6 hours and 54 hours) and three biological replicates. A total of 1.231 billion reads were obtained from all the samples (range 13.51 to 56.06 million) with a total length of 184 Gb. Adapters and low-quality reads (Q ≤ 20) were discarded by Trimmomatic, resulting in 1.204 billion (97.8%) high-quality reads. Out of the clean reads, 1.023 billion (84.9%) with an average of 21.31 million per sample were mapped on the rice reference genome using TopHat2 (Table S3). The normalized expression values (RPKM) were used to hierarchical clustering of the samples based on Euclidean distances and Average linkage method (Fig. 2a). The results showed that the samples are clustered based on two major groups of the roots and shoots. In the shoots, control and salinity samples separated into two subclusters. In control condition, the samples were separated according to the genotypes, and in each genotype, the most similarity was observed between the timepoints. However, in salinity condition, the samples separated by timepoints and there was the most similarity between the genotypes. In the roots, the expression profiles exhibited a more complex pattern, particularly in stress conditions. Unlike control condition and 6 hours salinity treatment, there was a remarkable difference between the sensitive IR28 and tolerant CSR28 genotypes at 54-hour timepoint. Figure 2b displays the expression patterns of the samples for the highly expressed genes in response to salinity stress. Most correlation between the genotypes was observed in the shoots and at 6-hour timepoint in the roots, whereas there was a pronounced difference between the roots of IR28 and CSR28 at 54-hour timepoint. As a result, the gene expression changes of the roots after long-term exposure of 150 mM NaCl treatment were the key factor to differentiating the genotypes in terms of salinity tolerance.
Identification of differentially expressed genes (DEGs) and transcription factors (TFs)
Out of 16 unique experimental samples, Cuffdiff identified a total of 15483 DEGs by 32 pairwise comparisons. The highest (7569 genes) and lowest (29 genes) number belonged to the comparison of the roots and shoots of IR28 at 6-hour timepoint under control condition (R-CT-6h-IR28 vs S-CT-6h-IR28) and the comparison of two timepoints in the shoots of CSR28 under control condition (S-CT-6h-CSR28 vs S-CT-54h-CSR28), respectively (Fig. S2). Figure 3 displays the comparisons between organs, treatments, genotypes and timepoints. Under salinity stress at 6-hour timepoint, two genotypes induced almost the same number of DEGs between the roots and shoots (Fig. 3a). However, there was a huge difference after 54 hours of treatment. CSR28 had interestingly 1283 genes overexpressed in the roots, while only 97 genes were induced in shoots. This result indicates that the higher salinity tolerance of CSR28 is predominantly due to altered gene expression in the roots. In response to salinity stress, a greater number of differentially expressed genes (DEGs) were observed in the shoots compared to the roots (Fig. 3b). Specifically, the roots of CSR28 exhibited a higher number of DEGs (1540 up-regulated genes and 848 down-regulated genes) as the duration of salinity exposure increased. Figure 3c illustrates the DEGs between CSR28 and IR28. Moreover, under 54 hours of salinity treatment, the roots of CSR28 displayed a significantly higher number of up-regulated genes (1373) compared to IR28 (521). Conversely, in both timepoints, the shoots of IR28 showed a greater number of up-regulated genes than CSR28. Furthermore, the number of genes that exhibited differential expression between the 6-hour and 54-hour timepoints is depicted in Fig. 3d. In both the roots and shoots, the number of DEGs increased in response to salinity stress. Specifically, in the roots, a higher number of genes were induced at 54-hour timepoint for CSR28, while IR28 showed more DEGs at 6-hour timepoint. In contrast, the shoots of IR28 had a significantly higher number of DEGs at 54-hour timepoint, and both timepoints in CSR28 exhibited almost the same level of gene induction. As a general conclusion, increasing the duration of salinity exposure from 6 to 54 hours was associated with an elevated number of expressed genes in the roots of salt-tolerant genotype and the shoots of salt-sensitive genotype. Venn diagram analysis was employed to identify specific and common genes between control condition and salinity stress for the genes which were differentially expressed between CSR28 and IR28. At 6-hour timepoint, 525 and 635 genes were specifically identified under salinity stress in the roots and shoots, respectively (Fig. 3e). We also identified 1472 and 606 genes at 54-hour timepoint, which were specifically expressed under salinity stress in the roots and shoots, respectively (Fig. 3f). We used these “salt-specific genes” for further analyses.
Using the PlantTFDB database, we identified the genes encoding transcription factors that responded to salinity by comparing control and salinity-treated samples. The most abundant transcription factors in response to salinity in different conditions belonged to the WRKY, MYB, bHLH, HB, and AP2-EREBP families (Fig. S3a). In the roots of the salt-tolerant CSR28, transcription factors belonging to the AP2-EREBP and MYB families were the most important at both timepoints, while in the salt-sensitive IR28, the WRKY family was the most abundant. In the shoots and in response to salinity at 6-hour timepoint, MYB played a greater role in both genotypes, and at 54-hour timepoint, the AP2-EREBP and MYB families were the most abundant in tolerant and sensitive genotypes, respectively. The expression level of genes belonging to the MYB family is shown as a heatmap in Fig. S3b. In response to salinity, the expression of various members of the MYB family increased, with the greatest increase observed in the roots of the salt-tolerant genotype at 54-hour timepoint. While some MYB members, like Os12g0564100 and Os04g0508500, increased in expression across all conditions in response to salinity, others like Os05g0350900 and Os07g0634900 decreased. Additionally, certain genes exhibited specific expression patterns in different conditions. For instance, Os11g0128500 was up-regulated in response to salinity specifically in the salt-tolerant genotype’ roots at 54-hour timepoint.
Functional classification of salt-specific genes by Gene Ontology (GO) and KEGG pathway enrichment analyses
In order to functional classification of the salt-specific genes at the three levels of molecular function (MF), cellular component (CC) and biological process (BP), the agriGO database was employed. Out of the 525 and 1472 salt-specific genes in the roots, 26 and 21 GO terms were identified at 6-hour and 54-hour timepoints, respectively (Fig. 4a, Table S4). At the molecular function (MF) level, the terms "binding" (GO:0005488) and "catalytic activity" (GO:0003824) were enriched with the highest number of genes in both timepoints. However, a larger number of genes were involved at 54-hour timepoint. The GO terms associated with important molecular responses of plants to environmental stresses, such as, "transcription regulator activity" (GO:0030528), "transporter activity" (GO:0005215), and "antioxidant activity" (GO:0016209), exhibited enrichment with more genes at 54-hour timepoint.
At the cellular component (CC) level, the terms "cell" (GO:0005623), "cell part" (GO:0044464), and "organelle" (GO:0043226) had the highest number of genes in both timepoints. The terms "organelle part" (GO:0044422) and "macromolecular complex" (GO:0032991) were significantly enriched with more genes at 54-hour timepoint, while the term "extracellular region" (GO:0005576) was specifically enriched at 6-hour timepoint.
In terms of biological process (BP), the terms "metabolic process" (GO:0008152) and "cellular process" (GO:0009987) were enriched with the largest number of genes in both timepoints. The functional term "response to stimulus" (GO:0050896) and the regulatory term "biological regulation" (GO:0065007) were significantly more enriched at 54-hour timepoint. On the other hand, terms related to development and reproduction, such as, "developmental process" (GO:0032502) and "reproduction" (GO:0000003), were specifically involved at 6-hour timepoint. Overall, it can be concluded that functional terms related to salinity stress in the roots at 54-hour timepoint played essential roles in differentiating the genotypes for salinity tolerance.
In the shoots, out of the 635 and 606 salt-specific genes, 22 and 27 GO terms were identified at 6-hour and 54-hour timepoints, respectively (Fig. 4a, Table S4). Similar to the roots, functional groups such as, "catalytic activity," "binding," "cell," "organelle," "metabolic process," and "cellular process" were more enriched than others in the shoots. However, unlike the roots, the GO terms in the shoots were more enriched at 6-hour timepoint than the 54-hour timepoint. Terms associated with stress response, such as, "response to stimulus" and "biological regulation," were more enriched at 6-hour timepoint, whereas terms like "antioxidant activity" at the MF level, "extracellular region" at the CC level, and "reproduction" and "multi-organism process" (GO:0051704) at the BP level were specifically involved at 54-hour timepoint. Fig. S4 illustrates a portion of the co-functional network of BP-related GO terms using salt-specific genes in the roots at 54-hour timepoint. This network revealed significant enrichment of functional groups such as, "response to stimulus" and "ion transport" (GO:0006811). Ion transporters are crucial proteins for sodium detoxification and maintenance of ion balance under salinity stress in roots.
In order to identify the pathways of molecular and metabolic processes, a pathway analysis of salt-specific genes was conducted using the KEGG database (Fig. 4b). The pathway of biosynthesis of secondary metabolites (ko01110), was found to be highly enriched in both organs and timepoints. We observed different functional pathways involved in the timepoints of the roots. Metabolic pathways (ko01100), phenylpropanoid biosynthesis (ko00940), amino sugar and nucleotide sugar metabolism (ko00520), fatty acid elongation (ko00062), and others were identified at 6-hour timepoint, while the pathways of ribosome (ko03010), porphyrin metabolism (ko00860), ABC transporters (ko02010), Carotenoid biosynthesis (ko00906), and others were involved at 54-hour timepoint. By contrast, the shoots represented almost same KEGG pathways in the both timepoints. In addition to biosynthesis of secondary metabolites, MAPK signaling pathway – plant (ko04016) and amino sugar and nucleotide sugar metabolism (ko00520) were commonly enriched at 6-hour and 54-hour timepoints, while plant hormone signal transduction (ko04075), plant-pathogen interaction (ko04626) and biotin metabolism (ko00780) specifically observed at 6-hour timepoint and metabolic pathways, phenylpropanoid biosynthesis and phenylalanine metabolism (ko00360) represented specific enrichment at 54-hour timepoint.
PPI-based network analysis revealed hub genes in salt-specific genes
Protein-protein interaction (PPI) network analysis was employed to explore the associations between the salt-responsive genes in the roots at 54-hour timepoint. Among 1130 up-regulated genes in CSR28, 50 hub genes were identified on the basis of Cytohubba plugin in Cytoscape, including OS02T0822600-01, OsJ_13095, OsJ_36392, RPS9, OS03T0122200-01, OsJ_22158, OS03T0265400-01, OsJ_03455, RPL5, OsJ_06142 and so on (Fig. 5a). Among 342 up-regulated genes in IR28, 25 hub genes with lower interconnected scores were identified, including OsJ_00707, CYCB2-2, OsJ_32048, OsJ_09263, KIN10A, OS01T0931200-01 and so on (Fig. 5b).
Functional classification of hub genes
The GO enrichment analysis of the hub genes using STRING revealed that GO terms such as, “translation”, “cellular nitrogen compound biosynthetic process” and “gene expression” in BP level, “rRNA binding” and “structural constituent of ribosome” in MF level and “chloroplast” and “ribosome” in CC level were significantly involved in CSR28 (Fig. S5, Table S5). In IR28, the GO terms such as, “cell division” and “microtubule-based process” in BP level, “microtubule binding” and “cytoskeletal protein binding” in MF level and “microtubule cytoskeleton”, “supramolecular complex” and “spindle” in CC level were involved with the highest level of significance (Fig. S5, Table S6). KEGG pathway analysis showed that “ribosome” was the only pathway enriched for the hub genes of CSR28. Out of 50 hub genes in CSR28, 21 genes were significantly involved in KEGG pathway of “ribosome” (Table S7). No KEGG pathway was involved for the hub genes of IR28’ roots at 54-hour timepoint.
Validation by qRT-PCR
The expression values of four selected genes in response to salinity stress were shown in Fig. S6a, indicating a similar pattern between RNA-seq and qRT-PCR approaches. The sequencing data was validated by a strong correlation (R2 = 0.85) with the expression values of qRT-PCR analysis (Fig. S6b).