Insight into the mechanisms of terminal HS-tolerance in wheat mutant with improved nutritional quality through de novo transcriptome sequencing

Terminal HS is one of the main bottle-neck in wheat yield and grain-quality. Here, we have developed wheat mutant (M 3 ) for HS-tolerance [parent-MP3054- C-306/CB.SPRING BW/CPAN2072 (Parentage)]. To elucidate the mechanism of thermotolerance in mutant, we performed de novo transcriptomic sequencing of mutant (M 3 ), parent (P 3 ), and mutant exposed to HS (M 3 H). We sequenced 6.5, 7.5, and 7.0 million reads in P 3 , M 3 and M 3 H and generated 3,05,537 genes and 5,88,788 transcripts with an N50 of 1,349 bp. We observed 6,120 upregulated and 4,428 downregulated transcripts (M 3 vs P 3 ), 11,354 upregulated and 12,408 downregulated genes (M 3 H vs P 3 ) and 4817 upregulated and 9085 downregulated genes (M 3 H vs M 3 ). Some of the highly upregulated genes observed were HSP20, SOD, ABC transporters, HSF, etc. and downregulated genes were starch synthase, sucrose synthase, etc. Gene Ontology analysis showed ‘ATP-binding’ to be most enriched category. Carbon metabolism pathway was observed most altered under HS. We identied 41940 SSRs, 1,10,772 SNPs and 2432 InDels. Potential markers were observed lying on HSP, SOD, STK, and starch synthase. Biochemical markers based characterization showed wheat mutant to be better in HS-tolerance and grain-quality, as compared to parent. 3 ), mutant (M 3 ) and mutant exposed to HS (M 3 H) were used for the de novo assembly using Trinity v 2.8.4 (http://trinityrnaseq.github.io) 22 . The contigs were clustered, and assembled using Trinity as per the paired-end information. All the assembled transcripts were searched against NCBI nr (NCBI non-redundant protein database). BLAST2GO (https://www.blast2go.com/) program was used for the Gene Ontology (GO) annotations of the assembled sequences 23 . and ADP-glucopyrophosphorylase (AGPase-LSU). Primers for qPCR analysis were designed manually using the Genesher2 primer designing software (supplementary Table S1). Trizol method was used for the total RNA isolation. expression (17.89-fold) in mutant exposed to HS. The transcript was observed lying on chr4D. Similarly, we identied putative HSFs (TRINITY_DN15834_c0_g2) of 1077 bp, which showed 17.3-fold upregulation in M 3 H, as compared to control. The TF was observed lying on the chr7A. We also identied 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 identied 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-,

Introduction modulating the tolerance level against various stresses. Some of the crops, wherein mutation breeding approach has been used to modulate different traits linked with tolerance are -sorghum 12 , Barley 13 , rice 14 , wheat 15 , etc. Different approaches has been used in the past to develop wheat mutants like using EMS or gamma irradiation in order to develop lines with very high tolerance level. The selected mutants were further used for the back cross breeding program in order to develop genotypes with high tolerance level 16 . The generation of desirable donor lines using the approach of mutation has helped the breeder in the development of climate smart crop in the past.
De novo sequencing has potential to identify the SAGs, single nucleotide polymorphism and their localization in the genes responsible for modulating the thermotolerance level in the crop plants 17 . The variations in the expression of stress-associated genes in response to varied environmental stresses modulate the tolerance of the plants adapting to adverse environmental conditions 15 . Heat stress causes abrupt increase in the expression of heat-responsive transcription factors and SAGs and is important for the adaptation of plant against HS 18 . Regulatory and structural genes showed signi cant variations in the expression under heat stress. The over-expression in response to HS causes accumulation of stress-associated proteins, chaperones such as heat shock proteins (HSPs), and various metabolites 19 . Although, identi cation and annotation of SAGs has been done in exhaustive manner in different crop species, still the number is not su cient to completely elucidate the mechanism underlying tolerance. Furthermore, isoforms of SAGs are reported to involve in different pathways showing temporal and spatial expression, and to be identi ed and characterized in order to understand the actual behaviour of different SAGs and SAPs involved in tolerance against different environmental stresses.
De novo transcriptome sequencing is best method for the identi cation and functional validation of genes and molecular mechanisms underlying different biological processes. It provides a robust option for creating an overview of transcript pro les in different samples under varied treatments in order to understand the mechanism and pathways involved.
In present investigation, a wheat mutant (M4) for HS-tolerance [parent-MP3054-C-306/CB.SPRING BW/CPAN2072 (Parentage)] was developed by Department of Atomic Energy (DAE), Bhabha Atomic Research Centre (BARC), Mumbai. The mutants were further, characterized for their HS-tolerance using RNA-seq in order to identify the possible SAGs, SSRs, SNPs and InDels responsible for modulating the HS-tolerance level of the wheat mutant with improved grain-quality. Even the parents were characterized along with the mutants for the comparative view.

Plant materials and treatments
A set of wheat mutants [parent-MP3054-C-306/CB.SPRING BW/CPAN2072 (Parentage)] were created by Department of Atomic Energy, Bhabha Atomic Research Center, Mumbai using the gamma irradiation (dose of 300 Gy) for heat stress-tolerance. The mutant population were screened at different location for thermotolerance and further the seeds of mutants (M4) were used for the molecular and biochemical characterization in order to elucidate the mechanism underlying HS-tolerance. The mutant seeds (in three lots) received from BARC, Mumbai were sown in separate pots (12'x12') inside regulated chambers at National Phytrotron Facility, IARI, New Delhi. The seeds were sterilized with 0.1% (w/v) sodium hypochlorite followed by washing with distilled water, and further germination was carried out on moistened sand for 20 days at 22±3 °C. The seedlings with uniform growth were transferred to pot lled with perlite and vermiculite in equal ratio with ve seedlings per pot. The pots were kept inside regulated growth chamber under a 16/8 h (day/night) photoperiod, with temperature of 22±3 °C (day-time) and 18±3°C (night-time), and 75% relative humidity. One group was exposed to HS [38±3°C day time and 28±3°C night time] in a sinusoidal mode starting from post-anthesis till harvesting. The samples for the transcriptome were collected during mealy-ripe stage (as per the Feekes scale) 20 . The samples were collected in triplicates and were further stored at -80°C for further molecular and biochemical research works. The experimental research and eld studies on wheat mutant complies with the rules and regulations of Indian Agricultural Research Institute (IARI), New Delhi.
De novo transcriptome sequencing Total RNA isolation and quality analysis The total RNA was isolated from the collected samples (root, stem, leaf and spikes) by RNeasy Plant Mini Kit (Qiagen, UK). The quality analysis was performed using the Nanodrop 2000c (Thermo sher scienti c, USA), and Agilent Tapestation. The total RNA having RIN value of more than 7.5 was used for the library preparation. High-quality RNA samples (RIN ≥ 6.5) were used for the construction of sequencing library after treating it with DNase I (TaKara) in order to remove the DNA contamination. Two biological replicates were used for each experiment.
Construction of RNA-Seq library and Illumina HiSeq4000 sequencing Transcriptome libraries for sequencing were constructed using the Illumina TruSeq stranded mRNA library construction kit as per the instructions given by the manufacturers (Part # 15008136; Rev. A, Illumina, USA). Brie y, 1 μg of pooled total RNA from the collected samples were used for the isolation of mRNA using oligo dT beads (TruSeq RNA Sample Preparation Kit, Illumina, USA). The puri ed mRNA was fragmented and reverse transcribed for the synthesis of rst strand using random hexamers and Superscript II Reverse transcriptase (Invitrogen, USA). Further, DNA polymerase-I was used for the second strand cDNA synthesis.
The cleaned cDNA was ligated with Illumina adapters after end repair. In order to enrich the adapter ligated fragments, the library generated was ampli ed using PCR. Quality analysis was done using Nanodrop and Bioanalyzer Chip (Agilent Technologies, California, USA). The denatured DNA from the libraries was used for the sequencing using Illumina HiSeq 4000, through the sequencing by synthesis method. The raw sequencing data in Fastq format were retrieved and subjected to quality analysis using SeqQC-V2.0 program (NGS data QC).

De novo assembly and sequence annotation
Pre-processing was carried out by removing the low-quality reads (Q-value < 20), adapter reads and ambiguous reads (>10% 'N' bases) in order to get the clean reads 21 . The clean reads obtained from the parent (P 3 ), mutant (M 3 ) and mutant exposed to HS (M 3 H) were used for the de novo assembly using Trinity v 2.8.4 (http://trinityrnaseq.github.io) 22 . The contigs were clustered, and assembled using Trinity as per the paired-end information. All the assembled transcripts were searched against NCBI nr (NCBI non-redundant protein database). BLAST2GO (https://www.blast2go.com/) program was used for the Gene Ontology (GO) annotations of the assembled sequences 23 .

GO enrichment analysis and biological pathway characterization
In order to characterise the involvement of assembled transcripts into different pathways in Kyoto Encyclopaedia of Genes and Genomes (KEGG), we have used KEGG Automatic Annotation Server (http://www.genome.jp/tools/kaas/). GO enrichment analyses of DEGs were performed using Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/expression.php) with GO annotation.

Expression analysis of differentially expressed transcripts (DETs)
The expression level of transcripts was analysed using the FPKM (fragments per kilobase of exon per million mapped fragments) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used for analysing the abundance of gene and isoforms. Comparative differential expression analysis was done based on the count of the transcripts in P3, M3 and M3H using EdgeR v 3.16.5 (Empirical analysis of Digital Gene Expression in R; http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) software (FDR < 0.05; log 2 FC ≥ 1).
The PCR program and calculation was followed as mentioned in our earlier publication 7 . The relative fold expression was calculated using the Pfa method 27 .
Biochemical characterization of wheat mutant for thermotolerance Antioxidant enzyme activity assay The samples (P 3 , M 3 and M 3 H) were used for the SOD activity assay using the protocol of Beauchamp and Fridovich 28 . SOD activity was measured by taking absorbance at 560 nm against the reagent blank (phosphate buffer). SOD activity was calculated using the formula: SOD measured = [Control-(Sample-Blank)]/ Control *100; 50% inhibition in O.D. was considered as 1 unit of SOD activity.

Total antioxidant capacity (TAC) analysis
The method of Benzie and Strain 29 was used for estimating the total antioxidant potential of collected samples. The absorbance of ferrous coloured form produced from the reduction of ferric-tripyridyltriazine complex was used for calculation of TAC (mM Fe/ g fresh weight).

Accumulation pattern of H 2 O 2
The accumulation of H 2 O 2 was estimated following the method of Loreto and Velikova 30 . Leaf samples (0.1 gm) were crushed in liquid N 2 and further 1 mL of Trichloroacetic acid (1%) was added and further centrifuged at 12000 rpm (4°C) for 15 min.
Further, a reaction mixture was prepared by adding 0.75 mL of supernatant with 0.75 mL of 10 mM phosphate buffer (pH 7.0) and 1.5 mL of potassium iodide (KI) and the OD was taken at 390 nm. The H 2 O 2 concentration was estimated using the extinction coe cient of H 2 O 2 (26.6 mM -1 cm-1 ).

In vivo detection of H 2 O 2 inside leaves
The in situ detection of hydrogen peroxide in the leaf was visualized following the protocol, as mentioned by Daudi et al. 31 . In brief, 10 mM Na 2 HPO 4 DAB staining solution was prepared by adding DAB (1 mg mL -1 ), in a separate test-tube with 2.5 mL of 200 mM Na 2 HPO 4 and 25 μL of 0.05% (v/v) Tween 20. We have selected very premature leaves (3 leaves per plant) from parent (P3), mutant (M3) and mutant exposed to HS (M3H) and was treated with 2 mL DAB staining solution in petriplate without exposing to light. For control treatment, we have used 2 mL of 10 mM Na 2 HPO 4. The petriplate was kept on standard laboratory shaker for 4-5 h at 80-100 rpm shaking speed. Further, bleaching solution (ethanol: acetic acid: glycerol = 3:1:1) was added by dispensing the DAB staining solution and was boiled (~90-95°C) for 15 min and the leaf was directly visualized for DAB staining.

Analysis of biochemical traits linked with grain quality
Activity assay of ADP-glucose pyrophosphorylase (AGPase), starch synthase and amylase enzymes The activity assay of AGPase was estimated using the method of Nivelle et al. 32 and Nakamura et al. 33 . The developing seeds were grind to ne powder using liq. nitrogen and 100 mg of ground powder was used for the analysis. The conversion of 1 µM of ADP-glucose into glucose-1-phosphate per min was expressed as one unit enzyme. The activity of SSS was assayed as per the protocol of Nakamura et al. 33 . The integrity assay of starch was carried out by analyzing the amylolytic activity in harvested seeds of P 3 , M 3 and M 3 H. The total amylolytic activity in the grains was estimated as per the modi ed protocol of Briggs 34 . In brief, 0.5 g of ne powder was resuspended in 2 mL of ice-cold calcium acetate buffer with pH 6.0 and supernatant was separated through centrifugation. Further, 100 μL of enzyme extract was added to reaction mixture containing 1 mL of substrate starch solution (1% w/v in reaction buffer) and 1 mL of calcium acetate buffer (pH 4.8) and was incubated at 24 °C for 10 min.
The absorbance was read at 620 nm against blank and standard curve of starch was used for the calculation of total amylase activity.

Estimation of Starch Content
The total starch in fractionated sample was estimated as per the protocol of Thayumanavan and Sadasivam 35 . In brief, the pellet extracted from 0.5 g sample was mixed with 6.5 mL of perchloric acid (52%) and 5.0 mL water followed by incubation at 0°C for 20 min. The supernatant obtained after centrifugation was pooled and the volume was make up to 100 mL. A reaction mixture (RM) was prepared by adding 0.1 mL aliquot, 4 mL anthrone reagent and the nal volume was made up to 5 mL followed by heating in a boiling water bath for 8 min. The absorbance was read at 630 nm and glucose standard was used for the estimation of starch by multiplying the value by a factor of 0.9.

Statistical analysis
For de novo transcriptome sequencing, we have used the pooled samples in duplicates, and for molecular and biochemical estimations, we have used three biological replicates. All data analyses were conducted using SPSS Statistics 20.

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 (P 3 ), 80444546 (M 3 ), and 73850514 (M 3 H) with more than 92% Q 30 bases. The GC content was ~48% and produced comparable data from the biological replicates. De novo assembly of ltered 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 Annotation and classi cation 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. Brie y, 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 signi cant 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,

Gene ontology (GO) annotation
The functions of the assembled transcripts were classi ed via GO analysis. The GO terms for all assembled transcripts were extracted wherever possible and BLAST2GO was used for the classi cation of groups into 3 major categories (supplementary Table S5). The total number of GO terms identi ed in biological processes, cellular component and molecular function processes are -3233, 761 and 1963.

DEGs in wheat mutant (M 3 ) under HS-treated condition
The expressed transcripts in the P 3 , M 3 , and M 3 H were identi ed 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 P 3 , M 3 , and M 3 H 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 P 3 , M 3 and M 3 H samples. The number of assembled transcripts considered for differential gene expression analysis was 588,788. DETs were ltered at log 2 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 cu inks package. Log 2 fold change cut-off 2 (p<=0.01 and 0.05) were separately used as cut-off for upregulated and downregulated genes and isoforms.  Table   S8). Some of the identi ed 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)

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 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 identi ed putative HSFs (TRINITY_DN15834_c0_g2) of 1077 bp, which showed 17.3-fold upregulation in M 3 H, as compared to control. The TF was observed lying on the chr7A. We also identi ed three different isoforms of chaperone clpB (TRINITY_DN25343_c0_g1,  TRINITY_DN11203_c0_g1, and TRINITY_DN20098_c0_g2) 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 identi ed 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

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 signi cant increase in the expression (3.74-fold) in M 3 , as compared to P 3 . Similarly, SOD showed 3.9-fold increase in the expression in M 3 compared with P 3 .
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 signi cant 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 (M 3 H) compared with parent (P 3 ). 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.4fold 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 signi cantly 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-signi cant variations was observed, as compared to parent. Expression analysis of M 3 H showed maximum downregulation of Suc Syn (-0.33-fold) followed by AGPase (-0.41-fold) and SPS (-0.44-fold), as compared to M 3 under control condition (Fig. 8b) 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 (r 2 =0.90) (Fig. 9). This con rms the reliability of the RNA-Seq data generated and used in present investigation.

Simple Sequence Repeat (SSRs) identi ed in wheat mutant
We have used the assembled de novo transcriptome data for the identi cation of putative Simple Sequence Repeats (SSRs) using Perl scripts of MISA (MIcroSAtellite identi cation tool). We identi ed total of 3388 SSRs in Parent (P 3 Table 5 and supplementary le 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 ltered at read depth of 5 and 10; Phred-score ≥30). The number of SNPs and InDels identi ed at RD5 (Read Depth5) and RD10 (Read Depth10) in the P 3 , M 3 and M 3 H are presented in Table 6 (supplementary le 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). In this study, we identi ed potential gene based InDel markers lying on differentially expressed transcripts identi ed in wheat mutant and parent. We identi ed a total of 2432 InDels at read depth of 5 and 1602 InDels at read depth of 10 (supplementary le S15). The InDels identi ed in parent were 798 (RD5) and 525 (RD10), mutants were 841 (RD5) and 548 (RD10) and mutant  (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 (M 3 H) 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 signi cant 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 signi cant 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 (M 3 H) during mealy-ripe stage.

Page 14/24
The accumulation pattern of hydrogen peroxide inside the leaf tissue of parent (P 3 ), mutant (M 3 ) and both exposed to HS were visually analysed though NBT assay test (Fig. S2). In case of parent (P 3 ) 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 (M 3 H), intense dark blotches appeared in the leaves showing the accumulation of H 2 O 2 .
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 signi cant (p<0.05) decrease in the AGPase activity in response to HS, though non-signi cant (p<0.05) variations were observed in the SSS and amylases activities (Fig. 11a) In present investigation, some of the highly upregulated genes identi ed were HSP20, SOD, HSF, PDI, ABC transporters, etc. Heat shock proteins plays dual role of protecting the nascent proteins from denaturation under HS, and also helps in protein folding in order to make the protein biologically active 42 . SOD has been observed to be at the pivotal position in antioxidant enzymes networking system neutralizing the harmful effect of free oxygen radicals and oxidative burst under HS 43 . Similarly, ATP-binding cassette transporters (ABC transporters) have been observed to regulate grain formation, other than ATP-powered translocation of many substrates across the membranes 44 . Nirmal et al. 45  In present investigation, the most enriched GO category observed was 'ATP binding' (35127), followed by 'nucleic acid binding' (15154), 'protein kinase activity', and 'serine threonine kinase activity'. Similarly, among the pathways -carbon metabolism, followed by ribosome, starch and sucrose metabolism, biosynthesis of amino acids, and photosynthesis was observed most altered. Zhang et al. 47