Identification of Candidate Genes for Cold Tolerance at Seedling Stage by GWAS in Rice (Oryza sativa L.)

1. Introduction

Rice is one of the most significant food crops, feeding more than half of the world's population and accounting for 30% of total grain output [1]. However, as the global climate changes, the frequency and intensity of severe temperatures are increasing, and very high and low temperatures have a significant detrimental influence on rice production, resulting in a decline in rice quality and output. China's rice planting region is vast, and low temperatures will affect practically all rice growing locations to variable degrees. As a result, low temperature has emerged as one of the primary abiotic stresses impacting rice growth and yield in China [2]. Therefore, the discovery of cold tolerance (CT) genes in rice, as well as the identification of signal transduction pathways and regulatory networks in response to low temperature, will provide theoretical guidance and application approaches for the study of molecular mechanisms of CT and rice genetic improvement.

More than 250 CT QTLs have been found by some domestic and international studies, and the research revealed that the majority of these QTLs came from japonica (geng) rice [3]. Additionally, more than 90 CT QTLs were identified at seedling stage, including 71 coming from geng rice, 20 from indicate (xian) rice, and 7 from wild rice [3]. According to data from the National Rice Data Center (https://www.ricedata.cn/), 113 CT-related genes had currently been cloned, mainly including the CTB4a [2], qLTG3-1 [4], HAN1 [5], and Ctb1 [6], etc. Zhang et al. (2017) [2] investigated a cold-tolerant QTL in NIL1913 and cloned the gene CTB4a, which encoded a conserved receptor-like kinase with a rich leucine repeat. Furthermore, the authors further found that CTB4a could interact with the β subunit AtpB of ATP synthase to positively regulate the content and activity of ATP under low temperature stress, thereby increasing pollen fertility, improving rice seed setting and yield, and consequently controlling rice tolerance to low temperature stress [2]. qLTG3-1 was located on chromosome 3 and encoded a protein product composed of 184 amino acids that was localized in the cytoplasm and nucleus [4]. During seed germination, qLTG3-1 was specifically expressed in the aleurone layer of the seed coat and the epiblast covering the coleoptile, which might improve the germination potential of seeds at low temperatures by regulating the cellular vacuolization of these tissues and thereby inducing the relaxation of these tissues [4]. Mao et al. (2019) [5] used a recombinant inbred line population constructed from rice cultivars 02428 and Teqing to clone a seedling CT QTL HAN1, namely LOC_Os11g29290 on chromosome 11, which encoded an oxidase that could catalyze the conversion of active JA-Il into inactive 12 OH-JA-Ile, and negatively regulated CT in rice. Saito et al. (2001) [7] identified two closely linked CT QTLs Ctb1 and Ctb2 related to anther length on chromosome 4 in the cold-resistant rice variety Norin-PL8. The Ctb1 was precisely located and the CT gene LOC_Os04g52830, encoding an F-box protein, was identified, which interacted with the E3 ubiquitin ligase component Skp1 and was implicated in the ubiquitin-proteasome pathway [6,7,8].

Genome-wide association study (GWAS) is a method for discovering genetic variants affecting complex traits by using single nucleotide polymorphism (SNP) in the genome [9]. CT associated loci can be quickly located in natural populations by GWAS. For instance, 67 QTLs for low temperature tolerance of seedlings were identified by GWAS in 295 rice RDP1 varieties, and the correlation between the OsRYH1 gene and CT was determined [10]. A total of 132 loci, affecting CT, were identified by GWAS using 529 rice materials [11]. Among them, the gene OsMYB2 related to CY was mapped, and the analysis of the gene haplotype revealed the xian-geng differentiation of this gene [11]. By GWAS, 53 QTLs related to seed germination under low temperatures were identified using a natural rice population of 187 accessions [12]. Researchers discovered that OsSAP16 was a potential gene influencing cold seed germination. OsSAP16 function loss decreased cold germination, but OsSAP16 over-expression increased low-temperature germination [12]. A large number of QTLs related to CT were discovered by GWAS, and there were a lot of candidate genes identified.

Based on DNA diversity, rice can be classified into 12 ecological types, including GJ-adm, GJ-subtropical, GJ-temperate, GJ-tropical, XI-1A, XI-1B, XI-2, XI-3, XI-adm, aus, basmati and admix [13,14,15]. Overall, the geng (GJ) rice subspecies, particularly GJ-temperate rice, displayed excellent CT while the xian (XI) rice subspecies was sensitive to low temperatures [5,16,17]. In this research, GWAS was conducted for the CT-related traits using 1992 rice germplasm accessions (including 473 geng and 1290 xian) from 3K Rice Genome Project (3K-RGP). This study clearly aimed to explore the genomic regions and candidate genes related to CT in the whole rice genome and screen breeding materials that could be used for molecular marker-assisted selection (MAS). Our results may facilitate rice breeding for CT to overcome the adverse effects of low temperature on rice yield.

2. Materials and Methods

2.1. Plant Materials

A total of 1992 accessions were chosen for this research from the 3K-RGP [15], including 473 geng, 1290 xian, 61 admix, 133 aus, and 35 basmati (Table S1).

2.2. Phenotypic Evaluation

The test materials were placed in an air oven set at 50°C for 72 hours to break dormancy, then soaked for 48 hours at 37°C, germinated for 24 hours at 37°C, and finally seeded in a seedling tray (60 cm length × 30 cm width × 5 cm height) filled with nutritious soil, with 21 accessions per seedling tray. Two replicates were set for the experiment, and the seeds were cultivated in the phytotron (13 h light at 28℃ and 11 h of darkness at 26℃) until the seedlings were at the three-leaf stage. After removing weak seedlings, cold treated was applied at 9℃ in an artificial climate chamber for 5 days, then moved out to recover in the phytotron for 7 days. Then, the survival rate (SR) was calculated as: the survival rate (%) = No. of surviving plants/total No. of plants × 100. Leaf score of cold tolerance (SCT) was recorded as follows: 1: only the tip of a leaf is yellow or no leaf damage; 3: one-third of the leaves on the third leaf are yellow or withered; 5: roughly two-thirds of the leaves on the second and third leaves are yellow and withered; 7: the second leaf is completely withered, the third leaf is mostly withered but the center is still intact; 9: the plant is completely dead and all of the leaves are wilted.

2.3. Genome-Wide Association Mapping

The SNP genotyping information of 1992 samples was extracted from the 4.8 M SNP dataset of the 3K-RGP [18] using PLINK software [19]. SNPs with the missing rate >20% and minor allele frequencies (MAF) <5% were filtered. The association analysis between the SNPs and CT was performed using a mixed linear model (MLM) in EMMAX [20]. The filtered SNPs were used to calculate the kinship matrix and principal components based on GCTA software, and the first three principal components as covariates to control the population structure [21]. The effective number of SNPs (N) was estimated utilizing GEC software [22]. The significant threshold (1/N) of the suggestive p-value was then determined using the Bonferroni correction method, setting at the whole population (p=1.97E-06), geng (p=5.27E-06) and xian (p=2.72E-06) subpopulations, respectively. The "qqman" package of the R software was used to create the manhattan and Q-Q plots [22]. Significant SNPs within 300-kb were classified as a locus based on the linkage disequilibrium (LD) decay distance of 3K-RGP [15], and the lead SNP in a locus was defined as the SNP with the lowest value of p, and the other significant SNPs within 150 kb on either side of the lead SNP were merged.

2.4. Haplotype Analysis for Candidate Genes

The following two steps were conducted to identify candidate genes: (1) We would match the genes with extremely significant as-sociations SNP (P-value＜1/N) in each characteristic (the downstream genes are selected if the SNP is located between the genes) as important candidate genes; (2) The funRiceGenes database [23] and the Nipponbare reference genome IRGSP 1.0 [24] were used to locate genes with abiotic stress-related functional annotations. The haplotype analysis was performed on all candidate genes using all SNPs within the gene coding sequence region [25]. Duncan’s multiple range post-hoc tests were used to compare phenotypic differences between haplotypes (n≥30 rice accessions).

3. Results

3.1. Phenotypic Variations in the Cold Tolerance (CT) of Rice at the Seedling Stage

Two traits including SR and SCT were measured for 1992 accessions, including 473 geng and 1290 xian accessions at the seedling stage (Table S1). We found geng subpopulation were significantly higher cold-tolerant than xian (Figure 1a,c). We observed that 51 accessions had SR of 100%, including five xian accessions and 46 geng accessions. And 613 accessions had SR of 0, including 557 xian accessions, 30 geng accessions, 7 admix accessions, 10 aus accessions and 9 bas accessions (Table S1). In addition, we found some accessions exhibited extreme phenotypes in the SCT trait, and 136 accessions had no any leaf damage in the CT test, including 11 xian accessions, 117 geng accessions, 4 admix accessions, 2 aus accessions and 2 bas accessions (Table S1), indicating that rice germplasms had wide differentiation in CT among different subpopulations.

There was obvious differentiation in CT even in the same subspecies. The mean SR in the geng subpopulation was 64.58%, ranging from 0 to 100% (Figure S1c). Similarly, its mean SCT was 4.56, with a range of 1.0 to 9.0 (Figure S1d). Among the four geng subgroups, GJ-tmp was significantly higher than the others, whereas most accessions with cold-sensitivity belonged to GJ-sbtrp (Figure 1b,d). The mean SR in the xian subpopulation was 15.54%, ranging from 0 to 100%, and the mean SCT was 7.67, with a range of 1.0 to 9.0 (Figure S1e,f). Among the five xian subgroups, XI-1A was significantly more resistant to cold than the others, indicating that subgroup XI-1A had stronger CT than other xian subgroups.

According to the correlation analysis among the CT traits in the whole population, SR was significantly negatively correlated with SCT (r = -0.84). In addition, correlations of CT-related traits in the subpopulations xian and geng were similar to that in the whole population.

3.2. GWAS for CT

We conducted GWAS based on the mixed linear model for the traits related to CT (Figure 1e,f; Figure 2). Using a Bonferroni correction based on the effective numbers, the genome-wide significant values of p thresholds were set at 1.97E-06, 5.27E-06 and 2.72E-06 for the whole population, geng subpopulation, and xian subpopulation, respectively. Total of 16, 2 and 12 QTLs significantly associated with CT at the seedling stage were identified in the whole population, geng subpopulation, and xian subpopulation, respectively (Table 1 and Table 2).

In the whole population, a total of 16 QTLs for two traits (SR and SCT) were identified, including 8 for SR and 8 for SCT (Table 1). Among them, qSR1.2a for SR and qSCT1.1a for SCT were mapped together in the region of 11.71-11.97 Mb on chromosome 1. The QTL qSR3.1a and qSCT3.1a were simultaneously identified in the region of 25.11-25.41 Mb on chromosome 3. The QTLs qSR4.1a and qSCT4.1a were identified together in the region of 3.58-3.78 Mb on chromosome 4.

A total of two QTLs were detected in the geng subpopulation on chromosomes 2 and 12, including qSR2.1g and qSCT12.1g (Table 2). The qSR2.1g was detected in the region of 19.64-19.94 Mb on chromosome 2 affecting the SR. Similarly, qSCT12.1g was detected in the regions of 11.01-11.31 Mb on chromosome 12 affecting the SCT.

For xian subpopulation, a total of 12 QTLs affecting CT-related traits were mapped on chromosomes 1, 4, 5, 6, 8, 11 and 12, including 6 QTLs for SR and 6 for SCT (Table 2). Among them, qSR1.1x and qSCT1.1x, qSR4.1x and qSCT4.1x, qSR11.1x and qSCT11.1x, and qSR12.1x and qSCT12.1x were commonly detected in the regions of 3.45-3.75 Mb on chromosome 1, 3.63-4.23 Mb on chromosome 4, 24.03-24.33 Mb on chromosome 11, and 5.71-7.25 Mb on chromosome 12, respectively.

Among the 30 QTLs detected for the two traits through GWAS in the three populations, four QTLs (qSR4.1a, qSCT4.1a, qSR4.1x, and qSCT4.1x) were co-identified in both the whole and xian populations, while two QTLs (qSR12.1a and qSR12.1x) were found exclusively in these two populations as well. However, no co-detected QTLs were observed between the whole and geng populations, nor were any QTLs detected simultaneously within the geng and xian subpopulations. These findings suggest a significant allele differentiation on CT between the xian and geng subspecies.

QTLs simultaneously detected by whole and xian could be used as important QTLs for excavating candidate genes associated with rice CT. Six important QTLs (qSR1.2a, qSR3.1a, qSR4.1a, qSR1.1x, qSR11.1x and qSR12.1a) were found in different traits by the whole and xian population. Among them, qSR1.2a and qSR1.1x were coincided with the previously cloned genes OsPLDα1 and OsLEA9 for CT in rice, respectively [26,27]. So, qSR3.1a, qSR4.1a, qSR11.1x and qSR12.1a were considered as newly identified QTLs for rice CT in this study, thus being used for further candidate gene analysis.

3.3. Haplotype Analyses of the Candidate Genes

Two candidate genes, LOC_Os04g06900 and LOC_Os04g07260, were identified using haplotype analysis and gene function analysis for qSR4.1a/qSCT4.1a/qSR4.1x/qSCT4.1x. Haplotype analysis indicated that LOC_Os04g06900 have three major haplotypes, Hap3 was determined as the favorable haplotype with significantly highest SR (Figure 3b, 3c). This haplotype was enriched in geng subpopulation (92.68%) accessions, with an average SR of 63% (Figure 3c, d). Hap1 and Hap2 as cold-sensitive haplotypes were significantly enriched for xian subgroup with SR values of 24% and 9%, respectively (Figure 3d). Furthermore, the proportion of whole accessions with LOC_Os04g06900^Hap3 was no difference between the landrace and modern varieties (Figure 3e). Candidate gene LOC_Os04g07260 was identified four major haplotypes (Figure 3f). Amont them, LOC_Os04g07260^Hap3 with significantly highest SR (65%) was considered as the favorable haplotype, which was significantly enriched in geng subpopulation (86.67%) (Figure 3g, 3h). LOC_Os04g07260^Hap2 and LOC_Os04g07260^Hap4 as cold-sensitive haplotypes were enriched in xian subpopulation with SR values of 23% and 11%, respectively (Figure 3h). The fraction of whole accessions containing LOC_Os04g07260^Hap3 increase from 4% in the landrace to 19% in the modern variety (Figure 5f).

Candidate gene LOC_Os12g10710 was identified at qSR12.1a/qSR12.1x/qSCT12.1x (Figure 4a). LOC_Os12g10710 was recognized two major haplotypes (Figure 4b). The favorable haplotype was identified to be Hap2 with considerably highest SR (59%) and significantly enriched in geng subpopulation accessions (Figure 4c, 4d). Hap1 as cold-sensitive haplotypes was all significantly enriched for subgroup xian with the SR values of 9% (Figure 4c, 4d). Furthermore, the fraction of the whole accessions containing LOC_Os12g10710^Hap2 increase from 40% in the landrace to 45% in the modern variety (Figure 4e).

For locus qSR3.1a/qSCT3.1a on chromosome 3, the candidate gene LOC_Os03g44760 was predicted by haplotype analysis and gene function analysis. LOC_Os03g44760 had a lead SNP rs3_25,249,852 (p=1.36E-07) for SR in the whole population (Figure 4f-h). Three major haplotypes was identified for LOC_Os03g44760 (Figure 4h), and Hap3 with significantly highest SR (49%) was determined as the favorable haplotype, which was significantly enriched in geng subpopulation (Figure 4i). Furthermore, the proportion of the whole accessions with LOC_Os03g44760^Hap3 drop from 8% in landrace to 6% in modern variety (Figure 4j).

The candidate gene LOC_Os11g40610 at qSR11.1x/qSCT11.1x on chromosome 11 (Figure 5a) was shown to be linked to stress tolerance, indicating that it is a relevant candidate gene (Figure 5b) [28]. The favorable haplotype among the three major haplotypes was found to be Hap3, which had a considerably highest SR (74%) (Figure 5c). For the whole population, the haplotype LOC_Os11g40610^Hap3 was highly enriched in the geng subpopulation accessions (Figure 5d). The LOC_Os11g40610^Hap3 increased from 3% in the landrace to 12% in the modern variety (Figure 5e).

3.4. Optimal Combination of CT-Haplotypes

Since SR is a key trait for CT, it was used to identify the favorable haplotype of each candidate gene. Based on haplotype analysis mentioned above, LOC_Os04g06900^Hap3, LOC_Os04g07260^Hap3, LOC_Os12g10710^Hap2, LOC_Os03g44760^Hap3 and LOC_Os11g40610^Hap3 were identified as favorable CT-haplotypes. Pyramiding effects of different haplotypes on CT was analyzed. There remained five groups comprising four candidate genes after removal of rare haplotype combinations (n≥10 accessions). The groups I, II and III was enriched in the geng subpopulation, whereas groups IV and V was mainly found in the xian subpopulation, which showed clear xian-geng differentiation on CT (Figure 6A). Compared to the other four groups, Group I exhibited the highest SR (80%) among all, as it contained the cold-tolerant haplotypes at LOC_Os04g07260, LOC_Os12g10710, LOC_Os03g44760 and LOC_Os11g40610 (Figure 6B). Group II also showed a relatively high SR (77%), possessing the cold-tolerant haplotypes at LOC_Os04g07260, LOC_Os12g10710, and LOC_Os11g40610 (Figure 6B). Both groups I and II had significantly higher SR than those groups III and IV which carried cold-tolerant haplotypes at two genes (LOC_Os03g44760 and LOC_Os12g10710) and one gene (LOC_Os12g10710), respectively (Figure 6B). These findings suggest that pyramiding more favorable haplotypes of these candidate genes enhances CT. After screening for favorable cold tolerant haplotypes of LOC_Os04g07260^Hap3, LOC_Os12g10710^Hap2, LOC_Os03g44760^Hap3 and LOC_Os11g40610^Hap3 in 3K-RGP materials, accessions with four favorable haplotypes (CX534) and three favorable haplotypes (B236 and IRIS_313-8565) were obtained. Furthermore, these results indicate that it is an effective strategy to improve rice CT at the seedling stage by pyramiding multiple favorable haplotypes/alleles.

4. Discussion

Rice is particularly vulnerable to low temperature stress, and the impact of low temperature on rice yield is sometimes severe in temperate rice areas [29]. Breeding for CT rice has always emphasized the evaluation of existing germplasm materials for their ability to withstand low temperatures, aiming to utilize those with strong CT for hybridization and ultimately develop new varieties that are tolerant to low temperatures. While this approach allows for more focused breeding efforts, it also limits the scope of available materials for CT breeding, resulting in a potential issue of limited genetic diversity in subsequent breeding work. The present study extensively evaluated the performance of CT at the seedling stage for 1992 accessions selected from the 3K-RGP, thereby providing a valuable dataset for identifying favorable alleles associated with these cloned CT genes. The identification of numerous accessions exhibiting extreme CT levels will significantly enhance genetic diversity and contribute to future rice breeding efforts targeting CT traits.

The difference of CT between the subspecies may be caused by the temperature difference of latitude and elevation and ecological habitats in different regions [30]. In this study, we observed differential responses to cold stress among various rice subgroups at the seedling stage. The geng subspecies exhibited higher CT compared to the xian subspecies; however, within the geng subgroup, there were variations in CT (Figure 1b, d). Notably, GJ-tmp from China and Europe displayed the highest level of CT with an average survival rate of 78.69% (Figure 1b, d). GJ-adm represented an intermediate type within the geng subspecies with a stronger CT (average survival rate: 69.64%). Moreover, GJ-trp from Southeast Asian island countries demonstrated better CT (58.68%) than GJ-sbtrp from Southeast Asian countries (44.97%) (Figure 1b, d). Similarly, in the xian subspecies, XI-A (44.97% for SR) from East Asia exhibited greater CT compared to other xian subspecies, namely XI-1B, XI-2, XI-3, and XI-adm (ranging from 10.10% to 16.29%). Notably, even within CT-sensitive xian subspecies, certain extreme CT accessions were identified such as Ha Goo, Bu Zhi Ming, Ai Da and Laozaogu of XI-1A with each displaying a SR of 100% from China, Taichun Sen Yu 214 of XI-1B with a SR of 100% from Taiwan, NCS 458 and Derawa of XI-2 with SRs of 97.62% and 96.46%, respectively, from India and Nepal, Gaset Bow of XI-3 with a SR of 94.72% from Nepal, and Padi Ladang Ase Polo Komek and ARC 10581 of XI-adm with SRs of 97.87% and 97.98%, respectively, from Indonesia and India. Therefore, to enhance the genetic diversity in rice breeding programs, it is strongly recommended to incorporate diverse CT varieties from xian subspecies as donor parents, alongside the utilization of traditional CT-landrace and temperate geng variety. Furthermore, previous research studies [31,40] have demonstrated the occurrence of transgressive performance in abiotic stress tolerance, including CT, surpassing that of the parental lines in a majority of BC populations for almost all abiotic stresses. It is noteworthy that identifying BC progeny with exceptional tolerances was a common observation [31,40]. Therefore, in rice CT breeding practices, it is feasible to incorporate more diverse varieties as donor parents irrespective of their individual performance when combined with stringent selection criteria. This approach will likely result in the broadest genetic variations in CT.

CT is a quantitative trait controlled by complex genetic networks in rice. GWAS is a technique for examining how phenotypic variation and genome-wide genetic variation are generally associated. As a result, GWAS analysis is currently a widely used tool for identifying candidate genes related to complex traits. Total of 16 loci associated with CT at the seedling stage were identified using the whole populations. By comparing the previously reported cloned genes for CT with the mapping results in this study, qSR1.2a and qSCT1.1a were located on chromosome 1 at 11.71-11.97 Mb and co-localize with the cold-tolerant gene OsLEA9 [27]. Overexpression of OsLEA9 significantly decreased the CT of rice during reproductive growth, and the CT of OsLEA9 knockout lines was significantly stronger than that of control line [27].The qSR1.1x, and qSCT1.1x were located on chromosome 1 at 3.45-3.75 Mb and co-localize with the cold-tolerant gene OsPLDα1, which played an important role in cold signal transduction in rice by producing phosphatidic acid (PA) and regulated the expression of OsDREB1 via OsMPK6, OsSIZ1, and other PA-binding proteins [26]. The qSCT2.1a was found in the 28.66-28.96 Mb of chromosome 2 and co-located with the transcription factor LGS1, which influenced grain size and could improve rice seedling CT and survival following cold stress treatment [35]. We detected some QTLs that had been previously located, indicating the accuracy of the mapping results by GWAS in CT-related traits in this study.

Bioinformatics was used to further infer candidate genes. We used MBKBASE's RNA-seq database to assess the expression patterns of candidate genes for each QTL to find potential candidate genes (Figure S2 and S3). Among the candidate gene, OsELF4a (LOC_Os11g40610) for qSR11.1x, combined with OsELF3-1 and OsLUX to constitute a terpolymer inhibitor complex OsEC1, which played a role in heading and stress tolerance [36]. Under short day, the heading time of mutant oself4a was delayed for 8 days, and all internodes were shortened with semi-dwarf, short ear, grain length, grain width and grain weight all decreased [28]. The survival rate and setting rate of oself4a decreased under salt stress [28]. So, LOC_Os11g40610 (Figure S2) was inferred as a most likely candidate gene affecting CT in rice. Another candidate gene of qSR3.1a, LOC_Os03g44760, which encoded a protein with a coiled-helix domain in an intermediate region, was expressed at the seedling, heading, grain milk, and endosperm filling stages, particularly in the ovules at the heading stage (Figure S3), and its expression level increased when exposed to abiotic stress (http://ipf.sustech.edu.cn/pub/ricerna/). The LOC_Os03g44760 gene was homologous to SWI1 in Arabidopsis and AM1 in maize, which was mostly found in prophase I of sex blast cells, and was essential for the development of the proper chromosomal shape during meiosis commencement [37]. However, the meiosis process was very vulnerable to environmental influences, particularly high and low temperature stress, resulting in pollen abortion and decreased seed setting rate [38,39]. LOC_Os03g44760 was a most likely candidate gene and worthy of further verification by gene-editing or transgenic approach. To identify CT of above candidate genes, we need to create transgenic materials to verify gene function.

To date, several CT genes have been cloned in rice, such as COLD1, bZIP73, and qPSR10, which enhance CT in cultivated rice derived from wild relatives [32,33,34]. Additionally, CTB4a and Ctb1 have been retained during the adaptation to cold climate conditions in temperate cultivars [2,6], while HAN1 has undergone mutations in temperate cultivars for improved cold adaptation [5]. However, it appears that these cloned CT genes have not yet found widespread utilization in rice breeding programs focused on enhancing CT or are already present in modern cultivars. Therefore, the identification of novel CT genes from germplasm resources remains a crucial task. Based on haplotype analysis of the important candidate genes identified in this study, the pyramiding of LOC_Os03g44760^Hap3 at qSR3.1a/qSCT3.1a, LOC_Os04g07260^Hap3 at qSR4.1a/qSCT4.1a/ qSR4.1x/qSCT4.1x, LOC_Os11g40610^Hap3 at qSR11.1x/qSCT11.1x, and LOC_Os12g10710^Hap2 at qSR12.1a/qSR12.1x/qSCT12.1x could significantly enhance CT at the rice seedling stage. Three accessions CX534, B236 and IRIS_313-8565 were identified, each having SR of 100% and SCT of 1. CX534 from China, carries four favorable haplotypes (LOC_Os03g44760^Hap3, LOC_Os04g07260^Hap3, LOC_Os11g40610^Hap and LOC_Os12g10710^Hap2) with a middle-late ripening cultivar with a flat blade leaf, oval grain type, red glume coloration along with glume tip. B236 from China carries three favorable haplotypes (LOC_Os04g07260^Hap3, LOC_Os11g40610^Hap, and LOC_Os12g10710^Hap2) with high alkali tolerance, fertilizer use efficiency, semidwarf and high seed setting rate; it also exhibits low temperature tolerance at the booting stage. IRIS_313-8565 from Thailand, harbors two favorable haplotypes (LOC_Os03g44760^Hap3 and LOC_Os12g10710^Hap2) with semidwarf, medium leaf size, high grain quality and lodging resistance. Our findings demonstrate that rice germplasms carrying more favorable haplotypes for candidate genes exhibit improved CT (Figure 6), suggesting that developing an optimal combination of haplotypes through pyramiding multiple favorable alleles can be an effective strategy to enhance CT at the seedling stage. Therefore, above three accessions (CX534, B236 and IRIS_313-8565) can serve as donor parents in rice breeding of CT through the introgression of diverse favorable alleles into elite varieties that are susceptible to cold stress by MAS.

5. Conclusions

Phenotypic analysis showed that geng accessions were more cold-tolerant than xian accessions and GJ-tmp subgroup had the strongest CT at the seedling stage. A total of 30 QTLs were identified associated with CT at seedling stage on the basis of two traits (SR and SCT) by GWAS analyses. Five candidate genes (LOC_Os03g44760, LOC_Os04g06900, LOC_Os04g07260, LOC_Os11g40610 and LOC_Os12g10710) were identified through gene function annotation and haplotype analysis. Our findings indicate that rice germplasms pyramiding more favorable haplotypes for candidate genes exhibit improved CT. Then, three accessions (CX534, B236 and IRIS_313-8565) carried different superior alleles were iidentified to be used for molecular breeding of CT in rice through marker-assisted selection (MAS). The result of this study provided resources for improving rice CT and laid the groundwork for future cloning of new CT genes.