RNA-seq and de novo transcriptome assembly
Total RNA extracted from the collected samples (leaf, stem, spike and root) were pooled and was used for the RNA-seq using Illumina HiSeq 4000 platform. The raw sequences obtained were subjected to different software’s for removing the adaptor sequences, low-quality reads (Q-value< 20) and ambiguous reads (containing more than 10% ‘N’ bases). The clean reads obtained in different replicates has been presented in supplementary Table S2. We acquired clean reads of ~65031642 (P3), 80444546 (M3), and 73850514 (M3H) with more than 92% Q30 bases. The GC content was ~48% and produced comparable data from the biological replicates. De novo assembly of filtered clean reads were performed using the Trinity v 2.8.4. The preliminary assembly generated 3,05,537 genes and 5,88,788 transcripts with an N50 of 1,349 bp. Among the assembled transcripts, the average contigs length was 836 bp (Table S3). We also predicted the isoforms of the genes predicted with total assembled bases of 18,44,83,645 and average contig length of 603.8 and N50 of 835. The assembled transcripts were plotted against the length of transcripts, as shown in figure 1. The transcript length in the range of 200-250 bp was observed most abundant in the assembled transcripts, followed by transcript length of 1000-1500 bp. GC distribution analysis in the assembled transcripts showed maximum GC percent in the range of 40-45%.
Annotation and classification of Triticum aestivum unigenes
All assembled transcripts were subjected to annotation using UniProt Plant Database, gene and protein annotation, organism annotation, gene ontology and pathway annotation. Briefly, the assembled transcriptome sequences were compared with UniPort Plant Database using the NCBI BlastX program with E-value cut-off of 10-5 and identity of 40% were retained for further annotation. Overall, we observed 3,62,502 assembled transcripts having at least one significant hit in UniProt Plant Protein Database. The annotation summary for all assembled transcripts is provided in the supplementary Table S4. The complete annotation is divided into the following categories: Total number of assembled transcripts for blastx was 5,88,788 and we observed 3,62,502 hit in uniprot plant database. Overall, the assembled transcripts showed maximum BLASTx similarities to gene sequences from Triticum aestivum (~2,20,000 transcripts), Triticum uratu (~41,000 transcripts), Aegilops tauchii (40,250 transcripts), etc. Other organisms showing maximum hit for the assembled transcripts are - Brachypodium distachyon, Oryza sativa, Cicer arietinum, etc. The top hits of each unigene and the organism name was extracted. The top 15 organisms are shown in figure 2. A summary of the top hit domains identified 1345 heat-responsive transcription factors involved in various biological processes. Similarly, we identified 579, 113 and 22 transcripts showing homology with HSPs, SODs and CDPKs.
Gene ontology (GO) annotation
The functions of the assembled transcripts were classified via GO analysis. The GO terms for all assembled transcripts were extracted wherever possible and BLAST2GO was used for the classification of groups into 3 major categories (supplementary Table S5). The total number of GO terms identified in biological processes, cellular component and molecular function processes are – 3233, 761 and 1963.
In total, assembled transcripts were categorized into 65 functional groups (Fig. 3), and further categorized into three major group i.e. ‘biological processes’, ‘cell component’, and ‘molecular function’. The dominant sub-categories of the classified genes included ‘regulation of transcription’ (10,680), ‘transcription’ (7,760) and ‘translation’ (5,463) in the ‘biological processes’ category; ‘integral component of membrane’ (61,163), ‘nucleus’ (15,877) and ‘cytoplasm’ (5,665) in the ‘cell component’ category; ‘ATP binding’ (35,127), ‘nucleic acid binding’ (15,154) and ‘protein kinase activity’ (12,473) in the ‘molecular function’ category.
DEGs in wheat mutant (M3) under HS-treated condition
The expressed transcripts in the P3, M3, and M3H were identified using the de novo-assembled transcriptome as a reference. A correlation analysis was carried out based on the FPKM (fragments per kilobase of gene per million mapped reads) values in order to verify the consistency among the samples. We observed very high correlation (Fig. 4; r > 0.91) between the replicates of P3, M3, and M3H samples, which shows very high consistency in the RNA-seq results. In order to identify the differentially expressed transcripts (DETs), a comparative expression analysis was carried out between P3, M3 and M3H samples. The number of assembled transcripts considered for differential gene expression analysis was 588,788. DETs were filtered at log2 fold change cut-off of +2/-2 (p<= 0.05). The differential expression analysis of the annotated transcripts was carried out using Cuffdiff program of cufflinks package. Log2 fold change cut-off 2 (p<=0.01 and 0.05) were separately used as cut-off for upregulated and downregulated genes and isoforms.
On comparative analysis of expression of DEGs in mutant (M3) over parent (P3), we observed 6,120 upregulated and 4,428 downregulated transcripts (p-value < 0.01) (Fig. 4a). Similarly, 11,354 upregulated and 12,408 downregulated genes were observed in HS-treated mutant (M3H) over parent (P3). The effect of HS on expression of DEGs in mutant showed 4817 upregulated and 9085 downregulated genes (M3H over M3). We also predicted the isoforms of identified differentially expressed genes in mutant (M3) and HS-treated mutant (M3H) over parents (P3) (Fig. 4b). We identified 6800 upregulated and 10,120 downregulated transcripts in M3H vs M3; 7196 upregulated and 7400 downregulated transcripts in M3 vs P3. Similarly, comparative analysis of DETs in mutant vs parent showed 14,228 upregulated and 15,485 downregulated transcripts.
We also plotted top 50 DE genes within each pairwise comparison based on the genes and isoforms. Comparative analysis of DE genes in M3 and P3 showed a stretch of transcripts showing upregulation in mutant, as compared to parent and vice versa (Fig. 5a). Similarly, more than 60% of the top 50 DE genes showed upregulation HS-treated mutant, as compared to mutant under control condition. The DE genes which showed upregulation in HS-treated mutant showed downregulation in parent and vice versa. Expression pattern of isoforms of top 50 DE genes within each pairwise comparison showed significant variations in the expression. Comparative analysis of mutant over parent showed a stretch of DETs showing significant upregulation in mutant and downregulation in parent. Similarly, a stretch of isoforms of DETs were identified in parent showing upregulation and downregulation in mutant (Fig. 5b). Comparative analysis of HS-treated mutant over mutant under control condition showed more than 75% of the isoforms of DETs showing upregulation under HS.
We analysed the average expression between the mean of the normalized counts of the samples used in the respective group through MA plot and volcano plot. The base 2 Logarithm of fold change in expression of M3H vs M3 showed much dispersed distribution of upregulated and downregulated transcripts above and below the mean line. Here, we observed more significant number of genes being upregulated, as high number of data points fall above the one threshold on the y-axis. Similarly, large number of genes was observed lying close to the threshold showing non-significant changes in the expression between the M3H and M3. We observed similar pattern of expression in M3 vs P3, though number of transcripts showing significant downregulation was observed more, as compared to upregulated genes. MA plot of M3H over P3 group showed very significant number of genes showing upregulation, as compared to downregulated genes with uniform distribution all over the mean average. We observed large chunk of transcripts lying on the threshold line showing non-significant changes in the Log2 fold expression in HS-treated mutant (M3H), as compared to parent (P3). We also generated the volcano plot analysis of mutant and parent for the visual identification of DEGs with statistically significant variations in the fold expression (Fig. S1).
Functional annotation of DEGs
The DEGs were identified in mutant as compared to parent (M3 vs P3), HS-treated mutant, as comparted to mutant (M3H vs M3) and HS-treated mutant as compared to parent (M3H vs P3). Comparative analysis of mutant over parent (M3 vs P3) showed 6118 upregulated and 4420 downregulated transcripts (Supplementary Table S6). Some of the highly upregulated transcripts identified were glutamate cysteine ligase (TRINITY_DN6079_c0_g1), Fructokinase-2 (TRINITY_DN9001_c0_g1), 70 kDa HSPs (TRINITY_DN437_c5_g1), Protein disulphide isomerase (TRINITY_DN3156_c1_g2), Kinases (TRINITY_DN68271_c0_g2), Peroxidases (TRINITY_DN9374_c1_g1), ABC transporters (TRINITY_DN60417_c0_g2), ATP synthase subunit alpha (TRINITY_DN7452_c0_g2), Superoxide dismutase [Cu-Zn] (TRINITY_DN17716_c0_g1), HSP20-like chaperone superfamily protein (TRINITY_DN32701_c0_g2), etc. (Table 1) Similarly, some of the downregulated transcripts identified in mutant over parent are - Xylose isomerase (TRINITY_DN19775_c0_g1), Fructose-bisphosphate aldolase (TRINITY_DN15530_c0_g1), Major Facilitator Superfamily protein (TRINITY_DN24775_c1_g1), RNA polymerase II-associated protein 1 (TRINITY_DN16713_c0_g1), Sucrose synthase (TRINITY_DN83591_c0_g1), Starch synthase family protein (TRINITY_DN34859_c0_g1), Sugar transporter (TRINITY_DN25168_c0_g1), UDP-sugar pyrophosphorylase (TRINITY_DN14149_c0_g2), Cysteine protease (TRINITY_DN48681_c0_g1), Zinc transporter (TRINITY_DN31037_c0_g1), etc. (Table 1).
Differentially expressed transcripts (DETs) identified in mutant treated with HS over parent (M3H vs M3) showed the presence of 9084 upregulated and 4816 downregulated transcripts (Supplementary table S7). Some of the potential upregulated transcripts identified are - putative chaperone clpb (TRINITY_DN13064_c0_g1), putative heat shock protein (TRINITY_DN21155_c0_g2), Peptidyl prolyl isomerase (TRINITY_DN16896_c0_g1), Calcium-binding EF-hand family protein-like (TRINITY_DN7556_c0_g1), Heat shock transcription factor (TRINITY_DN1205_c0_g1), Universal stress protein family protein (TRINITY_DN3408_c1_g1), Oxidative stress 3 (TRINITY_DN64554_c0_g1), Superoxide dismutase [Cu-Zn] (TRINITY_DN114794_c0_g1), etc. (Table 2). Similarly, some of the identified downregulated DETs are - Cysteine protease (TRINITY_DN14368_c1_g1), Pollen allergen Phl p 5 (TRINITY_DN16967_c0_g1), Late embryogenesis abundant protein (TRINITY_DN42728_c0_g1), L-ascorbate oxidase-like protein (TRINITY_DN4905_c0_g1), Kinase family protein (TRINITY_DN36483_c0_g1), Protein Zinc INDUCED FACILITATOR-LIKE 1 (TRINITY_DN26586_c0_g1), MYB-related transcription factor (TRINITY_DN10260_c0_g3), Pyruvate decarboxylase (TRINITY_DN2112_c1_g3), Alpha/beta-Hydrolases superfamily protein (TRINITY_DN15796_c0_g1), etc. (Table 2).
Further, the DETs were identified by comparing the transcripts of HS-treated mutant over the parent (M3H vs P3). We observed 12409 upregulated and 11533 downregulated transcripts in HS-treated mutant, as compared to parent (Supplementary Table S8). Some of the identified upregulated DETs are - Heat shock protein (TRINITY_DN8393_c7_g1), Acid phosphatase 1 (TRINITY_DN18979_c0_g1), ATP-dependent zinc metalloprotease FTSH protein (TRINITY_DN75467_c1_g1), putative chaperone clpb (TRINITY_DN25343_c0_g1), etc. Similarly, some of the downregulated DETs are - Sucrose synthase (TRINITY_DN4652_c1_g2), Hexosyltransferase (TRINITY_DN55003_c0_g1), Peptide transporter (TRINITY_DN46904_c0_g1), Auxin response factor (TRINITY_DN56713_c0_g1), α-1,4-glucan-protein synthase [UDP-forming] 1 (TRINITY_DN41347_c0_g1), etc. (Table 3).
Overall, we observed maximum transcripts of heat shock proteins showing Log2 fold expression of 17.0. Similarly, other genes showing high Log2 fold expression were HSFs, chaperones, peptidyl propyl isomerase, pentatricopeptide repeat-containing protein, universal stress proteins, sucrose phosphate synthase, etc. Similarly some of the downregulated transcripts identified were - Beta-galactosidase, Xylose isomerase, Fructose-bisphosphate aldolase, Protein disulfide-isomerase, Cytochrome P450 family protein, expressed Hexosyltransferase, Pyruvate dehydrogenase E1 component subunit alpha Histone H2B, Gibberellin-regulated protein 2, putative Sugar transporter, Phosphoglycerate mutase-like protein, etc.
Venn diagram analysis
Venn diagram analysis of DETs identified in M3 vs P3, M3H vs M3 and M3H vs P3 has been shown in figure 6. We observed 911 transcripts showing significant (P<0.01) upregulation in mutant, as compared to parent (Fig. 6a). Similarly, 2604 transcripts showed significant upregulation under HS in mutant. Similarly, 4170 transcripts showed significant upregulation in HS-treated mutant (M3H), as compared to parent (P3). We observed 2155 common transcripts showing upregulation in M3H vs M3 and M3H vs P3. Similarly, 5151 transcripts showed significant upregulation in M3H vs P3 and M3 vs P3. Overall, 57 transcripts showed significant upregulations (P<0.01) in M3H, M3 and P3.
We observed significant (P<0.01) downregulation of 2566 isoforms of DETs in M3 vs P3, 1131 transcripts in M3H vs M3, and 2729 transcripts in M3H vs P3 (Fig. 6b). We observed 1 common transcripts showing significant downregulation in M3 vs P3 and M3H vs M3. Similarly, 7824 transcripts were observed downregulated in M3H, as compared to M3 and P3; 1726 transcripts showed downregulation in M3H and M3, as compared to P3. Overall, we observed 128 transcripts showing significant (P<0.01) downregulation in HS-treated mutant (M3H), mutant (M3) and parent (P3).
Gene ontology enrichment analysis of DEGs
We have performed GO enrichment analysis using GOATOOLS (Fisher exact test, P-value ≤0.05) in order to get insight about the functional categories of the DEGs (Klopfenstein et al.36). The total numbers of GO terms identified under different categories are – molecular function (184815), biological processes (100390), and cellular component (114857). The most enriched GO category among these DEGs was ‘ATP binding’ (35127), followed by ‘nucleic acid binding’ (15154), ‘protein kinase activity’, ‘serine threonine kinase activity’, ‘regulation of transcription’ (11717)`, ‘translation’ (1726), ‘DNA integration’, ‘metabolic process’ (17015), ‘carbohydrate metabolic process’ (4412), ‘integral component of membrane’ (5717), ‘nucleus’ (15888), ‘chloroplast’ (5835), ‘cytoplasm’ (5776), and ‘ribosome’ (4808).
Biological pathway annotations
To identify the active biological pathways, functional analysis of proteins, cluster of orthologous groups (COG), orthologous protein-coding genes, protein domains, families and functional sites, the assembled 3,62,502 transcripts were characterized using KEGG, InterPro, RefSeq, eggnog, OrthoDB, Pfam, PROSITE, PANTHER, and TIGRFAMs. We observed 245 KEGG pathways based on the characterization of 9,815 KEGG-annotated unigenes (Supplementary Table S9). The pathways observed most altered under HS, and in mutant were - ‘carbon metabolism’, ‘ribosome’, ‘starch and sucrose metabolism’, ‘biosynthesis of amino acids’, photosynthesis, basal transcription factor, DNA replication, and calcium signalling pathways.
Identification of heat-responsive TFs and HSPs in wheat mutant
The annotated data was characterized and searched for the identification of putative HSFs and HSPs in wheat mutant along with their fold expression analysis in parent (P3), mutant (M3) and HS-treated mutant (M3H). We identified 232 HSFs and HSPs showing upregulation in wheat mutant under HS, as compared to control (Supplementary Table S10; Table 4). Heat shock protein (TRINITY_DN8393_c6_g1) showing homology with wheat HSP17 (TraesCS4D01G212500.1) showed maximum log fold expression (17.89-fold) in mutant exposed to HS. The transcript was observed lying on chr4D. Similarly, we identified putative HSFs (TRINITY_DN15834_c0_g2) of 1077 bp, which showed 17.3-fold upregulation in M3H, as compared to control. The TF was observed lying on the chr7A. We also identified three different isoforms of chaperone clpB (TRINITY_DN25343_c0_g1, TRINITY_DN11203_c0_g1, and TRINITY_DN20098_c0_g2) showing 15.5-, 14.7- and 14.4-fold increase in the expression in HS-treated mutant, as compared to control. We observed 39 different isoforms of high molecular weight heat shock proteins especially heat shock 70 kDa protein (TRINITY_DN14477_c0_g1, TRINITY_DN31360_c0_g1, TRINITY_DN6139_c0_g2, TRINITY_DN12920_c1_g1, TRINITY_DN437_c1_g1, TRINITY_DN13029_c0_g2, etc.) showing 11.6-, 9.83-, 8.55-, 8.5-, 8.37-fold upregulation in mutant exposed to HS. Similarly, we identified 6 putative Chaperone proteins DnaJ (TRINITY_DN18040_c0_g3, TRINITY_DN7672_c0_g1, TRINITY_DN23961_c0_g2, TRINITY_DN54525_c0_g1, TRINITY_DN2230_c0_g1, TRINITY_DN2883_c0_g1, TRINITY_DN58791_c0_g1) with digital fold expression of 13.08-, 7.38-, 6.05-, 5.93-, 5.78-, 5.75- and 4.78-fold. Overall, the expression of HSFs and HSPs were observed higher in wheat mutant M3H. Mutant and parent showed similar pattern of expression of HSFs and HSPs. We observed ~179 putative transcripts of HSPs showing upregulation in wheat mutant, as compared to control. Maximum HSPs were observed lying on chr3D followed by chr3B and chr4B.
Identification of transcripts involved in photosynthesis
We identified 68 transcripts showing homology with genes involved in photosynthesis (Supplementary Table S11). The expression of Photosynthetic NDH sub complex B3 (thylakoid transit peptide) was observed maximum (10.34-fold) in wheat mutant exposed to HS, as compared to control. Similarly, the expression of gene coding for L-ribulose-5-phosphate 3-epimerase (a metalloprotein) and involved in interconversion between D-ribulose 5-phosphate and D-xylulose 5-phosphate was observed highly upregulated (10.27-fold) in M3H compared with control. We identified 34 transcripts showing homology with photosystem I and II, for example - Photosystem II 10 kDa polypeptide family protein (10-fold), Photosystem II reaction center PsbP family protein (9.73-fold), Photosystem II CP47 reaction centre protein (9.65-fold), Photosystem I P700 chlorophyll a apoprotein A2 (9.31-fold), Photosystem II 10 kDa polypeptide family protein (9.17-fold), etc. We identified 16 putative transcripts of RuBisCo activase involved in conversion of inactive form of RuBisCo into active RuBisCo for further carbon fixation. Ribulose bisphosphate carboxylase/oxygenase activase (TRINITY_DN30380_c1_g1) showed 9.61-fold upregulation in wheat mutant, as compared to mutant grown under ambient condition. Other isoforms of RuBisCo activase identified were - TRINITY_DN21223_c0_g1 (7.93-fold), TRINITY_DN30380_c2_g1 (7.58-fold), TRINITY_DN30380_c2_g1 (7.49-fold), TRINITY_DN13450_c0_g1 (7.30-fold), etc. Similarly, we identified 4 sucrose synthase transcripts - TRINITY_DN44228_c0_g1 (4.29-fold), TRINITY_DN50041_c0_g1 (4.17-fold), TRINITY_DN40865_c0_g1 (2.91-fold), and TRINITY_DN92104_c0_g1 (2.84-fold) in wheat mutant.
Mining of transcriptome data for the identification of transcription factors (TFs)
We identified 217 TFs involved in various biological process linked with HS-tolerance in wheat mutant (Supplementary Table S12). Heat shock transcription factor family protein (TRINITY_DN15834_c0_g2) showed maximum expression (log2 FC 17.3) in M3H, as compared to M3 grown under control condition. We identified 36 transcripts showing homology with heat-responsive transcription factors, for example - TRINITY_DN133181_c4_g1 (log2 FC 13.3), TRINITY_DN301_c0_g1 (log2 FC 12.9), TRINITY_DN26865_c1_g1 (log2 FC 11.4), TRINITY_DN139419_c0_g1 (log2 FC 10.6), TRINITY_DN37726_c2_g1 (log2 FC 9.95), etc. Basic leucine zipper (bZIP) gene family plays very important role in modulating the defense mechanism of plants against abiotic stress. Here, we identified 4 putative bZIP TFs showing significant upregulation in HS-treated wheat mutant i.e. TRINITY_DN38923_c1_g2 (log2 FC 4.56), TRINITY_DN12902_c0_g4 (log2 FC 3.57), TRINITY_DN12902_c0_g3 (log2 FC 3.44), TRINITY_DN2410_c0_g1 (log2 FC 2.59), TRINITY_DN20899_c0_g1 (log2 FC 2.5), as compared to mutant under ambient condition. Similarly, MYB transcription factors (TFs) families, is involved in hormone signal transduction, and abiotic stress tolerance, etc. We identified 14 MYB TFs in wheat mutant with significantly higher expression, as compared to parent, for example - TRINITY_DN84781_c0_g1 (log2 FC 10.34), TRINITY_DN164804_c0_g1 (log2 FC 8.75), TRINITY_DN72943_c0_g2 (log2 FC 6.05), TRINITY_DN21219_c1_g1 (log2 FC 5.79), etc. WRKY transcription factors consist of 60-70 amino acid WRKY protein domains with a conserved motif and a zinc-finger region. Here, we observed 8 putative WRKY TFs showing significant upregulation in wheat mutant under HS, i.e. TRINITY_DN69699_c0_g1 (log2 FC 4.08), TRINITY_DN60592_c0_g1 (log2 FC 2.93), TRINITY_DN35554_c2_g1 (log2 FC 2.50), TRINITY_DN33330_c0_g1 (log2 FC 2.43), TRINITY_DN44136_c1_g1 (log2 FC 2.36), TRINITY_DN3286_c2_g1 (log2 FC 2.32), TRINITY_DN89270_c0_g1 (log2 FC 2.15), and TRINITY_DN7487_c0_g1 (log2 FC 2.13). Overall, we observed significant upregulation of identified heat-responsive TFs in wheat mutant.
Expression of genes involved in stress signalling in wheat mutant
We identified 128 transcripts showing homology with signalling molecule linked with HS-tolerance in wheat mutant M3H (Supplementary table S13). The expression of signal peptide peptidase-like protein was observed maximum (log2 FC 12.54) in M3H compared with M3. Similarly, sarcoplasmic/endoplasmic reticulum calcium ATPase 3, Calcium-binding EF-hand family protein-like, etc. showed significant upregulation in wheat mutant under HS. CDPK or Calcium-binding EF-hand family protein provides the first line of defense by triggering the signalling mechanism of the plant against HS. We observed 10 transcripts showing homology with CDPK or CBP and has significantly higher expression in wheat mutant under HS, as compared to mutant under control, for example - TRINITY_DN14736_c0_g1 (log2 FC 9.89), TRINITY_DN14736_c0_g2 (log2 FC 5.26), TRINITY_DN130772_c0_g3 (log2 FC 4.64), TRINITY_DN70876_c1_g1 (log2 FC 4.21), etc. Maximum CDPK/ CBP were observed lying on chr2D. We identified 4 MAP3Ks, 1 MAP2Ks and 2 MAPKs in wheat mutant showing significant upregulation under HS, for example – MAPK (TRINITY_DN62203_c0_g1, TRINITY_DN102354_c2_g1), MAP2Ks (TRINITY_DN37965_c0_g2, TRINITY_DN2937_c0_g1) and MAP3Ks (TRINITY_DN27062_c1_g1, TRINITY_DN27062_c1_g1, TRINITY_DN2937_c0_g1, TRINITY_DN2937_c0_g2).
Validation of DEGs using quantitative Real-Time PCR
In order to validate the accuracy and reproducibility of the RNA-Seq data, we have randomly selected 20 DEGs (10 up-regulated genes and 10 down-regulated genes in the dataset) and further, they were analysed for the expression using qRT-PCR (Fig. 7). The upregulated genes selected for the expression analysis were – calcium dependent protein kinase (CDPK), heat-shock protein (HSP70), small heat shock protein 17 (HSP17), heat-shock protein 101 (HSP101), heat shock factor (HSF), serine threonine kinase (STK), peptidyl-prolyl isomerases (PPIs), superoxide dismutase (SOD), oxygen evolving enhancer protein (OEEP), and ABC transporter (ABC-Tra). The expression analysis of these transcripts in mutant, as compared to parent showed maximum relative fold expression of ABC-Transporter (6.5-fold) followed by serine threonine kinases (STKs; 4.28-fold) (Fig. 7a). HSP17 showed significant increase in the expression (3.74-fold) in M3, as compared to P3. Similarly, SOD showed 3.9-fold increase in the expression in M3 compared with P3.
Expression analysis in mutant exposed to HS showed maximum expression of HSP17 (15.27-fold) followed by ABC-transporter (ABC-Tra; 7.81-fold), as compared to mutant under ambient condition (Fig. 7b). Expression analysis of signalling molecule showed maximum expression of CDPK (5.9-fold) in mutant exposed to HS. Similarly, we observed significant increase in the expression of HSPs and antioxidant enzymes in HS-treated mutant, as compared to control. Oxygen evolving enhancer protein (OEEP), being important enzyme involved in carbon assimilation showed 2.46-fold increase in the expression in mutant exposed to HS compared to control.
We also validated the differential expression in mutant exposed to HS (M3H) compared with parent (P3). The relative fold expression of small HSP-17 was observed maximum (19.2-fold), as compared to parent (Fig. 7c). Similarly, ABC-tra showed 9.4-fold increase in the expression in mutant exposed to HS, as compared to parent. Expression analysis of signalling molecule in mutant exposed to HS showed maximum expression of CDPK (5.27-fold) followed by STK (3.73-fold). Similarly, superoxide dismutase (SOD) showed 4.48-fold and OEEP 2.96-fold increase in the expression in HS-treated mutant, as compared to parents. Overall, HSP17, ABC-tra and CDPK showed significantly higher expression in mutant, as compared to parent and might be playing decisive role in modulating the thermotolerance of the mutant.
The downregulated genes selected for the expression analysis were – pyruvate decarboxylase (PDC), sucrose phosphate synthase (SPS), sucrose synthase (Suc Syn), sucrose transporter (SuT), soluble starch synthase (SSS), rubisco (small subunit; Rub-S), β-amylase (β-Amy), fructose bis-phosphate aldolase (FBA), debranching enzyme (DBE), ADP-gluco pyrophosphorylase (AGPase) (Fig. 8). The relative fold expression of AGPase was observed minimum (-0.58-fold) followed by SuT (-0.62-fold) in mutant, as compared to parent. RuBisCo (SSU) showed drastic downregulation in mutant (-0.77-fold), as compared to parents (Fig. 8a). Other genes involved in carbon assimilation especially SPS, Suc Syn, etc. showed downregulation in mutant; though non-significant variations was observed, as compared to parent.
Expression analysis of M3H showed maximum downregulation of Suc Syn (-0.33-fold) followed by AGPase (-0.41-fold) and SPS (-0.44-fold), as compared to M3 under control condition (Fig. 8b). The expression of sucrose transporter (SUT) was also observed downregulated (-0.61-fold) in M3H, which perturbed the transportation of photosynthates from the source (leaves) to the sink (endosperm). Even some of the genes linked with starch biosynthesis pathway like SSS and DBE showed downregulation in M3H, as compared to M3 under control condition. The gene of RuBisCo - most abundant enzyme – showed downregulation (-0.77-fold) in M3H.
Expression analysis of DEGs in M3 showed maximum downregulation of Suc Syn (-4.0-fold) followed by SPS (-0.46-fold), as compared to P3. Similarly, the genes associated with carbon assimilation like RuBisCo, SPS, AGPase, and SSS showed significant downregulation in M3, as compared to P3. Debranching enzyme, the product of DBE showed downregulation (-0.75-fold) in M3, as compared to P3.
Expression analysis in M3H showed significant downregulation of sucrose synthase (-0.33-fold) followed by sucrose phosphate synthase (-0.44-fold) and AGPase (-0.41-fold), as compared to M3 grown under control condition. The transcripts especially RuBisCo, sucrose transporter, linked with photosynthesis and SSS, DBE, etc. associated with starch biosynthesis pathway also showed significant downregulation in HS-treated mutant (M3H), as compared to control. Even catabolic enzyme β-amylase which acts on starch compromising its quality showed downregulation (-0.75-fold) in response to M3H. Expression analysis in M3H, as compared to P3 showed maximum downregulation in sucrose synthase (-0.4-fold) followed by sucrose phosphate synthase (-0.46-fold) and β-amylase (-0.62-fold), as compared to P3 (Fig. 8c). Even other important genes like pyruvate decarboxylase, RuBisCo and SSS showed downregulation under HS in mutant compared with parent.
Validation of RNA-seq data in wheat mutant under HS using qPCR
We have randomly selected 10 upregulated and 10 downregulated transcripts from the data-set for the validation of RNA-Seq data using the qPCR. The expression was analysed under ambient and HS-treated conditions. We established a positive correlation between the fold-change values (HS/control) of RNA-Seq and qPCR data (r2=0.90) (Fig. 9). This confirms the reliability of the RNA-Seq data generated and used in present investigation.
Simple Sequence Repeat (SSRs) identified in wheat mutant
We have used the assembled de novo transcriptome data for the identification of putative Simple Sequence Repeats (SSRs) using Perl scripts of MISA (MIcroSAtellite identification tool). We identified total of 3388 SSRs in Parent (P3), 4674 SSRs in mutant (M3) and 33878 SSRs in HS-treated mutant (M3H). A summary of different markers identified in parent, mutant and HS-treated mutant has been presented in Table 5 and supplementary file 14. We observed abundance of SSRs of trinucleotide type in parent (P3 – 1740), mutant (P3 -2345) and HS-treated mutant (P3 – 13140).
Variants detected in the DETs involved in modulating thermotolerance in wheat mutant
The assembled de novo transcriptome data were used for the prediction of variants using SAMTOOLS (quality filtered at read depth of 5 and 10; Phred-score ≥30). The number of SNPs and InDels identified at RD5 (Read Depth5) and RD10 (Read Depth10) in the P3, M3 and M3H are presented in Table 6 (supplementary file S15). The assembled transcripts were randomly mapped on to the 21 chromosomes of wheat. SNP calling predicted a total of 1,10,772 intervarietal SNPs at read depth 5 (RD5). We observed significant variations in the coding sequences of identified DETs in mutant over parent. Variant calling were analysed in the transcript sequences by comparing the DETs identified in mutant and parent. We observed potential SNPs in some of the keys DETs such as - Heat shock protein (TRINITY_DN49411_c0_g2_i2) lying on chr3B showed transversion of CC in parent to CG in mutant. Similarly, SOD (TRINITY_DN5409_c0_g1_i1) used as one of the potential markers for evaluating wheat germplasm for thermotolerance showed A/T allelic transversion in mutant. Serine threonine kinase (TRINITY_DN88200_c0_g3_i1) localized on chr3B is one of the important signalling molecules involved in triggering the defense response of the plants under HS showed C/G allelic transversion. ClpB chaperone – a part of stress induced chaperone system showed transversion of T/G allele. Myb TF, which is involved in hormonal signal transduction and abiotic stress tolerance in plants, showed transition of C/T allele in mutant over parent. bZIP Transcription Factor in mutant. The Basic Leucine Zipper Domain (bZIP domain) protein (TRINITY_DN29470_c0_g1_i3) localized on chr3D has been reported to modulate the tolerance of the plant to the maximum extent, and showed transition of T/C allele in wheat mutant. Universal Stress Protein (TRINITY_DN26896_c0_g1_i3) lying on chr2A showed transition of A/G allele in wheat mutant and has been reported to involved in modulating the stress tolerance of the plants. Similarly, 70 kDa HSP (TRINITY_DN865_c0_g2_i2) localized on chr5D showed transition of C/T allele in mutant. Similarly, Peptidyl prolyl isomerase (TRINITY_DN77799_c0_g1_i1) showed conversion of A/G allele in mutant compared to parent.
In this study, we identified potential gene based InDel markers lying on differentially expressed transcripts identified in wheat mutant and parent. We identified a total of 2432 InDels at read depth of 5 and 1602 InDels at read depth of 10 (supplementary file S15). The InDels identified in parent were 798 (RD5) and 525 (RD10), mutants were 841 (RD5) and 548 (RD10) and mutant exposed to HS were 793 (RD5) and 529 (RD10). Most of the InDels identified were observed lying on predicted genes, signalling molecules, chaperones and genes linked with antioxidant enzymes. For example, protein serine/threonine phosphatase (TRINITY_DN98416_c0_g1_i3) lying on chr3A showed the presence of InDels of alleles in mutant (GTCA/GTCA) and parent (G/G). Similarly, in starch synthase (TRINITY_DN7272_c0_g1_i3), we observed the presence of InDels TGCGC/TGCGC in mutant and TGC/TGC in parent. Being key enzyme of starch biosynthesis pathway, this InDel of SSS was observed lying on chr3A. The protein kinase gene (TRINITY_DN11281_c1_g1_i6) involved in regulating the activity of stress-associated proteins and lying on chr3B showed the presence of InDel GCAGAGTCAG/GCAGAGTCAG in wheat mutant and GCAGAGTCAG/GCAG in parent. Similarly, many InDels were observed lying on serine threonine kinase gene (STKs) which plays very important role as signalling molecule by triggering the defence mechanism against abiotic stress. For example, STK gene (TRINITY_DN12482_c0_g1_i3) lying on chr4B showed the presence of InDels alleles TTGCTGCTGC/TTGCTGCTGC in wheat mutant and TTGCTGCTGC/TTGCTGCTGC in parent.
Analysing the HS-tolerance level and grain-quality of wheat mutant
Characterizing the thermotolerance capacity of wheat mutant
The mutant [parent-MP3054- C-306/CB.SPRING BW/CPAN2072 (Parentage)] developed where further analysed for their HS-tolerance level and grain quality using different biochemical markers. The parent and mutant were grown under ambient condition [22±3°C (day time), 18±3°C (night time)] and HS-treated condition [38±3°C (day time), 30±3°C (night time)]. The flag leaf was used for analysing the tolerance level and the harvested grains were used for the quality analysis. SOD activity profiling in parent and mutant during pollination stage showed maximum activity in HS-treated parent (7.1 U/mg proteins) and minimum in mutant (4.25 U/mg Proteins) under control condition (Fig. 10a). During milky-ripe stage, the SOD activity was observed maximum in HS-treated mutant (6.95 U/mg proteins) and minimum in parent (2.85 U/mg proteins) under control condition. Further, we observed decrease in the SOD activity in both parent and mutant with increase in durations.
In case of guaiacol peroxidase, we observed maximum activity in HS-treated mutant (M3H) during pollination, milky-ripe and mealy-ripe stages of growth; maximum (63.5 U/mg proteins) was observed during mealy-ripe stage (Fig. 10b). The activity of GPX was observed minimum (10.8 U/mg proteins) in parent under control condition throughout the different stages of growth. Total antioxidant potential is recommended as one of the important biochemical markers for analysing the tolerance level of the plant under stress. Here, we observed significant increase in the TAC in parent (C and HS-treated) from pollination to milky-ripe stage, and further decrease in the TAC was observed during mealy-ripe stage. In case of mutant, we observed decrease in the TAC from pollination to milky-ripe stage, and further significant increase in the TAC was observed during mealy-ripe stage (Fig. 10c). The maximum TAC (16.1 mM/g FW) was observed in HS-treated mutant (M3H) during mealy-ripe stage.
The accumulation pattern of hydrogen peroxide inside the leaf tissue of parent (P3), mutant (M3) and both exposed to HS were visually analysed though NBT assay test (Fig. S2). In case of parent (P3) under control condition, we couldn’t observe any dark patch/spots on the bleached leaves, whereas prominent dark spots were observed in case of parent exposed to HS. Similarly, in case of mutant under control condition, the leaf was slightly darkish showing the minor accumulation of hydrogen peroxide. When the mutant was exposed to HS (M3H), intense dark blotches appeared in the leaves showing the accumulation of H2O2. Overall, the accumulation was observed maximum in the leaves of mutant exposed to HS.
Characterizing the grain-quality of wheat mutant and parent
The harvested grains were used for analysing the activities of starch synthesizing enzymes (AGPase, SSS), starch degrading enzymes (Amylase) and variations in starch composition under control and HS-treated conditions. In case of wheat parent under control and HS-treated condition, we observed significant (p<0.05) decrease in the AGPase activity in response to HS, though non-significant (p<0.05) variations were observed in the SSS and amylases activities (Fig. 11a). Similar pattern of decrease in the activities of AGPase, SSS and amylase was observed in case of mutant under HS; though percent decrease was observed less, as compared to parent. The variations in the SSS and amylase activity was observed non-significant (p<0.05) in mutant under HS. Overall, AGPase, SSS and amylase present in mutant were observed more thermostable under HS, as compared to parent.
We also characterize the effect of HS on starch composition, especially starch, amylose and amylopectin content in harvested grains of parent and mutant. The starch, amylose and amylopectin content was 78.2, 9.8 and 60.2 % in grains of parent under control condition, which further decreases to 55.1, 7.8, and 51.2 % under HS-treated conditions (Fig. 11b). Grain-quality analysis in mutant showed 75.1, 12.5, and 66 % starch amylose and amylopectin under control condition, which further decreases to 72.4, 10.1, and 58.2 % under HS-treated conditions. In case of parent, we observed ~30% decrease in the starch content under HS, whereas only 4% decrease was observed in mutant. Though, percent decrease was observed non-significant in case of amylose in both parent and mutant. Similarly, we observed 18% decrease in the amylopectin in parent under HS; though percent decrease was observed less (12%) in wheat mutant.