Navigating rice seedling cold resilience: QTL mapping in two inbred line populations and the search for genes

Due to global climate change resulting in extreme temperature fluctuations, it becomes increasingly necessary to explore the natural genetic variation in model crops such as rice to facilitate the breeding of climate-resilient cultivars. To uncover genomic regions in rice involved in managing cold stress tolerance responses and to identify associated cold tolerance genes, two inbred line populations developed from crosses between cold-tolerant and cold-sensitive parents were used for quantitative trait locus (QTL) mapping of two traits: degree of membrane damage after 1 week of cold exposure quantified as percent electrolyte leakage (EL) and percent low-temperature seedling survivability (LTSS) after 1 week of recovery growth. This revealed four EL QTL and 12 LTSS QTL, all overlapping with larger QTL regions previously uncovered by genome-wide association study (GWAS) mapping approaches. Within the QTL regions, 25 cold-tolerant candidate genes were identified based on genomic differences between the cold-tolerant and cold-sensitive parents. Of those genes, 20% coded for receptor-like kinases potentially involved in signal transduction of cold tolerance responses; 16% coded for transcription factors or factors potentially involved in regulating cold tolerance response effector genes; and 64% coded for protein chaperons or enzymes potentially serving as cold tolerance effector proteins. Most of the 25 genes were cold temperature regulated and had deleterious nucleotide variants in the cold-sensitive parent, which might contribute to its cold-sensitive phenotype.


Introduction
Rice is one of the most important crops and, due to its subtropical origin, is generally sensitive to even short exposures to cold stress at all stages of development (da Cruz et al., 2013;Li et al., 2022).Because of global climate change, recent weather extremes, such as colder winters affecting early spring rice planting, have become increasingly common (Johnson et al., 2018).Moreover, to either expand rice cultivation into higher latitudes or altitudes where air temperatures can fluctuate widely during the early growing season, it is important to elucidate rice cold stress tolerance response mechanisms that allow young rice seedlings to survive several days of exposure to chilling temperatures and resume growth and development during recovery from stress.Interestingly, because of the way rice was domesticated and spread by humans, this generally cold-sensitive crop has a wide range of temperature adaptation that is structured into subpopulations, and variation for cold tolerance within subpopulations may act through very different mechanisms (Li et al., 2022).These mechanisms can be studied genetically because there is considerable natural variation in cultivated self-pollinating plants such as rice.Previous studies by us and others showed the Japonica varietal group (VG), composed of temperate japonica (TEJ) and tropical japonica (TRJ) accessions, is considerably more cold tolerant than the Indica VG, composed of aus (AUS) and indica (IND) accessions, which allowed us to map various cold tolerance trait QTL using genome-wide association study (GWAS) mapping approaches (Schläppi et al., 2017;Shakiba et al., 2017;Shimoyama et al., 2020;Phan and Schläppi, 2021;reviewed in Li et al., 2022).Since many of these QTL cover relatively large genomic regions containing many genes, a recurrent challenge is to identify cold tolerance candidate genes and ultimately the causal genes and their alleles associated with those QTL.
The purpose of this study was to fine-map cold tolerance trait QTL using bi-parental mapping populations to identify smaller genomic regions that overlap with the larger QTL regions we previously identified through GWAS, thus allowing us to narrow down cold tolerance candidate genes within those smaller regions.Based on our previous results, we selected two robustly cold-tolerant and coldsensitive parents to generate recombinant inbred lines (RILs) and backcross inbred lines (BILs).The two populations were used to map agricultural traits as "quality control" to assess whether known heading date and plant height QTL and the associated genes could be uncovered.Young RIL and BIL seedlings were then exposed for 1 week to a chilling temperature of 10°C, after which the degree of membrane damage in the leaves of exposed plants was measured, while seedlings' survival after 1 week of recovery growth was determined.Our results validated the two mapping populations suitable for QTL mapping, and 16 cold tolerance trait QTL were uncovered, all of which overlapped with at least one of the QTL we found by GWAS mapping.The genomic regions delineated by those QTL contained at least 25 candidate genes with polymorphisms between the cold-tolerant and cold-sensitive parents, and their potential contributions to rice cold tolerance are discussed.

Rice materials and inbred line production
The three accessions used to generate inbred lines are part of the rice mini-core collection (RMC).Seeds were obtained from the Genetic Stocks-Oryza (GSOR) collection located at the USDA-ARS Dale Bumpers National Rice Research Center (DBNRRC; Stuttgart, Arkansas, USA; https://www.ars.usda.gov/GSOR).
To generate recombinant inbred lines (RILs), the cold-tolerant temperate japonica accession Krasnodarskij 3352 (GSOR311787; originally from the Krasnodar region of Western Russia) was crossed as female with the cold-sensitive aus accession Carolino 164 (GSRO311654; originally from the Chad region of Africa) to generate F 1 seeds.The F 2 seed produced from 22 F 1 plants was advanced by single seed descent to the F 8 generation, with 12 RILs only advanced to the F 7 , for a final population of 140 RILs.Leaf tissue was collected from the 140 F 7:8 RILs and the F 8:9 seed from 134 (or fewer, if limited seed) of these RILs was used to grow the seedlings for phenotyping.
To generate backcross inbred lines (BILs), the cold-tolerant temperate japonica accession WIR 911 (GSOR311685; originally from the Primorsky Krai region of Eastern Russia) was crossed as female with Carolino 164, and the resulting F 1 was backcrossed as female with Carolino 164, producing 93 BC 1 F 1 progeny that were advanced by single seed descent to the BC 1 F 5 generation, with three BILs only advanced to the BC 1 F 4 .Leaf tissue from these 93 BC 1 F 4:5 BILs was used for genotyping, and BC 1 F 5:6 seed from 92 (or fewer, if limited seed) of these BILs was used to grow the seedlings for phenotyping.

Germination and standard growth conditions
For phenotyping, seeds of the three parents, F 8:9 RILs and BC 1 F 5:6 BILs were germinated in the dark for 2 days at 37°C in deionized water containing 0.1% bleach to prevent bacterial contamination.Germinating seeds were transferred into PCR strips, placed into pipette tip boxes, and grown hydroponically in deionized water for 10 days in a growth chamber under 12-h light (approximately 150 µE photon flux)/12-h dark and 28°C day/25°C night cycles.To provide nutrients, on day 10, the water was replaced with ¼ Murashige-Skoog basal salt liquid medium.Each line was represented by up to eight plants per box in quadruplicates for up to 32 plants (four boxes of eight plants) per experiment.The four boxes were randomly arranged within the growth chamber.Each box contained 11 strips of RILs or BILs and one strip with the parents as controls.For the RILs, there were four Carolino 164 seedlings, the cold-sensitive control, and four Krasnodarskij 3352 seedlings, the cold-tolerant control.For the BILs, the cold-tolerant parent, WIR 911 was used as the control.

Chilling stress treatment
The four boxes containing 2-week-old seedlings at the 2-leaf stage were placed at random positions within a growth chamber set at a constant 10°C ± 1°C temperature and incubated for 7 days (12h light/12-h dark cycles).Watering was done every other day.At the end of the treatment, approximately 1 cm of the middle segments of the second leaf for three seedlings per line was collected for electrolyte leakage (EL) assays, and the boxes were returned to standard growth conditions for a recovery period of 7 days after which low-temperature seedling survivability (LTSS) was recorded (details below).

Phenotyping of BIL and RIL populations
Agronomic and cold tolerance trait data for the RIL population are summarized in Supplementary Table S1.Data for the BIL population are summarized in Supplementary Table S2.

Heading date
Seeds of the F 8:9 RIL and BC 1 F 5:6 BILs were germinated and grown under standard growth conditions to generate 14-day-old seedlings.During May 2021, three healthy seedlings per line were transplanted in "hills" into 1.22-m × 1.22-m raised bed paddies on a rooftop at Marquette University in Milwaukee, Wisconsin, 100 lines per paddy (10 rows × 10 columns; 10 cm ± 1 cm distance between seedling "hills").The heading date was recorded as the number of days from germination to panicle emergence.

Plant height
The height of the main shoot of each RIL or BIL was measured in cm from the soil level to the tip of the mature panicle.

Electrolyte leakage
At the end of the 7-day 10°C stress period, approximately 1 cm of the middle section of the second leaf from three individual seedlings per RIL, BIL, or control per box was collected and cut into three equally sized segments.The pieces were washed in deionized water and transferred into three different screw-cap glass tubes filled with 5 mL of deionized water (conductivity ≤ 2 µS/cm), then shaken at 200 rpm for 60 min at room temperature to release cellular electrolytes from low-temperature damaged tissues.The initial conductivity of the three replicates per box (a total of 12 replicates across the four randomly distributed boxes) was measured by taking out 120 µL of the solution and adding it into the well of a hand-held LAQUAtwin B-771 conductivity meter (Horiba Scientific, Japan).Leaf samples were boiled for 20 min after the initial reading to release total cellular electrolytes.Samples were shaken again at 200 rpm for 30 min after cooling to room temperature, and the final conductivity reading was done.The mean percent EL was determined as [(initial conductivity reading)/ (final conductivity reading)] × 100.

Low-temperature seedling survivability
At the end of the 7-day 10°C stress period, seedlings were allowed to recover at standard growth conditions for one week after which seedling survival was determined by visual inspection.Seedlings that were mostly green and stiff were scored as alive while seedlings that were mostly wilted and/or bleached and soft were scored as dead.The mean percent survivability was calculated as [(number of seedlings scored as alive)/(total number of stressed plants)] × 100.

Statistical analysis of traits
QTL analysis prefers phenotype values for every (RIL or BIL) line in the mapping population.Line means were calculated for the heading date and plant height traits.For the EL and percent LTSS traits, these values were calculated from a linear mixed model (LMM) to obtain the Best Linear Unbiased Predictions (BLUPs) for each line.The LMM to calculate BLUPs for EL and LTSS traits used an augmented design with a fixed effect variable (group) indicating whether the line is a control or RIL (or BIL), a random effect variable for the line nested in the group variable, and random effect variables for the set (growth chamber experiment date) and box.In addition to treating LTSS as a percent value, LTSS was used in a survival analysis with a binomial generalized mixed model (GLMM) and logit link function to predict the probability of survival under cold treatment for each line on the log-odds scale, because percent EL and LTSS data were not normally distributed (Supplementary Figures S1C, D).The binomial GLMM used the same fixed and random variables as the LMM used for EL and percent LTSS.

Genotyping of recombinant inbred lines and QTL mapping
Leaf tissue was harvested from the parents and an individual representative F 7:8 RIL and BC 1 F 4:5 plant and lyophilized.Genotyping was done using the Cornell-IR LD Rice Array (C7AIR) containing 7,098 SNP markers (Morales et al., 2020) that was commercially available as an Illumina Infinium array.For DNA extraction and genotyping, lyophilized leaf tissue was sent to Eurofins Diagnostics, Inc. (www.eurofinsgenomics.eu/en/genotyping-gene-expression/service-platforms/illuminaarrayplatforms/).

GWAS analysis of the RDP1 and RMC
An earlier GWAS study by Shimoyama et al. (2020) phenotyped 354 accessions included in the Rice Diversity Panel 1 (RDP1) for EL and LTSS.Taking advantage of the increased number of SNPs available through the imputation of the RDP1 genotypes (Wang et al., 2018), the phenotypic data from Shimoyama et al. (2020) was rerun with the imputed genotypes using the mixed linear model in Tassel V 5.0 (Bradbury et al., 2007).The traits examined were EL8, EL10, EL12, LTSS8, LTSS10, and LTSS12, with the numbers 8, 10, and 12 representing the temperature at which the plants were grown, 8°C, 10°C, and 12°C, respectively.To more effectively capture the genetic diversity, seven panels based on subpopulation(s) were used.All panels were filtered to exclude heterozygous markers and markers with a minor allele frequency of less than 5%.The first panel, "353," was composed of 353 accessions with 2,721,379 SNPs and contained all five rice subpopulations: indica (IND), aus (AUS), aromatic (ARO), tropical japonica (TRJ), temperate japonica (TEJ), and the admixtures.The subpopulation panels included IND (72 accessions) with 1,637,548 SNPs, AUS (51 accessions) with 1,675,068 SNPs, TRJ (89 accessions) with 921,375 SNPs, and TEJ (86 accessions) with 569,119 SNPs.The ARO subpopulation contained only a limited number of accessions and was included with the 353 set but not as its own panel.The Indica VG panel, INDAUS, included both IND and AUS (129 accessions) with 2,389,503 SNPs, and the Japonica VG panel, JAP, included TRJ and TEJ (204 accessions) with 1,090,261 SNPs.GWAS run conditions as well as PCA and kinship generation were the same, as summarized in Eizenga et al. (2023).Schläppi et al. (2017) conducted GWAS mapping of LTSS in the RMC utilizing the 157 markers (148 SSRs) genotypes that were available at that time.Subsequently, the RMC was resequenced, which significantly improved the marker density (Wang et al., 2016;Huggins et al., 2019).For this study, the LTSS phenotypes from five trials (Schläppi et al., 2017) were rerun utilizing the denser SNP genotypes.Eight panels were used to capture the genetic diversity present in specific subpopulations and groups.All panels were filtered to exclude heterozygous markers and markers with a minor allele frequency of less than 5%.The panel "All Lines" contained all 173 RMC accessions genotyped with 3,224,845 SNPs.For the subpopulation panels, IND had 58 accessions genotyped with 1,883,603 SNPs, AUS had 30 accessions genotyped with 1,753,773 SNPs, TRJ had 28 accessions genotyped with 1,067,722 SNPs, and TEJ had 30 accessions genotyped with 691,001 SNPs.The ARO subpopulation only contained 11 accessions, so it was not an individual panel but was included in the "All Lines" panel and in the JAPARO panel described below.Panels of combined subpopulations were Indica (INDAUS) composed of IND and AUS having 92 accessions genotyped with 2,477,397 SNPs, Japonica (JAP) composed of TRJ and TEJ having 61 accessions genotyped with 1,092,889 SNPs, and JAPARO composed of JAP and ARO having 77 accessions genotyped with 1,402,000 SNPs.A mixed linear model was used in Tassel V 5.0 (Bradbury et al., 2007) for GWAS.GWAS run conditions, PCA, and kinship matrix are summarized in Eizenga et al. (2023).
GWAS results were run through Perl scripts to identify associated chromosome regions from individual SNPs or groups of physically linked SNPs (Huggins et al., 2019).Chromosome regions were set to 50 kilobases (kb) in both directions around each individual significant SNP and were extended to include nearby significant SNPs occurring within 200 kb.A Perl script was used to designate a "Peak SNP" for each region, which corresponded to the SNP with the most significant p-value found within the region.

Description of RIL and BIL populations
Four overall cold-tolerant (CT) TEJ accessions and four overall cold-sensitive (CS) accessions from the Indica VG, including AUS and IND accessions, were selected from the RMC-based GWAS mapping of the RMC for cold tolerance (Schläppi et al., 2017).Crosses were attempted between four pairs of these CT and CS accessions.Figure 1 shows the phenotypes of a pair of CS and CT parents where 2-week-old seedlings were exposed for 1 week to constant chilling temperatures of 10°C, after which they were allowed to recover at warm temperatures for 1 week.In the example shown, most CS seedlings did not survive; they were bleached and started to droop during the recovery period.Conversely, most of the CT seedlings survived, remained green, and resumed growth and development during the recovery period.
Due to major flowering time and sterility issues with two of the pairs, only two crosses were selected for further population development.A recombinant inbred line (RIL) mapping population was developed from CT Krasnodarskij 3352 (TEJ) crossed with CS Carolino 164 (AUS), advanced by single-seed descent to the F 8:9 generation, for a final population of 140 RILs.To generate a backcross inbred line (BIL) mapping population with reduced sterility issues, WIR 911 (TEJ) was crossed with Carolino 164, and F 1 plants were backcrossed as female with Carolino 164.The backcross progeny was advanced by single-seed descent to the BC 1 F 4:5 , for a final population of 93 BILs.Leaf tissue from single plants of the 140 F 7:8 RILs and 93 BC 1 F 4:5 BILs was genotyped, and after curation, more than 2,500 SNPs were retained.For QTL mapping of the RIL population, 2,794 SNPs with an average of 233 SNPs per chromosome (range: 170-346 SNPs) were used (Supplementary Table S3), and for BIL population mapping, 2,590 SNPs with an average of 216 SNPs per chromosome (range: 154-324 SNPs) were used (Supplementary Table S4).For both populations, chromosome (chr.)9 had the fewest, and chr. 1 the most SNPs.

Mapping potential validation of RIL and BIL populations
To perform a "quality control" for the bi-parental mapping potential of the RIL and BIL populations, two-week-old chambergrown seedlings were transplanted into raised bed rooftop paddies at Marquette University (43°02′11.59″N/87°55′51.64″W).Heading date (HD) and plant height (PTHT) data indicated that the parents were at the extremes of each end of the normal HD distributions, with only a few "transgressive" plants, while the normal PTHT distributions showed increased transgressive variation (Supplementary Figures S1A, B; Tables 1, 2).QTL mapping of HD and PTHT QTL revealed eight qHD and 12 qPTHT with logarithm of odds (LOD) scores >3.0 (Table 3).The general locations of the RIL HD QTL, qHD3 and qHD7, and BIL HD QTL, qHD3-1, qHD3-2, qHD3-3, and qHD6 (Figure 2) matched previously obtained legacy HD QTL (Takahashi et al., 2001), which validated the RIL and BIL populations for QTL mapping.
This conclusion was further supported by plant height QTL mapping.Of the 10 PTHT QTL, two contained known genes affecting plant height (Supplementary Table S5).qPTHT1 and qPTHT5-2 uncovered the chitin-inducible gibberellin-responsive gene, plant height 1 (ph1), and the gibberellin 2-beta-dioxygenase gene, OsGAox4, respectively, and nonsynonymous amino acid   substitutions might contribute to plant height differences between the parents (Supplementary Table S6; Kovi et al., 2011;He et al., 2022).Taken together, heading date and plant height QTL mapping using the newly developed RIL and BIL populations described here yielded several QTL with high LOD scores containing known flowering time and plant height-regulating genes, thus validating their usefulness for mapping cold tolerance QTL and narrowing down associated cold tolerance genes.

Cold tolerance trait mapping
EL assays were done on 2-week-old seedlings after 1 week of exposure to 10°C as a measure of plasma membrane integrity right after the cold temperature exposure.After 1 week of recovery growth, LTSS was determined as a measure of the ability of coldstressed plants to maintain photosynthesis and resume growth and development.The chilling temperature and exposure time were chosen, because based on previous studies (Schläppi et al., 2017), the cold-tolerant parents had >90% survival while the cold-sensitive parent had approximately 10% survival, thus avoiding a scenario where most cold-sensitive RILs and BILs would have 0% survival.The RIL mapping yielded nine QTL with LOD scores of 2.7 or higher (three EL and five unique LTSS QTL; Table 4), and the BIL mapping also yielded nine QTL with LOD scores of 3.2 or higher (one EL and seven unique LTSS QTL; Table 4).Of the 16 unique cold tolerance trait QTL, 11 were found in clusters on chr.3, chr.6, and chr.8, while a single QTL was found on chr. 1, chr. 4, and chr.11 (Figure 2).Strikingly, all 16, in other words, 100% of the QTL (Supplementary Table S7) overlapped with seedling stage cold tolerance GWAS QTL discovered in the RMC or RDP1 (Schläppi et al., 2017;Shimoyama et al., 2020;Phan and Schläppi, 2021).Of note, for this study, GWAS QTL for the RMC was mapped more precisely using the phenotypes reported by Schläppi et al. (2017) and genotypes based on 3,224,845 SNPs (Wang et al., 2016;Huggins et al., 2019).Similarly, the precision of the GWAS mapping in the RDP1 was improved using the seedling stage cold tolerance phenotypes reported by Shimoyama et al. (2020) and genotypes based on 4,829,392 SNPs (Wang et al., 2018).These 16 overlapping QTL were analyzed to identify cold tolerance genes within the QTL regions.
The EL and LTSS QTL identified using the biparental RIL and BIL populations were smaller than those found by GWAS mapping in the RMC and RDP1, thus reducing the number of genes to investigate (Table 5; Schläppi et al., 2017;Shimoyama et al., 2020;Phan and Schläppi, 2021).Therefore, an analysis of genomic variation of nontransposable element genes between the CT and CS parents using a public database (Zhao et al., 2015) and literature searches were done to narrow down cold tolerance candidate genes within the QTL regions and 100,000 bp to the left and right marker SNPs, which uncovered 25 cold tolerance candidate genes in 15 QTL (Table 5).The SNP ID numbers, nucleotide changes, and SNP impacts in alleles of the cold-sensitive parent Carolino 164, such as nonsynonymous amino acid changes, frameshifts, premature stop, loss of start, loss of stop, and splice variants, are summarized in Table 6.
t h e m : t h e P -t y p e I I B C a 2 + A T P a s e g e n e O s A C A 6 (LOC_Os01g71240), the F-box protein gene OsFBL27 (LOC_Os06g06050), the serine esterase encoding gene LOC_Os06g06080, the mitogen-activated protein kinase gene OsMPK6 (LOC_Os06g06090), and the pentatricopeptide repeat type restorer-of-fertility gene RF6 (LOC_Os08g01870).OsACA6 was previously shown to enhance cold tolerance in transgenic tobacco (Huda et al., 2014), and due to nucleotide variations in the The RIL population had 134 progeny genotyped with 2,794 SNP markers and the BIL population had 92 progeny genotyped with 2,590 SNP markers.a Quantitative trait loci (QTL) are declared based on the Inclusive Composite Interval Mapping (ICIM) approach using the IciMapping software (Meng et al., 2015) and 0.5-step, 0.005-pin ICIM_Add settings unless otherwise noted.b A positive additive effect indicates that the Krasnodarskij 3352 or WIR 911 allele increases the phenotype for that trait; a negative additive effect indicates that the Carolino164 allele increases the trait.c PVE is the percentage of total phenotypic variation explained by an individual QTL, as estimated by R 2 values from ICIM analysis.d LOD score is close to the 90% confidence value of 2.774 (based on 1,000 permutations for the EL(BLUP) trait).e LOD score is >90% confidence value of 2.858 and <95% value of 3.205 (based on 1,000 permutations for the LTSS(logit) trait).
f LOD score is close to the 90% confidence value of 2.876 (based on 1,000 permutations for the LTSS(percent) trait).
g LOD score is close to the 95% confidence value of 3.384 (based on 1,000 permutations for the LTSS(percent) trait).Schläppi et al. 10.3389/fpls.2023.1303651Frontiers in Plant Science frontiersin.orgpromoter region, there might be differences in gene expression levels between CT and CS accessions.OsFBL27 activates OsMYB30, a negative regulator of cold tolerance, and OsFBL27 is downregulated by OsmiR528 (Tang and Thompson, 2019).Observed amino acid duplications and added nucleotides in the CS parent might enhance OsMYB30 activation due to increased OsFBL27 activity and/or lack of OsmiR528-mediated downregulation.OsMPK6 was previously shown to be cold inducible (Lee et al., 2008), and the 233 SNP variants between Carolino 164 and Krasnodarskij 3352 might affect kinase activity or regulation.Previously, other Ca 2+ ATPase, F-box protein, and mitogen-activated protein kinase genes were shown to be involved in cold stress tolerance response mechanisms leading to the production of compatible osmolytes that might stabilize biological membranes (reviewed in Li et al., 2022).The candidate genes uncovered here might have similar functions, which need to be tested in future studies.RF6 might be involved in nucleotide and nucleic acid metabolic processes in the mitochondrion (Huang et al., 2015), affecting reactive oxygen species (ROS) production, while the involvement of serine esterases in cold tolerance needs to be addressed in future studies.
The remaining 20 candidate genes were associated with 11 LTSS QTL, as follows: The two QTL on chr. 3 (qLTSS3-1 and qLTSS3-2) had four ankyrin-tetratricopeptide repeat (ANK-TPR) genes (OsANK; LOC_Os03g47460, LOC_Os03g47686, LOC_Os03g47693, LOC_Os03g47720) and the microsomal glutathione S-transferase gene OsGST3 (LOC_Os03g50130).The single QTL on chr. 4 (qLTSS4) had the two wall-associated kinase genes OsWAK38 (LOC_Os04g29680) and OsWAK40 (LOC_Os04g29790).The four QTL on chr.6 (qLTSS6-1, qLTSS6-2, qLTSS6-3, qLTSS6-4) had eight candidate genes: the receptor-like kinase gene OsRLK-803 (LOC_Os06g04370); the aminomethyltransferase gene LOC_Os06g04380; the histone chaperone gene OsNAPL2 (LOC_Os06g05660); the transferase domain-containing protein gene LOC_Os06g05790; the two peptidyl-prolyl cis/trans isomerase genes OsFKBP13 (LOC_Os06g45340) and CYP59A (LOC_Os06g45900); and the NAC transcription factor gene OsNAC011 (LOC_Os06g46270) and the glycosyl hydrolase family 31 gene LOC_Os06g46284.The three QTL on chr.8 (qLTSS8, qLTSS8-1, qLTSS8-2) had the transcription factor gene OsbHLH070 (LOC_Os08g08160), the protease inhibitor/ seed storage/LTP family protein precursor gene LTPL130 (LOC_Os08g27674), and the receptor-like cytoplasmic kinase gene OsRLCK253/CRINKLY4 (LOC_Os08g28710).The single QTL on chr.11 (qLTSS11) had the F-box genes OsFBX398 (LOC_Os11g09360) and the protein disulfide isomerase gene OsPDIL1-1 (LOC_Os11g09280).These 20 genes can be categorized into three general groups: (1) genes with kinase activity potentially involved in cold-mediated signal transduction; (2) genes that potentially regulate gene expression; and (3) genes that code for potential cold tolerance effectors such as chaperones or enzymes.Many of the 20 cold tolerance candidate genes were previously shown to be involved in stress tolerance responses.Specifically, in the kinase category, OsWAK38 and OsWAK40 upregulation was specifically shown to be associated with salt tolerance (Wang et al., 2021), and OsRLCK253 was previously shown to be the only candidate gene associated with QTL_19 (Le et al., 2021) because of its involvement in water-deficit improvement and salinity stress tolerance (Giri et al., 2011).In the gene expression category, OsNAC011 was specifically shown to regulate drought tolerance in rice (Fang et al., 2014), and OsFBX398 was previously identified as a candidate gene selected during the breeding of rice for adaptation to different environments in Vietnam (Higgins et al., 2022).In the potential cold tolerance effector category, chaperonetype ANK genes were shown to have numerous functions in plants (Huang et al., 2009), including roles in response to abiotic and biotic stresses, and histone chaperone NAP genes were previously shown to have a role in abiotic stress responses (Tripathi et al., 2016).For The candidate genes located in the QTL regions are listed.a The Mb position of the flanking markers for the RIL or BIL QTL interval is based on the Mb position of the markers that were closest to the cM position of the QTL region as defined by the QTL analysis and listed in Table 3. b Gene nomenclature followed the standardized nomenclature for rice genes used in Oryzabase (Yamazaki et al., 2010) enzyme-encoding genes, OsGST3 was specifically shown to respond to various stress hormones (Chen et al., 2023), and OsPDIL1-1 overexpression in transgenic rice was shown to enhance abiotic stress tolerance (Xia et al., 2018), possibly by controlling the production of ROS (Zhao et al., 2023).It was also shown that the aminomethyltransferase gene LOC_Os06g04380 was differentially expressed after arsenic treatment (Norton et al., 2008), that peptidyl-prolyl cis/trans isomerase genes have various functions in response to environmental stresses (Ahn et al., 2010), and that LTPL130 was previously selected as a candidate gene located on chr.8 at 16.86 Mb in a QTL, a QTL for several root and yield-associated traits under aerobic and irrigated conditions (Padmashree et al., 2023).Besides allelic differences between CT and CS rice accessions, another criterion for cold tolerance candidate genes is that their steady-state mRNA levels are up-or downregulated in response to cold temperatures.Strikingly, a search of publicly available databases determined that of the 24 candidate genes with expression data, 22, in other words, 91.7% were regulated by cold temperatures (Table 7).Of those 22, eight (i.e., 36.4%) were up-and 11 (i.e., 50.0%) were downregulated, while three (i.e., 13.6%) were either up-or downregulated in different data sources.Because genetic backgrounds, degree of cold, and lengths of cold exposures were different between data sets (Supplementary Table S8), it will be necessary to reexamine expression patterns of the cold tolerance candidate genes under the conditions we used for QTL mapping.However, it is conceivable that stably up-and downregulated genes might have positive and negative effects, respectively, in cold stress tolerance responses in rice.Such genes might be differentially regulated in CT and CS rice accessions, which, together with coding region SNP variants, would be another reflection of allelic differences between accessions with varying cold tolerance potentials.

Conclusions
In this study, we used two biparental mapping populations made from crosses between cold-tolerant and cold-sensitive parents (Figure 1) to fine-map cold tolerance QTL.We first performed "quality control" QTL analyses to demonstrate that heading date and plant height QTL contained known genes for those traits, such as Ehd4, Ghd7, Hd1, and OsGA3ox4, thus validating that the populations were suitable for cold tolerance trait QTL mapping (Table 3; Supplementary Table S5).All the 16 cold tolerance QTL identified here overlapped with the QTL we previously obtained through GWAS mapping approaches.This allowed us to narrow down the size of previously mapped QTL and identify 25 cold-tolerance candidate genes with mostly highimpact nucleotide variants between the cold-tolerant parents (Krasnodarskij 3352 and WIR 911) and the cold-sensitive parent (Carolino 164).Interestingly, most of those genes could be assigned to modules at the seedling stage that regulate osmolytes, ROS, and growth and development, as described in a recent rice cold tolerance review (Li et al., 2022).Specifically, five genes (20%) encoded receptor-like kinases that might be involved in the signal transduction of cold stress tolerance responses.That is, OsMPK6, OsRLK-803, OsRLCK253, OsWAK38, and OsWAK40 might be part of modules regulating osmolyte and/or ROS production, together with the two (8%) transcription factor genes OsbHLH070 and OsNAC011, and the two (8%) F-box genes OsFBL27 and OsFBX398 that might regulate transcription factor activities.The other 16 (64%) are potential cold stress tolerance effector genes that code for the following proteins: three (12%) chaperons that might help maintain biological macromolecule structure during cold stress; four (16%) ANK-TPR proteins that mediate interactions with protein partners of protein complexes with potential roles in cold stress tolerance responses; and nine (36%) enzymes involved in cellular metabolisms such as ROS production, maintenance of mitochondrial and/or chloroplast integrity, or functioning as ATPases, esterases, isomerases, hydrolases, or transferases.Of those, the two peptidyl-prolyl cis/trans isomerase genes OsFKBP13 and CYP59A might belong to a new module regulating the trade-off between stress response and growth and development, as previously shown for OsCYP20-2 (Ge et al., 2020).These candidate genes can be functionally analyzed in future studies using genomics, molecular genetics, and biochemical approaches, while the QTL regions containing those genes can be used for marker-assisted breeding of coldtolerant rice cultivars.
reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
FIGURE 1 Phenotypes of two rice varieties with varying cold tolerance potentials.(A) Two-week-old hydroponically grown rice seedlings.The cold-sensitive control Carolino 164 (aus) is shown on the left (red seedlings), and the cold-tolerant check Krasnodarskij 3352 (temperate japonica; green seedlings) is shown on the right.(B) Phenotype after 1 week of 10°C exposure and 1 week of warm temperature recovery growth.

TABLE 1
Summary statistics (overall mean, SE, range of progeny) of the four traits measured on the 134 Krasnodarskij 3352 × Carolino 164 progeny lines included in the quantitative trait locus analysis (Tables3, 4) and parents.
aThe SE of the mean was calculated as the SD divided by the square root of the number of entries.b Number of days from germination to panicle emergence.

TABLE 2
Summary statistics (overall mean, SE, range of progeny) of the four traits measured on the 92 (WIR 911 × Carolino 164) × Carolino 164 backcrossed progeny lines included in the quantitative trait locus analysis (Tables3, 4) and parents.
Trait distributions are shown in Supplementary FigureS1B.a The SE of the mean was calculated as the SD divided by the square root of the number of entries.b Number of days from germination to panicle emergence.
(Meng et al., 2015)had 134 progeny genotyped with 2,794 SNP markers and the BIL population had 92 progeny genotyped with 2,590 SNP markers.aQuantitativetrait loci (QTL) are declared based on the Inclusive Composite Interval Mapping (ICIM) approach using the IciMapping software(Meng et al., 2015)and 0.5-step, 0.005-pin ICIM_Add settings unless otherwise noted.b A positive additive effect indicates that the Krasnodarskij 3352 or WIR 911 allele increases the phenotype for that trait; a negative additive effect indicates that the Carolino164 allele increases the trait.c PVE is the percentage of total phenotypic variation explained by an individual QTL, as estimated by R 2 values from ICIM analysis.d LOD score is close to the 90% confidence value of 3.01 (based on 1,000 permutations for the PTHT trait).

TABLE 5 Continued
(Kawahara et al., 2013) Symbolization, Nomenclature, and Linkage gene symbol was used, if available.cThepseudomoleculeMbposition of the candidate gene in the Rice Genome Annotation Project Release 7(Kawahara et al., 2013).dRiceGenome Annotation Project locus identifier for the candidate gene(Kawahara et al., 2013).

TABLE 6
High impact SNP variants in cold-sensitive aus line Carolino 164 compared to cold-tolerant japonica line Krasnodarskij 3352.
a Variant in Krasnodarskij 3352; Carolino 164 same as the Nipponbare reference genome.b Total of 233 variants between Krasnodarskij 3352 and Carolino 164.

TABLE 7
Cold temperature-regulated gene expression of QTL-associated candidate genes.