Three rice genotypes Moroberekan, IR64 and N22 were used in this study. The details of the growth conditions followed in the analysis of Moroberekan, IR64 and N22 rice genotypes, as well as the relationship of the pollen grain stages on the progression to anthesis in rice to the stress regimes followed in this study, are shown in Figure 1 . Several independent sets of plants were transferred to growth chambers for the HS treatment to obtain plants with a variable time of flowering. For the ST_HS, anthers were harvested at ‘just before anthesis’ stage from plants subjected to one cycle of 42°C (4 h) HS ( Figure 1A ). For the LT_HS, anthers were harvested at a similar growth stage from plants subjected to 6 days of consecutive 38°C (6 h) HS ( Figure 1A ). In LT_HS treatment, the temperature regime followed was 28°C at 6 AM, gradually increased to 38°C by 9 AM, maintained at 38°C till 3 PM and gradually lowered to 28°C and maintained overnight. In the unstressed control (CON) condition, the temperature followed was: 28°C at 6 AM, gradually increased to 30°C and maintained till 6 PM and gradually brought down to 28°C overnight. Plants (15 to 35) were harvested (for each treatment), and flowering spikelets were collected on ice. The spikelets were dissected to separate the pollinated pistils and anthers. The anthers were harvested from only the middle portion of the whole spikelet to maintain uniformity in the development stage. Dissected anthers were collected in tubes containing liquid nitrogen and stored at –80°C until further use ( Figures 1B, C ).

Harvested anthers (~20 mg) were lysed and homogenized in TRIzol™ Reagent (Invitrogen™) to extract proteins for Western blotting. Extracted proteins were dissolved in 0.25% SDS and protein was quantified using Bradford reagent. For Western blotting, 10 μg of total soluble protein was loaded on the SDS-gel. Proteins were transferred to the nitrocellulose membrane (Amersham, UK). Hsp101 protein was probed using polyclonal anti-AtHsp101/ClpB primary antibodies (#AS08283, Agrisera, Sweden) and horseradish peroxidase (HRP) -conjugated anti-rabbit secondary antibodies (Sigma-Aldrich, USA), both at 1:10000 dilutions. The same membrane was probed with anti-tubulin alpha chain antibodies (# AS10680, Agrisera) for loading control. The blots were developed by enhanced chemiluminescence (ECL).

Proteins were extracted and quantified using a shotgun proteomics approach (Chaturvedi et al., 2013; Valledor and Weckwerth, 2014; Chaturvedi et al., 2015). Anther samples were freeze-dried and ground for 2 min in a shaking mill using steel balls (2 mm diameter). The homogenized samples were resuspended in 200 μL of protein extraction buffer (100 mM Tris- HCl, pH 8.0; 5% SDS, 10% glycerol; 10 mM DTT; 1% plant protease inhibitor cocktail (Sigma P9599) and incubated at room temperature for 5 min followed by incubation for 2.5 min at 95°C and centrifugation at 21000 × g for 5 min at room temperature. The supernatant was carefully transferred to a fresh tube. Two hundred microliter of 1.4 M sucrose was added to the supernatant. Proteins were extracted twice with 200 μL TE buffer-equilibrated phenol followed by counter extraction with 400 μL of 0.7 M sucrose. Phenol phases were combined and mixed with 2.5 volumes of 0.1 M ammonium acetate in methanol to precipitate proteins. After 16 h of incubation at -20°C, samples were centrifuged for 5 min at 5000 × g. The pellet was washed twice with 0.1 M ammonium acetate, once with acetone and air-dried at room temperature. The pellet was re-dissolved in 6 M urea and 5% SDS and protein concentration were determined using the bicinchoninic acid assay. Proteins were pre-fractionated by SDS-polyacrylamide gel electrophoresis. Forty micrograms of total protein were loaded onto a gel and run for 1.5 cm. Gels were fixed and stained with methanol: acetic acid: water: Coomassie Brilliant Blue R-250 (40:10:50:0.001). Gels were destained in methanol: water (40:60), and then each lane was divided into two fractions (Ghatak et al., 2021).

Gel pieces were destained, equilibrated and digested with trypsin (using Trypsin Sequencing Grade from Roche [11418475001]), desalted employing Bond-Elute C-18 SPEC plate (Agilent Technologies, Santa Clara, CA, USA) and concentrated in a Speed Vac concentrator (SCANVAC Cool Safe 110-4, Speed Vacuum concentrator, Labogene). Prior to mass spectrometric measurement, the tryptic peptide pellets were dissolved in 4% (v/v) acetonitrile, 0.1% (v/v) formic acid. One µg of the digested peptide (2 biological replicates from each condition) was loaded on a C18 reverse-phase column (Thermo Scientific, EASY‐Spray 500 mm, 2 µm particle size). The separation was achieved with a 90 min gradient from 98% solution A (0.1% formic acid) and 2% solution B (90% ACN and 0.1% formic acid) at 0 min to 40% solution B (90% ACN and 0.1% formic acid) at 90 min with a flow rate of 300 nl min⁻¹. nESI‐MS/MS measurements were performed on Orbitrap Elite (Thermo Fisher Scientific, Bremen, Germany) with the following settings: Full scan range 350–1800 m/z resolution 120000 max. 20 MS2 scans (activation type CID), repeat count 1, repeat duration 30 s, exclusion list size 500, exclusion duration 30 s, charge state screening enabled with the rejection of unassigned and +1 charge states, minimum signal threshold 500.

For protein identification, raw data were searched with the SEQUEST algorithm present in Proteome Discoverer version 1.3 (Thermo, Germany) as described in (Valledor and Weckwerth, 2014) using protein FASTA. We have used the following settings in Proteome Discoverer for data analysis which include: Peptide confidence: High, which is equivalent to 1% FDR, and Xcorr of 2, 3, 4, 5, 6 for peptides of charge 2, 3, 4, 5, 6. The variable modifications were set to acetylation of the N-terminus and oxidation of methionine, with a mass tolerance of 10 ppm for the parent ion and 0.8 Da for the fragment ion. The number of missed and/or non-specific cleavages permitted was 2. There were no fixed modifications, as dynamic modifications were used. The identified proteins were quantitated based on total ion count, followed by an NSAF normalization strategy (Paoletti et al., 2006).

The total number of spectra counts for the matching peptides from protein k (PSM) was divided by the protein length (L), then divided by the sum of PSM/L for all N proteins.

All the MS/MS spectra of the identified proteins and their meta-information were further uploaded to the PRIDE repository, Project accession: PXD035952.

A t-test was applied to find significantly altered proteins between control and stress condition in each genotype, and proteins with p-value < 0.05 were further considered. Differentially expressed proteins (DEPs), fold change (FC) were calculated for all identified proteins in HS compared to CON. Proteins with FC ≥ 2 (upregulation) or ≤ 0.5 (downregulation) were considered. Gene ontology was determined using RiceNetDB (http://bis.zju.edu.cn/ricenetdb/documentation.php). Only gene ontology (GO) terms with a p-value ≤ 0.01 were used for the analysis. Gene ontology network analysis for the identified DEPs was performed using Cytoscape software (Version 3.7.1) with ClueGO (Version 2.5.4) plug-in. ClueGo was used to visualize non-redundant biological terms for large clusters of genes in a functionally grouped network. The Oryza sativa (Japonica) marker list was used for searching biological processes, molecular function and cellular components related to proteins identified. P-value was calculated using a two-sided hypergeometric test and Benjamini-Hochberg adjustment for multiple test corrections. GO terms with P-value<0.05 were considered significant. The protein-protein interactions network in each set of up and down-regulated proteins was generated using the STRING database (Version 11.0). We restricted the parameter with the highest confidence to 0.9 and neglected the single molecules found in the interaction map to best predict the interactions.

Heat stress demonstrated phenotypic changes in N22 plants where mature spikelets were deformed, and pollen viability was reduced, especially after long-term heat stress compared to short-term heat stress ( Figure S1 ).

The detailed lists of the identified proteins in Moroberekan, IR64 and N22, are provided in supporting information ( Tables S1–S3 ). In total, proteins identified in Moroberekan were 2689 in control (referred to as MCT for Moroberekan control), 2529 in ST_HS (MST for Moroberekan ST_HS) and 2687 in LT_HS (MLT for Moroberekan LT_HS). Similarly, proteins identified in IR64 were 2739 in ICT (IR64 control), 2690 in IST (IR64 ST_HS) and 2973 in ILT (IR64 LT_HS) and in N22, 3072 proteins were identified in NCT (N22 control), 3009 in NST (N22 ST_HS) and 3045 in NLT (N22 LT_HS) as shown in Figure S2A . Thus, the number of identified proteins was in the order of N22>IR64>Moroberekan. The number of unique and common proteins in each rice genotypes is represented using the Venn diagram in Figure S2B . Between the control and ST_HS regimes, the number of overlapping proteins was higher in N22 compared to Moroberekan and IR64. More proteins were common in the control and LT_HS in Moroberekan than in the other two rice genotypes. Between ST_HS and LT_HS, overlapping protein numbers were higher in Moroberekan than in the other two rice genotypes. The identified proteome from three rice genotypes was compared with transcripts expressed during different developmental stages of rice plants. For this, co-expression analysis tools were used from the rice genome annotation project database. As an output, a z-score trend plot of the locus IDs of proteins was used. In all the rice genotypes, there was a clear dense peak just above the inflorescence P6 samples (mature pollen stage) in the plot generated using the inflorescence and seed developmental series (GSE6893) dataset and panicle 5 (heading panicle stage) and stamen (stamen 1 day before flowering) in the plot generated using tissue atlas from Minghui 63 rice (GSE19024) ( Figure S3 ). These results suggested that the proteins identified in the present work corresponded greatly with the transcripts expressed in anthers of the mature pollen stage.

The lists of the differentially expressed proteins (DEPs) in Moroberekan, IR64 and N22 rice types are provided in supporting information Table S4 . The total identified proteins were sorted into DEPs between the control and stressed samples by applying p-value<0.05 ( Figure 2A ). In Moroberekan, 734 DEPs were identified in MST and 254 in MLT. Of the 734 MST DEPs, 130 proteins were upregulated >2 fold, and of these, 103 proteins were exclusively identified in ST_HS. In the same dataset, 144 proteins were downregulated, and of these, 104 proteins were present only in the control but not in the HS condition. Of the 254 MLT DEPs, 74 proteins were upregulated, of which 64 were uniquely identified in LT_HS. In this analysis, 94 proteins were downregulated and of which 69 proteins were exclusively present in the control condition. In IR64, 243 DEPs were identified in IST, out of which 72 were upregulated, and 60 were noted exclusively in IST condition. In the same dataset, 85 DEPs were downregulated, of which 52 were exclusively present only in the control condition. Similarly, in ILT, 316 DEPs were identified, which include 131 upregulated, of which 105 were exclusively present in the LT_HS condition. In the same dataset, 79 DEPs were downregulated, and 42 were exclusively present in the control condition. In N22, 243 DEPs were identified in NST. At FC >2, 56 DEPs were upregulated and 33 were exclusively present in the NST. Under the NST regime, 82 DEPs were downregulated, and 71 were exclusively present in the control condition. Out of 312 DEPs identified in NLT, 78 were upregulated, of which 50 were unique under LT_HS. In the same dataset, 64 DEPs were downregulated, of which 29 were present in control but not in HS condition.

The analysis of the overlapping DEPs among ST_HS and LT_HS datasets showed that in Moroberekan, 30 upregulated proteins were common between MST and MLT treatments. In comparison, 55 downregulated proteins were common between MST and MLT datasets ( Figure 2B ). In IST and ILT, 22 DEPs were common. In NST and NLT datasets, 9 DEPs were commonly upregulated and 31 were downregulated ( Figure 2B ). The number of common DEPs between all three rice genotypes was smaller in ST_HS than in LT_HS treatments. Among the three rice genotypes, there were no common upregulated DEPs, while two common DEPs down-regulated were β-amylase (Os07g35880) and cytochrome P450 (Os02g38290) in the ST_HS treatment. In ST_HS, the numbers of uniquely upregulated DEPs were 127, 66 and 31 in Moroberekan, IR64 and N22, respectively ( Figure S4 ). In ST_HS, unique downregulated proteins were 129, 68 and 74 in Moroberekan, IR64 and N22, respectively. Overall, in LT_HS, unique upregulated proteins were 68, 126 and 70 in Moroberekan, IR64 and N22, respectively, and downregulated proteins were 88, 71 and 61 in Moroberekan, IR64 and N22, respectively ( Figure S4 ).

Cytosolic ClpB/HSP100 family proteins are molecular chaperones that promote the renaturation of protein aggregates and are required to develop acquired thermotolerance (Hong and Vierling, 2000). To validate the protein abundance of Hsp 101, a western blot analysis of control and heat-stressed anther tissues was performed. Moroberekan and IR64 showed higher levels of Hsp101 than N22 under ST_HS ( Figure 3 , upper panel). Similarly, the Hsp101 protein showed higher levels in Morobereken and IR64 than N22 ( Figure 3 , lower panel).

The relative abundance of the selected proteins differentially expressed in at least one of the rice genotypes and at least one condition related to diverse cellular pathways was used for constructing the heat maps for their relative expression ( Figure 4 ). In amino acid metabolism, most of the proteins were downregulated. S-adenosylmethionine synthetase was upregulated in all the rice genotypes except in MST. Amino methyltransferase, a glycine cleavage T family protein, was upregulated in both NST and NLT but downregulated in MST, IST and ILT. Phospho-2-dehydro-3-deoxyheptonate aldolase was downregulated in all the rice genotypes. GHMP kinase and ATP binding protein were specifically upregulated in NST and NLT. Several stress-associated genes were upregulated. These included Hsp20/alpha crystallin family protein, DnaK family protein, DnaJ homolog subfamily B member 11 and Hsp101. Genotype-specific Hsp20/alpha-crystallin family protein levels were identified in ST_HS and LT_HS of Moroberekan, IR64 and N22 ( Figure 4 ). On the other hand, the proteins in the biotic and drought/salt categories were down-regulated.

HVA22 (ABA and stress-inducible gene) was downregulated in MST, and no change was noted in NST and NLT. Endothelial differentiation-related factor 1 was upregulated in Moroberekan and downregulated in IR64 and N22. Most of the DEPs involved in fatty acid biosynthesis were unchanged or downregulated. Interestingly, most ribosomal proteins were downregulated in the three rice genotypes in all conditions ( Figure 4 ). However, a few ribosomal proteins, including 60S ribosomal protein subunits, 60S ribosomal protein L36-2 and ribosomal protein L7Ae were heat induced in NLT. Proteins related to post-translational modification were also heat inducible, especially in the case of NLT. Most photosynthesis-related proteins were downregulated except for a few chlorophyll-binding proteins, specifically in NLT.

GO enrichment analysis determined that three major biological processes (BP): cellular process, metabolism, and physiological process were enriched ( Figure 5 ). These terms were enriched in both up-and down-regulated DEPs of all three rice genotypes, except the cellular process category which was not represented in IST-UP (upregulated in IST of IR64). In response to the abiotic stimulus category, entries were identified in MST-UP and IST-UP, while no entries were found for it in NST-UP. Interestingly, the biological process of transport was enriched in MST-UP (with 26 proteins) and NST-UP (with 14 proteins) and ILT-DN. In the molecular function (MF), the transport category was enriched in MST-UP. In the cellular component, GO terms related to plastids, protein complexes, ribonucleoprotein, ribosomes, and thylakoids were enriched in MST-DN, NST-DN and MLT-DN. GO ontology of BP: protein biosynthesis and macromolecular biosynthesis were enriched in MST-DN only. This trend was reflected in enriched MF: structural molecule activity and RNA binding of MST-DN. The vacuole was enriched in MST-UP and MST-DN and, to a small extent, in NST-UP.

The catalytic activity was noted in all treatments except in protein regulation of LT-HS. To further interpret the biological relevance of DEPs, GO terms and pathways were functionally grouped in ClueGO ( Figure 6 ). ClueGO integrates GO terms and KEGG pathways, generating a functionally organized GO/pathway network. In MST-UP proteins, a network of proton export across the plasma membrane was identified. A major network in the N22 proteins NST-UP set was GTPase activator activity and vitamin E biosynthetic process. GTPases play an essential role in vesicular trafficking and act as molecular switches regulating plant developmental processes, including embryonic development and reproduction. Proteins in NLT-UP were enriched in unfolded protein binding process and cellular components cytosol, endoplasmic reticulum (ER) and chloroplast. Accumulation of misfolded proteins is at the core of HS in all cellular compartments, and in the ER, it activates a signaling unfolded protein response (UPR) pathway. The upregulation of unfolded protein binding process is in line with the reports showing components of the endoplasmic reticulum-folding machinery and UPR is upregulated at specific stages of pollen development, pollen germination, polar tube growth and fertilization (Fragkostefanakis et al., 2016; Chaturvedi et al., 2021).

The STRING database was used to identify the, interactions among proteins in three rice genotypes. Different functional clusters emerged among the interacting proteins. Only selected interactions with more than 0.9 confidence are shown ( Figure 7 ). In the MST-UP, nine distinct network clusters were noted (CI-CIX). These clusters included the respiratory metabolism cluster (CI), ribosomal proteins cluster (CII), and Hsps cluster (CIV). In the MST-DN, four distinct clusters were observed, including ribosomal protein (CI) and RNA modification (CII). The IST-UP revealed three clusters of interactions, including an Hsp cluster (CI). IST-DN category showed two clusters, one of which was the ribosomal proteins cluster. In both NST-UP and NST-DN, three clusters were identified. The NST-DN category contained a large cluster of ribosomal proteins, RNA splicing (CI) and starch metabolism (CIII). MLT-UP and MLT-DN two and four clusters were identified, respectively. The MLT-DN contained a cluster of ribosomal proteins (CI). ILT-UP and ILT-DN contained seven and four clusters of protein-protein interactions, respectively. ILT-UP had a proteasome proteins cluster (CI), whereas ILT-DN had a ribosomal protein interaction cluster (CI). In NLT-UP, two clusters were observed, both of which were of Hsps (CI and CII). In NLT-DN, four protein interaction clusters were noted, including one big ribosomal protein cluster (CI), one RNA splicing cluster (CIII) and one starch metabolism cluster (CIV). In the NLT-UP, cluster CI contained chaperone protein ClpB1 (Os05g44340) with 26.7 kDa sHsp (Os03g14180) and 18.0 kDa class II sHsp (Os01g08860), which further interacted with 24.1 kDa mitochondrial sHsp (Os02g52150). Similarly, 24.1 kDa sHsp was found in anther proteome of N22 in an earlier study, and its higher accumulation was linked to better thermotolerance of N22 (Jagadish et al., 2010).

In all three rice genotypes, different Hsps were highly induced. Hsps showed genotype - specific protein abundance pattern which can also be used as genotype-specific markers for the heat stress response. The abundance of Hsp20, chaperone protein DnaJ homolog subfamily B member Os05g06440.1 and Hsp101 (Os05g44340) were upregulated across all the three genotypes under stress condition. Hsp101 is an important protein in imparting thermotolerance in bacteria, yeast and Arabidopsis (Kumar et al., 2016; Kumar et al., 2020; Tiwari et al., 2021). The introduction of Arabidopsis thaliana Hsp101 (AtHsp101) into rice made the host variety more thermotolerant (Katiyar-Agarwal et al., 2003). Hsp20 provide thermotolerance to various levels of life forms such as bacteria and plants (Sarkar et al., 2009). Hsp20 co-aggregate with misfolded proteins (Sarkar et al., 2009; Kumar et al., 2015; Ungelenk et al., 2016; Sarkar et al., 2020). Furthermore, ATP-dependent Hsp70 and Hsp100 chaperones facilitate the solubilization and refolding of protein aggregates (Cherkasov et al., 2013; McLoughlin et al., 2016; Merret et al., 2017).

In NLT, upregulated proteins demonstrated a significant interaction with Hsps ( Figure 6 ). This interaction was not enriched in ILT and MLT. A major interaction with the ribosomal proteins suggests reduction in translation during heat stress. Other major interactions with proteasome and dehydrogenase also indicate that MST is involved in more protein degradation. Heat-induced ribosome inhibits pre-rRNA processing which triggers the imbalance in ribosomal profiles in Arabidopsis thaliana (Merret et al., 2015; Darriere et al., 2022).

Another category of stress proteins is related to carbohydrate metabolism. N22 was enriched in proteins such as trehalose-6-phosphate synthase, trehalose synthase, and haloacid dehalogenase-like hydrolase. The trehalose accumulation in rice protects plants from various stresses (Benaroudj et al., 2001). HS affects rice sugar partitioning. It appears that heat-tolerant plant types are better at preserving starch and sugar levels in pollen, which helps them germinate at higher temperatures (Pressman et al., 2002; Firon et al., 2006; Chaturvedi et al., 2015; Jegadeesan et al., 2018). Protein import into the nucleus is another process enriched in N22. Hsfs and other transcription factors are dynamically transported across the nucleus to induce the expression of Hsps and other stress-responsive genes (Scharf et al., 2012). This may indicate more efficient trafficking in N22 under heat stress conditions.

Thermotolerance involves different signaling pathways where phytohormones play important roles. HS leads to the abscission of reproductive organs due to increased abscisic acid and ethylene levels and reduced levels of auxins (Binder and Patterson, 2009). Also, the alterations in auxin biosynthesis and cytokinin content reduces pollen sterility and kernel filling, respectively, in cereals (Banowetz et al., 1999; Sakata et al., 2010). In the present study, DEPs involved in hormone metabolism were mostly downregulated under HS conditions. HVA22 was specifically induced in IST and ILT but downregulated in MST and MLT, and no change was noted in NST and NLT. This protein acts downstream of GAMyb and activates programmed cell death and other GA-mediated processes (Guo and Ho, 2008). Another protein, CYTOCHROME P450 51G1 (LOC_Os11g32240.1), was highly upregulated in IST, ILT and NLT. This protein participates in brassinosteroid biosynthesis. The application of brassinolide promotes panicle ripening in rice (Saka et al., 2015). Also, brassinolide is known for its induced resistance to abiotic stress in plants.

Our analysis revealed significant changes in major cellular processes and pathways associated with stress regimes and rice genotypes like amino acid metabolism, maintenance of cellular homeostasis and photosynthesis ( Figure 4 ). The amino acid metabolism is a major process affected by HS (Wang et al., 2018). In N22, higher number of proteins involved in amino acid synthesis were upregulated. S-adenosyl methionine (SAM) synthetase protein participates in various biological processes and methionine metabolic pathways. This protein was specifically downregulated in MST. Increased levels of this protein were noted in NST and NLT. This protein also plays a key role in pollen tube growth in the mutant of methionine adenosyltranferase3 (MAT3, one of the four kinds of SAM synthetase genes in Arabidopsis) (Chen et al., 2016). Some amino acid metabolism proteins like aminomethyltransferase, chorismate mutase, and cystathionine gamma synthase showed genotype- specific responses but not found in Moroberekan. This change in amino acid metabolism could be linked to the genetic variations in heat tolerance in rice.

Proteins upregulated in NLT were majorly enriched in biological processes like unfolded protein binding ( Figure 5 ). This process is critical for repairing unfolded or damaged proteins formed during HS (McLoughlin et al., 2019). Furthermore, proteins involved in carboxylic acid metabolism, including those involved in amino acid metabolism, were enriched under NLT-UP. This indicates that N22 anthers can efficiently maintain their energy metabolism under heat stress conditions. The pollen tissues at the mature stage are metabolically active and need energy and new proteins for the successful fertilization (Chaturvedi et al., 2013; Selinski and Scheibe, 2014). Proteins associated with aerobic respiration enrichment were found in MLT-UP, implicating an efficient TCA cycle for energy production. A major process noted in NST-UP was the GTPase activator. GTP metabolism is involved in protein trafficking (Hutagalung and Novick, 2011). This process may be important for maintaining cellular homeostasis during HS. It was also found that NST-UP was enriched in proteins involved in vitamin E biosynthesis ( Figure 5 ). Vitamin E protects photosystem II against heat and drought stresses (Havaux et al., 2005).

The amount of information available on photosynthesis-related proteins in rice anthers is limited. Several proteins related to light-harvesting complex II and photorespiration were upregulated in N22, especially in NLT. The chlorophyll a-b binding protein is key to balancing the excitation energy between the two photosystems and was found upregulated in NST. Considering the importance of ABC transporter in pollen exine formation and pollen-pistil interactions (Chang et al., 2016), ABC-type transporters found in both NST and NLT demonstrate relevance in the developmental process.