Genome-Wide Association Study Reveals the Genetic Basis of Chilling Tolerance in Rice at the Reproductive Stage

A genome-wide association study (GWAS) was used to investigate the genetic basis of chilling tolerance in a collection of 117 rice accessions, including 26 Korean landraces and 29 weedy rices, at the reproductive stage. To assess chilling tolerance at the early young microspore stage, plants were treated at 12 °C for 5 days, and tolerance was evaluated using seed set fertility. GWAS, together with principal component analysis and kinship matrix analysis, revealed five quantitative trait loci (QTLs) associated with chilling tolerance on chromosomes 3, 6, and 7. The percentage of phenotypic variation explained by the QTLs was 11–19%. The genomic region underlying the QTL on chromosome 3 overlapped with a previously reported QTL associated with spikelet fertility. Subsequent bioinformatic and haplotype analyses suggested three candidate chilling-tolerance genes within the QTL linkage disequilibrium block: Os03g0305700, encoding a protein similar to peptide chain release factor 2; Os06g0495700, encoding a beta tubulin, autoregulation binding-site-domain-containing protein; and Os07g0137800, encoding a protein kinase, core-domain-containing protein. Further analysis of the detected QTLs and the candidate chilling-tolerance genes will facilitate strategies for developing chilling-tolerant rice cultivars in breeding programs.


Introduction
Rice is a major crop that provides the staple diet of more than half of the global population, not only in Asia, but also in Australia, the USA, and Europe [1]. Rice is cultivated under various climates, and the effect of weather on the different developmental stages of rice growth has critical impacts on crop yields. The growth, development, and yield of rice are affected by unfavorable environmental conditions, such as drought, extreme temperatures, high salinity, and flooding [2]. Due to its origin in tropical and subtropical regions, rice is particularly sensitive to low nonfreezing temperature stress. In general, the optimum temperature range for rice cultivation is 25-30 • C [3], and damage from low nonfreezing temperature stress can impact the entire cultivation season, from germination to the grain-filling stage. In general, rice germination and seedling growth are adversely affected at temperatures lower than approximately 15 • C, whereas temperatures below approximately 17-20 • C induce chilling stress during the booting, flowering, and grain-filling stages [3].
Chilling stress in vegetative stages, including germination, causes leaf chlorosis, reduction in growth rate and tillering, and low seedling vigor, all of which lead to decreased yields [4,5]. Physiological processes are disrupted under chilling stress as a result of declining photosynthetic efficiency. Declines in chlorophyll content and fluorescence, alongside stomatal closure, contribute to the decline in photosynthetic rates [6,7]. During Illumina HiSeq 2500 Sequencing Systems Platform (Illumina, Inc., San Diego, CA, USA). Average genome coverage was 8×, and filtered reads were aligned to the rice reference genome sequence (IRGSP 1.0). Genomic sequences were filtered as follows, using VCFtools software v0.1.16 [26]: minor allele frequency (MAF) > 5%; max-missing, 0.75; minQ, 30; and minDP, 4. After filtering, 2,369,723 of a possible 6,243,699 SNP sites were retained.

Chilling-Tolerance Evaluation at the Reproductive Stage (Young Microspore Stage)
Individual plants were cultivated in the experimental paddy field at Konkuk University in 15 × 30 cm rows. Plants were cultivated with a conventional rice cultivation method in a natural rice paddy field until the early young microspore stage. Individual plants at the early young microspore stage, with auricle distance of −1 to +1 cm, were transferred to individual 15 cm diameter pots containing paddy field soil. Potted plants were then treated at 12 • C and 70% relative humidity under a 9 h day, with 25,000 Lux/15 h night cycle, for 5 days in a climate-controlled growth chamber (Hanbaek Sci. Bucheon, South Korea). After treatment, plants in a pot were placed in natural climate conditions in the experimental field and cultivated until fully matured. Seeds were collected from mature panicles and used for fertility assessment. Three biological replications were conducted.

Population Structure and Linkage Disequilibrium Decay Analysis
Population structure was analyzed based on SNPs in 117 rice accessions using AD-MIXTURE 1.3.0 software [27]. Cross-validation error was used for estimation of the number of subgroups. Admixture analysis was performed with the ideal number of ancestries (K) of 1-7, and K was estimated using cross-validation error from ADMIXTURE based on tenfold cross-validation. Results were visualized using Pophelper Structure Web APP v1.0.10 (http://pophelper.com, accessed on 8 June 2021) [28]. Plink software [29] was used to generate a principal component analysis (PCA) matrix, and the optimal number (5) of principal components (PCs) was used as a Q-matrix for GWAS correction. The PCA plot was visualized using R v4.0.3. Neighbor-joining (NJ) trees were generated using MEGA X [30], and the result was visualized using Tree of Life (iTOL) v6.1.2 (https://itol.embl.de, accessed on 15 April 2021) [31]. PopLDdecay software [32] was used to calculate linkage disequilibrium (LD). Pairwise LD was calculated for all SNPs and averaged across the whole genome. Chromosomal distance was estimated by LD block, which was half the average of the pairwise correlation coefficient (R 2 ).

GWAS
SNPs were filtered using VCFtools [33], with MAF < 5%; max-missing, 0.75; minQ, 30; and minDP, 4, and indels were removed. The fixed and random model Circulating Probability Unification (FarmCPU) [27], implemented using the rMVP FarmCPU package, was employed for GWAS. The kinship matrix was estimated using the mvp.data.kin function. The kinship matrix (K-matrix) was computed to allow FarmCPU to be used with Q-matrix and K-matrix. GWAS results and plots were automatically visualized in rMVP. The Benjamini-Hochberg FDR method was used to control for the multiple testing error rate [34]. Q-values of 0.10 were used for suggestive thresholds, calculated using the following Equation (1): Finally, the threshold was set as −log10(P) = 6.2873 for the identification of associated QTLs. SNP markers located at QTL peaks were designated as lead SNPs. LD decay analysis identified a 276 kb region on either side of the lead SNP as the candidate genome region for gene identification.

Candidate Gene Prediction and Haplotype Analysis
LD decay analysis was used to designate a region containing 276 kb either side of the lead SNP as the candidate region. Annotation of genes within candidate regions were Plants 2021, 10, 1722 4 of 12 derived from the Rice Annotation Project database (RAP-DB; https://rapdb.dna.affrc.go.jp, accessed on 23 May 2021). After the elimination of missing and heterozygous data, the full SNP marker set was used for haplotype analysis. Haplotypes (in at least five rice accessions) were used for comparative phenotypic analysis. One-way analysis of variance, followed by the least significant difference (LSD) test, was used to compare phenotypic differences between haplotypes.

Variability in Rice Chilling Tolerance during Meiosis in Anthesis Development
A total of 117 Korean rice accessions (KRICE-CORE) were used, comprising 26 Korean landraces, 29 weedy rice accessions collected in Korea, 50 rice germplasms from outside Korea, and 12 Korean rice varieties. In a previous study, these 117 accessions were categorized into five subgroups using 169 simple sequence repeat (SSR) marker genotypes: tropical japonica, temperate japonica, indica, aus, and admixed [35]. All of the 117 accessions showed normal fertility above 90% under the normal cultivated condition in the natural paddy field in Korea. All 26 Korean landraces and 29 of the weedy rice accessions are native to the Korean peninsula. Although the Korean peninsula is located in the temperate climate region, previous genotype analysis assigned Korean landraces and weedy rice accessions to four ecotype groups on the basis of their genetic relationships: aus, indica, tropical japonica, and temperature japonica [36,37]. Chilling tolerance at the early young microspore stage was evaluated with spikelet fertility. The distribution of chilling tolerance was skewed, with 57 accessions having <10% spikelet fertility ( Figure 1A). The average spikelet fertility of the 117 accessions was 26.82%. RWG-023 (Hanyangjo; Korean landrace) was the most chilling-tolerant accession, with 55% spikelet fertility. RWG-023 was grouped into the Aus ecotype, alongside RWG-122 (ChungdoHwayang12; Korean weedy rice), which had 45% spikelet fertility. RWG-021 (Saducho; Korean landrace) exhibited 41% spikelet fertility and was grouped into the indica ecotype. No statistical differences were detected in chilling tolerance between the five ecotype groups ( Figure 1B).
analysis identified a 276 kb region on either side of the lead SNP as the candidate genome region for gene identification.

Candidate Gene Prediction and Haplotype Analysis
LD decay analysis was used to designate a region containing 276 kb either side of the lead SNP as the candidate region. Annotation of genes within candidate regions were derived from the Rice Annotation Project database (RAP-DB; https://rapdb.dna.affrc.go.jp, accessed on 23 May 2021). After the elimination of missing and heterozygous data, the full SNP marker set was used for haplotype analysis. Haplotypes (in at least five rice accessions) were used for comparative phenotypic analysis. One-way analysis of variance, followed by the least significant difference (LSD) test, was used to compare phenotypic differences between haplotypes.

Variability in Rice Chilling Tolerance during Meiosis in Anthesis Development
A total of 117 Korean rice accessions (KRICE-CORE) were used, comprising 26 Korean landraces, 29 weedy rice accessions collected in Korea, 50 rice germplasms from outside Korea, and 12 Korean rice varieties. In a previous study, these 117 accessions were categorized into five subgroups using 169 simple sequence repeat (SSR) marker genotypes: tropical japonica, temperate japonica, indica, aus, and admixed [35]. All of the 117 accessions showed normal fertility above 90% under the normal cultivated condition in the natural paddy field in Korea. All 26 Korean landraces and 29 of the weedy rice accessions are native to the Korean peninsula. Although the Korean peninsula is located in the temperate climate region, previous genotype analysis assigned Korean landraces and weedy rice accessions to four ecotype groups on the basis of their genetic relationships: aus, indica, tropical japonica, and temperature japonica [36,37]. Chilling tolerance at the early young microspore stage was evaluated with spikelet fertility. The distribution of chilling tolerance was skewed, with 57 accessions having <10% spikelet fertility ( Figure 1A). The average spikelet fertility of the 117 accessions was 26.82%. RWG-023 (Hanyangjo; Korean landrace) was the most chilling-tolerant accession, with 55% spikelet fertility. RWG-023 was grouped into the Aus ecotype, alongside RWG-122 (ChungdoHwayang12; Korean weedy rice), which had 45% spikelet fertility. RWG-021 (Saducho; Korean landrace) exhibited 41% spikelet fertility and was grouped into the indica ecotype. No statistical differences were detected in chilling tolerance between the five ecotype groups ( Figure 1B analysis identified a 276 kb region on either side of the lead SNP as the candidate genome region for gene identification.

Candidate Gene Prediction and Haplotype Analysis
LD decay analysis was used to designate a region containing 276 kb either side of the lead SNP as the candidate region. Annotation of genes within candidate regions were derived from the Rice Annotation Project database (RAP-DB; https://rapdb.dna.affrc.go.jp, accessed on 23 May 2021). After the elimination of missing and heterozygous data, the full SNP marker set was used for haplotype analysis. Haplotypes (in at least five rice accessions) were used for comparative phenotypic analysis. One-way analysis of variance, followed by the least significant difference (LSD) test, was used to compare phenotypic differences between haplotypes.

Variability in Rice Chilling Tolerance during Meiosis in Anthesis Development
A total of 117 Korean rice accessions (KRICE-CORE) were used, comprising 26 Korean landraces, 29 weedy rice accessions collected in Korea, 50 rice germplasms from outside Korea, and 12 Korean rice varieties. In a previous study, these 117 accessions were categorized into five subgroups using 169 simple sequence repeat (SSR) marker genotypes: tropical japonica, temperate japonica, indica, aus, and admixed [35]. All of the 117 accessions showed normal fertility above 90% under the normal cultivated condition in the natural paddy field in Korea. All 26 Korean landraces and 29 of the weedy rice accessions are native to the Korean peninsula. Although the Korean peninsula is located in the temperate climate region, previous genotype analysis assigned Korean landraces and weedy rice accessions to four ecotype groups on the basis of their genetic relationships: aus, indica, tropical japonica, and temperature japonica [36,37]. Chilling tolerance at the early young microspore stage was evaluated with spikelet fertility. The distribution of chilling tolerance was skewed, with 57 accessions having <10% spikelet fertility ( Figure 1A). The average spikelet fertility of the 117 accessions was 26.82%. RWG-023 (Hanyangjo; Korean landrace) was the most chilling-tolerant accession, with 55% spikelet fertility. RWG-023 was grouped into the Aus ecotype, alongside RWG-122 (ChungdoHwayang12; Korean weedy rice), which had 45% spikelet fertility. RWG-021 (Saducho; Korean landrace) exhibited 41% spikelet fertility and was grouped into the indica ecotype. No statistical differences were detected in chilling tolerance between the five ecotype groups ( Figure 1B).

Population Structure and LD Decay
The genetic structure of 117 rice accessions was assessed using ADMIXTURE analyses ( Figure 2D). Cross-validation error values from ADMIXTURE analysis showed similar patterns with K values from 4-7. Although the CV error of K = 7 had the lowest value, clustering with K = 7 appeared noisy compared with clustering at K = 5 or K = 6 ( Figure 2A). Although previous research suggested five groups [37], in this study, K = 6 provided clearer clustering than K = 5. At K = 6, the accessions were divided into six groups that were mostly distinguished by their subspecies ( Figure 2D).
tions. Temperate japonica and tropical japonica were also distinguished from one another, but were closely related. Aus was distinguished from both indica and japonica, but was closely related to indica. Decay of LD along physical distance was computed for the full panel of rice accessions. The r 2 value declined with increasing physical distance. The threshold value for candidate regions was determined as half of the maximal r 2 value (0.276), producing a candidate genomic region of 276 kb.

GWAS of Chilling-Stress Tolerance during Anthesis Development
Manhattan plots of SNP markers significantly associated with chilling tolerance at the seedling stage are presented in Figure 3. SNPs with a −log10(P) score of 6.2873 were considered to be significantly associated with chilling tolerance. Nine significant SNPs associated with chilling tolerance were detected. Significant SNPs within a 552 kb range were considered as one association locus and, in total, five QTLs were mapped, on chromosomes 3, 6, and 7 ( Table 1). The percentage of phenotypic variation explained by QTLs was 11-19%, where qCTR7-1 showed the largest explained percentage of 19%. LD block Temperate japonica accessions were divided between cluster 2 and cluster 5. Tropical japonica was dominant in cluster 4, aus was dominant in cluster 1, and indica accessions were divided between cluster 3 and cluster 6. PCA was performed using 962,055 SNPs. The five previously described subgroups were well separated by the first two PCs ( Figure 2B), Plants 2021, 10, 1722 6 of 12 with PC1 and PC2 explaining 39.52% and 20.31% of the total variation in population structure, respectively.
Accessions grouped as temperate japonica, tropical japonica, indica, and aus formed significantly distinct clusters, whereas the admixed accessions exhibited ambiguous separation. In addition, an NJ tree constructed based on Nei's genetic distance revealed five clusters ( Figure 2C), consistent with the PCA separation among groups. Most accessions were clearly separated, with admixed accessions dispersed among the different clusters. PCA and NJ tree PCs ( Figure 2B,C) clearly differentiated the indica and japonica subpopulations. Temperate japonica and tropical japonica were also distinguished from one another, but were closely related. Aus was distinguished from both indica and japonica, but was closely related to indica. Decay of LD along physical distance was computed for the full panel of rice accessions. The r 2 value declined with increasing physical distance. The threshold value for candidate regions was determined as half of the maximal r 2 value (0.276), producing a candidate genomic region of 276 kb.

GWAS of Chilling-Stress Tolerance during Anthesis Development
Manhattan plots of SNP markers significantly associated with chilling tolerance at the seedling stage are presented in Figure 3. SNPs with a −log10(P) score of 6.2873 were considered to be significantly associated with chilling tolerance. Nine significant SNPs associated with chilling tolerance were detected. Significant SNPs within a 552 kb range were considered as one association locus and, in total, five QTLs were mapped, on chromosomes 3, 6, and 7 ( Table 1). The percentage of phenotypic variation explained by QTLs was 11-19%, where qCTR7-1 showed the largest explained percentage of 19%. LD block size was 552 kb surrounding the lead SNP, and the genomic segments corresponding to these five QTLs were defined and then compared with those of previously reported QTLs (Table  1). Overlaps were found between the loci detected by GWAS in this study and previously reported QTLs linked to various agronomic traits.  However, none of the previously reported QTLs were associated with chilling tolerance in the reproductive stage. QTLs associated with stress, fertility, or spikelet morphology were identified, and these traits may be indirectly associated with chilling tolerance during the reproductive stage. Comparative analysis revealed that qCTR3-1 on chromosome 3 overlapped with the AQFT003 QTL, which was associated with osmotic adjustment capacity [38]. On chromosome 3, qCTR3-2 overlapped with the AQCU106 QTL, which was associated with spikelet fertility [39]. On chromosome 6, qCTR6 overlapped with the CQAR19 QTL, which was associated with spikelet density [40]. On chromosome 7, qCTR7-1 overlapped with the CQAV4 QTL, which was also associated with osmotic adjustment capacity [41]. Moreover, on chromosome 7, qCTR7-2 overlapped with the AQEM004 QTL, which was associated with salt sensitivity [42].

Haplotype Analysis to Identify Chilling Tolerance Candidate Genes
In silico analysis was conducted to identify candidate genes responsible for chilling tolerance. The RAP-DB database (IRGSP 1.0) was used to identify annotated genes in the 552 kb regions encompassing the five lead SNPs. A total of 83 genes were found in the encompassing region of qCTR3-1, 81 genes for qCTR3-2, 41 genes for qCTR6, 70 genes for qCTR7-1, and 53 genes for qCTR7-2. Hypothetical protein genes were not used for subsequent haplotype analysis and, in total, 213 of the 328 identified genes in the five QTL regions were subjected to haplotype analysis. Phenotypic comparison was conducted for haplotypes containing at least five rice accessions. Finally, 26 candidate genes showing statistically significant differences between the haplotypes were detected (Supplementary  However, none of the previously reported QTLs were associated with chilling tolerance in the reproductive stage. QTLs associated with stress, fertility, or spikelet morphology were identified, and these traits may be indirectly associated with chilling tolerance during the reproductive stage. Comparative analysis revealed that qCTR3-1 on chromosome 3 overlapped with the AQFT003 QTL, which was associated with osmotic adjustment capacity [38]. On chromosome 3, qCTR3-2 overlapped with the AQCU106 QTL, which was associated with spikelet fertility [39]. On chromosome 6, qCTR6 overlapped with the CQAR19 QTL, which was associated with spikelet density [40]. On chromosome 7, qCTR7-1 overlapped with the CQAV4 QTL, which was also associated with osmotic adjustment capacity [41]. Moreover, on chromosome 7, qCTR7-2 overlapped with the AQEM004 QTL, which was associated with salt sensitivity [42].

Haplotype Analysis to Identify Chilling Tolerance Candidate Genes
In silico analysis was conducted to identify candidate genes responsible for chilling tolerance. The RAP-DB database (IRGSP 1.0) was used to identify annotated genes in the 552 kb regions encompassing the five lead SNPs. A total of 83 genes were found in the encompassing region of qCTR3-1, 81 genes for qCTR3-2, 41 genes for qCTR6, 70 genes for qCTR7-1, and 53 genes for qCTR7-2. Hypothetical protein genes were not used for subsequent haplotype analysis and, in total, 213 of the 328 identified genes in the five QTL regions were subjected to haplotype analysis. Phenotypic comparison was conducted for haplotypes containing at least five rice accessions. Finally, 26 candidate genes showing statistically significant differences between the haplotypes were detected (Supplementary  Table S2). Based on the phenotypic differences between haplotypes and the functional annotation of genes, three candidate genes were selected for further analysis. Haplotype analyses of the three candidate genes are presented in Figures 4-6. Heterozygous SNPs and SNPs with missing data were not included in the analysis. SNPs in exons were used for haplotype and haplotype variation analysis. Os03g0305700 (similar to peptide chain release factor 2) contained two SNPs within exons ( Figure 4A). The two SNPs were nonsynonymous and located in the coding region (T→A, Chr3_10850864, S→T substitution; G→C, Chr3_10852950, D→E substitution) of the dienelactone hydrolase family domain. The two SNPs formed three haplotypes, Hap1, Hap2, and Hap3 ( Figure 4B). The mean chilling tolerance of Hap2 (score, 17.6891; 28 accessions) was higher than that of Hap1 (score, 5.03673; 11 accessions). Os06g0495700 (beta tubulin, autoregulation binding-site-domaincontaining protein) contained four nonsynonymous SNPs in exons (C→T, Chr6_17275985, Y→H substitution; A→G, Chr6_17276385, R→H substitution, A→G, Chr6_17277077, A→T substitution; A→T, Chr6_17279115, N→Y substitution), forming four haplotypes (Hap1, Hap2, Hap3, and Hap4) (Figure 4). The mean chilling-tolerance score of Hap1 (23.4644, eight accessions) was higher than that of Hap3 (10.4898, three accessions) and that of Hap4 (11.7632, 70 accessions). Os07g0137800 (protein kinase, core-domain-containing protein) contained a single nonsynonymous SNP (C→G, Chr7_1993025, D→E substitution) that formed two haplotypes (Hap1 and Hap2) ( Figure 5B). The fertility of Hap1 differed significantly from that of Hap2. Chr6_17277077, A→T substitution; A→T, Chr6_17279115, N→Y substitution), forming four haplotypes (Hap1, Hap2, Hap3, and Hap4) (Figure 4). The mean chilling-tolerance score of Hap1 (23.4644, eight accessions) was higher than that of Hap3 (10.4898, three accessions) and that of Hap4 (11.7632, 70 accessions). Os07g0137800 (protein kinase, coredomain-containing protein) contained a single nonsynonymous SNP (C→G, Chr7_1993025, D→E substitution) that formed two haplotypes (Hap1 and Hap2) ( Figure  5B). The fertility of Hap1 differed significantly from that of Hap2.  four haplotypes (Hap1, Hap2, Hap3, and Hap4) (Figure 4). The mean chilling-tolerance score of Hap1 (23.4644, eight accessions) was higher than that of Hap3 (10.4898, three accessions) and that of Hap4 (11.7632, 70 accessions). Os07g0137800 (protein kinase, coredomain-containing protein) contained a single nonsynonymous SNP (C→G, Chr7_1993025, D→E substitution) that formed two haplotypes (Hap1 and Hap2) ( Figure  5B). The fertility of Hap1 differed significantly from that of Hap2.

Discussion
Chilling tolerance of rice during the reproductive stage is crucial for stable yield production. Chilling temperature exposure during this highly sensitive stage prevents devel-

Discussion
Chilling tolerance of rice during the reproductive stage is crucial for stable yield production. Chilling temperature exposure during this highly sensitive stage prevents development of normal pollen, resulting in significant yield loss. Chilling-tolerance distribution in the 117 rice accessions was skewed, with 47.86% of the accessions exhibiting <10% spikelet fertility. This reflected the high sensitivity of the early young microspore stage to chilling exposure. The temperate japonica rice ecotype is generally thought to be more chilling tolerant than the aus and indica ecotypes. However, no differences in chilling tolerance between the ecotypes were observed here. Moreover, accession RWG-023, which grouped as aus, was the most chilling-tolerant accession, with 55% spikelet fertility. RWG-122, which also grouped as aus, exhibited 45% spikelet fertility, while RWG-021, exhibiting 41% spikelet fertility, grouped as indica. Although these chilling-tolerant accessions were grouped as aus or indica, generally regarded as chilling-susceptible ecotypes, these lines belonged to the native germplasm of the Korean peninsula, which has a temperate climate. RWG-021 and RWG-023 are Korean landraces, and RWG-122 is a Korean weedy rice. The unexpected reproductive chilling tolerance of these accessions suggests that the chillingtolerance trait may not be associated with ecotype differentiation in the Korean peninsula. Although the dissociation between reproductive chilling tolerance and ecotype was not clearly elucidated in our study, these germplasms may be beneficial for introducing the chilling-tolerance trait into aus and indica varieties as a genetically closely related donor parent.
Previously reported QTLs overlapping the physical positions of the newly identified QTLs were investigated. Three stress-associated QTLs, one spikelet-fertility-associated QTL, and one panicle-associated QTL were found. The qCTR3-1 QTL overlapped the AQFT003 QTL, which was associated with osmotic adjustment capacity. AQFT003 was detected by QTL analysis of drought tolerance at the reproductive stage in DHLs derived from a cross between upland japonica (CT9993) and indica cultivars (IR62266) [38]. The qCTR3-2 QTL overlapped the AQCU106 QTL, which was associated with spikelet fertility. AQCU106 was identified in an RIL population derived from a cross between Lemont (japonica) and Teqing (indica) varieties [39]. The qCTR6 QTL overlapped the CQAR19 QTL, which was associated with spikelet density. CQAR19 was detected by QTL analysis of an F2 mapping population derived from indica Zhaiyeqing 8 and japonica Jingxi 17 [40]. The qCTR7-1 QTL overlapped the CQAV4 QTL, which was associated with osmotic adjustment capacity. CQAV4 was detected by QTL analysis of drought tolerance at the seedling stage, with RILs derived from a cross between an upland japonica cultivar (Moroberekan) and an indica cultivar (Co39) [41]. The qCTR7-2 QTL overlapped the AQEM004 QTL, which was associated with salt sensitivity. AQEM004 was detected by QTL analysis of F 2 and an equivalent F 3 population developed from a cross between a high-salt-tolerant indica variety, Nona Bokra, and a salt-susceptible elite japonica variety, Koshihikari, [42].
The rice genome database [43]. was used to identify 328 candidate genes in the LD blocks of the detected QTLs. After haplotype analysis and functional annotation of the genes, three candidate genes were selected for further consideration: Os03g0305700, Os06g0495700, and Os07g0138400. Expression profiling data for the three candidate genes were obtained from the Transcriptome Encyclopedia of Rice (TENOR) database (http://tenor.dna.affrc.go.jp/, accessed on 20 June 2021) [44]. Os03g0305700 is predicted to encode a peptide chain release factor 2 family protein. According to TENOR, expression of Os03g0305700 increased in seedlings under chilling conditions (Supplementary Figure S1). A release factor is a protein that allows for the termination of translation by recognizing the termination codon or stop codon in an mRNA sequence. Although it is widely believed that translation termination is a highly conserved process in eukaryotes, the role of translation termination in plant development is largely unknown. Eukaryotic release factor 1-2 is involved in hypersensitive responses to glucose and phytohormones during germination and early seedling development [45]. Os06g0495700 is predicted to encode a beta tubulin. According to TENOR, expression of Os06g0495700 increased under drought and flood conditions (Supplementary Figure S1). Beta tubulin microtubules have important roles in many cellular processes, including cell division and cell elongation in plants. The organization and stability of plant microtubules are affected by environmental stresses, such as dehydration, high salinity, low nonfreezing temperature, and aluminum exposure [46]. Post-meiotic radial microtubule arrays in Arabidopsis male gametes undergo depolymerization in response to chilling stress [47]. In Arabidopsis, chilling stress alters the formation of radial microtubule arrays at telophase II and consequently leads to defects in meiosis during pollen development [48]. Expression of beta tubulin also decreased in Arabidopsis leaves upon exposure to low nonfreezing temperatures [48]. Os07g0137800 is predicted to encode a protein kinase. According to TENOR, Os07g0137800 expression increased under high salinity, low phosphate, and high phosphate conditions (supplementary Figure S1). In addition, Ling et al. (2015) used qPCR analysis to confirm that Os07g0137800 was preferentially expressed in the mature anther [49]. Protein kinases regulate the biological activity of proteins by phosphorylation of specific amino acids, with ATP as the source of phosphate. Protein kinases are involved in various plant responses to environmental stresses such as drought, high salinity, chilling, and pathogen attack [50]. Increased expression of OsCDPK7, a calcium-dependent protein kinase, was detected upon exposure to chilling and salt stresses, and overexpression of OsCDPK7 in transgenic rice increased tolerance to chilling and salt stress during the vegetative stage, suggesting that OsCDPK7 acted as a positive regulator during the tolerance response to both stresses in rice [51]. Based on the function of protein kinase, the expression pattern of Os07g0137800, and the previously reported role of calcium-dependent protein kinase in chilling tolerance, Os07g0137800 is a strong candidate for further evaluation for its possible role in chilling tolerance during pollen development.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/plants10081722/s1, Figure S1: Expression profiles in rice seedling under the various environmental conditions. Table S1: A core Korean rice collection, 117 accessions;