Phenotypic Changes of Sugarcane Tillering Seedlings Under Drought Stress
“Guire 2” were planted in plastic pots and treated with drought at the tiller formation stage. The experiment applied drought stress by stop watering. The changes of soil water content (SWC) in the control group and drought treatment groups were shown in Fig. 1A. After 15 days of drought stress, the SWC in the treatment group decreased significantly from 22.25–7.12%. After 15 days of well watering, the SWC in the control group was maintained at about 21%~24%. From 5 days onwards, SWC was significantly different between both groups, consistently.
In the first 13 days of drought stress, experience group and watering control group did not have significant differences in tiller height or stem diameter. However, On the 15th day of drought stress, the height of tillering seedlings in the experimental group was significantly lower than that of the control group (P < 0.05), and meanwhile, the stem diameter was also significantly lower than that of the control group (P < 0.01) (Fig. 1B, C).
The leaves of “Guire 2” plant turn yellow after the 9th day of drought stress. After 15 days of drought treatment, “Guire 2” showed distinct wilting (Fig. 1D). Compared to the whole plant of sugarcane “Guire 2”, the tiller seedlings were less affected by drought stress. The phenotypes of tiller seedlings in the treatment and control groups were almost identical in the pre-drought (5th day) and mid-drought periods (9th day). And yet, the growth of “Guire 2” tiller seedlings becomes sluggish and tiller leaves begin to curl in late-drought (15th day). (Fig. 1E).
Base on the SWC and tiller phenotypic changes, we selected sugarcane that was drought stressed on day 5 (T1), day 9 (T2) and day 15 (T3) as the experimental groups, and sugarcane watered normally on day 5 (CK1), day 9 (CK2) and day 15 (CK3) as the control groups. The above six sets of tiller seedlings were further used for physiological assays and transcriptome sequencing.
Physiological Changes of “Guire 2” Tiller Seedlings Under Drought Stress
In order to study the physiological responses of sugarcane 'Guire 2' tiller seedlings to drought stress, the physiological indices such as superoxide dismutase (SOD), peroxide dismutase (POD), malondialdehyde (MDA), indoleacetic acid (IAA) and zeatin (ZT) were evaluated at the tiller formation stage.
In terms of oxidative stress, drought stress results in the production of reactive oxygenspecies (ROS), and excessive ROS would lead to oxidativestress, inhibit plant growth, and even cause cell death. Antioxidant enzymes such as SOD and POD form a defense system and play critical roles in removing the excessive ROS induced by drought stress. In our researcher, SOD activity increased from 482.82 U/g to 1042.90 U/g with the in creasing of drought treatment time. Compared with the treated group, the SOD activity of the control group was always maintained at a lower level (Fig. 2A). Results of POD activity showed that there was no significant difference between the drought stress and control on day 5 and day 9. Neverless, the drought stress group was significantly (p < 0.05) higher than the control group at 15 days of drought stress (Fig. 2B). Osmoregulatory substances play an equally important role in the response of plants to drought stress. In terms of osmoregulation, the result of proline (Pro) content showed no significant difference between the treatment and control groups on day 5 day 9 of drought stress. But, the Pro content of the control group increased significantly on the 15 day of drought stress (Fig. 2C). Soluble sugar (SS) and soluble protein (SP) contents, as well as starch content, remained consistent throughout the drought treatment and increased as drought stress increased. On day 5, there was no significant difference between the SS, SP and starch of treated and control groups. However, on the 15th day, in the treated group, SS, SP and starch content were significantly (P < 0.05) raised when compared with the control group (Fig. 2D,E and F). MDA is an important signal to evaluate plant cell damage. Results of MDA content showed that the drought-stressed group was consistently higher than the control group. But surprisingly, the MDA content did not increase with the increase of drought stress time and always maintained at a stable level (Fig. 2G). The data changes of MDA content imply that drought stress did not cause greater cellular damage to tiller seedlings within the experimental range.
Drought stress leads to changes in endogenous hormones from sugarcane tiller seedlings at each time point, as shown in Supplementary Fig. 1. A significant (P < 0.05) increase in ABA content was observed in all drought stress treatments over the control (Supplementary Fig. 1A ). Drought stress treatment significantly (P < 0.05) increased GA3 levels only on day 5 compared with the control, but in response to the increase of drought stress, GA3 content first decreases and then increases (Supplementary Fig. 1B ). Similarly, the SLs content of the drought stress group also showed the same trend (Supplementary Fig. 1E ). The IAA content exhibited a greater increase at 15 days than 5 and 9 days of drought stress. In addition, on day 9, the IAA content in the control group was significantly (P < 0.05) higher than that in the drought-stressed group (Supplementary Fig. 1C). Result of ZT content displayed that the differences between drought-stressed group and control group were significant at day 5 and 9, whereas on day 15, there were no longer significant differences between drought-stressed group and control group (Supplementary Fig. 1D).
RNA Sequencing Analysis of “Guire 2” Tiller Seedlings Under Drought Stress
To establish the key genes of sugarcane tiller seedlings in response to drought, the treated tiller seedlings of “Guire 2” were sequenced. A total of approximately 151.64 million raw reads were generated from the 18 ( 6 samples ×3 replications) cDNA libraries. After deleting 0.03 ~ 0.04% of adapter sequences, and filtering 0.31 ~ 0.41% of low-quality reads and 0.00% of n-containing reads, a total of approximately 151.01 million high-quality clean reads were finally confirmed (Supplementary table 1), with an average Q20 and Q30 value of 97.82% and 93.84%, respectively (Supplementary table 2). The mapping rate to the reference genome ranges from 82.91–88.31% (Supplementary table 3). These results showed that the transcriptome sequencing quality was sufficient for further analyses.
Normally, a stringent threshold absolute log2 FC ≥ 2 and FDR < 0.05 was used to screen out DEGs, GO enrichment and KEGG pathway analysis were applied to assay these DEGs to further understand their biological function. After 5 days of drought treatment, 18 up-regulated DEGs and 12 down-regulated DEGs were identified. The DEGs after 9 days of drought treatment, 1482 were up-regulated and 4381 were down-regulated. After 15 days of continuous drought treatment, the number of DEGs was the largest in the three treatment with 6634 up-regulated genes and 12084 down-regulated genes (Fig. 3A). On the whole, the number of down-regulated DEGs is higher than up-regulated DEGs, except for drought stress for 5 days. Together, the results revealed that the number of induced DEGs greatly increased with the continuation of drought stress time. Venn diagram shows one common DEG (MSTRG.65561) was consistently involved in CK1 vs T1, CK2 vs T2, and CK3 vs T3 (Fig. 3B), and the expression of this gene increased substantially with the enhancement of drought stress. However, MSTRG.65561 could not be annotated according to the KEGG and GO databases, and consequently its biological function could not be determined (Supplementary table 4).
In order to understand the function of these DEGs, GO enrichment analysis was performed for DEGs. GO enrichment analysis of CK1 vs T1 showed that the qvalue of GO enrichment was greater than 0.01, which was not analytically significant(Supplementary table 5). In CK2 and T2, DEGs were mainly enriched in 3 biological processes, include “photosynthesis, light reaction ”(GO:0019684), “photosynthesis” (GO:0015979) and “external encapsulating structure organization” (GO:0045229). Especially, the rich factor of GO:0019684 and GO:0015979 associated with photosynthesis were higher. Among the top 20 significantly enriched GO terms, 17 GO terms belong to cellular composition, they were “photosystem I” (GO:0009522), “photosystem” (GO:0009521), “plastid thylakoid” (GO:0031976), “photosynthetic membrane” (GO:0034357), “plastid thylakoid membrane” (GO:0055035), “chloroplast thylakoid membrane”(GO:0009535), “thylakoid part” (GO: 0044436), “chloroplast thylakoid” (GO:0009534), “thylakoid” (GO:0009579), “thylakoid membrane” (GO:0042651), “cell wall” (GO:0005618), “chloroplast part” (GO:0044434), “plastid envelope” (GO:0009526), “external encapsulation structure” (GO:0030312), “plastid part” (GO: 0044435), “chloroplast” (GO:0009507) and “plastid” (GO:0009536), respectively (Fig. 3C and supplementary table 6). Photosynthesis was seriously impacted in “Guire 2” tiller seedlings at 9 days of drought stress, as evidenced by GO enrichment analysis of CK2 vs T2.
The DEGs of CK3 vs T3 goup were subjected to a GO analysis comprising 7 biological processes, 8 cellular components, and 5 molecular functions among the top 20 significantly enriched GO terms (Fig. 3D and supplementary table 6). GO terms of molecular functions were increased, which were “catalytic activity” (GO:0003824), “hydrolase activity, acting on glycosyl” (GO:0016798), “tubulin binding” (GO:0015631), “motor activity” (GO:0003774), “hydrolase activity, hydrolyzing O-glycosyl compounds” (GO:0004553), respectively. As shown in Fig. 3D, most genes were down-regulated in GO terms of molecular function, thus we speculate that the hydrolase activity of tiller seedlings was reduced at this time due to drought stress. These GO enrichments of cellular components include “cell wall” (GO:0005618), “external encapsulating structure” (GO:0030312), “photosystem” (GO:0009521), “microtubule associated complex” (GO:0005875), “kinesin complex” (GO:0005871), “microtubule” (GO:0005874), “apoplast” (GO:0048046), “extracellular region” (GO:0005576). GO term for biological processes include: “external encapsulating structure organization” (GO:0045229), “cell wall organization or biogenesis” (GO:0071554), “cell wall organization” (GO:0071555), “plant-type cell wall organization or biogenesis” (GO:0071669), “phenylpropanoid metabolic process” (GO:0009698), “plant-type cell wall biogenesis occurrence” (GO:0009832), “DNA integration” (GO:0015074). We found that a number of significantly enriched GO items were associated with the development of cell wall tissue, and in particular, most of DEGs were down-regulated (Fig. 3D). The results indicate possibly that the cell wall development of tiller seedlings was negatively affected, thus corroborated our previous results that the height and stem diameter traits of the "Guire 2" treatment group were significantly lower than those of the control group (Fig. 1A, B).
The results of the KEGG enrichment analysis are shown in Fig. 3E and F with the first 20 top-ranking pathways indicated by the smallest significant qvalues. KEGG pathway enrichment analysis indicated that DEGs involved in photosynthesis, sugar metabolism and fatty acid metabolism pathways in both the mid ( CK2 vs T2 ) and late drought periods ( CK3 vs T3 ). For example, “Photosynthesis” (ko00195), “Photosynthesis - antenna proteins” (ko00196), “Carbon fixation in photosynthetic organisms” (ko00710) are related to photosynthesis. “Fructose and mannose metabolism” (ko01212), “Starch and sucrose metabolism” (ko00500) are related to sugar metabolism. “Carbon fixation in photosynthetic organisms”(ko00710)、”Fatty acid elongation”༈ko00062༉are related to fatty acid metabolism.
WGCNA of “Guire 2” Tiller Seedlings Under Drought Stress
The weighted gene co-expression network analysis (WGCNA) was used to analyze the connection between genes and physiological traits, discovering the hub genes associated with physiological and biological traits. Genes with a max FPKM of < 15were filtered out, and the 16415 selceted genes were assigned to 20 merged co-expression module with various colors (Supplementary Fig. 2). The weight value indicates the correlation of gene relationship pairs in the module, one relationship pair for every two genes. We selected the top 100 relationship pairs in terms of weight to draw the network regulatory map, and then selected the five genes with the highest degree of connectivity as the hub genes from them. It is commonly recognized that hub genes form the backbone of the network and tend to be critical in particular physiological events14.
As showed in Fig. 4A, we successfully identified five modules signicantly ( r > 0.7, p < 0.01) associated with physiological or biological traits of “Guire 2” tiller seedlings under drought stress. Based on the analysis results of the WGCNA co-expression module, the coexpression network of the darkgreen, grey60, skyblue, cyan and greenyellow modules were further analyzed. Additionally, the darkgreen module was positively correlated with SOD (r = 0.82, p = 4e-05), Pro (r = 0.82, p = 3e-05), sucrose (r = 0.82, p = 3e-05), protein (r = 0.78, p = 1e-04) and starch content (r = 0.83, p = 2e-05) under drought stress. The grey60 module was positively correlated with SOD content (r = 0.72, p = 8e-04) and also showed a good positive correlation with other physiological indicators. The skyblue module had a significant positive correlated with Pro (r = 0.74, p = 4e-04) and sucrose content (r = 0.72, p = 7e-04), while the cyan module was negatively correlated with ZT (r = -0.88, p = 0.003), SOD (r = -0.73, p = 6e-04), Pro (r = -0.7, p = 0.001) and sucrose (r = -0.79, p = 8e-05) content under drought stress. Particularly, the greenyellow module had a significant positive correlated with tiller stem diameter (r = 0.74, p = 5e-04).
We found that the darkgreen, grey60 and skyblue module were positively correlated with the physiological of "Guire 2" tiller seedlings (Fig. 4A). The correlation network of the darkgreen module is shown in Fig. 4B, Sspon.02G0045390-1B (LSG1-2), MSTRG.102651 (ERF1-2), Sspon.03G0003830-1A, Sspon.06G0000380-1A (SHKA), Sspon.03G0035940-2C were identified as candidate hub genes for this module(Supplementary table 8). In darkgreen module, eight significant pathways (qvalue < 0.05) are identified, and “Spliceosome” (ko03040) pathway was the most significant enrichment of genes (Fig. 4C). .
The correlation network of the grey60 module is shown in Fig. 5A, Sspon.04G0008400-1A (TIL), Sspon.01G0038620-2C (HSP18.1), Sspon.04G0002960-1A (HSP24.1), Sspon.01G0038620-3D (HSP 18.1), Sspon.03G0019770-1A (HSP16.1), MSTRG.111550 (HSFA6A) were regarded as module candidate hub genes(Supplementary table 8). Among these candidate hub genes, four genes belong to small molecule heat shock proteins (sHSPs), indicating that sHSP regulatory network may play a major role in SOD regulatory of drought stress. In grey60 module, KEGG pathway analysis showed that “Protein processing in endoplasmic reticulum” (ko04141), “Plant-pathogen interaction” (ko04626), “Endocytosis” (ko04144) and “Spliceosome” (ko03040) were most significantly (qvalue < 0.05) enriched by these genes (Fig. 5B).
As showed in Fig. 5C, the skyblue module could be signifiantly (qvalue < 0.05) enriched in “Phagosome” (ko04145), “Inositol phosphate metabolism” (ko00562), “Cyanoamino acid metabolism” (ko00460) and “Metabolic pathways” (ko01100). The correlation network of the skyblue module is shown in Fig. 5D, MSTRG.11350, MSTRG.102629, Sspon.06G0016280-4D, MSTRG.13435, Sspon.01G0017100-1A and MSTRG.84847 were identified as candidate hub genes for this module(Supplementary table 8).
The cyan module negatively correlated to physiological traits of tiller seedlings, and the genes in this module may have a suppressive effect on the physiological regulation of tiller seedlings. The correlation network shows that Sspon.01G0023790-2P (IAA30), Sspon.04G0004000-3D (MAP70.3), Sspon.02G0044960-1B (NEK2), MSTRG.92463 (CLASP) and Sspon.03G0023990-1P (At5g45910) were identified as candidate hub genes for the cyan module (Fig. 6A and Supplementary table 8). In additional, For KEGG pathway enrichment analysis as showed in Fig. 6B, the result showed that these genes were significantly (qvalue < 0.05) enriched in “Ubiquitin mediated proteolysis” (ko04120), “Plant hormone signal transduction” (ko04075), “Fatty acid elongation” (ko00062), “Glycerolipid metabolism” (ko00561).
Greenyellow module was positively correlated with stem diameter development of tiller seedlings. The correlation network of the skyblue module is shown in Fig. 6C, Sspon.01G0028080-4D (BHLH148), Sspon.01G0028080-1A (BHLH148), Sspon.04G0034710-2D (WAK3), Sspon.01G0028080-3C (BHLH148) and MSTRG.58338 (CPK4) were identified as candidate hub genes for this module(Supplementary table 8). Notably, members of the BHLH148 were important nodes with abundant connection points in the co-regulatory network, indicating that BHLH148 may be crucial to tiller seedlings development under drought stress. KEGG enrichment analysis for greenyellow module genes as showed in Fig. 6D, the most significantly (qvalue < 0.5) enriched pathways were “MAPK signaling pathway - plant” (ko04016), “Plant-pathogen interaction” (ko04626), “Linoleic acid metabolism” (ko00591), “alpha-Linolenic acid metabolism” (ko00592), “Cysteine and methionine metabolism” (ko00270), “Plant hormone signal transduction” (ko04075), “Phenylalanine metabolism” (ko00360) and “Phenylalanine, tyrosine and tryptophan biosynthesis” (ko00400).
Validation of Candidate Hub Gene Expression Changes by qRT - PCR
To verify the accuracy and reproducibility of the transcriptome analysis, 10 candidate hub genes were randomly selected to analysis the transcript abundance by quantitative real–time PCR (qRT-PCR). The results showed that most of genes the expression profiles detected by qRT-PCR were positively correlated (pearson correlation > 0.9) with the RNA-Seq results, only 3 candida te hub genes was inconsistent with the RNA-Seq results (Fig. 7). Overall, the reliability of the RNA-seq data was confirmed by the consistency between the qRT-PCR results and RNA-seq analyses.