Quantitative Trait Loci Mapping Analysis for Cold Tolerance Under Cold Stress and Brassinosteroid-Combined Cold Treatment at Germination and Bud Burst Stages in Rice

Low temperature is one of the major abiotic stresses limiting seed germination and early seedling growth in rice. Brassinosteroid (BR) application can improve cold tolerance in rice. However, the regulatory relationship between cold tolerance and BR in rice remains undefined. Here, we constructed a population of 140 backcross recombinant inbred lines (BRILs) derived from a cross between a wild rice (Dongxiang wild rice, DXWR) and a super rice (SN265). The low-temperature germination rate (LTG), survival rate (SR), plant height (PH), and first leaf length (FLL) were used as indices for assessing cold tolerance under cold stress and BR-combined cold treatment at seed germination and bud burst stages. A high-resolution SNP genetic map, covering 1,145 bin markers with a distance of 3188.33 cM onto 12 chromosomes, was constructed using the GBS technique. A total of 73 QTLs were detected, of which 49 QTLs were identified under cold stress and 24 QTLs under BR-combined cold treatment. Among these, intervals of 30 QTLs were pairwise coincident under cold stress and BR-combined cold treatment, as well as different traits including SR and FLL, and PH and FLL, respectively. A total of 14 candidate genes related to cold tolerance or the BR signaling pathway, such as CBF/DREB (LOC_Os08g43200), bHLH (LOC_Os07g08440 and LOC_Os07g08440), WRKY (LOC_Os06g30860), MYB (LOC_Os01g62410 and LOC_Os05g51160), and BRI1-associated receptor kinase 1 precursor (LOC_Os06g16300), were located. Among these, the transcript levels of 10 candidate genes were identified under cold stress and BR-combined cold treatment by qRT-PCR. These findings provided an important basis for further mining the genes related to cold tolerance or the BR signaling pathway and understanding the molecular mechanisms of cold tolerance in rice.


INTRODUCTION
Rice (Oryza sativa L.) is one of the most important crops for food production, which grows in tropical, subtropical, and temperate regions worldwide (Cheng et al., 2007). Cold stress is one of the most severe environmental factors limiting the growth, development, and yield formation of rice. It affects almost all growth stages of rice, including germination stage, bud burst stage, seedling stage, tillering stage, booting stage, flowering stage, and grain filling stage (Fujino et al., 2004;Xu et al., 2015). The germination rate and post-germination early seedling growth are two important traits that directly contribute to seedling vigor (Fujino et al., 2008). Cold stress at the germination and bud burst stages is the major limitation due to high sensitivity to cold at these stages, especially for the direct-seeded rice. It impairs the seed germination and early seedling growth of rice, which in turn leads to uneven stand establishment of seedling, delay of panicle development, spikelet sterility, and subsequent yield loss (Satoh et al., 2016;Najeeb et al., 2020). Therefore, improving the cold tolerance of germination and bud burst stages is an important objective in rice cultivation and breeding.
Brassinosteroids (BR) are essential plant steroid hormones that play important roles in growth, development, and tolerance to stresses in plants. Many studies have reported that exogenous BR application can obviously improve cold tolerance in plants (Chung et al., 2014;Sharma et al., 2017;Manghwar et al., 2022). In Arabidopsis thaliana, loss-of-function mutations of brassinosteroid insensitive 1 (BRI1) and brassinazole resistant 1 (BZR1) and overexpression of brassinosteroid insensitive 2 (BIN2) result in decreased cold tolerance, whereas overexpression of BRI1 and BZR1 enhances cold tolerance Manghwar et al., 2022). In addition, BR biosynthetic genes BR6ox2, DWF4, and CPD are rapidly downregulated under cold stress (Chung et al., 2014;Eremina et al., 2016). In rice, soaking with BR can effectively alleviate the damage of rice seeds under cold stress, and spraying with BR can also improve cold tolerance of rice seedlings (Hwang and Back, 2019;Wu et al., 2020;Sadura and Janeczko, 2021). However, few genes related to both cold tolerance and BR in rice have been reported to date. The regulatory relationship between cold tolerance and BR in rice remains undefined.
Cold tolerance is a very complex quantitative trait associated with a massive number of biochemical and physiological processes and environment factors, which is genetically controlled by multiple quantitative trait loci (QTLs) or genes (Zhang et al., 2014). QTL mapping is one of the most common methods to identify the QTLs and genes related to cold tolerance in plants. Many genomic regions on all 12 rice chromosomes have been reported to contain some QTLs related to cold tolerance at different development stages (Miura et al., 2001;Fujino et al., 2004Fujino et al., , 2008Yang et al., 2021). For instance, Ctb1 is the first cloned QTL for cold tolerance in rice. Ctb1 encodes a F-box protein, which improves cold tolerance by associating with a subunit of the E3 ubiquitin ligase (Skp1) at the booting stage of rice (Saito et al., 2010). CTB4a, a leucine-rich repeat receptor-like kinase, increases ATP synthase activity and ATP content by interacting with the β-subunit of ATP synthase, which enhances seed setting and improves yield under lowtemperature stress conditions . In addition, COLD1 functions as a GTPase-accelerating factor and regulates G-protein signaling under cold stress in rice. COLD1 interacts with regulator of GTPase-activating protein 1 (RGA1), rapidly activates inward current and Ca 2+ concentration, and ultimately enhances cold tolerance at the seedling stage in japonica rice (Ma et al., 2015). However, these are not sufficient to mine the useful genes and reveal the molecular mechanisms of cold tolerance in rice. With the rapid development of highthroughput sequencing technologies, faster and simpler ways to obtain single-nucleotide polymorphism (SNP) molecular markers based on sequencing technologies are increasingly used to perform QTL mapping, such as genotyping-by-sequencing (GBS), bulked segregant analysis (BSA), and bulked segregant RNA-seq (BSR). The introduction of these technologies has accelerated the identification of millions of SNPs across the whole genome (Jiang et al., 2016;Mielczarek and Szyda, 2016). In particular, GBS is becoming popular for QTL mapping, genetic diversity, and genomic selection, which has been successfully applied in QTL mapping for cold tolerance, salinity tolerance, aluminum tolerance, rice blast resistance, pericarp color, and some agronomic trait in rice (Spindel et al., 2013;Arbelaez et al., 2015;De Leon et al., 2016).
Wild rice (Oryza rufipogon Griff) is the relative ancestor of the cultivated rice, which can be used as a donor of novel and favorable alleles for rice breeding (Nakagahra et al., 1997;Atwell et al., 2014;Zhang et al., 2020). QTL mapping has been extensively used to identify the novel loci from wild rice, and these loci are expected to be used for improving the agronomic traits of rice (Koseki et al., 2010;Mao et al., 2015). Dongxiang wild rice (DXWR) possesses an extremely high innate tolerance to low-temperature stress at the seedling, booting, and flowering stages, and its underground stem can tolerate −12.8 • C (Li et al., 2010;Zhang et al., 2014). Some genes or QTLs related to cold tolerance have been isolated and characterized in DXWR (Liu et al., 2003;Tian et al., 2006;Li et al., 2010;Xiao et al., 2015;Deng et al., 2018;Liang et al., 2018Liang et al., , 2019Bai et al., 2021). Therefore, DXWR is an ideal germplasm for improving the cold tolerance of rice by hybridization, backcrossing, or genetic transformation.
In this study, we constructed a backcross recombinant inbred line (BRIL) population of 140 individuals derived from a cross between the DXWR and a super rice SN265 with excellent agronomic traits. On this basis, we used the GBS technique to construct a high-resolution genome-wide SNP genetic map for identification of cold-tolerant QTLs and candidate genes under cold stress and BR-combined cold treatment at the germination and bud burst stages. This study has important reference value for the discovery of cold tolerance genes and understanding of the regulation relationship between cold tolerance and BR in rice.

Plant Materials
The recipient parent SN265, a super rice variety with excellent agronomic traits from north China, was crossed with the donor parent DXWR, and the F 1 plants were backcrossed with SN265 four times to develop the BC 4 F 1 generation. The resulting BC 4 F 1 lines were selfed and advanced by using the single seed descent method to generate 140 BRIL individuals in the F 8 generation (BC 4 F 8 ).

Phenotypic Characterization for Cold Tolerance
The evaluation of cold tolerance at germination and bud burst stages was conducted based on previous studies with minor changes (Kim et al., 1999;Qiao et al., 2004;Akhtamov et al., 2020;Najeeb et al., 2020;Pan et al., 2020). Seeds of each accession from the BRILs were placed in a drying oven at 50 • C for 72 h to break dormancy, and surface-sterilized in 1% NaClO solution for 10 min, followed by washing three times in sterilized distilled water. Then, 50 seeds were soaked with sterilized distilled water (for cold treatment) and 0.1 µmol/L BR solution (for BRcombined cold treatment) in a petri dish (15 cm) at a constant temperature of 20 • C for 24 h, respectively. The 50 pre-soaked seeds were stressed at 10 • C for 10 days, and then moved to a greenhouse at 25 • C for recovery. The germinated seeds of each accession were counted daily and defined as the low-temperature germination rate (LTG). The germination rate was calculated as follows: Germination rate (%) = (number of germinated seeds/total number of seeds) × 100%. Cold tolerance based on the LTG was graded on a scale of 1-9 as follows: 1, more than 90%; 3, 80-90%; 5, 50-79%; 7, 1-49%; and 9, 0%.
The survival rate (SR) is the percentage of rice seedlings survived from the germination stage to the early seedling stage under cold stress. The pre-soaked seeds were germinated in a growth chamber at a constant temperature of 30 • C for 24 h. The germinated seeds with 5-mm coleoptiles were stressed at 10 • C for 7 days and then moved to a greenhouse at 25 • C for 10 days to allow early seedlings to resume normal growth. The SR was calculated as follows: Survival rate (%) = (number of survival seedlings/number of buds) × 100%. At least five pots of 30 plants were assessed for each line. Cold tolerance based on the SR was graded on a scale of 1-9 as follows: 1, more than 90%; 3, 70-90%; 5, 50-70%; 7, 10-49%; and 9, less than 10%.
Cold stress can directly influence the growth of coleoptile length, leaf length, and plant height at the bud burst stage. Therefore, phenotypic evaluation of cold tolerance at the bud burst stage is significant (Qiao et al., 2004;Najeeb et al., 2020;Pan et al., 2020). In the current study, cold tolerance was evaluated based on the phenotypic changes of plant height (PH) and first leaf length (FLL) at the bud burst stage. The PH and FLL of the early seedlings after recovering for 10 days were measured. The average of the three replicates of 10 seedlings for each treatment was analyzed. For cold treatment, the data obtained under normal temperature (25 • C) were used as control. Cold tolerance scores based on the reduction rate of the PH and FLL were calculated as follows: All experiments for evaluating cold tolerance at the bud burst stage were repeated three times under the same conditions, and the average cold tolerance scores from the three replicates were used for QTL mapping analysis. R software 3.5.1 was used to analyze frequency distributions and correlations of four phenotypic traits under their respective two treatments. 1 The "cor" function was used to calculate Pearson's correlation coefficient between any two treatments. The "hist" function computes a histogram of the average phenotypic scores of each individual under each treatment. The significance between any pair of treatment was classified as 0, 0.001, 0.01, 0.05, 0.1, and 1, which were denoted as " * * , " " * * , " " * , " "., " and "".

Genotyping and Single-Nucleotide Polymorphism Identification
Fresh young leaves were collected from 14-day-old seedlings of 140 BRILs along with their parental lines and subjected to DNA extraction using a DNA extraction kit (Aidlab, China). GBS technology was used to construct DNA fragment libraries from each BRILs and their parental lines. PstI and MspI (NEB) were used for digestion and T4 ligase (NEB) for ligation. The libraries were enriched by PCR amplification and sequenced using an Illumina HiSeq 4000 instrument (Huada Gene Technology Co. Ltd., Shenzhen, China). The data were analyzed by Tassel software (Glaubitz et al., 2014). The raw reads were then filtered and sorted according to indices, and the high-quality SNPs between parents were termed by alignment with Nipponbare reference genome.

Linkage Map Construction and Quantitative Trait Loci Mapping
The genotypic data from each BRILs with filtered SNP markers were used to construct linkage map using QTL IciMapping sofware v4.163. Furthermore, QTLs were mapped using the inclusive composite interval mapping of additive (ICIM-ADD) method. QTLs were computed by using a permutation test involving 1,000 runs at a significance level of p = 0.05. After completion of the permutation test, a window size of 10 cM and a walk speed of 1 cM were set to start analysis of composite interval mapping. The threshold for the logarithm of odds (LOD) scores was set to 3.0. QTL nomenclature was followed by the method of McCouch (2008).

Annotation and Validation of Candidate Genes
The QTL intervals were determined by the physical position of the SNP markers, flanking the respective QTLs. All candidate genes were predicted in the located QTL ranges based on the Rice Genome Annotation Project website. 2 The 3-day-old early seedlings with seeds from SN 265, DXWR, and three excellent cold-tolerant BRILs, and three non-excellent coldtolerant BRILs at normal temperature, 12-h under cold stress, and 12-h BR-combined cold treatment (mentioned before for SR) were subjected to total RNA extraction using an RNA extraction kit (Aidlab, China). Three biological replicates were collected, immediately frozen in liquid nitrogen, and then stored at −80 • C until further analysis. The real-time PCR (qRT-PCR) was performed to identify the expression patterns of 10 candidate genes. The qRT-PCR was carried out using a q225 Real-Time PCR System (Monad, China) under the following conditions: 95 • C for 5 min; 30 cycles of 95 • C for 10 s, 60 • C for 30 s, and 72 • C for 2 min; and 1 cycle of 72 • C for 10 min. Each sample was analyzed in three biological and three technical replicates. The relative expression levels were calculated using the Ct method (Livak and Schmittgen, 2002). The actin gene of rice (LOC_Os03g50885) was used as a reference gene. Ct was the difference between the Cts (threshold cycles) of the target gene and Cts of the reference gene.
Ct was calculated by subtracting Ct of the treatment group from Ct of the normal group. Fold change was calculated using the following formula: Fold change = 2 − Ct . All primer sequences are listed in Supplementary Table 1.

Phenotypic Characterization Under Cold Stress and Brassinosteroid-Combined Cold Treatment
The BRIL population was evaluated under cold stress and BRcombined cold treatment at the germination and bud burst stages. The BRILs showed varying levels of cold tolerance and significantly contrasting responses in LTG, SR, PH, and FLL (Figures 1, 2). For LTG, grade 5 had the highest number of BRILs under cold stress, followed by grades 3, 7, and 1. Under BR-combined cold treatment, the number of BRILs in grades 3 and 5 increased, whereas the number of BRILs in grade 7 decreased ( Figure 1B). For SR, grade 3 had the highest number of BRILs under cold stress, followed by grades 1, 5, and 7. After BR-combined cold treatment, the number of BRILs in grades 1 and 5 increased, whereas the number of BRILs in grade 7 decreased ( Figure 1C). Furthermore, the phenotypic changes of PH and FLL were analyzed under cold stress and BR-combined cold treatment (Figure 2). Under normal temperature, the range of PH was 5−8 cm, and most of BRILs were distributed in the range of 7−8 cm, while the range of FLL was 2−4 cm, and most of BRILs were distributed in the range of 3−4 cm. After cold stress and BR-combined cold treatment, the PH and FLL were significantly reduced in most BRILs, with the range of 4−7 cm for PH and 2−4 cm for FLL, respectively. Compared to cold stress, the number of BRILs with PH of a 6-cm range and FLL of both 2-cm and 3-cm ranges increased under BR-combined cold treatment (Figures 2B,C). These results indicated that the Correlation Analysis of Low-Temperature Germination Rate, Survival Rate, Plant Height, and First Leaf Length Pearson's correlation analysis revealed that all four traits in the BRIL populations showed continuous and approximately normal distributions. Under cold stress, highly significant positive correlations were observed between LTG and SR (r = 0.37, p ≤ 0.01), and between PH and FLL (r = 0.39, p ≤ 0.01), respectively. By contrast, there was a significant negative correlation between SR and PH (r = −0.21, p ≤ 0.05). Under BR-combined cold treatment, the positive correlations between LTG and SR (r = 0.37, p ≤ 0.01), as well as between PH and FLL (r = 0.39, p ≤ 0.01), and a negative correlation between LTG and SR (r = 0.37, p ≤ 0.01) were shown similar to those under cold stress. Additionally, under both cold stress and BRcombined cold treatment, all of LTG, SR, PH, and FLL showed highly significant positive correlations with r = 0.73 (p ≤ 0.01), r = 0.73 (p ≤ 0.01), r = 0.79 (p ≤ 0.01), and r = 0.73 (p ≤ 0.01), respectively (Figure 3).

Genetic Linkage Map of the Recombination Bin Markers
In total, 64.48 Gb of high-quality sequence data were obtained from 42.8 M raw reads via the GBS approach, and 96.19% of those reads were mapped to the Nipponbare reference genome. The ratio of Q20 for BRILs was above 90%, and the guaninecytosine (GC) content was 43.55%; thus, the quality of the data met requirements for further analysis. Finally, a total of 10,836 unfiltered SNPs were validated for the determination of recombinant events. After filtration followed by the ABH plugin, a total of 1,145 recombination bin markers were obtained to construct a recombination map for the BRIL population.
A genetic linkage map was developed by mapping the 1,145 bin marker to the 12 rice chromosomes. The 12 rice linkage groups varied in the number of markers and marker density. The total genetic distance of linkage map was 3188.33 cM, with the individual linkage groups ranging from 186.74 to 346.46 cM in length. The highest number of markers was 117 on chromosome 4, and the lowest number was 77 on chromosome 7. The averaged genetic distance between markers was 0.3 cM, and the interval of genetic distance among markers ranged from 1.9 to 146.3 cM ( Table 1).

Comprehensive Quantitative Trait Loci Mapping for Cold Tolerance in Backcross Recombinant Inbred Line Population
The 1,145 SNP markers were used in QTL mapping of cold tolerance based on LTG, SR, PH, and FLL, respectively. An LOD threshold of 3.0 was used to identify QTLs for each trait using the CIM analysis. A total of 73 QTLs were detected, of which 49  QTLs were identified under cold stress and 24 QTLs under BRcombined cold treatment. The phenotypic variation explained (PVE) by these QTLs ranged from 0.52 to 14.69%, while the range of LOD threshold was 3.02−9.38 (Figure 4 and Table 2).
A total of four QTLs related to LTG were localized on chromosomes 7, 8, and 10, respectively. Among them, qLTG7-1 with 3.93 LOD and 10.84% PVE was localized on chromosome 7 under cold stress, while qLTG7-2, qLTG8-1, and qLTG10-1, with the LOD ranging from 3.06 to 3.48 and PVE ranging from 4.74 to 12.41%, were localized on chromosomes 7, 8, and 10 under BR-combined cold treatment, respectively. The interval of qLTG7-1 localized under cold stress coincided with the interval of qLTG7-2 localized under BR-combined cold treatment (Table 2; Supplementary Figure 1).
For SR, a total of 43 QTLs were identified on all chromosomes, except chromosome 4, with LOD ranging from 3.05 to 9.38 and PVE ranging from 0.52 to 7.0%, 31 of which were localized   under cold stress and 12 were localized under BR-combined cold treatment, respectively. Intervals between some QTLs under both cold stress and BR-combined cold treatment were coincident, including qSR5-4 and qSR5-5, qSR6-4 and qSR6-3, qSR6-6 and qSR6-5, qSR7-1 and qSR7-2, qSR8-4 and qSR8-3, and qSR11-3 and qSR11-4, respectively (Table 2; Supplementary Figure 2).  Unlike the LTG and SR, the data analysis of PH and FLL for QTL mapping needed to set different treatments as the controls (see Section "Materials and Methods" for details). For PH, only one QTL (qPH1-2) was localized on chromosome 1 under cold stress. Under BR-combined cold treatment, three QTLs (qPH1-1, qPH4-1, and qPH7-1) were identified on chromosomes 1, 4, and 7 with the control of PH under the normal-temperature condition, while one QTL (qPH6-1) was localized on chromosome 6 with the control of PH under cold stress. The LOD values of these QTLs ranged from 3.63 to 8.58, and their PVE values ranged from 8.94 to 14.69% (Table 2; Supplementary Figure 3).
For FLL, a total of 21 QTLs with LOD ranging from 3.02 to 4.68 and PVE ranging from 1.27 to 11.12% were identified on all chromosomes, except chromosomes 3 and 10. Among these 21 QTLs, 16 were localized in cold stress and five were localized in BR-combined cold treatment. Under BR-combined cold treatment, two QTLs (qFLL1-2 and qFLL6-4) were identified on chromosomes 1 and 6 using PH under normal-temperature conditions as control, while three QTLs (qFLL6-2, qFLL8-1, and qFLL8-2) were localized on chromosomes 6 and 8 with PH under cold stress as control (Table 2; Supplementary Figure 4).

Identification of Candidate Genes
According to the Rice Genome Annotation Project, the coincident QTL regions and QTLs with a higher LOD score (more than 5.0) were chosen to identify candidate genes using flanking markers. A total of 121 genes were located in the QTL regions mentioned previously, of which 60 were annotated with known functions, while the remaining 61 genes were identified as expressed proteins, hypothetical proteins, transposon, and retrotransposon proteins (Supplementary Table 2). Among these annotated candidate genes, 14 genes were possibly associated with cold tolerance or the BR signaling pathway based on the previously reported studies ( Table 3). These genes included some transcription factors such as C-repeat-binding factor/DRE-binding factor (CBF/DREB, LOC_Os08g43200), basic helix-loop-helix (bHLH) transcription factor (LOC_Os07g08440 and LOC_Os05g50900), AP2 domaincontaining proteins (LOC_Os06g44750 and LOC_Os06g09390), WRKY (LOC_Os06g30860), and MYB (LOC_Os01g62410 and LOC_Os05g51160), ethylene-responsive transcription factor (ERF, LOC_Os05g34730), and some other proteins including leucine-rich repeat (LRR) family protein (LOC_Os03g01410), heat shock factor-binding protein (LOC_Os06g16270), zinc finger-containing protein (LOC_Os06g09310), and extra-large G-protein-related (LOC_Os05g50910), and brassinosteroid insensitive 1 (BRI1)-associated receptor kinase 1 precursor (LOC_Os06g16300). Surprisingly, BRI1-associated receptor kinase 1 precursor, a protein associated with the BR signaling pathway, was just identified under BR-combined cold treatment. Furthermore, we validated the expression levels of 10 candidate genes under cold stress and BR-combined cold treatment using qRT-PCR. The results showed that the transcript levels of all genes, except LOC_Os05g50910 and LOC_Os05g34730, were significantly higher under both cold stress and BR-combined cold treatment than under normal temperature in all materials tested, while the transcript levels of LOC_Os06g16300, LOC_Os08g43200, and LOC_Os05g50900 under BR-combined cold treatment were slightly higher than those under cold stress. For two parents, LOC_Os05g50900, LOC_Os06g16300, and LOC_Os03g01410 showed higher expression levels in DXWR than in SN265 under cold stress and BR-combined cold treatment, whereas the opposite expression pattern was observed for LOC_Os06g30860 and LOC_Os06g09310. For the representative BRILs, LOC_Os05g50900, LOC_Os06g16300, LOC_Os08g43200, LOC_Os05g50900, and LOC_Os06g09390 showed higher transcript levels in three excellent BRILs than three non-excellent BRILs both under cold stress and BR-combined cold treatment ( Figure 5; Supplementary Figure 5).

DISCUSSION
As a complex quantitative trait, cold tolerance in plants is controlled by various QTLs or genes. Thus, it is absolutely essential to mine more QTLs or genes to get a deeper understanding of molecular mechanisms for cold tolerance in plants (Zhang et al., 2014;Jha et al., 2017). Genetic diversity is the basis of the agronomic traits and stress resistance improvement in plants, and gene diversity is the most widely used index for measuring the level of genetic diversity. Wild species have high diversity of desirable genes for agronomic traits and resistance to stresses, and these allelic variations are important genetic resources for agronomic traits and stress resistance improvement (Koseki et al., 2010;Atwell et al., 2014). Common wild rice is the relative ancestor of the cultivated rice, which possesses abundant genetic diversity that can be used to improve agronomic traits and resistance to stress in cultivated rice. However, some gene resources might have been lost in cultivated rice after thousands of years of evolution and natural selection (Liu et al., 2003;Koseki et al., 2010;Liang et al., 2018Liang et al., , 2019Zhang et al., 2020).
As previously mentioned, DXWR has a strong cold tolerance at the seedling, booting, and flowering stages, which is an ideal material for identifying and offering genes related to cold tolerance for rice improvement (Zhang et al., 2006). To date, many rice populations have been constructed using DXWR as the paternal parent, while many QTLs and candidate genes related to cold tolerance have been identified in these populations (Mao et al., 2015;Liang et al., 2018Liang et al., , 2019). In the current study, SN265, a super rice variety with excellent agronomic traits, was crossed with the donor parent DXWR to obtain F1 populations. To have a uniform background, the material was propagated by backcrossing the F1 populations with female parent SN265 for additional four times (BC 4 F 1 ). Furthermore, the BC 4 F 1 populations were selfed for multiple generations to generate an advanced BRIL population (BC 4 F 8 ), and 140 BRIL individuals with abundant genetic diversity were chosen for cold tolerance evaluation and GBS sequencing. The backcrossing and selfing for multiple generations eliminated the background interference caused by different genetic bases, making the results of phenotypic evaluation and QTL mapping of cold tolerance more accurate and reliable.
Multiple studies have demonstrated that exogenous application of BR increased cold tolerance in plants. However, only several genes in BR biosynthesis and signaling pathways were reported to be involved in cold tolerance in Arabidopsis such as BRI1, BZR1, BIN2, BR6ox2, DWF4, and CPD (Chung et al., 2014;Eremina et al., 2016;Li et al., 2017;Sharma et al., 2017;Manghwar et al., 2022). In rice, few genes have been reported to be involved in both cold stress and BR biosynthesis or signaling pathways. FIGURE 5 | Expression patterns of the candidate genes. Biological triplicates were averaged and statistically analyzed via Student's t-test (*P < 0.05; **P < 0.01; ***textitP < 0.001). NT, normal temperature; C, cold stress; C + B, BR-combined cold treatment.
In the current study, the exogenous application of BR improved the cold tolerance of some BRILs at germination and bud burst stages under cold stress, which provided a basis for identifying the QTLs related to cold stress and BR pathways.
The data of LTG and SR could be derived directly from the calculation of corresponding proportion under different treatments, and the analysis process does not need the controls (see Section "Materials and Methods" for details). Unlike LTG and SR, however, the measurements of PH and FLL for QTL mapping need to have a control so that the data could embody the changes under different treatment. Thus, the measurements of PH and FLL under cold stress set the data under normal-temperature condition as the control, whereas the measurements of PH and FLL under BR-combined cold treatment set the data under normal-temperature condition and cold stress as the controls, respectively. This was done to obtain more comprehensive QTLs (Kim et al., 1999;Qiao et al., 2004;Akhtamov et al., 2020;Najeeb et al., 2020;Pan et al., 2020).
In QTL mapping, intervals between some QTLs were coincident under different treatments (cold stress and BRcombined cold treatment), as well as different traits including SR and FLL, and PH and FLL, respectively. These coincident QTLs were probably more accurate. Furthermore, many candidate genes were identified by using the coincident QTL regions and those QTLs with higher LOD scores (more than 5.0). Among these genes, there were some transcription factors related to cold tolerance such as CBF/DREB (LOC_Os08g43200), bHLH (LOC_Os07g08440 and LOC_Os07g08440), WRKY (LOC_Os06g30860), MYB (LOC_Os01g62410 and LOC_Os05g51160), and ERF (LOC_Os05g34730). In particular, DREB/CBF, bHLH, and WRKY transcription factors have been reported to be involved in BR signaling pathways. Under cold conditions, BR directs a bHLH transcription factor CESTA (CES) to contribute to the constitutive expression of the CBF regulators that control the expression of cold responsive genes in Arabidopsis (Eremina et al., 2016). BZR1, a key transcription factor in the BR signaling pathway, acts upstream of CBF1 and CBF2 to directly improve cold tolerance in Arabidopsis. Moreover, BZR1 also regulates WKRY6 to modulate plant response to cold stress . In the current study, qRT-PCR was applied to preliminarily characterize the transcript levels of all transcription factor genes under cold stress and BR-combined cold treatment. The transcript levels of these genes, except LOC_Os05g34730, were upregulated both under cold stress and BR-combined cold treatment. For the representative BRILs, LOC_Os05g50900 (bHLH), LOC_Os08g43200 (DREB/CBF), and LOC_Os05g50900 (bHLH) showed higher transcript levels in the excellent BRILs than in the non-excellent BRILs under cold stress and BR-combined cold treatment. In addition to the transcription factors, some genes related to cold tolerance were also identified such as LRR protein, heat shock factor-binding protein, and zinc finger-containing protein. Interestingly, a BRI1-associated receptor kinase 1 precursor (LOC_Os06g16300) was also identified under BR-combined cold stress. BRI1 is a leucine-rich repeat receptor-like kinase that functions as the cell surface receptor for BR, and loss-of-function mutation of BRI1 resulted in decreased cold tolerance in Arabidopsis (Eremina et al., 2016;Xu et al., 2022). Although the aforementioned candidate genes have been reported to be involved in cold tolerance or BR signaling pathways, few genes related to both cold stress and BR pathways have been studied in rice. In the current study, we initially demonstrated that the transcript levels of LOC_Os06g16300 (BRI1-associated receptor kinase 1), LOC_Os08g43200 (DREB/CBF), and LOC_Os05g50900 (bHLH) under BR-combined cold treatment were higher than those under cold stress. Further studies are needed to demonstrate whether these genes are involved in cold tolerance and BR pathways in rice. In conclusion, many QTLs and candidate genes were identified by GBS sequencing and mapping analysis, which provided an important basis for further mining the genes related to cold tolerance or BR pathways and studying the molecular mechanisms regulating cold tolerance in rice.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA836720.

AUTHOR CONTRIBUTIONS
ML and MZ designed and supervised the research work. ZG, MZ, and ML constructed the materials. ZG, HW, and JY determined phenotypic data, mapping analysis, and qRT-PCR experiments. YC and WZ performed the candidate gene association. JH and ZX contributed to revising the manuscript. All authors contributed to the article and approved the submitted version.