Phenotypic variation between CH891 and 02428 under Cd stress
Following exposure of the seedlings of rice strains CH891 and 02428 to Cd, large phenotypic variations were observed. The seedling lengths of CH891 and 02428 were shorter under Cd stress compared with the control condition (Fig. 1A-D). Furthermore, the seedling length of CH891 was significantly longer than 02428 in the control condition, while there were no significant differences between the two genotypes under Cd stress (Fig. 1E, F). This result indicated that CH891 was more sensitive to Cd than 02428. The experiment was carried out continuously for 7 days and we found that the two highest seedling length inhibition rates were on the 4th and 6th days (Fig. 1G). In addition, LOC_Os07 g15460 (OsNRAMP1), which was previously reported to be involved in Cd accumulation in rice, was detected on the 3rd and 5th day in CH891 and 02428 with a high expression level and it had significant differences between under stress and the control (Supplementary Table S1). This result indicated that this gene had higher expression on the 3rd and 5th days than on the other days. Therefore, the seedlings were harvested on the 3rd and 5th days, and their RNA was sequenced.
RNA sequencing of the seedling transcriptome of the two genotypes
A total of 1,070 million reads with an average of 46 million reads per sample were generated. The total number of valid reads was 1,057 million, and the average number of valid reads per sample was 46 million. The ratio of Q20 for each was above 99%, and the Q30 base percentage was above 95% (Supplementary Table S2). Therefore, the quality of the data was very high and met the requirements for further analysis.
To confirm the accuracy and reproducibility of the RNA-Seq results, 16 genes were selected for qRT-PCR using specific primers (Supplementary Table S3). The validation results for the 16 genes are shown in Fig. 5A, B. The qRT-PCR results were all consistent with the RNA-Seq data (Supplementary Table S3). To conclude, our transcriptome sequencing results were credible.
Identification of differentially expressed genes (DEGs)
By comparing samples of the same rice cultivar under different conditions (control and stress) and different rice cultivars (CH891 and 02428) under the same conditions at different stages, we assigned four comparison groups at the two stages, eight comparison groups in total. Then, DEGs in the different comparisons were selected by restricting pval≤0.05. These results indicated that not only the cultivar but also the treatments and time points affected the gene expression level. The number of upregulated and downregulated genes in each comparison group is shown in Figure 2A.
Functional enrichment analysis was performed for all of these DEGs (Supplementary Fig. S1). KEGG enrichment analysis showed that the plant-pathogen interaction, phenylpropanoid biosynthesis, tryptophan metabolism, flavonoid biosynthesis and diterpenoid biosynthesis pathways were enriched in Cd3_CH891 vs. Cd3_02428. Plant-pathogen interaction, tryptophan metabolism, phenylpropanoid biosynthesis, flavonoid biosynthesis, diterpenoid biosynthesis and amino sugar and nucleotide sugar metabolism pathways were enriched in CK3_CH891 vs. CK3_02428. Starch and sucrose metabolism, phenylpropanoid biosynthesis, glycolysis, carbon metabolism, alpha-linolenic acid metabolism, fatty acid degradation and brassinosteroid biosynthesis pathways were enriched in Cd3_CH891 vs. CK3_CH891. Plant-pathogen interaction and plant hormone signal transduction pathways were enriched in Cd3_02428 vs. CK3_02428. On the 5th day, phenylpropanoid biosynthesis, glutathione metabolism, flavonoid biosynthesis, brassinosteroid biosynthesis, amino sugar and nucleotide sugar metabolism and ABC transporter pathways were enriched in Cd5_CH891 vs. Cd5_02428. Tryptophan metabolism, sesquiterpenoid and triterpenoid biosynthesis and phenylpropanoid biosynthesis pathways were enriched in CK5_CH891 vs. CK5_02428. Phenylpropanoid biosynthesis, glutathione metabolism, diterpenoid biosynthesis, brassinosteroid biosynthesis, alpha-linolenic acid metabolism and ABC transporter pathways were enriched in Cd5_CH891 vs. CK5_CH891. Protein processing in the endoplasmic reticulum, porphyrin and chlorophyll metabolism, photosynthesis, phenylpropanoid biosynthesis and circadian rhythm-plant pathways were enriched in Cd5_02428 vs. CK5_02428.
All of these results revealed that there were more upregulated DEGs under Cd stress than under control conditions in CH891, while there were more downregulated DEGs under Cd stress than under control conditions in 02428 on either the 3rd or 5th day. Moreover, there was a higher number of DEGs obtained for Cd stress in CH891 (Cd-sensitive cultivar) than in 02428 (Cd-resistant cultivar).
Classification of differentially expressed genes (DEGs)
A total of 7204 unique DEGs were detected in the four groups on the 3rd day, and 6670 unique DEGs were obtained on the 5th day. These DEGs were divided into 15 subgroups (Fig. 2B, C). Excluding these DEGs of the groups (CK3_CH891 vs. CK3_02428 on the third day and the groups containing CK5_CH891 vs. CK5_02428 on the fifth day) irrelevant to Cd stress, the DEGs of the remaining seven subgroups were divided into three categories: genes from the sensitive variety with Cd-responsive (SCR), genes from the resistant variety with Cd-responsive (RCR), and common Cd-responsive (CCR) DEGs. There were 849, 676 and 770 DEGs in SCR, RCR and CCR on the 3rd day, respectively (Supplementary Table S4). In contrast, 959, 592 and 1516 DEGs were in the SCR, RCR and CCR, respectively, on the 5th day (Supplementary Table S4).
In addition, we selected important Cd-responsive (ICR) genes in the three categories with |log2 (fold change)|≥2 and FPKM≥2 in at least one group to select representative DEGs. On the 3rd day, 554 ICR genes were screened, including 194 for SCR, 154 for RCR and 206 for CCR (Supplementary Table S4). A total of 941 ICR genes were screened on the 5th day, of which 389 were screened for SCR, 216 for RCR and 336 for CCR (Supplementary Table S4). Subsequently, 17, 9 and 21 ICR genes were detected at the two stages in SCR, RCR and CCR, respectively (Supplementary Table S4).
Gene functional annotation analysis of SCR
Gene functional annotation analysis was conducted for the ICR of SCR. Cytoplasm and DNA binding were especially enriched in SCR3, while oxidation-reduction process and protein binding were particularly enriched in SCR5. Transcription factor activity, sequence-specific DNA binding regulation of transcription, DNA templating, integral components of the membrane and plasma membrane were enriched in both SCR3 and SCR5 (Supplementary Fig. S2).
To gain more biological information and regulatory networks, the Kyoto Encyclopedia of Genes and Genomes (KEEG) enrichment pathways were performed to understand the molecular mechanism of the Cd response in rice seedlings. KEGG analysis results showed that isoflavonoid biosynthesis, inositol phosphate metabolism, flavonoid biosynthesis, and anthocyanin biosynthesis pathways were enriched in SCR3 (Fig. 3A). The isoflavonoid biosynthesis and anthocyanin biosynthesis pathways were significantly enriched. One DEG involved in anthocyanin biosynthesis was downregulated, while one was upregulated in Cd3_CH891. Two DEGs involved in isoflavonoid biosynthesis were upregulated in Cd3_CH891. Phenylpropanoid biosynthesis, alpha-linolenic acid metabolism, linoleic acid metabolism, glutathione metabolism, anthocyanin biosynthesis, ABC transporters, brassinosteroid biosynthesis and diterpenoid biosynthesis pathways were enriched in SCR5 (Fig. 3B). However, phenylpropanoid biosynthesis, ABC transporters, brassinosteroid biosynthesis and diterpenoid biosynthesis pathways were the four extremely enriched pathways in SCR5. Eight DEGs involved in ABC transporters were all upregulated in Cd5_CH891. One DEG involved in phenylpropanoid biosynthesis was downregulated, while seventeen were upregulated in Cd5_CH891. Five DEGs involved in brassinosteroid biosynthesis were all upregulated, and six DEGs involved in diterpenoid biosynthesis were all upregulated in Cd5_CH891. Interestingly, most of the DEGs involved in the enriched pathways were upregulated in CH891 under Cd stress.
Gene functional annotation analysis of RCR
Plasma membrane, integral component of membrane and extracellular region were enriched in ICR of RCR3 in the GO enrichment analysis, while cytoplasm and chloroplast were enriched in the ICR of RCR5. Transcription factor activity, sequence-specific DNA binding protein binding, regulation of transcription and DNA templating were enriched in both RCR3 and RCR5 (Supplementary Fig. S2).
In terms of KEGG pathway enrichment analysis, ribosome, regulation of autophagy and the insulin resistance pathways were enriched in the ICR of RCR3 (Fig. 3C). Among these pathways, ribosome and regulation of autophagy pathways were the most enriched. One DEG involved in the regulation of autophagy was downregulated, and the other two were upregulated in Cd3_02428. One DEG involved in ribosomes was upregulated, and six were downregulated in Cd3_02428. Protein processing in the endoplasmic reticulum, plant hormone signal transduction, thiamine metabolism, glucosinolate biosynthesis and base excision repair pathways were enriched in the ICR of RCR5 (Fig. 3D). Protein processing in the endoplasmic reticulum and plant hormone signal transduction pathways were the most enriched. Four DEGs involved in protein processing in the endoplasmic reticulum were upregulated, and the other four were downregulated in Cd5_02428. One DEG involved in plant hormone signal transduction was upregulated, and the other ten were downregulated in Cd5_02428.
Gene functional annotation analysis of CCR
Regarding the ICR of CCR, chloroplast and zinc ion binding were enriched on the third day, while ATP binding was enriched on the fifth day. Protein binding regulation of transcription, DNA templating, cytoplasm, plasma membrane, integral component of membrane, transcription factor activity, sequence-specific DNA binding and DNA binding were enriched in both CCR3 and CCR5 (Supplementary Fig. S2).
In terms of KEGG pathway analysis, alpha-linolenic acid metabolism, phagosome, limonene and pinene degradation and fatty acid degradation pathways were enriched in CCR3 (Fig. 3E). Moreover, alpha-linolenic acid metabolism and fatty acid degradation pathways were markedly enriched in CCR3. One DEG involved in alpha-linolenic acid metabolism was downregulated, and two DEGs were upregulated. One DEG involved in fatty acid degradation was upregulated, while the other was downregulated. Plant-pathogen interaction, isoflavonoid biosynthesis, vitamin B6 metabolism, diterpenoid biosynthesis, and glutathione (GSH) metabolism pathways were enriched on the fifth day (Fig. 3F). Among these pathways, the plant-pathogen interaction, isoflavonoid biosynthesis, and glutathione (GSH) metabolism pathways were observably enriched. One DEG involved in isoflavonoid biosynthesis was downregulated, and two were upregulated. Two DEGs involved in glutathione metabolism were downregulated, and three were upregulated. Seven DEGs involved in plant−pathogen interactions were downregulated, and twenty were upregulated.
To understand the interaction of all of the enriched GO terms at all stages of SCR, RCR and CCR in general, we constructed a network of significantly enriched GO terms (biological process) (Fig. 4). The biological process was focused on transport, response to different stimuli, metabolism of hormones, gene expression, cellular activity, and secondary metabolic processes. This result suggested that different stress reactions were activated to respond to the Cd stress.
DEGs on both the 3rd and 5th days after stress (DAS)
In total, we obtained 46 common ICRs in three categories at two stages, and they were involved in several different pathways (Fig. 6A, B, C). Except for the pathways listed above, some other pathways were detected, such as inositol phosphate metabolism, cutin, suberine, and wax biosynthesis, biosynthesis of amino acids, starch and sucrose metabolism, ascorbate and aldarate metabolism, base excision repair, aminoacyl-tRNA biosynthesis, galactose metabolism, and circadian rhythm–plant. This result indicated that not only the enriched pathways above but also insignificantly enriched pathways interacted to regulate Cd stress in rice.