Phenotypic variation between CH891 and 02428 under Cd stress
Following exposure of seedlings of rice strains CH891 and 02428 to Cd, large phenotypic variations. 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 control condition, while there was no significant differences between the two genotypes with Cd stress (Fig. 1E, F). It indicated that CH891 was more sensitive to Cd than 02428. The experiment carried out continuously for 7 days and found that the top two seedling length inhibition rate were on 4th and 6th day (Fig. 1G). It indicated that genes had higher expression in the 3rd and the 5th day than other days. Therefore, seedlings were harvested at the 3rd and 5th day and then sequenced RNAs.
RNA sequencing of seedling transcriptome of the two genotypes
A total of 1,070 million reads with average of 46 million reads per sample were generated. While the total of valid reads is 1,057 million and the average valid read per sample is 46 million. The ratio of Q20 for each was above 99%, and the Q30 base percentage was above 95% (Supplementary Table S1). Therefore, the quality of the data was very high and meet the requirements for further analysis. To confirm the accuracy and reproducibility of the RNA-Seq results, 16 genes were selected for qRT-PCR using the specific primers (Supplementary Table S2). 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 S2). To conclude our transcriptome sequencing results were credible.
Identification of differentially expressed genes (DEGs)
By comparing samples of the same rice cultivar in different conditions (control and stress) and different rice cultivar (CH891 and 02428) in the same condition at different stages, we assigned four comparison groups at the two stages and eight comparison groups in total were obtained. And then selected DEGs in different comparision respectively by restricting pval ≤ 0.05. It indicated that not only the cultivar but also the treatments and time points affect the gene expression level. The number of up-regulated and down-regulated genes in each comparison group is shown in Fig. 3A.
Functional enrichment analysis was performed for all these DEGs (Supplementary Fig. S1). KEGG enrichment analysis showed that 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 transporters 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 transporters pathways were enriched in Cd5_CH891 vs. CK5_CH891. Protein processing in endoplasmic reticulum, porphyrin and chlorophyll metabolism, photosynthesis, phenylpropanoid biosynthesis and circadian rhythm-plant pathways were enriched in Cd5_02428 vs. CK5_02428.
All of the results revealed that there were more up-regulated DEGs under Cd stress than those under control condition in CH891, while there were more down-regulated DEGs under Cd stress than those under control condition 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 those in 02428 (Cd-resistant cultivar).
Classification of differentially expressed genes (DEGs)
A total of 7204 unique DEGs were detected in four groups on the 3rd day and 6670 unique DEGs obtained on the 5th day, these DEGs divided into 15 subgroups (Fig. 3B, C). Excluding these DEGs of the groups (CK3_CH891 vs. CK3_02428 at the third day and the groups contained CK5_CH891 vs. CK5_02428 at the fifth day) irrelevant to Cd stress, the DEGs of the rest seven subgroups divided into three categories: genes from sensitive variety with Cd-responsive (SCR), genes from resistant variety with Cd-responsive (RCR), and common Cd-responsive (CCR) DEGs. There were 849, 676 and 770 DEGs in SCR, RCR and CCR at the 3rd day, respectively (Supplementary Table S3). Whereas 959, 592 and 1516 DEGs were in SCR, RCR and CCR at the 5th day, respectively (Supplementary Table S3).
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. At the 3rd day, there were 554 ICR genes were screened, among which 194 for SCR, 154 for RCR and 206 for CCR (Supplementary Table S3). Whereas a total of 941 ICR genes were screened at the 5th day, of these 389 for SCR, 216 for RCR and 336 for CCR (Supplementary Table S3). Subsequently, there were 17, 9 and 21 ICR genes were detected at the two stages in SCR, RCR and CCR, respectively (Supplementary Table S3).
Gene functional annotation analysis of SCR
Gene functional annotation analysis was conducted for ICR of SCR. Cytoplasm and DNA binding were especially enriched in SCR3, oxidation-reduction process and protein binding were particularly enriched in SCR5. Transcription factor activity, sequence-specific DNA binding regulation of transcription, DNA-templated, integral component of membrane and plasma membrane were enriched in both SCR3 and SCR5 (Supplementary Fig. S2). In order to gain more biological information and regulatory network, the Kyoto Encyclopedia of Genes and Genomes (KEEG) enrichment pathways was performed for understanding the molecular mechanism of 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. 2A). While isoflavonoid biosynthesis and anthocyanin biosynthesis pathways were significantly enriched. One DEG involved in anthocyanin biosynthesis was down-regulated while one was up-regulated in Cd3_CH891. And two DEGs involved in isoflavonoid biosynthesis were up-regulated 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. 2B). However, phenylpropanoid biosynthesis, ABC transporters, brassinosteroid biosynthesis and diterpenoid biosynthesis pathways were four extremely enriched pathways in SCR5. Eight DEGs involved in ABC transporters were all up-regulated in Cd5_CH891. One DEG involved in phenylpropanoid biosynthesis was down-regulated while seventeen were up-regulated in Cd5_CH891. Five DEGs involved in brassinosteroid biosynthesis were all up-regulated, and six DEGs involved in diterpenoid biosynthesis were all up-regulated in Cd5_CH891. Interestingly, most of DEGs involved in enriched pathways were up-regulated 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 ICR of RCR5. Transcription factor activity, sequence-specific DNA binding protein binding, regulation of transcription and DNA-templated were enriched in both RCR3 and RCR5 (Supplementary Fig. S2).
In terms of KEGG pathway enrichment analysis, ribosome, regulation of autophagy and insulin resistance pathways were enriched in ICR of RCR3 (Fig. 2C). Among which ribosome and regulation of autophagy pathways were most enriched. One DEG involved in regulation of autophagy was down-regulated and the other two were up-regulated in Cd3_02428. One DEG involved in ribosome was up-regulated and six were down-regulated in Cd3_02428. Protein processing in endoplasmic reticulum, plant hormone signal transduction, thiamine metabolism, glucosinolate biosynthesis and base excision repair pathways were enriched in ICR of RCR5 (Fig. 2D). And protein processing in endoplasmic reticulum and plant hormone signal transduction pathways were most enriched. Four DEGs involved in protein processing in endoplasmic reticulum were up-regulated and other four were down-regulated in Cd5_02428. One DEG involved in plant hormone signal transduction was up-regulated and the other ten were down-regulated in Cd5_02428.
Gene functional annotation analysis of CCR
As far as ICR of CCR, chloroplast and zinc ion binding were enriched on the third day while ATP binding was enriched in the fifth day. And protein binding regulation of transcription, DNA-templated, 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. 2E). Moreover, alpha-linolenic acid metabolism and fatty acid degradation pathways were markedly enriched in CCR3. One DEG involved in alpha-linolenic acid metabolism was down-regulated and two DEGs were up-regulated. And one DEG involved in fatty acid degradation was up-regulated while the other one was down-regulated. Plant-pathogen interaction, isoflavonoid biosynthesis, vitamin B6 metabolism, diterpenoid biosynthesis, and glutathione (GSH) metabolism pathways were enriched in the fifth day (Fig. 2F). Among which plant-pathogen interaction, isoflavonoid biosynthesis, and glutathione (GSH) metabolism pathways were observably enriched. One DEG involved in isoflavonoid biosynthesis was down-regulated, two were up-regulated. Two DEGs involved in glutathione metabolism were down-regulated, three were up-regulated. Seven DEGs involved in plant − pathogen interaction were down-regulated and twenty were up-regulated.
To understand the interaction of all the enriched GO terms at all the stages of SCR, RCR and CCR in the general, we constructed a network of significantly enriched GO terms (biological process) (Fig. 4). It showed that the biological process was focus on transport, response to different stimulus, metabolic of hormone, gene express, cellular activity, and secondary metabolic process. It suggested that different stress reactions were activated to respond to Cd stress.
DEGs both in 3rd and 5th day after stress (DAS)
In total, we obtained 46 common ICR in three categories at two stages and they were involved in several different pathways (Fig. 5C, D, E). 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. It indicated that not only the enriched pathways above but also insignificantly enriched pathways had interaction to regulate Cd stress in rice.