Transcriptome and physiological analyses for revealing genes involved in wheat response to endoplasmic reticulum stress

Abstract

Background: Wheat production is largely restricted by adverse environmental stresses. Under many undesirable conditions, endoplasmic reticulum (ER) stress can be induced. However, the physiological and molecular responses of wheat to ER stress remain poorly understood. We used dithiothreitol (DTT) and tauroursodeoxycholic acid (TUDCA) to induce or suppress ER stress in wheat cells, respectively, with the aim to reveal the molecular background of ER stress responses using a combined approach of transcriptional profiling and morpho-physiological characterization.

Methods: To understand the mechanism of wheat response to ER stress, three wheat cultivars were used in our pre-experiments. Among them, the cultivar with a moderate stress tolerance, Yunong211 was used in the following experiments. We used DTT (7.5 mM) to induce ER stress and TUDCA (25 μg·mL^{- 1}) to suppress the stress. Under three treatment groups (Control, DTT and DTT + TUDCA), we firstly monitored the morphological, physiological and cytological changes of wheat seedlings. Then we collected leaf samples from each group for RNA extraction, library construction and RNA sequencing on an Illumina Hiseq platform. The sequencing data was then validated by qRT-PCR.

Results: Morpho-physiological results showed DTT significantly reduced plant height and biomass, decreased contents of chlorophyll and water, increased electrolyte leakage rate and antioxidant enzymes activity, and accelerated the cell death ratio, whereas these changes were all remarkably alleviated after TUDCA co-treatment. Therefore, RNA sequencing was performed to determine the genes involved in regulating wheat response to stress. Transcriptomic analysis revealed that 8204 genes were differentially expressed in three treatment groups. Among these genes, 158 photosynthesis-related genes, 42 antioxidant enzyme genes, 318 plant hormone-related genes and 457 transcription factors (TFs) may play vital roles in regulating wheat response to ER stress. Based on the comprehensive analysis, we propose a hypothetical model to elucidate possible mechanisms of how plants adapt to environmental stresses.

Conclusions: We identified several important genes that may play vital roles in wheat responding to ER stress. This work should lay the foundations of future studies in plant response to environmental stresses.

Keywords: Antioxidant enzymes; Chlorophyll; DTT; RNA-seq; TUDCA; Transcription factors.

Fig. 1

Morphological changes of wheat seedlings under different treatments after two days. (a) Whole view of wheat seedlings. (b, c) Seedling height and root length. (d, e) Fresh weight and dry weight. Different letters indicate significant difference among treatments at the 0.05 significance level based on Duncan’s multiple range tests. Bars represent the mean ± SD (n = 3)

Fig. 2

Physiological and biochemical changes under different treatments after two days. (a) Chlorophyll a content. (b) Chlorophyll b content. (c) Electrolyte leakage rate. (d) Water content. (e) SOD activity. (f) CAT activity. Different letters indicate significant difference among treatments at the 0.05 significance level based on Duncan’s multiple range tests. Bars represent the mean ± SD (n = 3)

Fig. 3

Comparison of wheat leaf and root under different treatments by trypan blue staining. (a) Trypan blue staining in leaf after 4-day’s treatment under microscope (X4). (b) Cell death ratio of leaf after 4-day’s treatment. (c, d) Trypan blue staining in seedling root after 1-day’s treatment under digital camera (c) Root system; (d) Root tip. (e) Root tip under microscope (X10). Bar = 500 μm in A and bar = 200 μm in E. Different letters of B indicate significant difference among treatments at the 0.05 significance level based on Duncan’s multiple range tests. Bars represent the mean ± SD (n = 3)

Fig. 4

Summary of RNA sequencing data. (a) Histogram of DEGs number under different groups. (b) Venn diagram showing the number of DEGs between every two groups and the number of joint DEGs. (c) Hierarchical clustering of DEGs, using the RNA sequencing data derived from three treatments based on log₁₀ (FPKM+ 1) values. The red bands indicate the higher expression, and the blue bands show the lower expression. (d) Gene expression pattern analysis of DEGs. The four subclusters obtained by K-means algorithm. Expression ratios are expressed as log2 values. The X-axis represents different treatments and the Y-axis represents the relative gene expression. C, control; D, DTT; T, DTT + TUDCA

Fig. 5

GO classification of DEGs under group “D vs. C”. The top 30 GO terms were determined by the corrected P-values. The X-axis indicates the number of genes, and the Y-axis is the enriched GO terms. Different colors are used to distinguish biological process, cell component, and molecular function, with “*” as the significantly enriched GO terms. C, control; D, DTT

Fig. 6

KEGG pathway classification of DEGs under group “D vs. C”. The Y-axis represents the pathway name, and the X-axis represents the rich factor. The size of the dot represents the number of DEGs in the pathway, and the color of the dot corresponds to a differently q-value (corrected P-value) range. The senior bubble was obtained by OmicShare tools, a free online platform for data analysis (www.omicshare.com/tools)

Fig. 7

DEGs relevant to the “protein processing in endoplasmic reticulum” pathway under group “D vs. C”. Red boxes refer to genes whose associated DEGs were un-regulated under DTT treatment, and green boxes refer to genes whose associated DEGs were down-regulated. Boxes with yellow color indicate genes that might have isoforms

Fig. 8

DEGs relevant to the “plant hormone signal transduction” pathway under group “D vs. C”. Red boxes refer to genes whose associated DEGs were un-regulated under DTT treatment, and green boxes refer to genes whose associated DEGs were down-regulated. Boxes with yellow color indicate genes that might have isoforms

Fig. 9

Cluster analysis of 371 DEGs under group “T vs. D” among three treatments. (a) Hierarchical clustering of 371 DEGs under group “T vs. D”, using the RNA sequencing data derived from three treatments based on log₁₀ (FPKM+ 1) values. The red bands indicate the higher expression, and the blue bands show the lower expression. (b) Gene expression pattern analysis of 371 DEGs between DTT and DTT + TUDCA under different treatments. The 6 subclusters obtained by h-cluster algorithm. Expression ratios are expressed as log2 values. The X-axis represents different treatments and the Y-axis represents the relative gene expression. C, control; D, DTT; T, DTT + TUDCA

Fig. 10

GO classification of DEGs under group ‘T vs. D’. The top 30 GO terms were determined by the corrected P-values. The X-axis indicates the number of genes, and the Y-axis is the enriched GO terms. Different colors are used to distinguish biological process, cell component, and molecular function, with “*” as the significantly enriched GO terms. D, DTT; T, DTT + TUDCA

Fig. 11

KEGG pathway classification of DEGs under group “T vs. D”. The Y-axis represents the pathway name, and the X-axis represents the rich factor. The size of the dot represents the number of DEGs in the pathway, and the color of the dot corresponds to a different q-value (a corrected P-value) range. The senior bubble was obtained by OmicShare tools, a free online platform for data analysis (www.omicshare.com/tools)

Fig. 12

Summary of transcription factor data. (a) Pie chart showing top 9 TF families which contain more than 50% of differentially expressed TFs. (b) Hierarchical clustering of TFs, using the RNA sequencing data derived from three treatments based on log₁₀ (FPKM+ 1) values. The red bands indicate the higher expression, and the blue bands show the lower expression. (c) Gene expression pattern analysis of TFs. The 4 subclusters obtained by K-means algorithm. Expression ratios are expressed as log2 values. The X-axis represents different treatments and the Y-axis represents the relative gene expression. C, control; D, DTT; T, DTT + TUDCA

Fig. 13

Validation of RNA sequencing data using qRT-PCR. The log₂ fold changes between group “D vs. C” or group “T vs. D” (Y-axis), were plotted against the log₂ fold changes of the same comparison, and determined through RNA-seq (X-axis). The function of the regression line and the R² are given

Fig. 14

A hypothetical model for wheat response to ER stress