A population of 128 DH lines was derived from a cross between the typical indica rice cultivar 93-11 and the japonica cultivar Milyang352 (Lee et al., 2020). In the summer of 2016, F₁ plants were generated, and their pollen was cultured on an N6Y2 medium to induce callus formation, from which regenerated plants were obtained. In the winter of 2016/2017, the regenerated A₀ generation plants were cultivated in a greenhouse. In the summer of 2017, 128 A₁ generation lines were selected to develop a DH population to assess the phenotypic data for cold tolerance during the seedling stage. Additionally, they were used to construct a molecular genetic map to identify QTLs that control survival under cold conditions during the seedling stage over 2 years (2022–2023).

In 2022 and 2023, cold water treatment of rice seedlings was conducted at the Chuncheon Substation of the National Institute of Crop Science, RDA (Chuncheon, Korea). The rice seedlings were cultivated in 72-cell trays (measuring 3.8 cm × 3.8 cm × 4.5 cm) filled with paddy soil to mimic field conditions. The seedlings were grown until the three-leaf stage before initiating cold tolerance treatment. As control varieties, ‘Satbyeolbyeo’ (cold-sensitive) and ‘Seolagbyeo’ (cold-tolerant) were used to validate the cold tolerance assessment under these conditions. The rice seedlings were placed in a cold water storage facility in Chuncheon, which measured 50 m in length and 8 m in width and could supply cold water (13°C ± 1°C) through outlets spaced 70 cm apart with an 8-mm diameter. During the experiment, the average temperature in Chuncheon was approximately 20.0°C in 2022 and 20.2°C in 2023, with relative humidity (RH) ranging from 60% and 80% RH during the same periods, reflecting typical conditions during the experimental period. After 2 weeks of growth under normal greenhouse conditions, the trays were placed in the cold water treatment facility and rinsed with cold water for 10 days. Subsequently, the seedlings were recovered for 7 days in the greenhouse (Jiang et al., 2011; Kim et al., 2021) ( Supplementary Figure S1 ). The phenotypes of the recovered seedlings were scored based on changes in leaf color, with scores ranging from 1 to 9, following the standard protocol provided by the International Rice Research Institute (IRRI), where 1 represents no change in leaf color, 5 indicates light yellow in half of the leaf, and 9 represents almost dead leaf (Donoso et al., 2015; IRRI, 2002).

Total genomic DNA was extracted from the fresh leaves of 4-week-old individual plants using the cetyltrimethylammonium bromide method, with certain modifications. DNA was quantified using a NanoDrop spectrophotometer (ND1000 spectrophotometer; Mettler Toledo, Greifensee, Switzerland). To construct the molecular map, Kompetitive allele-specific PCR (KASP) marker amplification and allelic discrimination were performed using the Nexar System (LGC Douglas Scientific, Alexandria, VA, USA) at the Seed Industry Promotion Center (Gimje, Korea) of the Foundation of Agri, Tech, Commercialization & Transfer in Korea. The KASP reaction mixture consisted of 0.8 μL of 2X master mix (LGC Genomics, London, UK), 0.02 μL of 72 KASP assay mix (LGC Genomics), and 5 ng genomic DNA template in a total volume of 1.6 μL in a 384-well array tape. KASP amplification was performed as described by Cheon et al. (2018). Fluidigm markers for single-nucleotide polymorphism (SNP) genotypes were determined using the BioMark™ HD system (Fluidigm, San Francisco, CA, USA) and 96.96 Dynamic Array IFC (96.96 IFC) chip, following the manufacturer’s protocol at the National Instrumentation Center for Environmental Management, Seoul National University (Pyeongchang, South Korea). Genotyping results were obtained using the Fluidigm SNP Genotyping Analysis software. A total of 240 SNP markers, including 106 KASP and 134 Fluidigm SNP markers, were used for genotyping the DH population. Sanger sequencing analysis was utilized to identify SNPs and insertion and deletion (InDel) regions that differ between the parental cultivars and explore candidate genes.

Linkage maps were constructed from genotype data using QTL IciMapping version 4.1 (Meng et al., 2015). Genetic distances were estimated using the Kosambi map function of the software (Kosambi, 1944). Additionally, over 2 years, QTL IciMapping was used for phenotypic data to analyze cold tolerance during the seedling stage. Inclusive composite interval mapping for additive QTL (ICIM-ADD) was used for analysis. The limit of detection (LOD) scores for ICIM-ADD were determined at an experiment-wise error rate of p > 0.05 by performing 1,000 permutations using QTL IciMapping.

The hybridization step followed the GeneTitan^® Multichannel Instrument User’s Manual using the KNU Axiom Oryza 580 K Genotyping Array. After ligation, the arrays were then stained and imaged using a GeneTitan MC Instrument (Kim et al., 2022). Image data were analyzed using the Affymetrix^® GeneChip^® Command Console^® software (Thermo Fisher Scientific, Waltham, MA, USA) to generate raw intensity values. Genotype calling was conducted using the Affymetrix Power Tools software. The BRLMM-P algorithm, which applies a Gaussian mixture model for one-dimensional clustering, was employed to ensure high accuracy in SNP genotype determination (BRLMM-P: a Genotype Calling Method for the SNP 5.0 Array). SNP markers with low call rates or poor clustering quality were excluded during quality control to ensure data reliability. For high-throughput SNP genotyping, the PolyHighResolution chip from the KNU Axiom Oryza 580 K Genotyping Array, which included 230,526 high-confidence SNP markers from 93-11 and Milyang352, was utilized. The identified QTL region was scanned to explore SNP markers exhibiting polymorphism between the parental lines. The fine-mapping analysis focused on these polymorphic SNP markers to narrow down the candidate region and improve the resolution of QTL localization. This approach allowed for the identification of specific SNP markers within the QTL region, which are critical for MAS and for identifying candidate genes underlying cold tolerance in rice.

Total RNA was isolated from the leaves during the seedling stage using an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). DNase treatment and cDNA synthesis used a QuantiTect Reverse Transcription kit (Qiagen). The reactions were triplicated, with actin expression (OsActin1 and Os03g0718100) as an internal baseline control. Samples were kept on ice until qRT-PCR analysis using QuantStudio 1 (Thermo Fisher Scientific, Boston, MA, USA). The run parameters followed a standard PCR protocol: reverse transcription at 50°C for 15 min, denaturation at 95°C for 2 min, followed by 50 amplification cycles with denaturation at 95°C for 15 s and annealing and extension at 60°C for 1 min, during which fluorescence was acquired. Data were analyzed using the QuantStudio software version 1.5.2 (Thermo Fisher Scientific, Boston). The cycle threshold was set to 0.05, and the baseline was automatically chosen. Two-tailed t-tests were performed using the GraphPad software to compare the samples. Statistical significance was set at p < 0.01.

For the haplotype analysis of the 137 KRICE_CORE population, genomic data were generated with an average coverage of approximately eightfold using the Illumina HiSeq 2500 Sequencing Systems Platform (Illumina Inc., San Diego, CA, USA). Raw reads were aligned to the rice reference genome (IRGSP 1.0) to generate genotype calls, ensuring high accuracy in the subsequent analyses. To create a reliable genotype dataset, haplotype analysis was performed using the parameters missing values <1%, minor allele frequency >5%, and heterozygosity ratio <5%, as implemented in the PLINK software. These thresholds were chosen to minimize data noise and enhance the accuracy of the haplotype associations. Phenotypic data for each variety were averaged to determine cold tolerance scores, which were subsequently used to identify haplotypes significantly associated with cold tolerance. The 137 KRICE_CORE population was also subjected to cold water treatment during the seedling stage, and the resulting cold tolerance scores were assessed following standardized protocols. A neighbor-joining tree (NJ-Tree) was constructed using the MEGA X software to visualize the genetic relationships among the varieties. The NJ-Tree was further processed in the Newick format and visualized using the Interactive Tree of Life (iTOL) tool for enhanced interpretability (Kumar et al., 2018; Letunic and Bork, 2011).

The cold tolerances of the rice varieties 93-11, Milyang352, and 128 DH lines were assessed during the seedling stage under cold water storage conditions over 2 years. Phenotypic analysis facilitated the classification of parent cultivars and DH populations into nine distinct types based on their cold tolerance scores. Across both years, the indica parent 93-11 consistently exhibited higher susceptibility to cold, with scores ranging from 6 to 7, whereas the japonica parent Milyang352 exhibited strong cold tolerance, with scores ranging from 2 to 3 ( Figures 1A–D ). In 2022, 23 DH lines were identified as cold-tolerant during the seedling stage (scores of 1–3), whereas 34 DH lines were classified as cold-susceptible (scores of 7–9) ( Figure 1C ). In 2023, a similar pattern was observed, with 14 DH lines exhibiting cold tolerance (scores 1–3) and 30 DH lines exhibiting susceptibility (scores 7–9) ( Figure 1D ). Figure 1E illustrates the variability among the DH populations across the 2 years. Despite this variability, a significant positive correlation was observed between the cold tolerance scores in 2022 and 2023 (r = 0.67, p-value < 0.01), indicating a consistent response to cold stress in the DH population across different years.

To identify the genetic loci associated with cold tolerance during the seedling stage in the DH population, a linkage map was constructed, and QTL analysis was performed using ICIM analysis. Among the KASP and Fluidigm markers tested, 229 SNPs exhibited polymorphic patterns between parents. These markers were used to construct linkage maps from genotype data using QTL IciMapping version 4.1 ( Supplementary Figure S2 ). Additionally, QTL IciMapping was used for QTL analysis using the phenotypic data, specifically through ICIM-ADD. The LOD threshold was obtained based on a permutation test (1,000 permutations, p = 0.05) for each trait.

In the QTL analysis, three QTLs for cold tolerance during the seedling stage were detected on chromosomes 3, 10, and 11 in 2022, and one QTL was detected on chromosome 10 in 2023. In 2022, the QTL qCTS3²² (LOD score = 5.09) was mapped on chromosome 3 (from 3,524 to 5,676 kb), accounting for 16.90% of the total phenotypic variance. On chromosome 10, qCTS10²² (LOD score = 5.11) was mapped from 10,730 to 15,128 kb, accounting for 16.06% of the total phenotypic variance. On chromosome 11, qCTS11²² (LOD score = 3.88) was mapped from 680 to 2,455 kb, accounting for 11.38% of the total phenotypic variance. In 2023, qCTS10²³ (LOD score = 16.70) was located on chromosome 10, at a position similar to that of qCTS10²² from 2022, accounting for 40.55% of the phenotypic variance. The allele qCTS11²² was confirmed to have been derived from 93-11, and the remaining QTLs were confirmed to have been derived from Milyang352 ( Figure 2 , Table 1 ). Based on the QTL analysis conducted over 2 years, the analysis focused on QTLs detected at the same loci (qCTS10²² and qCTS10²³ ).

To narrow down the candidate region associated with cold tolerance during the seedling stage, a fine-mapping analysis of the allele qCTS10^22/23 on chromosome 10 was conducted using a 580 K SNP array (Kim et al., 2022). Polymorphic SNPs within the 580 K SNP array were identified by comparing the sequences of the parent cultivars (93-11 and Milyang352) between markers id10002842 and KJ10_019 on chromosome 10. This analysis confirmed the presence of 44 SNPs within this allele. These 44 SNPs were subsequently used for QTL dissection, facilitating further refinement of the candidate region using phenotypic data collected over 2 years. Therefore, the interval narrowed significantly from 5 Mb to 300 kb. The fine-mapping focused on four primary SNPs between AX-155116904 and AX-295867315 that are crucial for reducing the candidate region. The refined target region now encompasses a 170-kb interval with six predicted genes (Os10g0408020, Os10g0408360, Os10g0408700, Os10g0409366, Os10g0409383, and Os10g0409400) according to The Rice Annotation Project Database (https://rapdb.dna.affrc.go.jp/) ( Figure 3 ). Among the six candidate genes of qCTS10²² and qCTS10²³ , three candidate genes (Os10g0408020, Os10g0408360, and Os10g0408700) were monomorphic between 93-11 and Milyang352. The remaining three genes (Os10g0409366, Os10g0409383, and Os10g0409400) were sequenced to identify their respective regions. The Os10g0409366 gene showed six SNPs, Os10g0409383 had one SNP, and Os10g0409400 exhibited 23 SNPs and two InDel regions between the parents ( Table 2 , Figure 4 ; Supplementary Figure S3 ).

To identify genes associated with cold tolerance at the seedling stage, real-time PCR was used to assess gene expression in leaves after 0, 12, and 24 h of cold treatment at the seedling stage. The expression levels of the candidate genes Os10g0409366, Os10g0409383, and Os10g0409400, located in the fine-mapped region associated with cold tolerance, were assessed using a primer set designed based on their coding sequences. The results demonstrated that Os10g0409400 was significantly upregulated in the Milyang352 cultivar compared to that in the 93-11 cultivar, specifically after 12 h and 24 h of cold treatment (p < 0.01). At 12 h, Milyang352 exhibited expression levels approximately 12 times higher than those of 93-11, indicating a robust transcriptional response to cold stress. This increased expression persisted at 24 h, further supporting the association of Os10g0409400 with cold tolerance in Milyang352.

In contrast, Os10g0409366 and Os10g0409383 showed no significant differences in expression levels between the two cultivars ( Figure 5 ). Additionally, 93-11 exhibited relatively low and stable expression levels of Os10g0409400 across all time points, indicating a weaker response to cold stress. These findings indicate that Os10g0409400 plays a crucial role in the cold tolerance observed in Milyang352 during the seedling stage. The significant differential expression between the two cultivars, specifically under prolonged cold stress, highlights the potential of Os10g0409400 as a primary genetic determinant for enhancing cold tolerance in rice.

To explore the genetic diversity and potential associations with cold tolerance within the KRICE_CORE rice population, we performed a haplotype analysis of the Os10g0409400 gene. We selected this gene because of its location within a QTL region associated with cold tolerance during the seedling stage. We identified five distinct haplotypes of Os10g0409400 within the 137 KRICE_CORE population based on 35 SNPs and four InDels, with specific variations highlighted in blue and mustard to represent japonica and indica types, respectively ( Figure 6A ). The haplotype analysis revealed that among the five identified haplotypes, Hap1, Hap2, and Hap3 were predominantly associated with the indica type, whereas Hap4 and Hap5 were associated with the japonica type. Specifically, the indica parent 93-11 obtained from the DH population was classified under Hap2, and the japonica parent Milyang352 was classified under Hap4. The cold tolerance scores for Hap1, Hap2, and Hap3 were 7.90, 8.01, and 7.40, respectively, indicating susceptibility to cold stress during the seedling stage. In contrast, Hap4 and Hap5 exhibited lower scores, indicating a higher tolerance to cold stress ( Figure 6B ). Neighbor-joining tree analysis further categorized the 137 rice varieties into two distinct subgroups. The first subgroup included Hap1, Hap2, and Hap3, representing indica types, whereas the second subgroup included Hap4 and Hap5, representing japonica types ( Figure 6C ). The haplotype analysis demonstrated that Hap1 and Hap2 included 18 varieties each, with Hap1 comprising 3 and 15 Aus and indica types, respectively. Hap2 had entirely indica types, whereas Hap3 included 16 varieties with five Aus types, one aromatic, and 10 indica types. Hap4, the most diverse haplotype, also included 79 varieties with two aromatic types, one admixture type, 57 temperate japonica types, and 19 tropical japonica types. Hap5 contained five varieties, all temperate japonica types ( Figure 6D ).