The anther transcriptome of cold-tolerant rice cultivars is largely insensitive to temperature changes

The anther transcriptome of cold-tolerant rice cultivars is largely insensitive to temperature changes | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article The anther transcriptome of cold-tolerant rice cultivars is largely insensitive to temperature changes Koichi Yamamori, Seiya Ishiguro, Kei Ogasawara, Kayyis Lubba, and 4 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4399503/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Many studies of stress tolerance in plants have characterized genes that show differences among a small number of lines with clearly distinct tolerance or sensitivity to the given stress. From the few cloned genes, it is difficult to genetically interpret intermediate tolerance or susceptibility levels and explain the complexity of stress responses and tolerance. In this study, we explored the changes in the transcriptome of anthers from 13 rice lines with different cold tolerance grown under control conditions or exposed to 4 days of cold stress to look for correlations between cold tolerance at the booting stage and expression levels. When examining the overall expression patterns in anthers at low temperature, the cold-tolerant lines tended to have relatively few highly expressed genes, and the expression levels of ribosome-related genes tended to be lower in cold-tolerant lines than in cold-sensitive lines. Importantly, we observed these different expression patterns between the cold-tolerant and -sensitive lines regardless of whether cold stress had been applied. Minimal expression changes under cold stress tended to be characteristic of the cold-tolerant lines, especially in repetitive sequences. We also identified unknown genes whose expression was cold responsive and common to all the lines studied. We conclude that rice lines whose transcriptome remains constant or insensitive in response to cold stress are more tolerant to low-temperature exposure during the booting stage than rice lines with more widespread expression changes. Anther booting stage cold stress flowering stage microarray pollen RNA-seq sterility transcriptome analysis Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Background Rice ( Oryza sativa ) is a tropical crop that is susceptible to low temperatures. Low-temperature stress, especially during the booting stage, causes irreversible pollen sterility, resulting in lower yields (Hayase et al. 1969 ; Satake 1976 ). Rice, native to the tropics and subtropics, can now be grown at high latitudes above 40°N (Sagehashi et al. 2022 ). Acquisition of cold tolerance traits during the booting stage is important for rice acclimation to low temperatures at high latitudes (Li et al. 2021 ; Liu et al. 2018 ; Zhang et al. 2017 ). Genetic analyses have been conducted in various rice cultivars to identify factors that confer cold tolerance at booting stage, and many quantitative trait loci (QTLs) have been identified (Andaya and Mackill 2003 ; Dai et al. 2004 ; Kuroki et al. 2007 ; Li et al. 1997 ; Li et al. 2021 ; Oh et al. 2004 ; Saito et al. 2010 ; Saito et al. 2001 ; Shimono et al. 2016 ; Shirasawa et al. 2012 ; Suh et al. 2010 ; Sun et al. 2019 ; Takeuchi et al. 2001 ; Xu et al. 2008 ; Zhang et al. 2017 ; Zhu et al. 2015 ). Typically, these analyses use two lines with contrasting cold tolerance at the booting stage. However, of the 30 QTLs identified to date, only four genes have been cloned: Cold tolerance at the booting stage 1 ( CTB1 ), CTB2 , CTB4a , and qCTB7 (Li et al. 2021 ; Saito et al. 2010 ; Yang et al. 2023 ; Zhang et al. 2017 ). Numerous genes are involved in abiotic stress responses, including low temperature, relying on multiple pathways that act independently or interact with one another (Kidokoro et al. 2021 ; Kidokoro et al. 2022 ; Kuroki et al. 2007 ; Liu et al. 2018 ; Lou et al. 2022 ; Nakashima et al. 2007 ; Park et al. 2010 ; Zhang et al. 2022b ). Therefore, we hypothesize that many genes participate in cold tolerance at the booting stage. While the cloning of specific genes is crucial, it would not allow a global understanding of the genetic mechanisms behind cold tolerance at this developmental stage that would apply across many rice genotypes. In a previous study, we investigated the relationship between genome-wide cold responsiveness and cold tolerance in anthers at the booting stage by microarray analysis, focusing on repetitive sequences (RSs) (Ishiguro et al. 2014 ). Expression of transposable elements (TEs), a major component of RSs, is repressed by epigenetic regulation, thus maintaining genome stability (Choi and Lee 2020 ; Parker et al. 2022 ). However, TE repression is released under stress conditions (Hu et al. 2012 ; Ito et al. 2011 ; Zhou et al. 2017 ). For example, among five rice cultivars differing in their cold tolerance, more tolerant cultivars tended to show fewer changes in RS expression in response to cold stress than more sensitive cultivars (2014). This result suggested that cold-tolerant lines at the booting stage tend to suppress TE activation at lower temperatures and maintain genome stability. To our knowledge, the idea that genome-wide expression changes, rather than those of specific genes, are related to stress response and tolerance has not been previously proposed in stress tolerance research. However, the difficulty of explaining the evocation of cold tolerance at the booting stage from individual loci or common QTLs lends some validity to this claim. In this study, to explore rice cold tolerance at the booting stage, we analyzed the relationship between genome-wide responsiveness to cold stress and tolerance using 13 rice lines representing a continuous gradient of cold tolerance. We examined anther transcript levels during the booting stage with 60K microarrays containing probes for about 38,000 genes and 23,000 RSs. We calculated the pollen fertility index (PFI) of each genotype and looked for correlations with genome-wide expression levels, identified differentially expressed genes (DEGs) and RSs (DERs) whose expression was altered by low temperature, and characterized genes whose change in expression correlated with PFI. We propose that genome-wide expression patterns can well explain the differences for cold tolerance at the booting stage among rice lines. In support of this idea, previously cloned cold tolerance genes did not always show high correlations with cold tolerance among rice lines. Results Pollen fertility of 13 rice lines under cold stress and control conditions The 13 rice lines used in this study comprised eight temperate japonica , one tropical japonica , one indica , and three recombinant inbred lines (RILs) derived from a cross between a temperate japonica cultivar (‘A58’) and an indica cultivar (‘I33’) (Table 1). These lines are genetically diverse and expected to represent various degrees of cold tolerance. We grew all plants under greenhouse conditions (day: 25°C, night: 19°C) before exposing them to continuous low temperature (12°C) for 4 days at the booting stage in the growth chamber. To examine pollen fertility from these plants, we collected the third, fourth, and fifth flowers from the top of each primary rachis branch after heading (Yamamori et al. 2021 ). Pollen fertility for the 13 genotypes ranged from 59.7–92.6% for control plants maintained under 25°C/19°C day/night temperature cycles and from 2.2–85.6% in the cold-treated plants (Fig. 1 ). The line with the lowest pollen fertility under the cold treatment was the RIL ‘R75’ and the highest was the tropical japonica ‘Lambayeque1’. We calculated the PFI in response to cold treatment as the ratio between the average pollen fertility under the cold treatment and the average pollen fertility under control conditions. Genotypes with a higher PFI are more cold tolerant. We prepared another set of materials to verify the results obtained by microarray analysis and to perform a transcriptome deep-sequencing (RNA-seq) analysis with 4 of the 13 genotypes, choosing Lambayeque1 and ‘PL9’ as cold-tolerant lines and ‘Sasanishiki’ and ‘Nipponbare’ as sensitive lines. For this RNA-seq analysis, we independently grew the four lines under the same conditions as the 13 genotypes used for microarray analysis and collected samples in triplicate. The pollen fertility of the cold-treated materials was 3.9% (Nipponbare), 24.8% (Sasanishiki), 72.1% (Lambayeque1), and 85.8%, (PL9) (Supplemental Fig. S1 ). Although the exact order of the four lines differed between the microarray and RNA-seq analyses, Nipponbare and Sasanishiki were clearly more sensitive to cold treatment than Lambayeque1 and PL9 in both datasets. Anther microarray analyses with 13 lines and RNA-seq analyses with 4 lines We collected anthers at the booting stage from control and cold-treated plants for each of the 13 lines, as two biological replicates, to examine their transcriptomes by microarray analysis. We used the 60K microarray designed with 37,955 probes from protein-coding genes and 23,013 probes from RSs (Supplemental Table S1 ). This microarray allowed us to comprehensively survey nearly all rice genes and RSs across the 13 lines. We converted the raw fluorescence intensity data of each probe on the microarray by logarithmic transformation (log 2 ) and normalization by the quantile method (Bolstad et al. 2003 ). We plotted the normalized expression values for each probe in control and cold-treated plants to assess the extent of correlation. In the resulting scatterplots (Fig. 2 ), each rice line showed a high coefficient of determination ( R 2 ), ranging from 0.82 (A58) to 0.95 (Sasanishiki). Notably, probes falling outside of the diagonal (defined as y = x ± 1) varied among the lines (Fig. 2 ). To validate the microarray analysis, we performed an RNA-seq analysis of anthers from the Nipponbare, Sasanishiki, Lambayeque1, and PL9 lines. These plants were grown under the same conditions as for the microarray analysis but were independently collected. We obtained expression values for 37,860 genes using the Nipponbare reference genome (IRGSP-1.0 2020-12-02). When we plotted the RNA-seq data between the control and cold-treated anther samples in each line as a scatterplot, we obtained R 2 values of 0.89, 0.88, 0.88, and 0.89 for Nipponbare, Sasanishiki, Lambayeque1, and PL9, respectively (Supplemental Fig. S2 ). We also plotted the expression values obtained through microarray and RNA-seq analyses for each of these four lines for control or cold-treated anthers. We calculated Pearson’s correlation coefficients ( r ) for each comparison, resulting in r values of 0.71 or above, indicating positive correlations. Specifically, the r values for the control and cold-treated samples were 0.82 and 0.80 for Nipponbare, 0.77 and 0.76 for Sasanishiki, 0.73 and 0.71 for Lambayeque1, and 0.75 and 0.77 for PL9, respectively (Supplemental Fig. S3). The r values of each line were comparable between the control and cold-treated samples and positive, suggesting sufficient equivalence between the microarray and RNA-seq datasets. Correlation between expression levels and pollen fertility As the 13 lines showed varying levels of cold tolerance, based on their PFI values, we wished to examine the possible correlation between gene expression levels in anthers and cold tolerance. Before conducting this analysis, we selected reliable and reproducible probes from the microarrays. In general, microarray probes with weaker fluorescence are less accurate. Indeed, probes with log 2 expression < 5 were associated with larger coefficients of variation between replicates (Supplemental Fig. S4). We therefore applied a cutoff of log 2 expression = 5, only considering probes for genes and RSs with higher values for analysis. From 60,168 original probes, we retained 38,060 and 37,807 probes with log 2 expression ≥ 5 in common among the 13 lines for control or cold-treated anthers, respectively. To examine the relationship between the PFI in the 13 lines and the expression of each gene and RS detected by the microarrays, we calculated the r values obtained from a scatterplot of expression (on the y -axis) as a function of the PFI (on the x -axis). In the control samples, we detected a significant and positive correlation ( P < 0.05, r ≥ 0.55) for 668 probes (367 genes and 301 RSs) and a significant and negative correlation ( P < 0.05 and r ≤ − 0.55) for 724 probes (633 genes and 91 RSs) among the 38,060 probes (Fig. 3 A, B). A similar analysis conducted on the 37,807 reliably expressed probes (with log 2 expression ≥ 5) in the cold-treated samples yielded 234 probes (169 genes and 65 RSs) with a significant and positive correlation and 530 probes (479 genes and 51 RSs) with a significant and negative correlation (Fig. 3 C, D). In both cases, we identified more probes with a negative correlation to the PFI than with a positive correlation. This difference was significant in the cold-treated samples, based on Fisher's exact test of departure from a 1:1 ratio between probe types ( P < 0.001) (Fig. 3 D). The number of probes from genes and RSs whose expression levels were negatively correlated with the PFI was already slightly higher than that for positively correlated probes under normal growth conditions; this difference became more pronounced after cold treatment (Fig. 3 B, D). We extracted probes (for microarrays) or genes (for RNA-seq) showing a significantly different expression level between the cold-tolerant lines (Lambayeque1 and PL9) and the cold-sensitive lines (Nipponbare and Sasanishiki) with a false discovery rate (FDR) < 0.05 and an absolute log 2 (FC) ≥ 1. From the microarrays, we identified 23 probes with higher expression levels in the tolerant lines than in the sensitive lines when grown under control conditions and another 229 probes with higher expression levels in the sensitive lines than in the tolerant lines (Supplemental Fig. S5A). In the cold-treated samples, 19 probes had higher signals in the tolerant lines than in the sensitive lines, while another 167 probes yielding higher signal intensity in the sensitive lines than in the tolerant lines (Supplemental Fig. S5A). We observed similar patterns in the RNA-seq data for both growth conditions, as more genes showed a higher expression in the sensitive lines than in the cold-tolerant lines ( P < 0.001). Indeed, there were 918 genes with higher expression in the sensitive lines under control conditions, compared to 530 more highly expressed genes in the tolerant lines under the same conditions; likewise, in the cold-treated samples, 2,453 genes were more highly expressed in the cold-sensitive lines, compared to only 1,223 more highly expressed genes in the cold-tolerant lines (Supplemental Fig. S5B). Although the exact numbers of probes or genes differed between the microarray and RNA-seq datasets, both analyses support the notion that the sensitive lines have significantly more up-regulated genes than the cold-tolerant lines under both growth conditions. In summary, the expression levels of genes and RSs in anthers tended to be less responsive to cold treatment in the cold-tolerant lines relative to the cold-sensitive lines. Notably, this pattern was already present in the anthers collected from plants grown under control conditions. Thus, the cold tolerance of these lines may be latent before cold treatment. Identification of differentially expressed genes and repetitive sequences in response to cold stress We identified DEGs and RSs (DERs) between the control and cold treatment samples using the microarray datasets from the 13 lines, yielding 27,908 DEGs and DERs with at least a 2-fold difference in expression (in either direction) between the control and cold-treated samples in at least one line. Among these, there were 1,897–4,906 up-regulated genes and RSs (Ure) with higher expression (log 2 [FC) ≥ 1) in the cold-treated samples relative to control samples, and 1,942–5,625 down-regulated genes and RSs (Dre) with lower expression (log 2 [FC] ≤ − 1) upon cold treatment (Fig. 4 A). We obtained more DEGs than DERs within the Ure and Dre lists (Fig. 4 A). We then compared the number of DEGs and DERs among the 13 lines. We observed fewer common DERs compared to common DEGs as more lines were considered; for instance, the number of DEGs common to seven lines was 975 for Ure and 1,039 for Dre, compared to 10 Ure RSs and 104 Dre RSs (Fig. 4 B). While the numbers of Dre and Ure genes were comparable, the number of Dre RSs exceeded that of Ure RSs, suggesting that down-regulation of RSs is a distinct cold response. As more lines were considered, the numbers of common Ure and Dre genes decreased; nonetheless, there were about 100 common genes each for Ure and Dre with the same expression pattern in 12 lines (Fig. 4 B). By contrast, fewer RSs showed a similar expression pattern across many lines, with fewer than 100 Ure and Dre RSs when considering five and eight lines, respectively (Fig. 4 B). We identified 20 Ure genes and 24 Dre genes showing the same expression profile in all 13 lines (Tables 2, 3, Supplemental Fig. S6), which may be characteristic of the intrinsic transcriptome responses to low temperature at the booting stage. With the exception of the Ure gene Os06g0246500 ( r = 0.66), we observed no significant correlation between the log 2 (FC) of the 43 remaining DEGs and the PFI values (Supplemental Fig. S7). Of these 44 common DEGs, we identified six as Ure genes and nine as Dre genes in the RNA-seq data (Tables 2, 3). Among these six Ure genes, we noticed Os04g0568700, which encodes a heat shock factor (HSF) related to multiple stress responses including cold (Chauhan et al. 2011 ; Guo et al. 2016 ; Zhang et al. 2022b ) (Supplemental Fig. S6A, Table 2). Of the nine Dre genes, the most significantly down-regulated gene was Os11g0150400, encoding a stress-responsive alpha/beta barrel domain–containing protein (Supplemental Fig. S6B, Table 3). Correlation between expression change and pollen fertility by cold treatment Based on the microarray data with the 13 lines, we analyzed the relationship between expression change and pollen fertility (as measured by the PFI) in response to cold treatment. To this end, we calculated the r values for each probe from a scatterplot between the PFI and the log 2 (FC) values for all 13 lines. Of the 27,908 DEGs and DERs, 905 showed a significant negative correlation ( r ≤ − 0.55, P < 0.05) with the PFI, and another 514 had a significant positive correlation ( r ≥ 0.55, P < 0.05) with the PFI (Fig. 5 A). We further divided the DEGs and DERs based on the median log 2 (FC) values (designated expression change index [ECI] hereafter) across the 13 lines. A positive ECI indicates an Ure probe in response to cold treatment, while a negative ECI indicates a Dre probe under the same condition. With these parameters, we divided the 1,419 probes above into four types by combining r values and ECI values (Fig. 5 B). Type 1 indicated that 600 Ure probes (137 genes and 463 RSs) tend to have weaker up-regulation in cold-tolerant lines than in cold-sensitive lines; 233 Ure probes (152 genes and 81 RSs) in type 2 have stronger up-regulation in the tolerant lines than in the sensitive lines; 305 Dre probes (172 genes and 133 RSs) in type 3 have stronger down-regulation in the tolerant lines than in the sensitive lines; 281 Dre probes (155 genes and 126 RSs) in type 4 have weaker down-regulation in the tolerant lines than in the sensitive lines. The typical expression patterns with highest r values in types 1–4 are illustrated by Os01g0699400, Os04g0604000, Os01g0191200, and Os01g0730600, respectively (Fig. 5 C). The number of genes was comparable across types, ranging from 137 to 172. Notably, while the number of RSs was similar among types 2–4, type 1 was characterized by far more (463) RSs (Fig. 5 B). This result suggested that most RSs related to pollen fertility are up-regulated in the sensitive lines relative to the tolerant lines. Consistent with this finding, Ishiguro et al. ( 2014 ) previously reported that the expression levels of RSs in cold-sensitive rice lines changed more in anthers in response to low temperatures than in cold-tolerant lines. The expression of RSs in tolerant lines may therefore be less responsive to cold treatment (Fig. 5 B). Since we only used four lines for RNA-seq analysis, we opted for different selection criteria to divide DEGs into four types equivalent to those obtained by the analysis of the microarray data. Specifically, we defined type 1 genes as those with a positive ECI value in the four lines and with a lower average log 2 (FC) for the two tolerant lines (Lambayeque1 and PL9) than for the two sensitive lines (Nipponbare and Sasanishiki); type 2 genes had a positive ECI value and with a higher average log 2 (FC) for the two tolerant lines than for the two sensitive lines; type 3 genes had negative ECI values and a lower mean log 2 (FC) for the two tolerant lines than for the two sensitive lines; and type 4 genes had a negative ECI value and a greater mean log 2 (FC) value for the two tolerant lines than for the two sensitive lines. The above analysis returned 3,167 genes, of which 295 belonged to the same type in the microarray and RNA-seq datasets, with 52 genes for type 1, 48 genes for type 2, 154 genes for type 3, and 41 genes for type 4 (Supplemental Fig. S8). Overlapping genes were relatively more abundant for type 3 DEGs than for the other types of DEGs and corresponded to genes whose expression levels decreased upon cold treatment in the cold-tolerant lines. Overall, the correlations between expression change and PFI clearly showed the tendency of numerous RSs with activated expressions in the sensitive lines and more genes with reduced expression in the tolerant lines. GO term enrichment analysis of anther transcriptomes From the transcriptome analysis of the 13 rice lines above, we obtained three sets of information: (1) the correlation between expression levels and the PFI, (2) a list of common DEGs across the lines, and (3) the extent of correlation between the ECI and PFI values. In each case, we generated lists of genes with characteristic expression patterns, which might be related directly or indirectly to the pollen sterility caused by cold exposure at the booting stage. To explore the function of these genes, we performed a Gene Ontology (GO) term enrichment analysis on each gene list (Supplemental Tables S2–S8), focusing on terms associated with fewer than 100 genes across the entire rice genome and exhibiting significant (FDR < 0.001) enrichment (Fig. 6 ). We first characterized the genes showing a negative or positive correlation in their expression levels under control growth conditions with the PFI values. From genes whose expression levels are negatively correlated with the PFI, we identified 15 significantly enriched GO terms, most of which were related to “cytosolic small ribosomal subunit” (GO:0022627) and “small ribosomal subunit” (GO:0015935) (Fig. 6 A right panel). Ribosome biosynthesis is known to be related to pollen development and the stress response in Arabidopsis ( Arabidopsis thaliana ) (Renák et al. 2014 ). Structural changes in ribosomes contribute to stress acclimation by forming stress-responsive ribosomes under stress conditions (Dias-Fields and Adamala 2022 ). This result also suggests that cold-sensitive lines may exhibit higher expression levels for ribosome-related genes than cold-tolerant lines. The cold tolerance of plants might thus be predicted from the expression levels of specific ribosome-related genes under normal growth conditions. For genes whose expression levels were positively correlated with the PFI in control samples, we obtained six significantly enriched GO terms (Fig. 6 A left panel). These GO terms were jasmonic acid (JA)- and fatty acid–related GO terms: “jasmonic acid mediated signaling pathway” (GO:0009867), “cellular response to jasmonic acid stimulus” (GO:0071395), “response to jasmonic acid” (GO:0009753), “cellular response to fatty acid” (GO:0071395), and “response to fatty acid” (GO:0070542) (Fig. 6 A left panel). JA is a phytohormone related to multiple stress responses including cold stress (Raza et al. 2021 ). The most common unsaturated fatty acids in plants, 18-carbon species, can be modified into bioactive molecules, including JA (He and Ding 2020 ). Therefore, we speculate that the expression of genes involved in JA signal transduction is activated by cold stress. The expression pattern of these JA-related genes remained higher in the cold-tolerant lines than in the sensitive lines under normal growth conditions. We also looked for GO terms enriched in Ure or Dre genes common to 11 lines or more. With fewer common DEGs across more lines, we failed to identify significant GO terms for DEGs shared by 12 or more lines (Supplemental Fig. S9). We detected five GO terms for DEGs shared by 11 lines (Fig. 6 B left). The GO term enriched among Ure genes was “positive regulation of response to salt stress” (GO:1901002), which would be consistent with a shared signaling pathway between salinity and cold stress responses (Zhang et al. 2022b ; Zhang and Xia 2023 ). The four GO terms enriched among Dre genes were related to DNA metabolism such as “DNA packaging” (GO:0044815), “nucleosome” (GO:0000786), “protein-DNA complex” (GO:0032993), and “structural constituent of chromatin” (GO:0030527) (Fig. 6 B right). These GO terms suggest that in most of the lines, DNA metabolism is reduced due to cold stress during the booting stage. We also obtained the list of GO terms associated with each of the four types of genes defined by correlations between ECI and PFI values (Supplemental Table S8). Genes showing expression patterns related to low temperature at the booting stage Here, we highlighted known and unknown genes related to cold tolerance based on the GO term enrichment analysis of the microarray and RNA-seq data. From the microarray data, we identified 43 and 49 genes whose expression levels were positively or negatively correlated with the PFI ( r ≥ 0.7 or r ≤ − 0.7), respectively, under control conditions (Supplemental Tables S9, S10). Among the 43 genes with a positive correlation, Os04g0623300 ( POLYAMINE OXIDASE3 , OsPAO3 ) showed the highest correlation with the PFI ( r = 0.85) and was previously reported to be up-regulated in response to cold stress in rice seedlings (Sagor et al. 2021 ). For the remaining 91 genes, an association with cold stress has not been reported until now, although we noticed three ribosome-related genes among them: Os06g0550000, Os01g0962600, and Os08g0234000. The expression levels of these three ribosome-related genes were negatively correlated ( r ≤ − 0.7) with the PFI under control conditions (Fig. 6 A), in agreement with the GO term enrichment analysis above. For cold-treated plants, the expression levels of 16 and 28 genes were significantly and positively or negatively correlated ( r ≥ 0.7 or r ≤ − 0.7), respectively, with the PFI (Supplemental Tables S9, S10). However, none of these genes have been reported to have a role in cold stress responses, except for MITOGEN-ACTIVATED PROTEIN KINASE KINASE KINASE70 ( OsMKKK70 ) (Mei et al. 2022 ). OsMKKK70 , whose expression was negatively correlated ( r = − 0.63) with the PFI (Supplemental Fig. S10), was reported to be associated with cold tolerance at the booting stage (Mei et al. 2022 ). Other genes whose expression levels were positively correlated with the PFI values included Os01g0962700 ( PEROXIDASE20 , PRX20 ), Os06g0156400 ( GDSL ESTERASE/LIPASE76 , OsGELP76 ), and Os03g0682200 ( ARGONAUTE12 , OsAGO12 ) (Supplemental Table S9). PRX20 , which encodes a peroxidase, was highly expressed in cold-tolerant lines ( r = 0.73 for the expression–PFI relationship) and was previously reported to participate in the response to salinity stress (Kim and Kim 2023 ). As with samples collected from control plants, of the 28 genes whose expression levels were negatively correlated ( r ≤ − 0.7) with the PFI, we noticed a ribosome-related gene, Os03g0284400 ( r = − 0.73), suggesting that ribosome status in rice plants under normal or cold stress conditions potentially influences pollen fertility upon cold exposure. We took a closer look at the 20 and 24 common Ure and Dre genes shared by all 13 lines (Tables 2, 3, Supplemental Fig. S6). The expression of these genes responds to cold stress in a similar manner regardless of their innate cold tolerance or sensitivity. The 20 common Ure genes included NAC-type transcription factor genes such as OsNAC5 and ONAC088 (Table S2 , Supplemental Fig. S6). We did not detect known stress response genes among the 24 common Dre genes but were intrigued by the presence of several histone genes, Os08g0427700 ( H2A ), Os04g0583600 ( H4 ), and Os08g0490900 ( H2B ), whose functional association with cold stress might be cell division (Table 3, Supplemental Fig. S6B). These 44 common DEGs did not include any known genes associated with cold tolerance at the booting stage. However, three genes known to be associated with cold tolerance were among common Ure or Dre genes when considering nine or more of the 13 lines: LATE EMBRYOGENESIS ABUNDANT9 ( OsLEA9 ), OsMKKK70 , and OsMYB4 (Table 4) (Lou et al. 2022 ; Mei et al. 2022 ; Park et al. 2010 ). Returning to genes from types 1–4 defined based on their correlation between the log 2 (FC) values and PFI values (Fig. 5 ), we focused on genes with pronounced changes in expression. Using a selection criterion of ECI > 1 or ECI < − 1, we obtained 86 genes (Supplemental Table S13), consisting of 20 type 1 genes, 13 type 2 genes, 25 type 3 genes, and 28 type 4 genes. Five of the 86 genes were also present among the DEGs classified into four types using the ECI and the log 2 (FC) values from the RNA-seq data. Two of these genes are known to participate in the cold tolerance of rice. Os02g0579000, a NAC-type transcription factor gene induced at low temperatures (Fang et al. 2014 ), was among the type 1 genes, which tend to be more highly expressed in the cold-sensitive lines. This gene showed an increased expression in 11 lines in response to cold stress, reaching higher levels in the cold-sensitive lines than in the tolerant lines (Supplemental Table S13). The other NAC-type transcription factor gene, Os03g0133000, known to respond to some abiotic stresses, including cold (Hong et al. 2016 ), was a type 2 gene, more highly expressed in cold-tolerant lines and exhibiting increased expression in 10 lines in response to cold stress. The identification of several known genes and their expression changes in response to low temperatures support the validity of the microarray-based analysis in this study. Discussion Expression levels and cold tolerance: negative correlation between pollen fertility and expression of ribosome-related genes We defined genes whose expression levels were negatively or positively correlated with the PFI, a measure of cold tolerance at the booting stage. There were significantly more genes with low expression in the cold-tolerant lines and with high expression in the cold-sensitive lines in response to cold stress than genes with high expression in the cold-tolerant lines and with low expression in the cold-sensitive lines (Fig. 3 C, D). We corroborated this result with our RNA-seq data obtained from two cold-tolerant liens and two cold-sensitive lines (Supplemental Fig. S5B). Among genes whose expression levels were negatively correlated with the PFI under cold conditions, we observed a significant enrichment for ribosome-related GO terms (Supplemental Table S5). When plants experience abiotic stresses, such as low temperatures, they produce inoperative ribosomes (Dias-Fields and Adamala 2022 ). This phenomenon is thought to function to conserve energy in response to stress by slowing down protein translation and growth (Bechtold and Field 2018 ; Zandalinas et al. 2022 ; Zhang et al. 2022a ). The lower induction of expression of these ribosome-related genes in cold-tolerant lines than in sensitive lines may reflect energy conservation, which may affect pollen development (Table 5). The expression of ribosome-related genes in anthers of the cold-tolerant lines tended to be lower than that of the sensitive lines under both cold treatment and control conditions (Fig. 6 A). Therefore, the expression of ribosome-related genes in the anthers of cold-tolerant lines may be lower than in cold-sensitive lines at all times. A GO analysis of genes whose expression levels under control conditions were significantly and positively correlated with the PFI identified JA-related genes, which are known stress-responsive genes, indicating that these genes maintain higher expression levels in cold- tolerant lines than in cold- sensitive lines (Fig. 6 A, Supplemental Table S2 ) (Devireddy et al. 2021 ; Kim et al. 2021 ). This finding suggests that the high expression level of JA-related genes may reflect the degree of cold tolerance in each line. Low temperature–induced changes in expression: common DEGs among lines In each line, we identified 4,117 to 12,191 genes on the microarrays whose expression levels changed in either direction by more than 2-fold upon cold treatment, defining a set of DEGs and DERs. The number of up-regulated and down-regulated (Ure and Dre) probes in each line was 1,897–4,906 and 1,942–5,625, respectively, with no remarkable difference between the tolerant and sensitive lines. The low-temperature tolerance of rice emerged from the analysis of DEGs: We detected four genes annotated “regulation of response to salt stress” (GO:1901000) among Ure genes in more than 11 lines (Fig. 6 B, Supplemental Table S6), possibly reflecting a common pathway involved in responses to salinity and cold stresses (Devireddy et al. 2021 ; Hu et al. 2008 ; Kim et al. 2021 ; Nakashima et al. 2012 ; Ren et al. 2023 ). We also identified eight genes related to chromosome segregation, such as “DNA packaging complex” (GO:0044815), among Dre genes; they may be related to diminished cell division due to low temperatures (Supplemental Table S7) (Qi and Zhang 2020 ). On a gene-by-gene basis, we identified 20 and 24 Ure and Dre genes, respectively, common to all 13 lines (Tables 2, 3, Supplemental Fig. S6). The 20 Ure genes included OsNAC5 and ONAC088 , NAC-type transcription factor genes involved in stress responses including low temperature (Table 2, Supplemental Fig. S6A) (Nakashima et al. 2012 ; Sun et al. 2015 ; Takasaki et al. 2010 ). The inclusion of known low temperature–responsive genes among the Ure genes validates our detection of cold stress responses in anthers. The 24 Dre genes did not contain any known low temperature–responsive genes but did contain genes encoding multiple histones: Os08g0427700 ( H2A ), Os04g0583600 ( H4 ), and Os08g0490900 ( H2B ) (Table 3, Supplemental Fig. S6B). In addition to known DEGs in response to low temperature, we also discovered genes whose expression levels respond to low temperature in a similar manner among lines with different genetic backgrounds. One of the 20 Ure genes, Os06g0246500, whose expression positively correlated with the PFI, encodes a pyruvate dehydrogenase E1 alpha subunit-like protein, which has not been reported to be associated with cold tolerance (Supplemental Fig. S7). Complexity of cold tolerance: expression patterns of known cold tolerance genes across multiple genetic backgrounds We looked at the expression levels of 12 genes reported to be associated with cold tolerance at the booting stage in rice (Table 4). OsLEA9 , OsMKKK70 , and OsMYB4 all showed increased expression in at least nine lines, confirming their association with cold responsiveness reported in previous studies (Lou et al. 2022 ; Mei et al. 2022 ; Park et al. 2010 ). Of the 12 genes, OsMYB4 and OsAPX1 promote cold tolerance when overexpressed (Park et al. 2010 ; Sato et al. 2011 ). However, the expression levels or changes in expression between control and cold stress conditions for OsMYB4 and OsAPX1 showed no significant correlation with the PFI in this study. This result indicates the complexity of cold stress responses during the booting stage. Most previous studies on stress response have focused on a few loci, revealed from comparisons among a small number of lines with contrasting cold tolerance or on a few genes conferring a strong stress response (Mehrotra et al. 2020 ; Raj and Nadarajah 2023 ). Thus, the extent to which genes involved in stress tolerance can explain differences in diverse degrees of cold tolerance among a large number of lines with distinct genetic backgrounds has not been tested. Here, the expression levels for most known cold tolerance genes did not correlate with the PFI. Thus, individual genes cannot universally explain differences in cold tolerance during the booting stage of rice because the factors contributing to cold tolerance vary from line to line (Table 5). Cold tolerance and changes in expression of repetitive sequences In this study, we performed a correlation analysis between the extent of change in expression for each probe and the PFI; we then classified genes into one of four types based on their ECI and r - values (Fig. 5 B). While we generally identified more DEGs than DERs, this was not the case for type I genes, with 137 genes and 463 RSs (Fig. 5 B). This result indicates a genome-wide trend of increased expression of more RSs in the cold-sensitive lines (Table 5). TEs, which constitute the majority of RSs, induce mutations in the host genome through transposition when they are activated, and their expression is normally suppressed by epigenetic mechanisms such as DNA methylation and histone modifications (Choi and Lee 2020 ; Parker et al. 2022 ). However, plants exposed to abiotic stresses such as high or low temperatures experience a weakened repression of RSs, leading to their increased expression (Hu et al. 2012 ; Ito et al. 2011 ; Zhou et al. 2017 ). Therefore, the above trend suggests that the expression of RSs is more suppressed upon cold stress conditions in cold-tolerant lines than in cold-sensitive lines. A previous transcriptome analysis of cold-treated anthers at the booting stage in rice showed that fewer RSs fluctuated in their expression under cold stress in cold-tolerant lines than in cold-sensitive lines (Ishiguro et al. 2014 ). The present study and the previous transcriptome study (Ishiguro et al., 2014 ) thus describe a genome-wide trend whereby cold-tolerant lines are less responsive to cold stress (Table 5). The expression levels of known cold tolerance genes did not correlate with the cold tolerance of the 13 lines tested in this study. Conventional stress-related studies, which focus on a small number of genes from pairwise comparisons of specific lines, therefore cannot fully explain the diversity in cold tolerance responses observed among multiple lines. By contrast, this study focused on genome-wide responses to low-temperature stress and revealed that the transcriptome of more cold-tolerant lines is relatively insensitive to low-temperature stress. Focusing on this feature should be a useful new approach to comprehensively elucidate the differences in stress tolerance intensity among lines. In this study, we focused on genome-wide responsiveness to stress and successfully explained the difference between stress-tolerant and -sensitive plants. Materials & Methods Plant materials Thirteen rice lines were used in this study: eight temperate japonica lines, ‘Hokkai PL9’ (PL9), ‘Kokushokuine 2 go’ (A58), ‘Koshihikari’, ‘Hokkai 287’, Nipponbare, ‘Kirara 397’, Sasanishiki, and ‘Fukoku’; one tropical japonica line (Lambayeque1); one indica line (‘Kasalath’); and three recombinant inbred lines (RILs) ‘R45’, R75, and R85, which were derived from a cross between the japonica line A58 and the indica line ‘Surjamkhi’ (I33) (Table 1). Materials were grown with 20 individuals sown in one Wagner pot (1/5,000 100 m 2 ) in a greenhouse under 25°C (day)/19°C (night) temperature cycles with only the main stem maintained by cutting off the offshoots. For the cold treatment, plants were exposed to 12°C for 4 days during the booting stage, when the microspores are between the late tetrad stage and the uninucleated pollen stage after meiosis. The booting stage was determined by auricle distance (AD) as in our previous study, with AD = − 2 to − 4 cm for Hokkaido and AD = 0–2 cm for other lines (Yamamori et al. 2021 ). Anthers were collected from control plants and plants exposed to cold treatment at the same growth stage and immediately stored at − 80°C. Anthers were collected from the third, fourth, and fifth spikelets of primary rachis branches. Samples for microarrays were anthers collected from multiple individuals for each line and treatment, with two samples per line. Samples for RNA-seq were collected in triplicates using anthers from several individuals for each replicate. To investigate pollen fertility, some individuals were returned to the greenhouse after the cold treatment and allowed to continue growing. Observation of pollen fertility Anthers were collected from the third, fourth, and fifth anthers from the tip of the primary branches during the flowering period. Pollen grains were stained with Lugol's iodine solution [KI–I 2 : 1.5% (w/v) KI, 0.15% (w/v) I 2 ]. Pollen fertility was calculated as the ratio between the number of fertile pollen grains and the total number of pollen grains for each line and under each growth condition. Pollen fertility was determined by sampling at least nine flowers per line, and the average of these was used as the pollen fertility for each line and each growth condition. Microarray probe design In this experiment, a 60K microarray was designed with 23,103 probes derived from repetitive sequences (RSs) and 37,955 derived from protein-coding genes and pre-miRNAs (Supplemental Table S1 ). The probes derived from RSs were selected from the 44K microarray designed by Ishiguro et al. ( 2014 ), taking probes with a fluorescence signal value of 100 or higher. All probes derived from protein-coding genes and pre-miRNAs were from the Agilent catalog array (G2519F) (Agilent technologies, California, USA). Microarray analysis Total RNA used for microarray analysis was extracted from the above samples and purified using a TRIzol Plus RNA purification kit (Life Technologies, California, USA) after freezing and pulverization in liquid nitrogen. RNA concentration was measured with a Nano Drop ND-2000 (Thermo Fisher Scientific, Massachusetts, USA) and quality-checked with an Agilent 2100 Bioanalyzer (Agilent Technologies, California, USA). The microarray analysis was performed with One-Color Spike-Mix containing cyanine-3 (Cy3)-CTP (Agilent) added to the RNA samples, which were labeled using a Quick Amp Labeling Kit (Agilent). After fragmentation of 600 ng complementary RNA (cRNA) synthesized from RNA samples, hybridization to the microarray was performed at 65°C for 17 h using the Agilent Gene Expression Hybridization Kit (Agilent). The slides were then washed, and fluorescence was measured using an Agilent Technologies C version scanner. Agilent Feature Extraction software (Agilent) was used to quantify the fluorescence signals. The data from all samples were log 2 -transformed using R (4.1.0) and normalized by the quantile normalization method using the limma (3.48.3) package in R (Ritchie et al. 2015 ). RNA-seq Total RNA was extracted from anthers as described for microarray analysis; three replicates were used per genotype. The extracted RNA was subjected to strand-specific library preparation (dUTP method) using the NEBNext Ultra™ ll Directional RNA Library Prep Kit (New England Biolabs, Massachusetts, USA). After library preparation, libraries were sequenced as paired-end 150-bp reads on a NovaSeq 6000 instrument (Illumina, California, USA). Library preparation and sequencing were performed by Rhelixa (Tokyo, Japan). From the data obtained by sequencing, kallisto (0.46.2) was used to obtain expression levels of each gene (Bray et al. 2016 ). IRGSP-1.0 (2021-05-10), published by RAP-DB, was used as a reference for RNA-seq analysis (Kawahara et al. 2013 ; Sakai et al. 2013 ). The bootstrap value was set to 100, and the expression levels were normalized to transcripts per million (TPM) using kallisto. Sleuth (0.30.0) was used to compare expression levels between cold-treated and control plants for each line and between lines (Pimentel et al. 2017 ). Differences in expression levels between treatments in each line were statistically tested using the likelihood ratio test. Gene annotation and GO analysis The genes selected in this study were annotated based on IRGSP-1.0 (2021-05-10). GO analysis of the selected genes was performed using PANTHER ( http://geneontology.org/ ) provided by The Gene Ontology Consortium (Mi et al. 2019 ). Fisher's exact test was used to identify GO terms that were significantly (FDR < 0.05) more abundant in the selected genes than across the genome. Conclusion We investigated the relationship between genome-wide response to low-temperature stress and cold stress tolerance based on transcriptome analysis of 13 rice lines with varying degrees of cold tolerance. We interpreted the results of these analyses from three perspectives: gene expression levels and cold tolerance, changes in expression in response to cold, and changes in expression as a function of cold tolerance. The relationship between expression levels and cold tolerance revealed that, compared to cold-sensitive lines, the cold-tolerant lines tested here tended to have more genes with lower expression than with higher expression, regardless of growth conditions. Some of the genes with lower expression in cold-tolerant lines were related to ribosomes, suggesting that transcription and translation of the genes are stabilized in the cold-tolerant lines. Looking at genes with altered expression in response to cold, we identified genes with common differential expression across the 13 lines tested. We thus identified new genes that are universally responsive to cold in rice anthers. Many of these genes had not previously been reported to be associated with cold stress. We finally revealed that the RSs whose changes in expression levels showed significant correlations with cold tolerance most often exhibited increased expression in the cold-sensitive lines but remained largely constant in cold-tolerant lines. This result is consistent with our previous study (Ishiguro et al. 2014 ) and reaffirms the genomic stability under low-temperature stress displayed by cold-tolerant lines. In summary, this study provides new insights into cold stress responses during the booting stage by focusing on the relationship between the various cold tolerance intensities exhibited by the lines and genome-wide expression changes. Abbreviations DEG differentially expressed gene DER differentially expressed repetitive sequences Dre down-regulated ECI expression change index PFI pollen fertility index RS repetitive sequence TPM transcripts per million Ure up-regulated Declarations Ethical Approval and Consent to participate This study complied with the ethical standards of Japan, Hokkaido University, where this research was conducted. Consent for publication All authors have consented to the publication of this manuscript. Availability of supporting data Our microarray data were recorded as GSExxx in the National Center for Biotechnology Information (NCBI) provided by http://www. ncbi.nlm.nih.gov/yyyy. Our RNA-seq data were recorded as GSExxx in the NCBI provided by http://www. ncbi.nlm.nih.gov/yyyy . Competing interests The authors declare that they have no competing interests. Funding This work was supported by JSPS KAKENHI grants (19H00937 and 23H02180) to Y.Ki and Scientific Technique Research Promotion Program for Agriculture, Forestry, Fisheries and Food Industry to Y.S. Author contributions K.Y., Y.S. and Y.K. designed the research plans; S.I. and K.O. performed the cold stress evaluations and the microarray experiments; K.Y. and K.M.L. performed the cold stress evaluations and the RNA-seq experiments; K.Y. and K.F. performed the data analysis; K.F. and Y.S. provided instruments for the experiments; K.O. produced RIL lines; K.Y. and Y.K. contributed to manuscript preparation; K.Y. and Y.K. wrote the article with contributions from all the authors. Acknowledgments We gratefully acknowledge Dr. Y. Koide (Research Faculty of Agriculture, Hokkaido University) for his valuable suggestions concerning this study. References Andaya VC, Mackill DJ (2003) QTLs conferring cold tolerance at the booting stage of rice using recombinant inbred lines from a japonica x indica cross. Theor Appl Genet 106:1084-1090 Bechtold U, Field B (2018) Molecular mechanisms controlling plant growth during abiotic stress. J Exp Bot 69:2753-2758 Bolstad BM, Irizarry RA, Astrand M, Speed TP (2003) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19:185-193 Bray NL, Pimentel H, Melsted P, Pachter L (2016) Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 34:525-527 Chauhan H, Khurana N, Agarwal P, Khurana P (2011) Heat shock factors in rice (Oryza sativa L.): genome-wide expression analysis during reproductive development and abiotic stress. Mol Genet Genomics 286:171-187 Choi JY, Lee YCG (2020) Double-edged sword: The evolutionary consequences of the epigenetic silencing of transposable elements. PLoS Genet 16:e1008872 Dai LY, Lin XH, Ye CR, Ise KZ, Saito K, Kato A, Xu FR, Yu TQ, Zhang DP (2004) Identification of quantitative trait loci controlling cold tolerance at the reproductive stage in Yunnan landrace of rice, Kunmingxiaobaigu. Breeding Science 54:253-258 Devireddy AR, Zandalinas SI, Fichman Y, Mittler R (2021) Integration of reactive oxygen species and hormone signaling during abiotic stress. Plant J 105:459-476 Dias-Fields L, Adamala KP (2022) Engineering ribosomes to alleviate abiotic stress in plants: A perspective. Plants (Basel) 11:2097 Fang Y, Xie K, Xiong L (2014) Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot 65:2119-2135 Guo M, Liu JH, Ma X, Luo DX, Gong ZH, Lu MH (2016) The plant heat stress transcription factors (HSFs): Structure, regulation, and function in response to abiotic stresses. Front Plant Sci 7:114 Hayase H, Satake T, Nishiyama I, Ito N (1969) Male sterility caused by cooling treatment at the meiotic stage in rice plants: 2. The most sensitive stage to cooling and the fertilizing ability of pistils. Japanese Journal of Crop Science 38:706-711 He M, Ding NZ (2020) Plant unsaturated fatty acids: multiple roles in stress response. Front Plant Sci 11:562785 Hong Y, Zhang H, Huang L, Li D, Song F (2016) Overexpression of a stress-responsive NAC transcription factor gene ONAC022 improves drought and salt tolerance in rice. Front Plant Sci 7:4 Hu H, You J, Fang Y, Zhu X, Qi Z, Xiong L (2008) Characterization of transcription factor gene SNAC2 conferring cold and salt tolerance in rice. Plant Mol Biol 67:169-181 Hu Y, Zhang L, He S, Huang M, Tan J, Zhao L, Yan S, Li H, Zhou K, Liang Y, Li L (2012) Cold stress selectively unsilences tandem repeats in heterochromatin associated with accumulation of H3K9ac. Plant Cell Environ 35:2130-2142 Ishiguro S, Ogasawara K, Fujino K, Sato Y, Kishima Y (2014) Low temperature-responsive changes in the anther transcriptome's repeat sequences are indicative of stress sensitivity and pollen sterility in rice strains. Plant Physiol 164:671-682 Ito H, Gaubert H, Bucher E, Mirouze M, Vaillant I, Paszkowski J (2011) An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress. Nature 472:115-119 Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S, Childs KL, Davidson RM, Lin H, Quesada-Ocampo L, Vaillancourt B, Sakai H, Lee SS, Kim J, Numa H, Itoh T, Buell CR, Matsumoto T (2013) Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice (N Y) 6:4 Kidokoro S, Hayashi K, Haraguchi H, Ishikawa T, Soma F, Konoura I, Toda S, Mizoi J, Suzuki T, Shinozaki K, Yamaguchi-Shinozaki K (2021) Posttranslational regulation of multiple clock-related transcription factors triggers cold-inducible gene expression in Arabidopsis. Proc Natl Acad Sci U S A 118:e2021048118 Kidokoro S, Shinozaki K, Yamaguchi-Shinozaki K (2022) Transcriptional regulatory network of plant cold-stress responses. Trends Plant Sci 27:922-935 Kim H, Seomun S, Yoon Y, Jang G (2021) Jasmonic acid in plant abiotic stress tolerance and interaction with abscisic acid. Agronomy-Basel 11:1886 Kim TH, Kim SM (2023) Identification of candidate genes for salt tolerance at the seedling stage using integrated genome-wide association study and transcriptome analysis in rice. Plants-Basel 12:1401 Kuroki M, Saito K, Matsuba S, Yokogami N, Shimizu H, Ando I, Sato Y (2007) A quantitative trait locus for cold tolerance at the booting stage on rice chromosome 8. Theor Appl Genet 115:593-600 Li HB, Wang J, Liu AM, Liu KD, Zhang Q, Zou JS (1997) Genetic basis of low-temperature-sensitive sterility in indica-japonica hybrids of rice as determined by RFLP analysis. Theor Appl Genet 95:1092-1097 Li J, Zeng Y, Pan Y, Zhou L, Zhang Z, Guo H, Lou Q, Shui G, Huang H, Tian H, Guo Y, Yuan P, Yang H, Pan G, Wang R, Zhang H, Yang S, Guo Y, Ge S, Li J, Li Z (2021) Stepwise selection of natural variations at CTB2 and CTB4a improves cold adaptation during domestication of japonica rice. New Phytol 231:1056-1072 Liu C, Ou S, Mao B, Tang J, Wang W, Wang H, Cao S, Schlappi MR, Zhao B, Xiao G, Wang X, Chu C (2018) Early selection of bZIP73 facilitated adaptation of japonica rice to cold climates. Nat Commun 9:3302 Lou Q, Guo H, Li J, Han S, Khan NU, Gu Y, Zhao W, Zhang Z, Zhang H, Li Z, Li J (2022) Cold-adaptive evolution at the reproductive stage in Geng/japonica subspecies reveals the role of OsMAPK3 and OsLEA9. Plant J 111:1032-1051 Mehrotra S, Verma S, Kumar S, Kumari S, Mishra BN (2020) Transcriptional regulation and signalling of cold stress response in plants: An overview of current understanding. Environmental and Experimental Botany 180 Mei EY, Tang JQ, He ML, Liu ZQ, Tian XJ, Bu QY (2022) OsMKKK70 negatively regulates cold tolerance at booting stage in rice. IJMS 23:14472 Mi HY, Muruganujan A, Ebert D, Huang XS, Thomas PD (2019) PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res 47:D419-D426 Nakashima K, Takasaki H, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K (2012) NAC transcription factors in plant abiotic stress responses. Bba-Gene Regul Mech 1819:97-103 Nakashima K, Tran LSP, Van Nguyen D, Fujita M, Maruyama K, Todaka D, Ito Y, Hayashi N, Shinozaki K, Yamaguchi-Shinozaki K (2007) Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J 51:617-630 Oh CS, Choi YH, Lee SJ, Yoon DB, Moon HP, Ahn SN (2004) Mapping of quantitative trait loci for cold tolerance in weedy rice. Breeding Science 54:373-380 Park MR, Yun KY, Mohanty B, Herath V, Xu FY, Wijaya E, Bajic VB, Yun SJ, De Los Reyes BG (2010) Supra-optimal expression of the cold-regulated transcription factor in transgenic rice changes the complexity of transcriptional network with major effects on stress tolerance and panicle development. Plant Cell Environ 33:2209-2230 Parker AH, Wilkinson SW, Ton J (2022) Epigenetics: a catalyst of plant immunity against pathogens. New Phytol 233:66-83 Pimentel H, Bray NL, Puente S, Melsted P, Pachter L (2017) Differential analysis of RNA-seq incorporating quantification uncertainty. Nature Methods 14:687-690 Qi FF, Zhang FX (2020) Cell cycle regulation in the plant response to stress. Front Plant Sci 10:1765 Raj SRG, Nadarajah K (2023) QTL and candidate genes: Techniques and advancement in abiotic stress resistance breeding of major cereals. Int J Mol Sci 24:6 Raza A, Charagh S, Zahid Z, Mubarik MS, Javed R, Siddiqui MH, Hasanuzzaman M (2021) Jasmonic acid: a key frontier in conferring abiotic stress tolerance in plants. Plant Cell Rep 40:1513-1541 Ren HM, Zhang YT, Zhong MY, Hussian J, Tang YT, Liu SK, Qi GN (2023) Calcium signaling-mediated transcriptional reprogramming during abiotic stress response in plants. Theor Appl Genet 136 Renák D, Gibalová A, Solcová K, Honys D (2014) A new link between stress response and nucleolar function during pollen development in mediated by AtREN1 protein. Plant Cell Environ 37:670-683 Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015) limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43:e47 Sagehashi Y, Ikegaya T, Fujino K (2022) Integration of genetic engineering into conventional rice breeding programs for the next generation. Euphytica 218:145 Sagor GHM, Inoue M, Kusano T, Berberich T (2021) Expression profile of seven polyamine oxidase genes in rice (Oryza sativa) in response to abiotic stresses, phytohormones and polyamines. Physiol Mol Biol Plants 27:1353-1359 Saito K, Hayano-Saito Y, Kuroki M, Sato Y (2010) Map-based cloning of the rice cold tolerance gene. Plant Science 179:97-102 Saito K, Miura K, Nagano K, Hayano-Saito Y, Araki H, Kato A (2001) Identification of two closely linked quantitative trait loci for cold tolerance on chromosome 4 of rice and their association with anther length. Theor Appl Genet 103:862-868 Sakai H, Lee SS, Tanaka T, Numa H, Kim J, Kawahara Y, Wakimoto H, Yang CC, Iwamoto M, Abe T, Yamada Y, Muto A, Inokuchi H, Ikemura T, Matsumoto T, Sasaki T, Itoh T (2013) Rice Annotation Project Database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol 54:e6 Satake T (1976) Determination of the most sensitive stage to sterile type cool injury in rice plants. Research bulletin of the Hokkaido National Agricultural Experiment Station 113:1-43 Sato Y, Masuta Y, Saito K, Murayama S, Ozawa K (2011) Enhanced chilling tolerance at the booting stage in rice by transgenic overexpression of the ascorbate peroxidase gene, OsAPXa. Plant Cell Rep 30:399-406 Shimono H, Abe A, Aoki N, Koumoto T, Sato M, Yokoi S, Kuroda E, Endo T, Saeki KI, Nagano K (2016) Combining mapping of physiological quantitative trait loci and transcriptome for cold tolerance for counteracting male sterility induced by low temperatures during reproductive stage in rice. Physiol Plant 157:175-192 Shirasawa S, Endo T, Nakagomi K, Yamaguchi M, Nishio T (2012) Delimitation of a QTL region controlling cold tolerance at booting stage of a cultivar, 'Lijiangxintuanheigu', in rice, Oryza sativa L. Theor Appl Genet 124:937-946 Suh JP, Jeung JU, Lee JI, Choi YH, Yea JD, Virk PS, Mackill DJ, Jena KK (2010) Identification and analysis of QTLs controlling cold tolerance at the reproductive stage and validation of effective QTLs in cold-tolerant genotypes of rice (Oryza sativa L.). Theor Appl Genet 120:985-995 Sun L, Huang L, Hong Y, Zhang H, Song F, Li D (2015) Comprehensive analysis suggests overlapping expression of rice ONAC transcription factors in abiotic and biotic stress responses. Int J Mol Sci 16:4306-4326 Sun ZH, Du J, Pu XY, Ali MK, Yang XM, Duan CL, Ren MR, Li X, Zeng YW (2019) Near-isogenic lines of rice revealed new QTLs for cold tolerance at booting stage. Agronomy-Basel 9:40 Takasaki H, Maruyama K, Kidokoro S, Ito Y, Fujita Y, Shinozaki K, Yamaguchi-Shinozaki K, Nakashima K (2010) The abiotic stress-responsive NAC-type transcription factor OsNAC5 regulates stress-inducible genes and stress tolerance in rice. Mol Genet Genomics 284:173-183 Takeuchi Y, Hayasaka H, Chiba B, Tanaka I, Shimano T, Yamagishi M, Nagano K, Sasaki T, Yano M (2001) Mapping quantitative trait loci controlling cool-temperature tolerance at booting stage in temperate rice. Breeding Science 51:191-197 Xu LM, Zhou L, Zeng YW, Wang FM, Zhang HL, Shen SQ, Li ZC (2008) Identification and mapping of quantitative trait loci for cold tolerance at the booting stage in a rice near-isogenic line. Plant Science 174:340-347 Yamamori K, Ogasawara K, Ishiguro S, Koide Y, Takamure I, Fujino K, Sato Y, Kishima Y (2021) Revision of the relationship between anther morphology and pollen sterility by cold stress at the booting stage in rice. Annals of Botany 128:559-575 Yang L, Lei L, Wang J, Zheng H, Xin W, Liu H, Zou D (2023) qCTB7 positively regulates cold tolerance at booting stage in rice. Theor Appl Genet 136:135 Zandalinas SI, Balfagon D, Gomez-Cadenas A, Mittler R (2022) Plant responses to climate change: metabolic changes under combined abiotic stresses. J Exp Bot 73:3339-3354 Zhang H, Zhu J, Gong Z, Zhu JK (2022a) Abiotic stress responses in plants. Nat Rev Genet 23:104-119 Zhang M, Zhao R, Huang K, Huang S, Wang H, Wei Z, Li Z, Bian M, Jiang W, Wu T, Du X (2022b) The OsWRKY63-OsWRKY76-OsDREB1B module regulates chilling tolerance in rice. Plant J 112:383-398 Zhang Y, Xia P (2023) The DREB transcription factor, a biomacromolecule, responds to abiotic stress by regulating the expression of stress-related genes. Int J Biol Macromol 243:125231 Zhang Z, Li J, Pan Y, Li J, Zhou L, Shi H, Zeng Y, Guo H, Yang S, Zheng W, Yu J, Sun X, Li G, Ding Y, Ma L, Shen S, Dai L, Zhang H, Yang S, Guo Y, Li Z (2017) Natural variation in CTB4a enhances rice adaptation to cold habitats. Nat Commun 8:14788 Zhou H, Hirata M, Osawa R, Fujino K, Kishima Y (2017) Detainment of Tam3 transposase at plasma membrane by its BED-Zinc finger domain. Plant Physiol 173:1492-1501 Zhu Y, Chen K, Mi X, Chen T, Ali J, Ye G, Xu J, Li Z (2015) Identification and fine mapping of a stably expressed QTL for cold tolerance at the booting stage using an interconnected breeding population in rice. PLoS One 10:e0145704 Tables Tables 1 to 5 are available in the Supplementary Files section. Additional Declarations No competing interests reported. Supplementary Files Tables.docx YamamorietalSupplementaryFigures2024.4.17.pptx YamamorietalSupplementaryTables2024.4.17.xlsx SupplementalFigLegends.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4399503","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":304420299,"identity":"e09fe226-9c70-4070-a3a6-ab71e190cdb2","order_by":0,"name":"Koichi Yamamori","email":"","orcid":"","institution":"Hokkaido University","correspondingAuthor":false,"prefix":"","firstName":"Koichi","middleName":"","lastName":"Yamamori","suffix":""},{"id":304420300,"identity":"6557dfc3-7dcd-4493-8985-cf30efeff016","order_by":1,"name":"Seiya Ishiguro","email":"","orcid":"","institution":"Hokkaido University","correspondingAuthor":false,"prefix":"","firstName":"Seiya","middleName":"","lastName":"Ishiguro","suffix":""},{"id":304420301,"identity":"ff2c8eea-a581-4515-b5d3-3530e154da85","order_by":2,"name":"Kei Ogasawara","email":"","orcid":"","institution":"Hokkaido University","correspondingAuthor":false,"prefix":"","firstName":"Kei","middleName":"","lastName":"Ogasawara","suffix":""},{"id":304420302,"identity":"11ff828c-c262-4f13-82ca-c08578e186bf","order_by":3,"name":"Kayyis Lubba","email":"","orcid":"","institution":"Hokkaido University","correspondingAuthor":false,"prefix":"","firstName":"Kayyis","middleName":"","lastName":"Lubba","suffix":""},{"id":304420303,"identity":"e98b372e-20ec-4fb2-b4b2-ed67365de0eb","order_by":4,"name":"Kaien Fujino","email":"","orcid":"","institution":"Hokkaido University","correspondingAuthor":false,"prefix":"","firstName":"Kaien","middleName":"","lastName":"Fujino","suffix":""},{"id":304420304,"identity":"5643f356-418e-4a4c-9a3f-3c6f9f502ef3","order_by":5,"name":"Kazumitsu Onishi","email":"","orcid":"","institution":"Obihiro University of Agriculture and Veterinary Medicine","correspondingAuthor":false,"prefix":"","firstName":"Kazumitsu","middleName":"","lastName":"Onishi","suffix":""},{"id":304420305,"identity":"b65374b9-2248-4605-bf6c-f1894ba93d4e","order_by":6,"name":"Yutaka Sato","email":"","orcid":"","institution":"National Agriculture and Food Research Organization","correspondingAuthor":false,"prefix":"","firstName":"Yutaka","middleName":"","lastName":"Sato","suffix":""},{"id":304420306,"identity":"84637057-d71c-4cd0-b529-746a2f450da9","order_by":7,"name":"Yuji Kishima","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABEElEQVRIiWNgGAWjYBACCWYGNoYEBgY5EIeZAcQ8wNgA5eDXYkyCFgagFiBIbEBoIeAwyXb2Zw8e5til9/evfcBcUJGW2Hf8cAPDjxoGdnMcWqSZecwNErcl58648dyAecaZnMSZZxIbGHuOMTBbNmDXIsfMwyaRuI05t+EGUBVvW0XihgNAR/I2MDAb4HChHDP7M6CW+nR5uJbzDxsY/+LRIs3MYAbUcjjB4HwbSEtO4oYbiQ3M+GyRbOYBaTluuPEGG8NhnjNpxjNvPGw4LHNMAqdfJM4ffyb5c1u1vNz5Y4yPeSqSZfvOpz98+KbGJhlXiCFpTkDECJAhkWxAUAs/mtPtCGsZBaNgFIyCEQIAJhpb5DeIXjkAAAAASUVORK5CYII=","orcid":"","institution":"Hokkaido University","correspondingAuthor":true,"prefix":"","firstName":"Yuji","middleName":"","lastName":"Kishima","suffix":""}],"badges":[],"createdAt":"2024-05-10 09:01:22","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4399503/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4399503/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":56936849,"identity":"f7863970-32d5-4705-929c-091e45b377f5","added_by":"auto","created_at":"2024-05-22 11:20:11","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":215435,"visible":true,"origin":"","legend":"\u003cp\u003eEvaluation of pollen fertility in 13 rice lines.\u003c/p\u003e\n\u003cp\u003ePollen fertility was calculated as the ratio between the number of fertile pollen grains and the total number of pollen grains for each line and under each growth condition.Values are means ± standard error of the mean (SEM); \u003cem\u003en\u003c/em\u003e = 9 flowers per line.\u003c/p\u003e","description":"","filename":"YamamorietalFigure1.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/ba54baa78a8914642976d297.png"},{"id":56937642,"identity":"845140f6-1acd-4c47-bfcc-ca1c73bbc561","added_by":"auto","created_at":"2024-05-22 11:28:12","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":595209,"visible":true,"origin":"","legend":"\u003cp\u003eExtent of correlation between the anther microarray data from 13 rice lines subjected to a 12°C low-temperature treatment and control plants.\u003c/p\u003e\n\u003cp\u003eFor each line, the \u003cem\u003ex\u003c/em\u003e-axis shows the log\u003csub\u003e2\u003c/sub\u003e-normalized expression values from control plants; the \u003cem\u003ey\u003c/em\u003e-axis shows the log\u003csub\u003e2\u003c/sub\u003e-normalized expression values from cold-treated plants. The coefficient of determination (\u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e) value in each plot indicates the degree of similarity between the anther transcriptomes from the cold-treated and control plants. The dotted lines indicate \u003cem\u003ey\u003c/em\u003e = \u003cem\u003ex\u003c/em\u003e ± 1.\u003c/p\u003e","description":"","filename":"YamamorietalFigure2.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/1954e826e109d9c121868067.png"},{"id":56938342,"identity":"a467e3d4-8281-41cf-857f-3bd31cae5711","added_by":"auto","created_at":"2024-05-22 11:36:11","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":534664,"visible":true,"origin":"","legend":"\u003cp\u003eCorrelation between cold tolerance and expression levels based on microarrays for the 13 lines.\u003c/p\u003e\n\u003cp\u003e(A) Distribution of correlation values between the extent of cold tolerance, estimated by the pollen fertility index (PFI), and expression levels of plants from each of the 13 lines grown under control conditions. (B) Heatmap representation of probe expression levels (genes and repetitive sequences (RSs)) in plants grown under control conditions. The rice lines were sorted according to increasing PFI values. Under control conditions, the expression levels of 367 genes and 301 RSs were positively correlated (\u003cem\u003er\u003c/em\u003e ≥ 0.55) with PFI; the expression levels of 633 genes and 91 RSs were negatively correlated (\u003cem\u003er\u003c/em\u003e ≤ −0.55) with PFI. The\u003cem\u003e P\u003c/em\u003e-value was calculated by Fisher's test. (C) Distribution of correlation values between the extent of tolerance, estimated by the PFI, and expression levels of plants from each of the 13 lines grown in cold-treated conditions. (D) Heatmap representation of probe expression levels in plants grown under cold-treated conditions. Under cold treatment, the expression levels of 169 genes and 65 RSs were positively correlated (\u003cem\u003er\u003c/em\u003e ≥ 0.55) with PFI; the expression levels of 479 genes and 51 RSs were negatively correlated (\u003cem\u003er\u003c/em\u003e ≤ −0.55) with PFI. The\u003cem\u003e P\u003c/em\u003e-value was calculated by Fisher's test.\u003c/p\u003e","description":"","filename":"YamamorietalFigure3.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/11641cc01033b289faf52520.png"},{"id":56938343,"identity":"629274ae-e221-4c4e-9ff9-65bcd2fb2cbf","added_by":"auto","created_at":"2024-05-22 11:36:12","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":135096,"visible":true,"origin":"","legend":"\u003cp\u003eIdentification of cold-responsive differentially expressed genes and repetitive sequences in the 13 rice lines.\u003c/p\u003e\n\u003cp\u003e(A) Number of probes identified as up-regulated (Ure) or down-regulated (Dre) in each of the 13 lines. Probes corresponding to genes are shown in dark colors; probes corresponding to repetitive sequences (RSs) are shown in lighter colors. (B) Number of genes (left) and RSs (right) with differential expression in common between 2 to 13 lines. The \u003cem\u003ey\u003c/em\u003e-axis represents the number of lines (2 to 13) with the indicated numbers of common expression changes. The insets are magnified views for genes or RSs identified as differentially expressed in 10 or more lines.\u003c/p\u003e","description":"","filename":"YamamorietalFigure4.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/72bc2fb0bd610a93e969a99b.png"},{"id":56936857,"identity":"b1ced001-6d4f-4a05-a24a-3fd0f2ec3e78","added_by":"auto","created_at":"2024-05-22 11:20:12","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":359965,"visible":true,"origin":"","legend":"\u003cp\u003eCorrelation between changes in expression levels measured by microarrays and cold tolerance of the 13 lines.\u003c/p\u003e\n\u003cp\u003e(A) Distribution of correlation values between cold tolerance, estimated by the pollen fertility index (PFI) of each line, and the change in expression of each probe due to cold stress across the 13 lines. Red bars indicate significant positive correlations (\u003cem\u003er\u003c/em\u003e ≥ 0.55, \u003cem\u003eP\u003c/em\u003e \u0026lt; 0.05). Blue bars indicate significant negative correlations (\u003cem\u003er\u003c/em\u003e ≤ −0.55, \u003cem\u003eP\u003c/em\u003e \u0026lt; 0.05). (B) Scatterplots showing the relationships between the expression change index (ECI) and their correlation values (as shown in A) for genes (left) or RSs (right). The four types of probes were selected based on the following criteria: type 1, negative correlation (\u003cem\u003er\u003c/em\u003e ≤ −0.55) and ECI \u0026gt; 0; type 2, positive correlation (\u003cem\u003er\u003c/em\u003e ≥ 0.55) and ECI \u0026gt; 0; type 3, negative correlation (\u003cem\u003er\u003c/em\u003e ≤ −0.55) and ECI \u0026lt; 0; and type 4, positive correlation (\u003cem\u003er\u003c/em\u003e ≥ 0.55) and ECI \u0026lt; 0. The number of each type of probe is indicated in parentheses. (C) Examples of genes from each probe type showing the relationship between the PFI and expression changes (log\u003csub\u003e2\u003c/sub\u003e [FC]).\u003c/p\u003e","description":"","filename":"YamamorietalFigure5.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/a98a6578714ad9fa79dda370.png"},{"id":56936859,"identity":"b5249f3d-6015-43aa-9e50-3e9af9aa3a59","added_by":"auto","created_at":"2024-05-22 11:20:12","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":258840,"visible":true,"origin":"","legend":"\u003cp\u003eSignificant Gene Ontology terms for genes whose expression is correlated with cold tolerance or differentially expressed in most lines.\u003c/p\u003e\n\u003cp\u003e(A) Gene Ontology (GO) term enrichment analysis for genes showing positive (left) or negative (right) correlations between their expression levels in control conditions and pollen fertility index (PFI). (B) GO term enrichment analysis for the genes detected as up-regulated (Ure, left) or down-regulated (Dre, right) in more than 11 lines. The extracted GO terms were formed with fewer than 100 rice genes and indicated with FDR \u0026lt; 0.001 (dashed lines) (Supplemental Tables S2–S8). Green, GO term “biological process”; gray, “molecular function”; blue, “cellular component.”\u003c/p\u003e","description":"","filename":"YamamorietalFigure6.png","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/f97bf9f8368d2113ce39be9c.png"},{"id":58559587,"identity":"bbf17561-493a-4734-b740-f9ed031944cc","added_by":"auto","created_at":"2024-06-18 08:47:43","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2886500,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/cfa05637-d4ff-49f2-925b-af573bb75041.pdf"},{"id":56936850,"identity":"96843053-d638-4554-ae3e-0227352507eb","added_by":"auto","created_at":"2024-05-22 11:20:11","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":777691,"visible":true,"origin":"","legend":"","description":"","filename":"Tables.docx","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/c97a102e87c4a2e5667e32af.docx"},{"id":56936858,"identity":"4891a055-a5a4-4b54-b29e-dc46fa7154b2","added_by":"auto","created_at":"2024-05-22 11:20:12","extension":"pptx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":4041210,"visible":true,"origin":"","legend":"","description":"","filename":"YamamorietalSupplementaryFigures2024.4.17.pptx","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/b476f75704869848b6de7835.pptx"},{"id":56936853,"identity":"c56ba67d-5daa-4070-bf2c-a9cda6d01b09","added_by":"auto","created_at":"2024-05-22 11:20:12","extension":"xlsx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":56317,"visible":true,"origin":"","legend":"","description":"","filename":"YamamorietalSupplementaryTables2024.4.17.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/ca2fad9b0d661618146e646d.xlsx"},{"id":56937640,"identity":"acbc32af-33d4-49fa-aa7f-5d5e0206da38","added_by":"auto","created_at":"2024-05-22 11:28:11","extension":"docx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":15875,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementalFigLegends.docx","url":"https://assets-eu.researchsquare.com/files/rs-4399503/v1/396f6a901b483fb1319aa79e.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"The anther transcriptome of cold-tolerant rice cultivars is largely insensitive to temperature changes","fulltext":[{"header":"Background","content":"\u003cp\u003eRice (\u003cem\u003eOryza sativa\u003c/em\u003e) is a tropical crop that is susceptible to low temperatures. Low-temperature stress, especially during the booting stage, causes irreversible pollen sterility, resulting in lower yields (Hayase et al. \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e1969\u003c/span\u003e; Satake \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e1976\u003c/span\u003e). Rice, native to the tropics and subtropics, can now be grown at high latitudes above 40\u0026deg;N (Sagehashi et al. \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Acquisition of cold tolerance traits during the booting stage is important for rice acclimation to low temperatures at high latitudes (Li et al. \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Liu et al. \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Genetic analyses have been conducted in various rice cultivars to identify factors that confer cold tolerance at booting stage, and many quantitative trait loci (QTLs) have been identified (Andaya and Mackill \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2003\u003c/span\u003e; Dai et al. \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2004\u003c/span\u003e; Kuroki et al. \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Li et al. \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e1997\u003c/span\u003e; Li et al. \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Oh et al. \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2004\u003c/span\u003e; Saito et al. \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Saito et al. \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Shimono et al. \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Shirasawa et al. \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Suh et al. \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Sun et al. \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Takeuchi et al. \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Xu et al. \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Zhu et al. \u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Typically, these analyses use two lines with contrasting cold tolerance at the booting stage. However, of the 30 QTLs identified to date, only four genes have been cloned: \u003cem\u003eCold tolerance at the booting stage 1\u003c/em\u003e (\u003cem\u003eCTB1\u003c/em\u003e), \u003cem\u003eCTB2\u003c/em\u003e, \u003cem\u003eCTB4a\u003c/em\u003e, and \u003cem\u003eqCTB7\u003c/em\u003e (Li et al. \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Saito et al. \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Yang et al. \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2017\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eNumerous genes are involved in abiotic stress responses, including low temperature, relying on multiple pathways that act independently or interact with one another (Kidokoro et al. \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Kidokoro et al. \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Kuroki et al. \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Liu et al. \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Lou et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Nakashima et al. \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Park et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2022b\u003c/span\u003e). Therefore, we hypothesize that many genes participate in cold tolerance at the booting stage. While the cloning of specific genes is crucial, it would not allow a global understanding of the genetic mechanisms behind cold tolerance at this developmental stage that would apply across many rice genotypes.\u003c/p\u003e \u003cp\u003eIn a previous study, we investigated the relationship between genome-wide cold responsiveness and cold tolerance in anthers at the booting stage by microarray analysis, focusing on repetitive sequences (RSs) (Ishiguro et al. \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Expression of transposable elements (TEs), a major component of RSs, is repressed by epigenetic regulation, thus maintaining genome stability (Choi and Lee \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Parker et al. \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). However, TE repression is released under stress conditions (Hu et al. \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Ito et al. \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Zhou et al. \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). For example, among five rice cultivars differing in their cold tolerance, more tolerant cultivars tended to show fewer changes in RS expression in response to cold stress than more sensitive cultivars (2014). This result suggested that cold-tolerant lines at the booting stage tend to suppress TE activation at lower temperatures and maintain genome stability. To our knowledge, the idea that genome-wide expression changes, rather than those of specific genes, are related to stress response and tolerance has not been previously proposed in stress tolerance research. However, the difficulty of explaining the evocation of cold tolerance at the booting stage from individual loci or common QTLs lends some validity to this claim.\u003c/p\u003e \u003cp\u003eIn this study, to explore rice cold tolerance at the booting stage, we analyzed the relationship between genome-wide responsiveness to cold stress and tolerance using 13 rice lines representing a continuous gradient of cold tolerance. We examined anther transcript levels during the booting stage with 60K microarrays containing probes for about 38,000 genes and 23,000 RSs. We calculated the pollen fertility index (PFI) of each genotype and looked for correlations with genome-wide expression levels, identified differentially expressed genes (DEGs) and RSs (DERs) whose expression was altered by low temperature, and characterized genes whose change in expression correlated with PFI. We propose that genome-wide expression patterns can well explain the differences for cold tolerance at the booting stage among rice lines. In support of this idea, previously cloned cold tolerance genes did not always show high correlations with cold tolerance among rice lines.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003ePollen fertility of 13 rice lines under cold stress and control conditions\u003c/h2\u003e \u003cp\u003eThe 13 rice lines used in this study comprised eight temperate \u003cem\u003ejaponica\u003c/em\u003e, one tropical \u003cem\u003ejaponica\u003c/em\u003e, one \u003cem\u003eindica\u003c/em\u003e, and three recombinant inbred lines (RILs) derived from a cross between a temperate \u003cem\u003ejaponica\u003c/em\u003e cultivar (\u0026lsquo;A58\u0026rsquo;) and an \u003cem\u003eindica\u003c/em\u003e cultivar (\u0026lsquo;I33\u0026rsquo;) (Table\u0026nbsp;1). These lines are genetically diverse and expected to represent various degrees of cold tolerance. We grew all plants under greenhouse conditions (day: 25\u0026deg;C, night: 19\u0026deg;C) before exposing them to continuous low temperature (12\u0026deg;C) for 4 days at the booting stage in the growth chamber. To examine pollen fertility from these plants, we collected the third, fourth, and fifth flowers from the top of each primary rachis branch after heading (Yamamori et al. \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Pollen fertility for the 13 genotypes ranged from 59.7\u0026ndash;92.6% for control plants maintained under 25\u0026deg;C/19\u0026deg;C day/night temperature cycles and from 2.2\u0026ndash;85.6% in the cold-treated plants (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). The line with the lowest pollen fertility under the cold treatment was the RIL \u0026lsquo;R75\u0026rsquo; and the highest was the tropical \u003cem\u003ejaponica\u003c/em\u003e \u0026lsquo;Lambayeque1\u0026rsquo;. We calculated the PFI in response to cold treatment as the ratio between the average pollen fertility under the cold treatment and the average pollen fertility under control conditions. Genotypes with a higher PFI are more cold tolerant.\u003c/p\u003e \u003cp\u003eWe prepared another set of materials to verify the results obtained by microarray analysis and to perform a transcriptome deep-sequencing (RNA-seq) analysis with 4 of the 13 genotypes, choosing Lambayeque1 and \u0026lsquo;PL9\u0026rsquo; as cold-tolerant lines and \u0026lsquo;Sasanishiki\u0026rsquo; and \u0026lsquo;Nipponbare\u0026rsquo; as sensitive lines. For this RNA-seq analysis, we independently grew the four lines under the same conditions as the 13 genotypes used for microarray analysis and collected samples in triplicate. The pollen fertility of the cold-treated materials was 3.9% (Nipponbare), 24.8% (Sasanishiki), 72.1% (Lambayeque1), and 85.8%, (PL9) (Supplemental Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). Although the exact order of the four lines differed between the microarray and RNA-seq analyses, Nipponbare and Sasanishiki were clearly more sensitive to cold treatment than Lambayeque1 and PL9 in both datasets.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eAnther microarray analyses with 13 lines and RNA-seq analyses with 4 lines\u003c/h2\u003e \u003cp\u003eWe collected anthers at the booting stage from control and cold-treated plants for each of the 13 lines, as two biological replicates, to examine their transcriptomes by microarray analysis. We used the 60K microarray designed with 37,955 probes from protein-coding genes and 23,013 probes from RSs (Supplemental Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). This microarray allowed us to comprehensively survey nearly all rice genes and RSs across the 13 lines. We converted the raw fluorescence intensity data of each probe on the microarray by logarithmic transformation (log\u003csub\u003e2\u003c/sub\u003e) and normalization by the quantile method (Bolstad et al. \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2003\u003c/span\u003e). We plotted the normalized expression values for each probe in control and cold-treated plants to assess the extent of correlation. In the resulting scatterplots (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), each rice line showed a high coefficient of determination (\u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e), ranging from 0.82 (A58) to 0.95 (Sasanishiki). Notably, probes falling outside of the diagonal (defined as \u003cem\u003ey\u003c/em\u003e\u0026thinsp;=\u0026thinsp;\u003cem\u003ex\u003c/em\u003e\u0026thinsp;\u0026plusmn;\u0026thinsp;1) varied among the lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eTo validate the microarray analysis, we performed an RNA-seq analysis of anthers from the Nipponbare, Sasanishiki, Lambayeque1, and PL9 lines. These plants were grown under the same conditions as for the microarray analysis but were independently collected. We obtained expression values for 37,860 genes using the Nipponbare reference genome (IRGSP-1.0 2020-12-02). When we plotted the RNA-seq data between the control and cold-treated anther samples in each line as a scatterplot, we obtained \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e values of 0.89, 0.88, 0.88, and 0.89 for Nipponbare, Sasanishiki, Lambayeque1, and PL9, respectively (Supplemental Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). We also plotted the expression values obtained through microarray and RNA-seq analyses for each of these four lines for control or cold-treated anthers. We calculated Pearson\u0026rsquo;s correlation coefficients (\u003cem\u003er\u003c/em\u003e) for each comparison, resulting in \u003cem\u003er\u003c/em\u003e values of 0.71 or above, indicating positive correlations. Specifically, the \u003cem\u003er\u003c/em\u003e values for the control and cold-treated samples were 0.82 and 0.80 for Nipponbare, 0.77 and 0.76 for Sasanishiki, 0.73 and 0.71 for Lambayeque1, and 0.75 and 0.77 for PL9, respectively (Supplemental Fig. S3). The \u003cem\u003er\u003c/em\u003e values of each line were comparable between the control and cold-treated samples and positive, suggesting sufficient equivalence between the microarray and RNA-seq datasets.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eCorrelation between expression levels and pollen fertility\u003c/h2\u003e \u003cp\u003eAs the 13 lines showed varying levels of cold tolerance, based on their PFI values, we wished to examine the possible correlation between gene expression levels in anthers and cold tolerance. Before conducting this analysis, we selected reliable and reproducible probes from the microarrays. In general, microarray probes with weaker fluorescence are less accurate. Indeed, probes with log\u003csub\u003e2\u003c/sub\u003e expression\u0026thinsp;\u0026lt;\u0026thinsp;5 were associated with larger coefficients of variation between replicates (Supplemental Fig. S4). We therefore applied a cutoff of log\u003csub\u003e2\u003c/sub\u003e expression\u0026thinsp;=\u0026thinsp;5, only considering probes for genes and RSs with higher values for analysis. From 60,168 original probes, we retained 38,060 and 37,807 probes with log\u003csub\u003e2\u003c/sub\u003e expression\u0026thinsp;\u0026ge;\u0026thinsp;5 in common among the 13 lines for control or cold-treated anthers, respectively.\u003c/p\u003e \u003cp\u003eTo examine the relationship between the PFI in the 13 lines and the expression of each gene and RS detected by the microarrays, we calculated the \u003cem\u003er\u003c/em\u003e values obtained from a scatterplot of expression (on the \u003cem\u003ey\u003c/em\u003e-axis) as a function of the PFI (on the \u003cem\u003ex\u003c/em\u003e-axis). In the control samples, we detected a significant and positive correlation (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05, \u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026ge;\u0026thinsp;0.55) for 668 probes (367 genes and 301 RSs) and a significant and negative correlation (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and \u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.55) for 724 probes (633 genes and 91 RSs) among the 38,060 probes (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA, B). A similar analysis conducted on the 37,807 reliably expressed probes (with log\u003csub\u003e2\u003c/sub\u003e expression\u0026thinsp;\u0026ge;\u0026thinsp;5) in the cold-treated samples yielded 234 probes (169 genes and 65 RSs) with a significant and positive correlation and 530 probes (479 genes and 51 RSs) with a significant and negative correlation (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC, D). In both cases, we identified more probes with a negative correlation to the PFI than with a positive correlation. This difference was significant in the cold-treated samples, based on Fisher's exact test of departure from a 1:1 ratio between probe types (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD). The number of probes from genes and RSs whose expression levels were negatively correlated with the PFI was already slightly higher than that for positively correlated probes under normal growth conditions; this difference became more pronounced after cold treatment (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eB, D).\u003c/p\u003e \u003cp\u003eWe extracted probes (for microarrays) or genes (for RNA-seq) showing a significantly different expression level between the cold-tolerant lines (Lambayeque1 and PL9) and the cold-sensitive lines (Nipponbare and Sasanishiki) with a false discovery rate (FDR)\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and an absolute log\u003csub\u003e2\u003c/sub\u003e(FC)\u0026thinsp;\u0026ge;\u0026thinsp;1. From the microarrays, we identified 23 probes with higher expression levels in the tolerant lines than in the sensitive lines when grown under control conditions and another 229 probes with higher expression levels in the sensitive lines than in the tolerant lines (Supplemental Fig. S5A). In the cold-treated samples, 19 probes had higher signals in the tolerant lines than in the sensitive lines, while another 167 probes yielding higher signal intensity in the sensitive lines than in the tolerant lines (Supplemental Fig. S5A). We observed similar patterns in the RNA-seq data for both growth conditions, as more genes showed a higher expression in the sensitive lines than in the cold-tolerant lines (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001). Indeed, there were 918 genes with higher expression in the sensitive lines under control conditions, compared to 530 more highly expressed genes in the tolerant lines under the same conditions; likewise, in the cold-treated samples, 2,453 genes were more highly expressed in the cold-sensitive lines, compared to only 1,223 more highly expressed genes in the cold-tolerant lines (Supplemental Fig. S5B). Although the exact numbers of probes or genes differed between the microarray and RNA-seq datasets, both analyses support the notion that the sensitive lines have significantly more up-regulated genes than the cold-tolerant lines under both growth conditions.\u003c/p\u003e \u003cp\u003eIn summary, the expression levels of genes and RSs in anthers tended to be less responsive to cold treatment in the cold-tolerant lines relative to the cold-sensitive lines. Notably, this pattern was already present in the anthers collected from plants grown under control conditions. Thus, the cold tolerance of these lines may be latent before cold treatment.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eIdentification of differentially expressed genes and repetitive sequences in response to cold stress\u003c/h2\u003e \u003cp\u003eWe identified DEGs and RSs (DERs) between the control and cold treatment samples using the microarray datasets from the 13 lines, yielding 27,908 DEGs and DERs with at least a 2-fold difference in expression (in either direction) between the control and cold-treated samples in at least one line. Among these, there were 1,897\u0026ndash;4,906 up-regulated genes and RSs (Ure) with higher expression (log\u003csub\u003e2\u003c/sub\u003e[FC)\u0026thinsp;\u0026ge;\u0026thinsp;1) in the cold-treated samples relative to control samples, and 1,942\u0026ndash;5,625 down-regulated genes and RSs (Dre) with lower expression (log\u003csub\u003e2\u003c/sub\u003e[FC]\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;1) upon cold treatment (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA). We obtained more DEGs than DERs within the Ure and Dre lists (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA).\u003c/p\u003e \u003cp\u003eWe then compared the number of DEGs and DERs among the 13 lines. We observed fewer common DERs compared to common DEGs as more lines were considered; for instance, the number of DEGs common to seven lines was 975 for Ure and 1,039 for Dre, compared to 10 Ure RSs and 104 Dre RSs (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB). While the numbers of Dre and Ure genes were comparable, the number of Dre RSs exceeded that of Ure RSs, suggesting that down-regulation of RSs is a distinct cold response. As more lines were considered, the numbers of common Ure and Dre genes decreased; nonetheless, there were about 100 common genes each for Ure and Dre with the same expression pattern in 12 lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB). By contrast, fewer RSs showed a similar expression pattern across many lines, with fewer than 100 Ure and Dre RSs when considering five and eight lines, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB).\u003c/p\u003e \u003cp\u003eWe identified 20 Ure genes and 24 Dre genes showing the same expression profile in all 13 lines (Tables\u0026nbsp;2, 3, Supplemental Fig. S6), which may be characteristic of the intrinsic transcriptome responses to low temperature at the booting stage. With the exception of the Ure gene Os06g0246500 (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.66), we observed no significant correlation between the log\u003csub\u003e2\u003c/sub\u003e(FC) of the 43 remaining DEGs and the PFI values (Supplemental Fig. S7). Of these 44 common DEGs, we identified six as Ure genes and nine as Dre genes in the RNA-seq data (Tables\u0026nbsp;2, 3). Among these six Ure genes, we noticed Os04g0568700, which encodes a heat shock factor (HSF) related to multiple stress responses including cold (Chauhan et al. \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Guo et al. \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2022b\u003c/span\u003e) (Supplemental Fig. S6A, Table\u0026nbsp;2). Of the nine Dre genes, the most significantly down-regulated gene was Os11g0150400, encoding a stress-responsive alpha/beta barrel domain\u0026ndash;containing protein (Supplemental Fig. S6B, Table\u0026nbsp;3).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eCorrelation between expression change and pollen fertility by cold treatment\u003c/h2\u003e \u003cp\u003eBased on the microarray data with the 13 lines, we analyzed the relationship between expression change and pollen fertility (as measured by the PFI) in response to cold treatment. To this end, we calculated the \u003cem\u003er\u003c/em\u003e values for each probe from a scatterplot between the PFI and the log\u003csub\u003e2\u003c/sub\u003e(FC) values for all 13 lines. Of the 27,908 DEGs and DERs, 905 showed a significant negative correlation (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.55, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) with the PFI, and another 514 had a significant positive correlation (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026ge;\u0026thinsp;0.55, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) with the PFI (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA). We further divided the DEGs and DERs based on the median log\u003csub\u003e2\u003c/sub\u003e(FC) values (designated expression change index [ECI] hereafter) across the 13 lines. A positive ECI indicates an Ure probe in response to cold treatment, while a negative ECI indicates a Dre probe under the same condition. With these parameters, we divided the 1,419 probes above into four types by combining \u003cem\u003er\u003c/em\u003e values and ECI values (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). Type 1 indicated that 600 Ure probes (137 genes and 463 RSs) tend to have weaker up-regulation in cold-tolerant lines than in cold-sensitive lines; 233 Ure probes (152 genes and 81 RSs) in type 2 have stronger up-regulation in the tolerant lines than in the sensitive lines; 305 Dre probes (172 genes and 133 RSs) in type 3 have stronger down-regulation in the tolerant lines than in the sensitive lines; 281 Dre probes (155 genes and 126 RSs) in type 4 have weaker down-regulation in the tolerant lines than in the sensitive lines. The typical expression patterns with highest \u003cem\u003er\u003c/em\u003e values in types 1\u0026ndash;4 are illustrated by Os01g0699400, Os04g0604000, Os01g0191200, and Os01g0730600, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC). The number of genes was comparable across types, ranging from 137 to 172. Notably, while the number of RSs was similar among types 2\u0026ndash;4, type 1 was characterized by far more (463) RSs (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). This result suggested that most RSs related to pollen fertility are up-regulated in the sensitive lines relative to the tolerant lines. Consistent with this finding, Ishiguro et al. (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) previously reported that the expression levels of RSs in cold-sensitive rice lines changed more in anthers in response to low temperatures than in cold-tolerant lines. The expression of RSs in tolerant lines may therefore be less responsive to cold treatment (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB).\u003c/p\u003e \u003cp\u003eSince we only used four lines for RNA-seq analysis, we opted for different selection criteria to divide DEGs into four types equivalent to those obtained by the analysis of the microarray data. Specifically, we defined type 1 genes as those with a positive ECI value in the four lines and with a lower average log\u003csub\u003e2\u003c/sub\u003e(FC) for the two tolerant lines (Lambayeque1 and PL9) than for the two sensitive lines (Nipponbare and Sasanishiki); type 2 genes had a positive ECI value and with a higher average log\u003csub\u003e2\u003c/sub\u003e(FC) for the two tolerant lines than for the two sensitive lines; type 3 genes had negative ECI values and a lower mean log\u003csub\u003e2\u003c/sub\u003e(FC) for the two tolerant lines than for the two sensitive lines; and type 4 genes had a negative ECI value and a greater mean log\u003csub\u003e2\u003c/sub\u003e(FC) value for the two tolerant lines than for the two sensitive lines. The above analysis returned 3,167 genes, of which 295 belonged to the same type in the microarray and RNA-seq datasets, with 52 genes for type 1, 48 genes for type 2, 154 genes for type 3, and 41 genes for type 4 (Supplemental Fig. S8). Overlapping genes were relatively more abundant for type 3 DEGs than for the other types of DEGs and corresponded to genes whose expression levels decreased upon cold treatment in the cold-tolerant lines.\u003c/p\u003e \u003cp\u003eOverall, the correlations between expression change and PFI clearly showed the tendency of numerous RSs with activated expressions in the sensitive lines and more genes with reduced expression in the tolerant lines.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eGO term enrichment analysis of anther transcriptomes\u003c/h2\u003e \u003cp\u003eFrom the transcriptome analysis of the 13 rice lines above, we obtained three sets of information: (1) the correlation between expression levels and the PFI, (2) a list of common DEGs across the lines, and (3) the extent of correlation between the ECI and PFI values. In each case, we generated lists of genes with characteristic expression patterns, which might be related directly or indirectly to the pollen sterility caused by cold exposure at the booting stage. To explore the function of these genes, we performed a Gene Ontology (GO) term enrichment analysis on each gene list (Supplemental Tables S2\u0026ndash;S8), focusing on terms associated with fewer than 100 genes across the entire rice genome and exhibiting significant (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.001) enrichment (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eWe first characterized the genes showing a negative or positive correlation in their expression levels under control growth conditions with the PFI values. From genes whose expression levels are negatively correlated with the PFI, we identified 15 significantly enriched GO terms, most of which were related to \u0026ldquo;cytosolic small ribosomal subunit\u0026rdquo; (GO:0022627) and \u0026ldquo;small ribosomal subunit\u0026rdquo; (GO:0015935) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA right panel). Ribosome biosynthesis is known to be related to pollen development and the stress response in Arabidopsis (\u003cem\u003eArabidopsis thaliana\u003c/em\u003e) (Ren\u0026aacute;k et al. \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Structural changes in ribosomes contribute to stress acclimation by forming stress-responsive ribosomes under stress conditions (Dias-Fields and Adamala \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). This result also suggests that cold-sensitive lines may exhibit higher expression levels for ribosome-related genes than cold-tolerant lines. The cold tolerance of plants might thus be predicted from the expression levels of specific ribosome-related genes under normal growth conditions.\u003c/p\u003e \u003cp\u003eFor genes whose expression levels were positively correlated with the PFI in control samples, we obtained six significantly enriched GO terms (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA left panel). These GO terms were jasmonic acid (JA)- and fatty acid\u0026ndash;related GO terms: \u0026ldquo;jasmonic acid mediated signaling pathway\u0026rdquo; (GO:0009867), \u0026ldquo;cellular response to jasmonic acid stimulus\u0026rdquo; (GO:0071395), \u0026ldquo;response to jasmonic acid\u0026rdquo; (GO:0009753), \u0026ldquo;cellular response to fatty acid\u0026rdquo; (GO:0071395), and \u0026ldquo;response to fatty acid\u0026rdquo; (GO:0070542) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA left panel). JA is a phytohormone related to multiple stress responses including cold stress (Raza et al. \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). The most common unsaturated fatty acids in plants, 18-carbon species, can be modified into bioactive molecules, including JA (He and Ding \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Therefore, we speculate that the expression of genes involved in JA signal transduction is activated by cold stress. The expression pattern of these JA-related genes remained higher in the cold-tolerant lines than in the sensitive lines under normal growth conditions.\u003c/p\u003e \u003cp\u003eWe also looked for GO terms enriched in Ure or Dre genes common to 11 lines or more. With fewer common DEGs across more lines, we failed to identify significant GO terms for DEGs shared by 12 or more lines (Supplemental Fig. S9). We detected five GO terms for DEGs shared by 11 lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB left). The GO term enriched among Ure genes was \u0026ldquo;positive regulation of response to salt stress\u0026rdquo; (GO:1901002), which would be consistent with a shared signaling pathway between salinity and cold stress responses (Zhang et al. \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2022b\u003c/span\u003e; Zhang and Xia \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). The four GO terms enriched among Dre genes were related to DNA metabolism such as \u0026ldquo;DNA packaging\u0026rdquo; (GO:0044815), \u0026ldquo;nucleosome\u0026rdquo; (GO:0000786), \u0026ldquo;protein-DNA complex\u0026rdquo; (GO:0032993), and \u0026ldquo;structural constituent of chromatin\u0026rdquo; (GO:0030527) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB right). These GO terms suggest that in most of the lines, DNA metabolism is reduced due to cold stress during the booting stage. We also obtained the list of GO terms associated with each of the four types of genes defined by correlations between ECI and PFI values (Supplemental Table S8).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eGenes showing expression patterns related to low temperature at the booting stage\u003c/h2\u003e \u003cp\u003eHere, we highlighted known and unknown genes related to cold tolerance based on the GO term enrichment analysis of the microarray and RNA-seq data.\u003c/p\u003e \u003cp\u003eFrom the microarray data, we identified 43 and 49 genes whose expression levels were positively or negatively correlated with the PFI (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026ge;\u0026thinsp;0.7 or \u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.7), respectively, under control conditions (Supplemental Tables S9, S10). Among the 43 genes with a positive correlation, Os04g0623300 (\u003cem\u003ePOLYAMINE OXIDASE3\u003c/em\u003e, \u003cem\u003eOsPAO3\u003c/em\u003e) showed the highest correlation with the PFI (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.85) and was previously reported to be up-regulated in response to cold stress in rice seedlings (Sagor et al. \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). For the remaining 91 genes, an association with cold stress has not been reported until now, although we noticed three ribosome-related genes among them: Os06g0550000, Os01g0962600, and Os08g0234000. The expression levels of these three ribosome-related genes were negatively correlated (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.7) with the PFI under control conditions (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA), in agreement with the GO term enrichment analysis above.\u003c/p\u003e \u003cp\u003eFor cold-treated plants, the expression levels of 16 and 28 genes were significantly and positively or negatively correlated (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026ge;\u0026thinsp;0.7 or \u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.7), respectively, with the PFI (Supplemental Tables S9, S10). However, none of these genes have been reported to have a role in cold stress responses, except for \u003cem\u003eMITOGEN-ACTIVATED PROTEIN KINASE KINASE KINASE70\u003c/em\u003e (\u003cem\u003eOsMKKK70\u003c/em\u003e) (Mei et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). \u003cem\u003eOsMKKK70\u003c/em\u003e, whose expression was negatively correlated (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.63) with the PFI (Supplemental Fig. S10), was reported to be associated with cold tolerance at the booting stage (Mei et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Other genes whose expression levels were positively correlated with the PFI values included Os01g0962700 (\u003cem\u003ePEROXIDASE20\u003c/em\u003e, \u003cem\u003ePRX20\u003c/em\u003e), Os06g0156400 (\u003cem\u003eGDSL ESTERASE/LIPASE76\u003c/em\u003e, \u003cem\u003eOsGELP76\u003c/em\u003e), and Os03g0682200 (\u003cem\u003eARGONAUTE12\u003c/em\u003e, \u003cem\u003eOsAGO12\u003c/em\u003e) (Supplemental Table S9). \u003cem\u003ePRX20\u003c/em\u003e, which encodes a peroxidase, was highly expressed in cold-tolerant lines (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.73 for the expression\u0026ndash;PFI relationship) and was previously reported to participate in the response to salinity stress (Kim and Kim \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). As with samples collected from control plants, of the 28 genes whose expression levels were negatively correlated (\u003cem\u003er\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;\u0026minus;\u0026thinsp;0.7) with the PFI, we noticed a ribosome-related gene, Os03g0284400 (\u003cem\u003er\u003c/em\u003e\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.73), suggesting that ribosome status in rice plants under normal or cold stress conditions potentially influences pollen fertility upon cold exposure.\u003c/p\u003e \u003cp\u003eWe took a closer look at the 20 and 24 common Ure and Dre genes shared by all 13 lines (Tables\u0026nbsp;2, 3, Supplemental Fig. S6). The expression of these genes responds to cold stress in a similar manner regardless of their innate cold tolerance or sensitivity. The 20 common Ure genes included NAC-type transcription factor genes such as \u003cem\u003eOsNAC5\u003c/em\u003e and \u003cem\u003eONAC088\u003c/em\u003e (Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e, Supplemental Fig. S6). We did not detect known stress response genes among the 24 common Dre genes but were intrigued by the presence of several histone genes, Os08g0427700 (\u003cem\u003eH2A\u003c/em\u003e), Os04g0583600 (\u003cem\u003eH4\u003c/em\u003e), and Os08g0490900 (\u003cem\u003eH2B\u003c/em\u003e), whose functional association with cold stress might be cell division (Table\u0026nbsp;3, Supplemental Fig. S6B). These 44 common DEGs did not include any known genes associated with cold tolerance at the booting stage. However, three genes known to be associated with cold tolerance were among common Ure or Dre genes when considering nine or more of the 13 lines: \u003cem\u003eLATE EMBRYOGENESIS ABUNDANT9\u003c/em\u003e (\u003cem\u003eOsLEA9\u003c/em\u003e), \u003cem\u003eOsMKKK70\u003c/em\u003e, and \u003cem\u003eOsMYB4\u003c/em\u003e (Table\u0026nbsp;4) (Lou et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Mei et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Park et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2010\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eReturning to genes from types 1\u0026ndash;4 defined based on their correlation between the log\u003csub\u003e2\u003c/sub\u003e(FC) values and PFI values (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e), we focused on genes with pronounced changes in expression. Using a selection criterion of ECI\u0026thinsp;\u0026gt;\u0026thinsp;1 or ECI\u0026thinsp;\u0026lt;\u0026thinsp;\u0026minus;\u0026thinsp;1, we obtained 86 genes (Supplemental Table S13), consisting of 20 type 1 genes, 13 type 2 genes, 25 type 3 genes, and 28 type 4 genes. Five of the 86 genes were also present among the DEGs classified into four types using the ECI and the log\u003csub\u003e2\u003c/sub\u003e(FC) values from the RNA-seq data. Two of these genes are known to participate in the cold tolerance of rice. Os02g0579000, a NAC-type transcription factor gene induced at low temperatures (Fang et al. \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2014\u003c/span\u003e), was among the type 1 genes, which tend to be more highly expressed in the cold-sensitive lines. This gene showed an increased expression in 11 lines in response to cold stress, reaching higher levels in the cold-sensitive lines than in the tolerant lines (Supplemental Table S13). The other NAC-type transcription factor gene, Os03g0133000, known to respond to some abiotic stresses, including cold (Hong et al. \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2016\u003c/span\u003e), was a type 2 gene, more highly expressed in cold-tolerant lines and exhibiting increased expression in 10 lines in response to cold stress. The identification of several known genes and their expression changes in response to low temperatures support the validity of the microarray-based analysis in this study.\u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eExpression levels and cold tolerance: negative correlation between pollen fertility and expression of ribosome-related genes\u003c/h2\u003e \u003cp\u003eWe defined genes whose expression levels were negatively or positively correlated with the PFI, a measure of cold tolerance at the booting stage. There were significantly more genes with low expression in the cold-tolerant lines and with high expression in the cold-sensitive lines in response to cold stress than genes with high expression in the cold-tolerant lines and with low expression in the cold-sensitive lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC, D). We corroborated this result with our RNA-seq data obtained from two cold-tolerant liens and two cold-sensitive lines (Supplemental Fig. S5B). Among genes whose expression levels were negatively correlated with the PFI under cold conditions, we observed a significant enrichment for ribosome-related GO terms (Supplemental Table S5). When plants experience abiotic stresses, such as low temperatures, they produce inoperative ribosomes (Dias-Fields and Adamala \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). This phenomenon is thought to function to conserve energy in response to stress by slowing down protein translation and growth (Bechtold and Field \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Zandalinas et al. \u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Zhang et al. \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2022a\u003c/span\u003e). The lower induction of expression of these ribosome-related genes in cold-tolerant lines than in sensitive lines may reflect energy conservation, which may affect pollen development (Table\u0026nbsp;5). The expression of ribosome-related genes in anthers of the cold-tolerant lines tended to be lower than that of the sensitive lines under both cold treatment and control conditions (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). Therefore, the expression of ribosome-related genes in the anthers of cold-tolerant lines may be lower than in cold-sensitive lines at all times.\u003c/p\u003e \u003cp\u003eA GO analysis of genes whose expression levels under control conditions were significantly and positively correlated with the PFI identified JA-related genes, which are known stress-responsive genes, indicating that these genes maintain higher expression levels in cold- tolerant lines than in cold- sensitive lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA, Supplemental Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e) (Devireddy et al. \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Kim et al. \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). This finding suggests that the high expression level of JA-related genes may reflect the degree of cold tolerance in each line.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eLow temperature\u0026ndash;induced changes in expression: common DEGs among lines\u003c/h2\u003e \u003cp\u003eIn each line, we identified 4,117 to 12,191 genes on the microarrays whose expression levels changed in either direction by more than 2-fold upon cold treatment, defining a set of DEGs and DERs. The number of up-regulated and down-regulated (Ure and Dre) probes in each line was 1,897\u0026ndash;4,906 and 1,942\u0026ndash;5,625, respectively, with no remarkable difference between the tolerant and sensitive lines. The low-temperature tolerance of rice emerged from the analysis of DEGs: We detected four genes annotated \u0026ldquo;regulation of response to salt stress\u0026rdquo; (GO:1901000) among Ure genes in more than 11 lines (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB, Supplemental Table S6), possibly reflecting a common pathway involved in responses to salinity and cold stresses (Devireddy et al. \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Hu et al. \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Kim et al. \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Nakashima et al. \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Ren et al. \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). We also identified eight genes related to chromosome segregation, such as \u0026ldquo;DNA packaging complex\u0026rdquo; (GO:0044815), among Dre genes; they may be related to diminished cell division due to low temperatures (Supplemental Table S7) (Qi and Zhang \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). On a gene-by-gene basis, we identified 20 and 24 Ure and Dre genes, respectively, common to all 13 lines (Tables\u0026nbsp;2, 3, Supplemental Fig. S6). The 20 Ure genes included \u003cem\u003eOsNAC5\u003c/em\u003e and \u003cem\u003eONAC088\u003c/em\u003e, NAC-type transcription factor genes involved in stress responses including low temperature (Table\u0026nbsp;2, Supplemental Fig. S6A) (Nakashima et al. \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Sun et al. \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Takasaki et al. \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). The inclusion of known low temperature\u0026ndash;responsive genes among the Ure genes validates our detection of cold stress responses in anthers. The 24 Dre genes did not contain any known low temperature\u0026ndash;responsive genes but did contain genes encoding multiple histones: Os08g0427700 (\u003cem\u003eH2A\u003c/em\u003e), Os04g0583600 (\u003cem\u003eH4\u003c/em\u003e), and Os08g0490900 (\u003cem\u003eH2B\u003c/em\u003e) (Table\u0026nbsp;3, Supplemental Fig. S6B). In addition to known DEGs in response to low temperature, we also discovered genes whose expression levels respond to low temperature in a similar manner among lines with different genetic backgrounds. One of the 20 Ure genes, Os06g0246500, whose expression positively correlated with the PFI, encodes a pyruvate dehydrogenase E1 alpha subunit-like protein, which has not been reported to be associated with cold tolerance (Supplemental Fig. S7).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eComplexity of cold tolerance: expression patterns of known cold tolerance genes across multiple genetic backgrounds\u003c/h2\u003e \u003cp\u003eWe looked at the expression levels of 12 genes reported to be associated with cold tolerance at the booting stage in rice (Table\u0026nbsp;4). \u003cem\u003eOsLEA9\u003c/em\u003e, \u003cem\u003eOsMKKK70\u003c/em\u003e, and \u003cem\u003eOsMYB4\u003c/em\u003e all showed increased expression in at least nine lines, confirming their association with cold responsiveness reported in previous studies (Lou et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Mei et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Park et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Of the 12 genes, \u003cem\u003eOsMYB4\u003c/em\u003e and \u003cem\u003eOsAPX1\u003c/em\u003e promote cold tolerance when overexpressed (Park et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Sato et al. \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). However, the expression levels or changes in expression between control and cold stress conditions for \u003cem\u003eOsMYB4\u003c/em\u003e and \u003cem\u003eOsAPX1\u003c/em\u003e showed no significant correlation with the PFI in this study. This result indicates the complexity of cold stress responses during the booting stage. Most previous studies on stress response have focused on a few loci, revealed from comparisons among a small number of lines with contrasting cold tolerance or on a few genes conferring a strong stress response (Mehrotra et al. \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Raj and Nadarajah \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Thus, the extent to which genes involved in stress tolerance can explain differences in diverse degrees of cold tolerance among a large number of lines with distinct genetic backgrounds has not been tested. Here, the expression levels for most known cold tolerance genes did not correlate with the PFI. Thus, individual genes cannot universally explain differences in cold tolerance during the booting stage of rice because the factors contributing to cold tolerance vary from line to line (Table\u0026nbsp;5).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eCold tolerance and changes in expression of repetitive sequences\u003c/h2\u003e \u003cp\u003eIn this study, we performed a correlation analysis between the extent of change in expression for each probe and the PFI; we then classified genes into one of four types based on their ECI and \u003cem\u003er\u003c/em\u003e- values (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). While we generally identified more DEGs than DERs, this was not the case for type I genes, with 137 genes and 463 RSs (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). This result indicates a genome-wide trend of increased expression of more RSs in the cold-sensitive lines (Table\u0026nbsp;5). TEs, which constitute the majority of RSs, induce mutations in the host genome through transposition when they are activated, and their expression is normally suppressed by epigenetic mechanisms such as DNA methylation and histone modifications (Choi and Lee \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Parker et al. \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). However, plants exposed to abiotic stresses such as high or low temperatures experience a weakened repression of RSs, leading to their increased expression (Hu et al. \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Ito et al. \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Zhou et al. \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Therefore, the above trend suggests that the expression of RSs is more suppressed upon cold stress conditions in cold-tolerant lines than in cold-sensitive lines. A previous transcriptome analysis of cold-treated anthers at the booting stage in rice showed that fewer RSs fluctuated in their expression under cold stress in cold-tolerant lines than in cold-sensitive lines (Ishiguro et al. \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). The present study and the previous transcriptome study (Ishiguro et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) thus describe a genome-wide trend whereby cold-tolerant lines are less responsive to cold stress (Table\u0026nbsp;5).\u003c/p\u003e \u003cp\u003eThe expression levels of known cold tolerance genes did not correlate with the cold tolerance of the 13 lines tested in this study. Conventional stress-related studies, which focus on a small number of genes from pairwise comparisons of specific lines, therefore cannot fully explain the diversity in cold tolerance responses observed among multiple lines. By contrast, this study focused on genome-wide responses to low-temperature stress and revealed that the transcriptome of more cold-tolerant lines is relatively insensitive to low-temperature stress. Focusing on this feature should be a useful new approach to comprehensively elucidate the differences in stress tolerance intensity among lines. In this study, we focused on genome-wide responsiveness to stress and successfully explained the difference between stress-tolerant and -sensitive plants.\u003c/p\u003e \u003c/div\u003e "},{"header":"Materials \u0026 Methods","content":"\u003cdiv id=\"Sec15\" type=\"MaterialsAndMethods\" class=\"Section2\"\u003e \u003cdiv id=\"Sec16\" class=\"Section3\"\u003e \u003ch2\u003ePlant materials\u003c/h2\u003e \u003cp\u003eThirteen rice lines were used in this study: eight temperate \u003cem\u003ejaponica\u003c/em\u003e lines, \u0026lsquo;Hokkai PL9\u0026rsquo; (PL9), \u0026lsquo;Kokushokuine 2 go\u0026rsquo; (A58), \u0026lsquo;Koshihikari\u0026rsquo;, \u0026lsquo;Hokkai 287\u0026rsquo;, Nipponbare, \u0026lsquo;Kirara 397\u0026rsquo;, Sasanishiki, and \u0026lsquo;Fukoku\u0026rsquo;; one tropical \u003cem\u003ejaponica\u003c/em\u003e line (Lambayeque1); one \u003cem\u003eindica\u003c/em\u003e line (\u0026lsquo;Kasalath\u0026rsquo;); and three recombinant inbred lines (RILs) \u0026lsquo;R45\u0026rsquo;, R75, and R85, which were derived from a cross between the \u003cem\u003ejaponica\u003c/em\u003e line A58 and the \u003cem\u003eindica\u003c/em\u003e line \u0026lsquo;Surjamkhi\u0026rsquo; (I33) (Table\u0026nbsp;1). Materials were grown with 20 individuals sown in one Wagner pot (1/5,000 100 m\u003csup\u003e2\u003c/sup\u003e) in a greenhouse under 25\u0026deg;C (day)/19\u0026deg;C (night) temperature cycles with only the main stem maintained by cutting off the offshoots. For the cold treatment, plants were exposed to 12\u0026deg;C for 4 days during the booting stage, when the microspores are between the late tetrad stage and the uninucleated pollen stage after meiosis. The booting stage was determined by auricle distance (AD) as in our previous study, with AD\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;2 to \u0026minus;\u0026thinsp;4 cm for Hokkaido and AD\u0026thinsp;=\u0026thinsp;0\u0026ndash;2 cm for other lines (Yamamori et al. \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Anthers were collected from control plants and plants exposed to cold treatment at the same growth stage and immediately stored at \u0026minus;\u0026thinsp;80\u0026deg;C. Anthers were collected from the third, fourth, and fifth spikelets of primary rachis branches. Samples for microarrays were anthers collected from multiple individuals for each line and treatment, with two samples per line. Samples for RNA-seq were collected in triplicates using anthers from several individuals for each replicate. To investigate pollen fertility, some individuals were returned to the greenhouse after the cold treatment and allowed to continue growing.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eObservation of pollen fertility\u003c/h2\u003e \u003cp\u003eAnthers were collected from the third, fourth, and fifth anthers from the tip of the primary branches during the flowering period. Pollen grains were stained with Lugol's iodine solution [KI\u0026ndash;I\u003csub\u003e2\u003c/sub\u003e: 1.5% (w/v) KI, 0.15% (w/v) I\u003csub\u003e2\u003c/sub\u003e]. Pollen fertility was calculated as the ratio between the number of fertile pollen grains and the total number of pollen grains for each line and under each growth condition. Pollen fertility was determined by sampling at least nine flowers per line, and the average of these was used as the pollen fertility for each line and each growth condition.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eMicroarray probe design\u003c/h2\u003e \u003cp\u003eIn this experiment, a 60K microarray was designed with 23,103 probes derived from repetitive sequences (RSs) and 37,955 derived from protein-coding genes and pre-miRNAs (Supplemental Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). The probes derived from RSs were selected from the 44K microarray designed by Ishiguro et al. (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e), taking probes with a fluorescence signal value of 100 or higher. All probes derived from protein-coding genes and pre-miRNAs were from the Agilent catalog array (G2519F) (Agilent technologies, California, USA).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eMicroarray analysis\u003c/h2\u003e \u003cp\u003eTotal RNA used for microarray analysis was extracted from the above samples and purified using a TRIzol Plus RNA purification kit (Life Technologies, California, USA) after freezing and pulverization in liquid nitrogen. RNA concentration was measured with a Nano Drop ND-2000 (Thermo Fisher Scientific, Massachusetts, USA) and quality-checked with an Agilent 2100 Bioanalyzer (Agilent Technologies, California, USA). The microarray analysis was performed with One-Color Spike-Mix containing cyanine-3 (Cy3)-CTP (Agilent) added to the RNA samples, which were labeled using a Quick Amp Labeling Kit (Agilent). After fragmentation of 600 ng complementary RNA (cRNA) synthesized from RNA samples, hybridization to the microarray was performed at 65\u0026deg;C for 17 h using the Agilent Gene Expression Hybridization Kit (Agilent). The slides were then washed, and fluorescence was measured using an Agilent Technologies C version scanner. Agilent Feature Extraction software (Agilent) was used to quantify the fluorescence signals. The data from all samples were log\u003csub\u003e2\u003c/sub\u003e-transformed using R (4.1.0) and normalized by the quantile normalization method using the limma (3.48.3) package in R (Ritchie et al. \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2015\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eRNA-seq\u003c/h2\u003e \u003cp\u003eTotal RNA was extracted from anthers as described for microarray analysis; three replicates were used per genotype. The extracted RNA was subjected to strand-specific library preparation (dUTP method) using the NEBNext Ultra\u0026trade; ll Directional RNA Library Prep Kit (New England Biolabs, Massachusetts, USA). After library preparation, libraries were sequenced as paired-end 150-bp reads on a NovaSeq 6000 instrument (Illumina, California, USA). Library preparation and sequencing were performed by Rhelixa (Tokyo, Japan). From the data obtained by sequencing, kallisto (0.46.2) was used to obtain expression levels of each gene (Bray et al. \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). IRGSP-1.0 (2021-05-10), published by RAP-DB, was used as a reference for RNA-seq analysis (Kawahara et al. \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2013\u003c/span\u003e; Sakai et al. \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). The bootstrap value was set to 100, and the expression levels were normalized to transcripts per million (TPM) using kallisto. Sleuth (0.30.0) was used to compare expression levels between cold-treated and control plants for each line and between lines (Pimentel et al. \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Differences in expression levels between treatments in each line were statistically tested using the likelihood ratio test.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003eGene annotation and GO analysis\u003c/h2\u003e \u003cp\u003eThe genes selected in this study were annotated based on IRGSP-1.0 (2021-05-10). GO analysis of the selected genes was performed using PANTHER (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://geneontology.org/\u003c/span\u003e\u003cspan address=\"http://geneontology.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) provided by The Gene Ontology Consortium (Mi et al. \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Fisher's exact test was used to identify GO terms that were significantly (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05) more abundant in the selected genes than across the genome.\u003c/p\u003e \u003c/div\u003e"},{"header":"Conclusion","content":"\u003cp\u003eWe investigated the relationship between genome-wide response to low-temperature stress and cold stress tolerance based on transcriptome analysis of 13 rice lines with varying degrees of cold tolerance. We interpreted the results of these analyses from three perspectives: gene expression levels and cold tolerance, changes in expression in response to cold, and changes in expression as a function of cold tolerance. The relationship between expression levels and cold tolerance revealed that, compared to cold-sensitive lines, the cold-tolerant lines tested here tended to have more genes with lower expression than with higher expression, regardless of growth conditions. Some of the genes with lower expression in cold-tolerant lines were related to ribosomes, suggesting that transcription and translation of the genes are stabilized in the cold-tolerant lines. Looking at genes with altered expression in response to cold, we identified genes with common differential expression across the 13 lines tested. We thus identified new genes that are universally responsive to cold in rice anthers. Many of these genes had not previously been reported to be associated with cold stress. We finally revealed that the RSs whose changes in expression levels showed significant correlations with cold tolerance most often exhibited increased expression in the cold-sensitive lines but remained largely constant in cold-tolerant lines. This result is consistent with our previous study (Ishiguro et al. \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) and reaffirms the genomic stability under low-temperature stress displayed by cold-tolerant lines. In summary, this study provides new insights into cold stress responses during the booting stage by focusing on the relationship between the various cold tolerance intensities exhibited by the lines and genome-wide expression changes.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cdiv class=\"DefinitionList\"\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eDEG\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003edifferentially expressed gene\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eDER\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003edifferentially expressed repetitive sequences\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eDre\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003edown-regulated\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eECI\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eexpression change index\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003ePFI\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003epollen fertility index\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eRS\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003erepetitive sequence\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eTPM\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003etranscripts per million\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003e\u003cb\u003eUre\u003c/b\u003e\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eup-regulated\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthical Approval and Consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis study complied with the ethical standards of Japan, Hokkaido University, where this research was conducted.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll authors have consented to the publication of this manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of supporting data\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eOur microarray data were recorded as GSExxx in the National Center for Biotechnology Information (NCBI) provided by http://www. ncbi.nlm.nih.gov/yyyy. Our RNA-seq data were recorded as GSExxx in the NCBI provided by http://www. ncbi.nlm.nih.gov/yyyy .\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by JSPS KAKENHI grants (19H00937 and 23H02180) to Y.Ki and Scientific Technique Research Promotion Program for Agriculture, Forestry, Fisheries and Food Industry to Y.S.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eK.Y., Y.S. and Y.K. designed the research plans; S.I. and K.O. performed the cold stress evaluations and the microarray experiments; K.Y. and K.M.L. performed the cold stress evaluations and the RNA-seq experiments; K.Y. and K.F. performed the data analysis; K.F. and Y.S. provided instruments for the experiments; K.O. produced RIL lines; K.Y. and Y.K. contributed to manuscript preparation; K.Y. and Y.K. wrote the article with contributions from all the authors.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe gratefully acknowledge Dr. Y. Koide (Research Faculty of Agriculture, Hokkaido University) for his valuable suggestions concerning this study.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAndaya VC, Mackill DJ (2003) QTLs conferring cold tolerance at the booting stage of rice using recombinant inbred lines from a japonica x indica cross. Theor Appl Genet 106:1084-1090\u003c/li\u003e\n\u003cli\u003eBechtold U, Field B (2018) Molecular mechanisms controlling plant growth during abiotic stress. J Exp Bot 69:2753-2758\u003c/li\u003e\n\u003cli\u003eBolstad BM, Irizarry RA, Astrand M, Speed TP (2003) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19:185-193\u003c/li\u003e\n\u003cli\u003eBray NL, Pimentel H, Melsted P, Pachter L (2016) Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 34:525-527\u003c/li\u003e\n\u003cli\u003eChauhan H, Khurana N, Agarwal P, Khurana P (2011) Heat shock factors in rice (Oryza sativa L.): genome-wide expression analysis during reproductive development and abiotic stress. Mol Genet Genomics 286:171-187\u003c/li\u003e\n\u003cli\u003eChoi JY, Lee YCG (2020) Double-edged sword: The evolutionary consequences of the epigenetic silencing of transposable elements. PLoS Genet 16:e1008872\u003c/li\u003e\n\u003cli\u003eDai LY, Lin XH, Ye CR, Ise KZ, Saito K, Kato A, Xu FR, Yu TQ, Zhang DP (2004) Identification of quantitative trait loci controlling cold tolerance at the reproductive stage in Yunnan landrace of rice, Kunmingxiaobaigu. Breeding Science 54:253-258\u003c/li\u003e\n\u003cli\u003eDevireddy AR, Zandalinas SI, Fichman Y, Mittler R (2021) Integration of reactive oxygen species and hormone signaling during abiotic stress. Plant J 105:459-476\u003c/li\u003e\n\u003cli\u003eDias-Fields L, Adamala KP (2022) Engineering ribosomes to alleviate abiotic stress in plants: A perspective. Plants (Basel) 11:2097\u003c/li\u003e\n\u003cli\u003eFang Y, Xie K, Xiong L (2014) Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot 65:2119-2135\u003c/li\u003e\n\u003cli\u003eGuo M, Liu JH, Ma X, Luo DX, Gong ZH, Lu MH (2016) The plant heat stress transcription factors (HSFs): Structure, regulation, and function in response to abiotic stresses. Front Plant Sci 7:114\u003c/li\u003e\n\u003cli\u003eHayase H, Satake T, Nishiyama I, Ito N (1969) Male sterility caused by cooling treatment at the meiotic stage in rice plants: 2. The most sensitive stage to cooling and the fertilizing ability of pistils. Japanese Journal of Crop Science 38:706-711\u003c/li\u003e\n\u003cli\u003eHe M, Ding NZ (2020) Plant unsaturated fatty acids: multiple roles in stress response. Front Plant Sci 11:562785\u003c/li\u003e\n\u003cli\u003eHong Y, Zhang H, Huang L, Li D, Song F (2016) Overexpression of a stress-responsive NAC transcription factor gene ONAC022 improves drought and salt tolerance in rice. Front Plant Sci 7:4\u003c/li\u003e\n\u003cli\u003eHu H, You J, Fang Y, Zhu X, Qi Z, Xiong L (2008) Characterization of transcription factor gene SNAC2 conferring cold and salt tolerance in rice. Plant Mol Biol 67:169-181\u003c/li\u003e\n\u003cli\u003eHu Y, Zhang L, He S, Huang M, Tan J, Zhao L, Yan S, Li H, Zhou K, Liang Y, Li L (2012) Cold stress selectively unsilences tandem repeats in heterochromatin associated with accumulation of H3K9ac. Plant Cell Environ 35:2130-2142\u003c/li\u003e\n\u003cli\u003eIshiguro S, Ogasawara K, Fujino K, Sato Y, Kishima Y (2014) Low temperature-responsive changes in the anther transcriptome\u0026apos;s repeat sequences are indicative of stress sensitivity and pollen sterility in rice strains. Plant Physiol 164:671-682\u003c/li\u003e\n\u003cli\u003eIto H, Gaubert H, Bucher E, Mirouze M, Vaillant I, Paszkowski J (2011) An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress. Nature 472:115-119\u003c/li\u003e\n\u003cli\u003eKawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S, Childs KL, Davidson RM, Lin H, Quesada-Ocampo L, Vaillancourt B, Sakai H, Lee SS, Kim J, Numa H, Itoh T, Buell CR, Matsumoto T (2013) Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice (N Y) 6:4\u003c/li\u003e\n\u003cli\u003eKidokoro S, Hayashi K, Haraguchi H, Ishikawa T, Soma F, Konoura I, Toda S, Mizoi J, Suzuki T, Shinozaki K, Yamaguchi-Shinozaki K (2021) Posttranslational regulation of multiple clock-related transcription factors triggers cold-inducible gene expression in Arabidopsis. Proc Natl Acad Sci U S A 118:e2021048118\u003c/li\u003e\n\u003cli\u003eKidokoro S, Shinozaki K, Yamaguchi-Shinozaki K (2022) Transcriptional regulatory network of plant cold-stress responses. Trends Plant Sci 27:922-935\u003c/li\u003e\n\u003cli\u003eKim H, Seomun S, Yoon Y, Jang G (2021) Jasmonic acid in plant abiotic stress tolerance and interaction with abscisic acid. Agronomy-Basel 11:1886\u003c/li\u003e\n\u003cli\u003eKim TH, Kim SM (2023) Identification of candidate genes for salt tolerance at the seedling stage using integrated genome-wide association study and transcriptome analysis in rice. Plants-Basel 12:1401\u003c/li\u003e\n\u003cli\u003eKuroki M, Saito K, Matsuba S, Yokogami N, Shimizu H, Ando I, Sato Y (2007) A quantitative trait locus for cold tolerance at the booting stage on rice chromosome 8. Theor Appl Genet 115:593-600\u003c/li\u003e\n\u003cli\u003eLi HB, Wang J, Liu AM, Liu KD, Zhang Q, Zou JS (1997) Genetic basis of low-temperature-sensitive sterility in indica-japonica hybrids of rice as determined by RFLP analysis. Theor Appl Genet 95:1092-1097\u003c/li\u003e\n\u003cli\u003eLi J, Zeng Y, Pan Y, Zhou L, Zhang Z, Guo H, Lou Q, Shui G, Huang H, Tian H, Guo Y, Yuan P, Yang H, Pan G, Wang R, Zhang H, Yang S, Guo Y, Ge S, Li J, Li Z (2021) Stepwise selection of natural variations at CTB2 and CTB4a improves cold adaptation during domestication of japonica rice. New Phytol 231:1056-1072\u003c/li\u003e\n\u003cli\u003eLiu C, Ou S, Mao B, Tang J, Wang W, Wang H, Cao S, Schlappi MR, Zhao B, Xiao G, Wang X, Chu C (2018) Early selection of bZIP73 facilitated adaptation of japonica rice to cold climates. Nat Commun 9:3302\u003c/li\u003e\n\u003cli\u003eLou Q, Guo H, Li J, Han S, Khan NU, Gu Y, Zhao W, Zhang Z, Zhang H, Li Z, Li J (2022) Cold-adaptive evolution at the reproductive stage in Geng/japonica subspecies reveals the role of OsMAPK3 and OsLEA9. Plant J 111:1032-1051\u003c/li\u003e\n\u003cli\u003eMehrotra S, Verma S, Kumar S, Kumari S, Mishra BN (2020) Transcriptional regulation and signalling of cold stress response in plants: An overview of current understanding. Environmental and Experimental Botany 180\u003c/li\u003e\n\u003cli\u003eMei EY, Tang JQ, He ML, Liu ZQ, Tian XJ, Bu QY (2022) OsMKKK70 negatively regulates cold tolerance at booting stage in rice. IJMS 23:14472\u003c/li\u003e\n\u003cli\u003eMi HY, Muruganujan A, Ebert D, Huang XS, Thomas PD (2019) PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res 47:D419-D426\u003c/li\u003e\n\u003cli\u003eNakashima K, Takasaki H, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K (2012) NAC transcription factors in plant abiotic stress responses. Bba-Gene Regul Mech 1819:97-103\u003c/li\u003e\n\u003cli\u003eNakashima K, Tran LSP, Van Nguyen D, Fujita M, Maruyama K, Todaka D, Ito Y, Hayashi N, Shinozaki K, Yamaguchi-Shinozaki K (2007) Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J 51:617-630\u003c/li\u003e\n\u003cli\u003eOh CS, Choi YH, Lee SJ, Yoon DB, Moon HP, Ahn SN (2004) Mapping of quantitative trait loci for cold tolerance in weedy rice. Breeding Science 54:373-380\u003c/li\u003e\n\u003cli\u003ePark MR, Yun KY, Mohanty B, Herath V, Xu FY, Wijaya E, Bajic VB, Yun SJ, De Los Reyes BG (2010) Supra-optimal expression of the cold-regulated transcription factor in transgenic rice changes the complexity of transcriptional network with major effects on stress tolerance and panicle development. Plant Cell Environ 33:2209-2230\u003c/li\u003e\n\u003cli\u003eParker AH, Wilkinson SW, Ton J (2022) Epigenetics: a catalyst of plant immunity against pathogens. New Phytol 233:66-83\u003c/li\u003e\n\u003cli\u003ePimentel H, Bray NL, Puente S, Melsted P, Pachter L (2017) Differential analysis of RNA-seq incorporating quantification uncertainty. Nature Methods 14:687-690\u003c/li\u003e\n\u003cli\u003eQi FF, Zhang FX (2020) Cell cycle regulation in the plant response to stress. Front Plant Sci 10:1765\u003c/li\u003e\n\u003cli\u003eRaj SRG, Nadarajah K (2023) QTL and candidate genes: Techniques and advancement in abiotic stress resistance breeding of major cereals. Int J Mol Sci 24:6\u003c/li\u003e\n\u003cli\u003eRaza A, Charagh S, Zahid Z, Mubarik MS, Javed R, Siddiqui MH, Hasanuzzaman M (2021) Jasmonic acid: a key frontier in conferring abiotic stress tolerance in plants. Plant Cell Rep 40:1513-1541\u003c/li\u003e\n\u003cli\u003eRen HM, Zhang YT, Zhong MY, Hussian J, Tang YT, Liu SK, Qi GN (2023) Calcium signaling-mediated transcriptional reprogramming during abiotic stress response in plants. Theor Appl Genet 136\u003c/li\u003e\n\u003cli\u003eRen\u0026aacute;k D, Gibalov\u0026aacute; A, Solcov\u0026aacute; K, Honys D (2014) A new link between stress response and nucleolar function during pollen development in mediated by AtREN1 protein. Plant Cell Environ 37:670-683\u003c/li\u003e\n\u003cli\u003eRitchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015) limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43:e47\u003c/li\u003e\n\u003cli\u003eSagehashi Y, Ikegaya T, Fujino K (2022) Integration of genetic engineering into conventional rice breeding programs for the next generation. Euphytica 218:145\u003c/li\u003e\n\u003cli\u003eSagor GHM, Inoue M, Kusano T, Berberich T (2021) Expression profile of seven polyamine oxidase genes in rice (Oryza sativa) in response to abiotic stresses, phytohormones and polyamines. Physiol Mol Biol Plants 27:1353-1359\u003c/li\u003e\n\u003cli\u003eSaito K, Hayano-Saito Y, Kuroki M, Sato Y (2010) Map-based cloning of the rice cold tolerance gene. Plant Science 179:97-102\u003c/li\u003e\n\u003cli\u003eSaito K, Miura K, Nagano K, Hayano-Saito Y, Araki H, Kato A (2001) Identification of two closely linked quantitative trait loci for cold tolerance on chromosome 4 of rice and their association with anther length. Theor Appl Genet 103:862-868\u003c/li\u003e\n\u003cli\u003eSakai H, Lee SS, Tanaka T, Numa H, Kim J, Kawahara Y, Wakimoto H, Yang CC, Iwamoto M, Abe T, Yamada Y, Muto A, Inokuchi H, Ikemura T, Matsumoto T, Sasaki T, Itoh T (2013) Rice Annotation Project Database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol 54:e6\u003c/li\u003e\n\u003cli\u003eSatake T (1976) Determination of the most sensitive stage to sterile type cool injury in rice plants. Research bulletin of the Hokkaido National Agricultural Experiment Station 113:1-43\u003c/li\u003e\n\u003cli\u003eSato Y, Masuta Y, Saito K, Murayama S, Ozawa K (2011) Enhanced chilling tolerance at the booting stage in rice by transgenic overexpression of the ascorbate peroxidase gene, OsAPXa. Plant Cell Rep 30:399-406\u003c/li\u003e\n\u003cli\u003eShimono H, Abe A, Aoki N, Koumoto T, Sato M, Yokoi S, Kuroda E, Endo T, Saeki KI, Nagano K (2016) Combining mapping of physiological quantitative trait loci and transcriptome for cold tolerance for counteracting male sterility induced by low temperatures during reproductive stage in rice. Physiol Plant 157:175-192\u003c/li\u003e\n\u003cli\u003eShirasawa S, Endo T, Nakagomi K, Yamaguchi M, Nishio T (2012) Delimitation of a QTL region controlling cold tolerance at booting stage of a cultivar, \u0026apos;Lijiangxintuanheigu\u0026apos;, in rice, Oryza sativa L. Theor Appl Genet 124:937-946\u003c/li\u003e\n\u003cli\u003eSuh JP, Jeung JU, Lee JI, Choi YH, Yea JD, Virk PS, Mackill DJ, Jena KK (2010) Identification and analysis of QTLs controlling cold tolerance at the reproductive stage and validation of effective QTLs in cold-tolerant genotypes of rice (Oryza sativa L.). Theor Appl Genet 120:985-995\u003c/li\u003e\n\u003cli\u003eSun L, Huang L, Hong Y, Zhang H, Song F, Li D (2015) Comprehensive analysis suggests overlapping expression of rice ONAC transcription factors in abiotic and biotic stress responses. Int J Mol Sci 16:4306-4326\u003c/li\u003e\n\u003cli\u003eSun ZH, Du J, Pu XY, Ali MK, Yang XM, Duan CL, Ren MR, Li X, Zeng YW (2019) Near-isogenic lines of rice revealed new QTLs for cold tolerance at booting stage. Agronomy-Basel 9:40\u003c/li\u003e\n\u003cli\u003eTakasaki H, Maruyama K, Kidokoro S, Ito Y, Fujita Y, Shinozaki K, Yamaguchi-Shinozaki K, Nakashima K (2010) The abiotic stress-responsive NAC-type transcription factor OsNAC5 regulates stress-inducible genes and stress tolerance in rice. Mol Genet Genomics 284:173-183\u003c/li\u003e\n\u003cli\u003eTakeuchi Y, Hayasaka H, Chiba B, Tanaka I, Shimano T, Yamagishi M, Nagano K, Sasaki T, Yano M (2001) Mapping quantitative trait loci controlling cool-temperature tolerance at booting stage in temperate rice. Breeding Science 51:191-197\u003c/li\u003e\n\u003cli\u003eXu LM, Zhou L, Zeng YW, Wang FM, Zhang HL, Shen SQ, Li ZC (2008) Identification and mapping of quantitative trait loci for cold tolerance at the booting stage in a rice near-isogenic line. Plant Science 174:340-347\u003c/li\u003e\n\u003cli\u003eYamamori K, Ogasawara K, Ishiguro S, Koide Y, Takamure I, Fujino K, Sato Y, Kishima Y (2021) Revision of the relationship between anther morphology and pollen sterility by cold stress at the booting stage in rice. Annals of Botany 128:559-575\u003c/li\u003e\n\u003cli\u003eYang L, Lei L, Wang J, Zheng H, Xin W, Liu H, Zou D (2023) qCTB7 positively regulates cold tolerance at booting stage in rice. Theor Appl Genet 136:135\u003c/li\u003e\n\u003cli\u003eZandalinas SI, Balfagon D, Gomez-Cadenas A, Mittler R (2022) Plant responses to climate change: metabolic changes under combined abiotic stresses. J Exp Bot 73:3339-3354\u003c/li\u003e\n\u003cli\u003eZhang H, Zhu J, Gong Z, Zhu JK (2022a) Abiotic stress responses in plants. Nat Rev Genet 23:104-119\u003c/li\u003e\n\u003cli\u003eZhang M, Zhao R, Huang K, Huang S, Wang H, Wei Z, Li Z, Bian M, Jiang W, Wu T, Du X (2022b) The OsWRKY63-OsWRKY76-OsDREB1B module regulates chilling tolerance in rice. Plant J 112:383-398\u003c/li\u003e\n\u003cli\u003eZhang Y, Xia P (2023) The DREB transcription factor, a biomacromolecule, responds to abiotic stress by regulating the expression of stress-related genes. Int J Biol Macromol 243:125231\u003c/li\u003e\n\u003cli\u003eZhang Z, Li J, Pan Y, Li J, Zhou L, Shi H, Zeng Y, Guo H, Yang S, Zheng W, Yu J, Sun X, Li G, Ding Y, Ma L, Shen S, Dai L, Zhang H, Yang S, Guo Y, Li Z (2017) Natural variation in CTB4a enhances rice adaptation to cold habitats. Nat Commun 8:14788\u003c/li\u003e\n\u003cli\u003eZhou H, Hirata M, Osawa R, Fujino K, Kishima Y (2017) Detainment of Tam3 transposase at plasma membrane by its BED-Zinc finger domain. Plant Physiol 173:1492-1501\u003c/li\u003e\n\u003cli\u003eZhu Y, Chen K, Mi X, Chen T, Ali J, Ye G, Xu J, Li Z (2015) Identification and fine mapping of a stably expressed QTL for cold tolerance at the booting stage using an interconnected breeding population in rice. PLoS One 10:e0145704\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Tables","content":"\u003cp\u003eTables 1 to 5 are available in the Supplementary Files section.\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Anther, booting stage, cold stress, flowering stage, microarray, pollen, RNA-seq, sterility, transcriptome analysis","lastPublishedDoi":"10.21203/rs.3.rs-4399503/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4399503/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eMany studies of stress tolerance in plants have characterized genes that show differences among a small number of lines with clearly distinct tolerance or sensitivity to the given stress. From the few cloned genes, it is difficult to genetically interpret intermediate tolerance or susceptibility levels and explain the complexity of stress responses and tolerance. In this study, we explored the changes in the transcriptome of anthers from 13 rice lines with different cold tolerance grown under control conditions or exposed to 4 days of cold stress to look for correlations between cold tolerance at the booting stage and expression levels. When examining the overall expression patterns in anthers at low temperature, the cold-tolerant lines tended to have relatively few highly expressed genes, and the expression levels of ribosome-related genes tended to be lower in cold-tolerant lines than in cold-sensitive lines. Importantly, we observed these different expression patterns between the cold-tolerant and -sensitive lines regardless of whether cold stress had been applied. Minimal expression changes under cold stress tended to be characteristic of the cold-tolerant lines, especially in repetitive sequences. We also identified unknown genes whose expression was cold responsive and common to all the lines studied. We conclude that rice lines whose transcriptome remains constant or insensitive in response to cold stress are more tolerant to low-temperature exposure during the booting stage than rice lines with more widespread expression changes.\u003c/p\u003e","manuscriptTitle":"The anther transcriptome of cold-tolerant rice cultivars is largely insensitive to temperature changes","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-05-22 11:20:07","doi":"10.21203/rs.3.rs-4399503/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"c6487f25-ce01-411a-9f59-7aa15253361c","owner":[],"postedDate":"May 22nd, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-06-18T08:39:35+00:00","versionOfRecord":[],"versionCreatedAt":"2024-05-22 11:20:07","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4399503","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4399503","identity":"rs-4399503","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}