Rice (Oryza sativa L.), one of the most important staple crops, feeds more than 2.7 billion people worldwide (Fairhurst and Dobermann 2002). In general, rice is exposed to various environmental conditions during its life cycle. Especially, rice is more sensitive to cold stress than are other crops, such as wheat (Triticum aestivum L.) and barley (Hordeum vulgare L.), because rice originates in tropical and subtropical regions. Low temperatures in reproductive stages critically to affect grain fertility and quality and reduce yield in high-latitude or high-altitude regions of China, Japan, Korea, and other parts of the world (Jena et al. 2012). For example, 750,000 kg of rice production was reduced by cold stress at the seedling stage in the Hangzhou region in 1980 (Xiong 1984). Rice production of South Korea is also affected by cold stress every 10 to 15 years because of its geographical location and 20% of agricultural land faces possible cold stress in every year (Yea et al. 2008).

Cold stress can be classified as chilling (0-15°C) or freezing (< 0°C) stress. Cold stress at the seedling stage causes the delayed of growth and defective root growth. Especially, defective roots caused by cold stress cause suppress early-stage growth and tiller growth (Ham et al. 2021). The typical symptoms of cold stress are discoloration, yellowing, and withering, including reduced chlorophyll contents and fluorescence, which hinders the photosynthesis of rice (Kanneganti and Gupta 2008; Kim et al. 2009). The cold response also produces reactive oxygen species and malondialdehyde, which contributes to damaged cells and impaired cellular metabolism (Nakashima et al. 2007; Xie et al. 2009). On cold-tolerance responses of rice, to maintain the proper functioning of rice cells, proline must be accumulated in a cell to stabilize protein synthesis (Kandpal and Rao 1985). For cultivation method, phosphorus fertilizer and deep irrigation may be applied to prevent cold-stress damages, but doing so is difficult, because weather forecasts can’t perfectly forecast low temperatures and the damage appearance has various patterns and causes (Ham et al. 2021).

In general, the most efficient method for preventing cold-stress damage is a genetic modification to develop tolerant varieties with the aid of genetic studies, many of which have tried to find QTLs and have identified genes related to cold tolerance. A Genome-wide association study (GWAS) is the analysis that finds the linkage between a single nucleotide polymorphism (SNP) marker and phenotype through the association analysis. Recently, a GWAS has been widely used to discover the genetic control of cold tolerance in rice seedlings (Shakiba et al. 2017). Ham et al. (2021), evaluated the cold tolerance in a temperature-controlled growth chamber, which did not ideally represent the natural field conditions or the status of roots. Here, we evaluated cold tolerance by using cold-water irrigation to provide a natural rice-growing environment and to direct the cold treatment to root area in the controlled growth chamber, and natural cold regions. With 136 rice germplasms in the Korean Rice Core Set (Kim et al. 2007; Kim et al. 2016), we conducted a GWAS to survey the QTLs related to cold tolerance in the seedling stage.

Plant materials and genotype

The KRICE_CORE of 136 accessions (Kim et al. 2007; Kim et al. 2016) consisted of 62 Temperate japonica, 19 Tropical japonica, 42 Indica, 8 Aus, 3 Aromatic and 2 Admixture accessions. We implemented sequencing analysis of 136 accessions from the Illumina HiSeq 2500 Sequencing Systems Platform (Illumina Inc., San Diego, CA, USA). In brief, average genome coverage was 8× and we aligned filtered reads were aligned to the rice reference genome (IRGSP 1.0). The filtering parameters for genotypes were as follows; minor allele frequency (MAF) > 5%, missing data < 1%, and heterozygosis ratio < 5%, as implemented using Plink software v1.90 (Kim et al. 2016).

Cold water treatment and cold tolerance evaluation

The cold-water treatment for rice seedling were conducted in 2020 at the Chuncheon Substation of the National Institute of Crop Science, RDA (Chunchon, Korea, 127’2’E, 37 16’N, 36 m a.s.l). Rice seedlings were placed in the facility of Chuncheon’s cold-water storage unit which is 50 m long and 8 m wide. This facility can supply 13 ± 1°C of cold water through the cold-water outlet (70 cm gap, 8 mm diameter). Germinated rice seeds were placed in 72-cells trays (each cell size was 3.8 cm × 3.8 cm × 4.5 cm). After two weeks of growing in the normal greenhouse condition, the trays were placed in the cold-water treatment facility and were irrigated 13 ± 1°C cold water for ten days. After the cold treatment, the seedlings recovered for seven days in the greenhouse. The phenotypes of the recovered seedlings were scored. The score standard is ranged from 1 to 9, based on the is evaluated change in leaf color (Fig. 1).

Figure 1. Representation of cold tolerant score stress at seedling stage. 1: Resistance-9: Highly sensitive.

The score of 1 represents symptoms of no change in leaf color, 5 represents a change to light yellow in half of the leaf, and 9 represents an almost entirely dead leaf.

Genome wide association study

In the first step, we removed SNPs displaying a minimal allele frequency (frequency of the minor allele) lower than 5% were removed. In the GWAS analysis, we used the Mixed Linear Model (MLM) implemented in the GAPIT R package (version 3.0) were used. For MLM with Q-matrix and K-matrix, we computed the kinship (K matrix) and the first three principal components (PCs) of a PCA of genomic data as the Q matrix. The PCA matrix was generated by Plink software (Purcell et al. 2007). We did the GWAS and plot visualization using the GAPIT package (version 3.0) in R (Yang et al. 2014). We calculated the association threshold by the formula: ‒log10 (1/number of independent SNPs). Finally, we set the threshold as ‒log10(P) = 4.6 for identification of association loci. We designated the SNP markers located at locus peaks as lead SNPs for the detected loci. LD decay analysis identified 250 kb around the lead SNPs as candidate genomic regions for gene identification.

Haplotype analysis of candidate genes

We extracted candidate genes ware extracted in the range of ± 250 kb on the lead SNP position (region of QTL) based on RAP-DB (http://rapdb.dna.affrc.go.jp/). We obtained the expression patterns of the candidate genes inside the QTLs from the Rice Expression Profile Database (Rice XPro: https://ricexpro.dna.affrc.go.jp/). We used haplotype analysis for the selected candidate genes, using all SNP markers from the gene, but excluded missing and heterozygote data. Among the groups of haplotypes, we did phenotypic comparisons with a one-way ANOVA followed by LSD through SAS 9.4.

Phenotypic variations in cold tolerance

We used six rice subgroups consisting of 136 rice accessions were used for evaluating cold tolerance under the cold-water irrigating. The six subgroups were Temperate japonica, Tropical japonica, Indica, Aus, Aromatic, and Admixture. The distribution of the cold-tolerantce score skewed to the right, showing the highest frequency in score 9 (Fig. 2A).

Figure 2. The distribution of cold tolerant score for the 136 rice accessions. (A) Histogram of the cold tolerant score for 136 accessions. The dotted line is the moving average. Different colors indicated the cold tolerance score: blue 1, orange 2, gray 3, yellow 4, sky blue 5, green 6, black 7, pale gray 8, and dark gray 9. (B) Cold tolerant score grouped by box plot according to their ecotype.

The average score of the 136 lines was 5.8, where 39 accessions showed a cold-tolerance score of 9. We compared the variation between the subspecies groups (Fig. 2B). The average score of Indica (42 accessions), Temperate japonica (62 accessions), Tropical japonica (19 accessions), Admixture (2 accessions), Aromatic (3 accessions), Aus (8 accessions) was 7.67, 5.8, 5.53, 7.00, 9.00, and 7.75 (Fig. 2B), respectively. The Japonica group showed stronger tolerance than did the Indica and Aus groups.

GWAS of cold-stress tolerance in the rice seedling stage

The Manhattan plots for the markers significantly associated with cold tolerance at the seedling stage are shown in Fig. 3

Figure 3. Manhattan plots and Q-Q plot for the cold tolerant scores of 136 rice accessions. (A) Manhattan plots (B) Q-Q plot. It indicated the ‒log₁₀(P) value on y axis and the SNP position of each chromosome on x axis. The horizontal blue line indicated thresholds (‒log₁₀(P) = 4,60).

. Associations showing higher than the threshold of ‒log10(P) = 4.6 were detected as a significant association. The significantly associated SNP makers located within ± 250 kb on the lead SNP position was designated to one QTL locus of the lead SNP markers. Seven association loci were detected (Fig. 3, Table 1)

Table 1 . The locations of QTL of detected in GWAS and previously reported QTL.

QTL	Lead SNP	Chr	‒log10(P)	PVE(%)	Reported QTL	Reference of previously reported QTLs
					Related trait
qCWS1-1	7,239,225	1	4.63	8.63	Drought tolerance	Wan et al. 2003
qCWS1-2	24,297,420	1	4.78	8.96	Drought tolerance	Li et al. 2005
qCWS3	1,094,354	3	4.6	8.55	Drought tolerance	Hemamalini et al. 2000
qCWS6	24,641,646	6	4.74	8.87	Drought tolerance	Zhang et al. 2001
qCWS7	20,640,777	7	4.51	8.34	Cold tolerance	Hou et al. 2004
qCWS10	18,189,402	10	4.92	9.29	Drought tolerance	Ali et al. 2000
qCWS12	22,021,391	12	5.02	9.51	Blast resistance	Bagali et al. 1998

. Two associated QTLs were detected on chromosome 1, designated as qCWS1-1 and qCWS1-2, explaining 8.63% and 8.96% of the total phenotypic variations, respectively. One associated QTL was detected on chromosome 3, designated as qCWS3, explaining 8.55% of total phenotypic variations. One associated QTL was detected on chromosome 6, designated as qCWS6, explaining 8.87% of total phenotypic variations. One associated QTL was detected on chromosome 7, designated as qCWS7, explaining 8.34% of total phenotypic variations. One associated QTL was detected on chromosome 10, designated as qCWS10, explaining 9.29% of total phenotypic variations. One associated QTL was detected on chromosome 12, designated as qCWS12, explaining 9.51% of total phenotypic variations. For segments of the seven QTLs, we compared previously reported QTLs (Table 1). The location of qCWS7 was overlapped with the previously reported QTL for low temperature germination (Hou et al. 2004). Five QTLs (qCWS1-1, qCWS1-2, qCWS3, qCWS6, and qCWS10) were overlapped with previously reported QTLs for drought tolerance (Ali et al. 2000; Hemamalini et al. 2000; Zhang et al. 2001; Wan et al. 2003; Li et al. 2005).

Haplotype analysis of candidate genes

To identify candidate genes responsible for cold tolerance under the cold-water irrigating, we extracted all annotated genes located within 500 kb of the QTL regions were extracted based on the RAP-DB (IRGSP 1.0). There were 966 genes were in the seven QTL regions where 286 genes were located in chromosome 1, 176 genes in chromosome 3, 132 genes in chromosome 6, 160 genes in chromosome 7, 139 genes in chromosome 10, and 73 genes in chromosome 12. We looked for the expression patterns of the 966 genes in the previously reported RNA-seq database (Kawahara et al. 2016). Considering our cold treating method, we focused on the gene expression in the roots and shoots at the seedling stage. We have no doubt that other genes not showing specific gene expression patterns in roots might be involved in the cold-tolerant response. However, there are many candidate genes in the regions of the seven QTLs; so we focused only on the genes showing specific expression patterns in the roots and shoots. We detected an increase of 70 genes in the seedling roots and subjected them to haplotype analysis. Given those results and the previously known functions of the genes, we arbitrarily selected five genes showing significant differences in the haplotype analysis as candidate genes for cold tolerance under cold-water irrigation. The haplotype analyses of the five candidate genes are presented in Fig. 4-8.

Figure 4. Haplotype analysis result of Os01g0228600. (A) Gene structure and SNPs positions on Os01g0228600. Black box and line indicated exons and introns, respectively. Red marks indicated SNPs. (B) Haplotypes and their average score of cold tolerance.

We used only the SNPs in exon regions for analyzing haplotypes and haplotype variations, where heterozygous SNPs and SNPs with missing data were filtered. The candidate gene of Os01g0228600 encodes a cytosolic hydroxypyruvate reductase. It contained six non-synonymous SNPs (C→T, Chr1_7100465, A→V substitution; G→C, Chr1_7100466, A→P substitution; C→T, Chr1_7100621, A→V sub-stitution; G→A, Chr1_7102701, D→N substitution; T→G, Chr1_7102745, L→R substitution; T→C, Chr1_7102844, L→P substitution) that formed five haplotypes (Fig. 4). The Hap5 of Os01g0357800 was the superior genotype in cold tolerance. The candidate gene of Os03g0115000 encodes a cupredoxin domain containing protein. It contained two non-synonymous SNPs SNPs (G→T, Chr3_851240, A→S substitution; T→C, Chr3_851346, L→P substitution) that formed three haplotypes (Fig. 5). The cold-tolerance score of Hap1 differed from that of Hap2 and Hap3.

Figure 5. Haplotype analysis result of Os03g0115000. (A) Gene structure and SNPs positions on Os03g0115000. Black box and line indicated exons and introns, respectively. Red marks indicated SNPs. (B) Haplotypes and their average score of cold tolerance.

The candidate gene of Os06g0612800 encodes a stress associated protein (SAP) gene family. It contained a non-synonymous SNPs SNPs (C→T, Chr6_851240, M→T substitution) that formed two haplotypes (Fig. 6). The cold tolerance score of Hap1 differed significantly from that of Hap2. The candidate gene of Os10g0482900 encodes a thioredoxin fold domain containing protein. It contained three non-synonymous SNPs (T→G, Chr10_18232294, S→R substitution; A→G, Chr10_18232294, K→E sub-stitution; C→T, Chr10_18232851, A→V substitution) that formed three haplotypes (Fig. 7). The cold tolerance score of Hap1 differed from that of Hap2 and Hap3. The candidate gene of Os12g0552500 encode a universal stress protein domain containing protein. It contained two non-synonymous SNPs (T→G, Chr12_22430675, S→A substitution; G→T, Chr12_22431727, A→S substitution) that formed two haplotypes (Fig. 8). The Hap2 differed significantly from that of Hap1.

Figure 6. Haplotype analysis result of Os06g0612800. (A) Gene structure and SNPs positions on Os06g0612800. Black box and line indicated exons and introns, respectively. Red marks indicated SNPs. (B) Haplotypes and their average score of cold tolerance.

Figure 7. Haplotype analysis result of Os10g0482900. (A) Gene structure and SNPs positions on Os10g0482900. Black box and line indicated exons and introns, respectively. Red marks indicated SNPs. (B) Haplotypes and their average score of cold tolerance.

Figure 8. Haplotype analysis result of Os12g0482900. (A) Gene structure and SNPs positions on Os12g0482900. Black box and line indicated exons and introns, respectively. Red marks indicated SNPs. (B) Haplotypes and their average score of cold tolerance.

In general, rice cold tolerance in seedling rice has been evaluated in the temperature-controlled growth chamber, where the rice shoot is mainly affected by cold air. In this study, we evaluated the phenotypic variations as a cold-tolerantce score for rice seedling under cold- water irrigating, where the irrigated cold water continuously affected mainly the roots and basal shoots. The treated seedlings recovered from the cold treatment by replacing them in the greenhouse for seven days. Therefore, the cold-tolerance score represented both the damage from the cold treatment and the recovering ability. The overall distribution of cold-tolerance scores was skewed to the right (cold sensitive), implyingd that the cold treatment in this study was effective.

Among the detected seven QTLs, five were overlapped with previously reported QTLs for unfavorite environ-mental conditions. Interestingly, the location of qCWS7 was overlapped with the previously reported QTL for low temperature germination(Hou et al. 2004). Even though the QTLs were detected from the RIL population where the cold tolerant parent was the japonica variety USSR5 and the target traits was germinating ability under the low temperature, the overlapped detection of qCWS7 in this study suggested that the possible role of cold tolerance in this area. Despite of the haplotype analysis, we did not select a candidate gene in this region. Since we focused on the genes expressed in roots and shoots thus further investigation with other approaching for selecting candidate genes needs to be applied in these QTL regions. Besides the QTL for cold tolerance, five QTLs (qCWS1-1, qCWS1-2, qCWS3, qCWS6, and qCWS10) overlapped with previously reported QTLs for drought tolerance, suggesting that these QTLs may play a role in the general response to abiotic stresses.

Based on the gene expression patterns, haplotype analysis, and previously known function of the genes, we suggested five candidate genes as for cold tolerance. Os01g0228600 encodes a cytosolic hydroxypyruvate reductase (HPR) which is responsible for the conversion of hydroxypyruvate into glycerate in the photorespiration cycle (Shi et al. 2021). Photorespiration helps maintain normal physiological metabolism in higher plants. In barely, the important involve of photorespiration in drought response and the involve of HPR induction (WINGLER et al. 1999). Os06g0612800 encodes a stress associated protein (SAP) gene family and Os12g0552500 encodes a universal stress protein (USP). From the annotating and the functional study, these two stress related genes were known to be involved in stress response. For the SAP, overexpression in transgenic rice and tobacco conferred tolerance to salt, drought, cold stress (Kanneganti and Gupta 2008). For the USP, 44 USP genes have been detected in rice and analysis of the expression data in multiple stress conditions including -heat, salt, drought, dehydration, and submergence- revealed that Os12g0552500 was upregulated in those stress conditions (Arabia et al. 2021). Os03g0115000 encodes a cupredoxin domain containing protein that has a type I copper centers and is involved in inter-molecular electron transfer reaction (Dennison 2005). In maize, cupredoxin domain-containing protein was reported to be involved in modulation of redox homeostasis in kernel development (Zhang et al. 2021). Os10g0482900 encodes a thioredoxin fold domain containing protein, which is an antioxidants in reactive oxygen species (ROS) scavenging in plant responses against stresses (Park et al. 2019). The increase of Os10g0482900 under the nitrogen starving stress condition was reported (Shin et al. 2018). These two genes of Os03g0115000 and Os10g0482900 are involved in redox homeostasis in plant cells. Various redox-active compounds, such as ascorbate, glutathione, NAD(P)H, and proteins from the thioredoxin superfamily, contribute to the general redox homeostasis in the plant cells (Kapoor et al. 2015). Accumulating evidence has shown that the redox state of the cell affects its proliferation/differentiation program and redox poising, which involves the movement of electrons from multiple sites in respiratory and photosynthetic electron-transport chains to molecular oxygen to triggers the formation of superoxide to produce other ROS (Huang et al. 2019). With the redox homeostasis gene in plant cells, the cells percepts abiotic stress and triggers physiological responses where antioxidants play a significant role as a signaling compounds, because they also possess a vibrant metabolic interface (Huang et al. 2019).

From previously reported studies and the known functions of these five candidate genes, the direct and indirect evidence suggested their possible role in cold tolerance in rice; so it is highly worthwhile to conduct future evaluation of the molecular mechanisms of these genes in stress response in cold tolerance and other abiotic stress as well.

This research was funded by Korea Institute of Planning and Evaluation for Technology in Food, Agriculture and Forestry (IPET) through Agri-Food Export Business Model Development Program, funded by Ministry of Agriculture, Food and Rural Affairs (MAFRA) (320105-3).

This research was funded by Rural Development Administration, grant number PJ01579403.