Identification of QTL under Brassinosteroid-Combined Cold Treatment at Seedling Stage in Rice Using Genotyping-by-Sequencing (GBS)

Cold stress is a major threat to the sustainability of rice yield. Brassinosteroids (BR) application can enhance cold tolerance in rice. However, the regulatory mechanism related to cold tolerance and the BR signaling pathway in rice has not been clarified. In the current study, the seedling shoot length (SSL), seedling root length (SRL), seedling dry weight (SDW), and seedling wet weight (SWW) were used as the indices for identifying cold tolerance under cold stress and BR-combined cold treatment in a backcross recombinant inbred lines (BRIL) population. According to the phenotypic characterization for cold tolerance and a high-resolution SNP genetic map obtained from the GBS technique, a total of 114 QTLs were identified, of which 27 QTLs were detected under cold stress and 87 QTLs under BR-combined cold treatment. Among them, the intervals of many QTLs were coincident under different treatments, as well as different traits. A total of 13 candidate genes associated with cold tolerance or BR pathway, such as BRASSINAZOLE RESISTANT1 (OsBZR1), OsWRKY77, AP2 domain-containing protein, zinc finger proteins, basic helix-loop-helix (bHLH) protein, and auxin-induced protein, were predicted. Among these, the expression levels of 10 candidate genes were identified under different treatments in the parents and representative BRIL individuals. These results were helpful in understanding the regulation relationship between cold tolerance and BR pathway in rice.


Introduction
Rice (Oryza sativa L.) is one of the most important food crops that grows in tropical to temperate regions worldwide [1]. Low temperature is a severe environmental restriction that most strongly impacts rice growth and development, especially at the early seedling and reproductive stages [2]. Cold stress at the early vegetative stage can result in stunted growth and increased seedling mortality, which in turn, leads to uneven seedling stand establishment, delayed panicle development, spikelet sterility, and subsequently, decreased yield of rice [3,4]. Therefore, it is a crucial object in rice cultivation and breeding to improve cold tolerance at seedling stage.
Cold tolerance of rice is a complex quantitative trait controlled by multiple quantitative trait loci (QTL), which is directly related to a large amount of physiological and biochemical processes and environmental factors [3,4]. To date, numerous backcross inbred lines, recombinant inbred lines, and near-isogenic lines have been developed using many cold tolerance varieties as donors, and a massive number of cold tolerance-related QTLs

Phenotypic Characterization under Cold Stress and BR-Combined Cold Treatment
Cold tolerance of each BRIL individual was evaluated under cold stress and BRcombined cold treatment at the seedling stage. The BRIL population showed varying levels of cold tolerance and significant contrasting response in SSL, SRL, SDW, and SWW ( Figure 1). SSL was in the range of 20-40 cm and the majority of BRIL individuals were distributed in the range of 30-35 cm under normal temperature. Meanwhile, SRL was in the range of 3-9 cm and the majority of BRIL individuals were concentrated in the range of 4-7 cm. After cold stress and BR-combined cold treatment, the SSL of the majority of BRIL individuals was significantly reduced. SSL was in the range of 10-35 cm and the majority of BRIL individuals were distributed in the range of 20-30 cm. Compared to cold stress, the number of BRIL individuals with SSL in the range of 25-30 cm significantly increased under the BR-combined cold treatment (Figure 1b). After cold stress and BR-combined cold treatment, the SRL of the majority of BRIL individuals also reduced; however, it was less pronounced than the SSL. The number of BRIL individuals with SRL in the range of 4-5 cm, 6-7 cm, and 8-9 cm decreased, respectively. In addition, the SRL difference between cold stress and BR-combined cold treatment was not significant (Figure 1c). under the BR-combined cold treatment (Figure 1b). After cold stress and BR-combined cold treatment, the SRL of the majority of BRIL individuals also reduced; however, it was less pronounced than the SSL. The number of BRIL individuals with SRL in the range of 4-5 cm, 6-7 cm, and 8-9 cm decreased, respectively. In addition, the SRL difference between cold stress and BR-combined cold treatment was not significant (Figure 1c). For SDW and SWW, both cold stress and BR-combined cold treatment resulted in a significant decrease in SDW and SWW of the majority of BRIL individuals. The change in the number of BRILs with SDW in the range of 0.01-0.03 g and SWW in the range of 0.1-0.2 g were remarkable, respectively. Compared to cold stress, the number of BRILs with SDW in the range of 0.02-0.03 g and SWW in the range of 0.15-0.2 g increased under the BR-combined cold treatment, respectively (Figure 1d,e).
The correlation among all four traits with different treatments in the BRIL population was analyzed. The results revealed that all traits showed continuous and approximately normal distributions. Most of the traits showed a significant positive correlation regardless of cold stress or BR-combined cold treatment ( Figure 2). For example, the significant positive correlations were observed between SSL and SRL (r = 1.00, p ≤ 0.01), SSL and SDW (r = 0.29, p ≤ 0.01), SSL and SWW (r = 0.31, p ≤ 0.01), SRL and SDW (r = 0.29, p ≤ 0.01), SRL and SWW (r = 0.31, p ≤ 0.01), and SDW and SWW (r = 0.85, p ≤ 0.01) under cold stress, For SDW and SWW, both cold stress and BR-combined cold treatment resulted in a significant decrease in SDW and SWW of the majority of BRIL individuals. The change in the number of BRILs with SDW in the range of 0.01-0.03 g and SWW in the range of 0.1-0.2 g were remarkable, respectively. Compared to cold stress, the number of BRILs with SDW in the range of 0.02-0.03 g and SWW in the range of 0.15-0.2 g increased under the BR-combined cold treatment, respectively (Figure 1d,e).
The correlation among all four traits with different treatments in the BRIL population was analyzed. The results revealed that all traits showed continuous and approximately normal distributions. Most of the traits showed a significant positive correlation regardless of cold stress or BR-combined cold treatment ( Figure 2). For example, the significant positive correlations were observed between SSL and SRL (r = 1.00, p ≤ 0.01), SSL and SDW (r = 0.29, p ≤ 0.01), SSL and SWW (r = 0.31, p ≤ 0.01), SRL and SDW (r = 0.29, p ≤ 0.01), SRL and SWW (r = 0.31, p ≤ 0.01), and SDW and SWW (r = 0.85, p ≤ 0.01) under cold stress, respectively. Under BR-combined cold treatment, the positive correlations between SSL and SRL (r = 1.00, p ≤ 0.01) as well as between SDW and SWW (r = 0.57, p ≤ 0.01) were shown similar to under cold stress. Additionally, under both cold stress and BR-combined cold treatment, all of the SSL, SRL, SDW, and SWW showed highly significant positive correlations with r = 0.19 (p ≤ 0.1), r = 0.19 (p ≤ 0.1), r = 0.43 (p ≤ 0.01), and r = 0.76 (p ≤ 0.01), respectively. These results suggested that all four traits showed continuous single-peak pattern distributions, which met the requirements for QTL mapping. Moreover, the treatment of BR improved the SSL, SDW, and SWW of some BRILs under cold stress.
respectively. Under BR-combined cold treatment, the positive correlations between SSL and SRL (r = 1.00, p ≤ 0.01) as well as between SDW and SWW (r = 0.57, p ≤ 0.01) were shown similar to under cold stress. Additionally, under both cold stress and BR-combined cold treatment, all of the SSL, SRL, SDW, and SWW showed highly significant positive correlations with r = 0.19 (p ≤ 0.1), r = 0.19 (p ≤ 0.1), r = 0.43 (p ≤ 0.01), and r = 0.76 (p ≤ 0.01), respectively. These results suggested that all four traits showed continuous single-peak pattern distributions, which met the requirements for QTL mapping. Moreover, the treatment of BR improved the SSL, SDW, and SWW of some BRILs under cold stress.

Construction of the Linkage Maps
Through the GBS approach, a total of 64.48 Gb high-quality filtered reads were filtered and 96.19% of these reads were mapped to the rice reference genome. Furthermore, a total of 10,836 SNPs were validated for the determination of recombinant events. The ratio of Q20 for BRILs was above 90% and the guanine-cytosine (GC) content was 43.55%, thus the quality of the data met the requirements for further analysis. A total of 10,836 unfiltered SNPs were validated for the determination of recombinant events. Finally, after filtration, followed by the ABH plugin, a total of 1145 bin markers were retained and used in the linkage map construction using the polymorphic markers to map all 12 rice chromosomes. All 12 linkage groups varied in the number and density of markers, with a total of 3188.33 cM in genetic distance and the individual linkage groups ranged from 186.74 to 346.46 cM in length. The interval of genetic distance of each polymorphic marker ranged from 1.9 to 146.3 cM and the average genetic distance between markers was 0.3 cM. The highest number of markers was 117 on chromosome 4, while the lowest number was 77 on chromosome 7 (Table 1).

Construction of the Linkage Maps
Through the GBS approach, a total of 64.48 Gb high-quality filtered reads were filtered and 96.19% of these reads were mapped to the rice reference genome. Furthermore, a total of 10,836 SNPs were validated for the determination of recombinant events. The ratio of Q20 for BRILs was above 90% and the guanine-cytosine (GC) content was 43.55%, thus the quality of the data met the requirements for further analysis. A total of 10,836 unfiltered SNPs were validated for the determination of recombinant events. Finally, after filtration, followed by the ABH plugin, a total of 1145 bin markers were retained and used in the linkage map construction using the polymorphic markers to map all 12 rice chromosomes. All 12 linkage groups varied in the number and density of markers, with a total of 3188.33 cM in genetic distance and the individual linkage groups ranged from 186.74 to 346.46 cM in length. The interval of genetic distance of each polymorphic marker ranged from 1.9 to 146.3 cM and the average genetic distance between markers was 0.3 cM. The highest number of markers was 117 on chromosome 4, while the lowest number was 77 on chromosome 7 (Table 1).

Prediction of Candidate Genes
The QTL regions with coincident intervals and higher LOD (more than 4.0) were used to predict candidate genes based on the Rice Genome Annotation Project website. A total of 156 candidate genes were predicted, of which 65 were annotated with known function, while the other genes were identified as retrotransposon and transposon proteins (40), expressed proteins (44), and hypothetical proteins (7), respectively (Table S2).

Validation of Candidate Genes
To further validate functions of the candidate genes, qRT-PCR was used to determine the expression levels of the 10 candidate genes under cold stress and BR-combined cold treatment in the parents and representative BRIL individuals (three excellent cold tolerance and three non-excellent cold tolerance).
The expression levels of all transcription factor genes, including LOC_Os01g40260 (OsWRKY77), LOC_Os11g45740 (MYB), LOC_Os07g05805 (OsBZR1), LOC_Os07g08440 (bHLH), LOC_Os06g44750 (AP2), and LOC_Os05g39380 (zinc finger), were higher under both cold stress and BR-combined cold treatment than those under normal temperature in all of the materials mentioned above. For different treatments, the expression levels of LOC_Os07g05805 (OsBZR1) and LOC_Os07g08440 (bHLH) under BR-combined cold treatment were higher than those under cold stress. For the maternal and paternal parents, the expression levels of LOC_Os11g45740 (MYB) and LOC_Os07g08440 (bHLH) were higher in DXWR than in SN265 under both cold stress and BR-combined cold treatment, whereas the opposite expression pattern was monitored for LOC_Os07g05805 (OsBZR1). For the representative BRIL individuals, the expression levels of LOC_Os11g45740 (MYB), LOC_Os07g05805 (OsBZR1), and LOC_Os07g08440 (bHLH) in three excellent BRIL individuals were higher than those in three non-excellent BRIL individuals under cold stress as well as under BR-combined cold treatment, whereas the opposite expression pattern was detected for LOC_Os01g40260 (OsWRKY77) (Figure 4). In addition to the candidate genes mentioned above, the expression levels of other genes LOC_Os06g01966 (auxin-induced protein), LOC_Os12g41820 (heat shock protein), LOC_Os07g01480 (oxygen evolving enhancer protein), and LOC_Os09g27660 (F-box protein) were not affected by cold stress and BR-combined cold treatment ( Figure S5). as well as under BR-combined cold treatment, whereas the opposite expression pattern was detected for LOC_Os01g40260 (OsWRKY77) (Figure 4). In addition to the candidate genes mentioned above, the expression levels of other genes LOC_Os06g01966 (auxininduced protein), LOC_Os12g41820 (heat shock protein), LOC_Os07g01480 (oxygen evolving enhancer protein), and LOC_Os09g27660 (F-box protein) were not affected by cold stress and BR-combined cold treatment ( Figure S5).

Discussion
The rich genetic diversity is an important basis for improving the agronomic trait and stresses tolerance in plants. Previous studies regarding evolution and genetic divergence showed that cultivated rice was domesticated from wild rice species. However, during the domestication process with thousands of years of evolution and natural selection, the diversity of many morphological traits was reduced, and many gene resources were lost in cultivated rice [14,[31][32][33]. As the relative ancestor of the cultivated rice, common wild rice has high diversity of desirable genes, which are very important genetic resources for

Discussion
The rich genetic diversity is an important basis for improving the agronomic trait and stresses tolerance in plants. Previous studies regarding evolution and genetic divergence showed that cultivated rice was domesticated from wild rice species. However, during the domestication process with thousands of years of evolution and natural selection, the diversity of many morphological traits was reduced, and many gene resources were lost in cultivated rice [14,[31][32][33]. As the relative ancestor of the cultivated rice, common wild rice has high diversity of desirable genes, which are very important genetic resources for improving the agronomic trait and stresses tolerance in cultivated rice [34][35][36]. DXWR (O. rufipogon) was discovered in Dongxiang County, Jiangxi Province, China, which has the northernmost distribution in all species of wild rice (28 • 14 N). Since DWXR plays important roles in basic rice research and industrial development, it is known as "the panda of wild plants" [17,18,37]. DXWR possesses abundant genetic resources associated with wide cross-compatibility, fertility restoration, cytoplasmic male sterility, high grain yield, and resistance to a range of biotic and abiotic stresses, especially low-temperature tolerance.
Although wild rice and cultivated rice have the same genome type, they exhibit many differences in their genome sequences. If the temporary genetic populations, such as F 2 or BC 1 populations and the permanent primary genetic populations, such as recombinant inbred line were constructed and used for QTL mapping, they could show not only a low QTL detection efficiency, but also poor stability. Therefore, a stable and reliable mapping population is essential for utilizing wild rice [32,33,38]. To date, multiple types of mapping populations have been constructed using DXWR as paternal parent, while a large number of QTLs associated with cold tolerance at the germination, seedling, booting, and flowering stages have been mapped to all chromosomes [16][17][18][19]30]. In this study, the female parent SN265 was used to backcross for an additional four times after a cross between the female parent SN265 and female parent DXWR (BC 4 F 1 ). Furthermore, the BC 4 F 1 population was self-crossed for eight times to obtain advanced BRIL populations (BC 4 F 8 ). The back-and self-crossing for multiple generations reduced the genetic differences caused by different genetic bases between cultivated rice and wild rice, which made the result of QTL mapping more stable and reliable.
It is well known that the low-temperature tolerance of plants is a complicated quantitative trait controlled by polygenes. To cope with low-temperature stress, plants modify their physiology, metabolism, and growth by means of signal transduction and expression regulation of many genes associated with cold tolerance when plants are exposed to lowtemperature stress [5]. A substantial number of genes that facilitate cold signaling and control the expression of cold regulons have been identified in many plants. A cascade signaling pathway, ICE-CBF-COR, is the most intensively studied, and is thought to be pretty important for cold tolerance in plants. This pathway contains the core components CBF/DREB (C-repeat binding factors/dehydration-responsive element-binding proteins) transcriptional factor, ICE (inducer of CBF expression) activated factor, and diverse downstream functional proteins called cold-regulated (COR) proteins. The CBFs belong to the AP2/ERF (apetala 2/ethylene response factor) family of transcription factors, which are the most important in the ICE-CBF-COR pathway. The CBFs control the expression of COR genes in response to cold stress, while the ICE acts as a positive upstream regulator of CBF genes [39][40][41].
In rice, 10 CBF/DREB homologous genes (OsDREB1A-OsDREB1J) have been identified. Among them, OsDREB1A, OsDREB1B, OsDREB1D, OsDREB1F, and OsDREB1G enhanced cold tolerance in Arabidopsis or rice [42,43]. In the case of QTL mapping studies, OsDREB1G and OsDREB1J were identified as the candidate genes related to cold tolerance for three QTLs, qLOP2, qPSR2-1 and qSR8-3, using the DXWR populations as the donor parent, respectively [16,30]. In addition, an ICE1-like gene, OrbHLH001, was isolated from DXWR, which enhanced cold tolerance when expressed in transgenic Arabidopsis [44]. Apart from these core components, many other activators or repressors are related to this pathway that directly or indirectly affect the low-temperature tolerance of plants. A MYB transcriptional factor, MYB15, interacted with ICE1 and was bound to MYB recognition sequences in the promoters of CBF genes, which repressed the expression of CBF genes. ICE1 activated the transcription of CBF3 by binding to MYC recognition elements in the promoter [39,45]. OsbHLH002, a homolog of ICE1 in rice, positively regulated cold tolerance by promoting the expression of OsTPP1, which encoded a key enzyme for trehalose biosynthesis [46]. Two F-box proteins, EIN3-binding F-box 1/2 (EBF1/2), positively regulated the expression of CBF genes by regulating the degradation of EIN3 and PIF3 [47,48]. In the current study, we identified some genes similar to those described above, such as bHLH transcription factor (LOC_Os07g08440), F-box and another domaincontaining protein (LOC_Os09g27660 and LOC_Os10g04590), MYB family transcription factor (LOC_Os11g45740), and AP2 domain-containing protein (LOC_Os06g44750). The transcript levels of these genes except for LOC_Os09g27660 and LOC_Os10g04590 (F-box) were upregulated both under cold stress and BR-combined cold treatment. In the representative BRILs, LOC_Os07g08440 (bHLH), LOC_Os11g45740 (MYB), and LOC_Os06g44750 (AP2) showed higher transcript levels in the excellent BRILs than the non-excellent BRILs under cold stress and BR-combined cold treatment.
Previous studies have shown that the exogenous application of BR can improve lowtemperature tolerance at different growth stages in various plants, including rice [27,28], maize [49], wheat [50], winter rye [51], and Arabidopsis [24,26], etc. To date, the molecular mechanisms regarding both cold stress and BR signaling pathway were only intensively studied in Arabidopsis. Moreover, the cold tolerance in BR signaling pathway was involved in ICE-CBF-COR signaling pathway. Under cold stress, BR directed a bHLH transcription factor CESTA (CES) to regulate the expression of the CBF and downstream COR genes in Arabidopsis [45]. BZR1 was a very important transcription factor in BR signaling pathway, acting upstream of CBF1 and CBF2 to directly enhance low-temperature tolerance in Arabidopsis. Interestingly, OsBZR1 (LOC_Os07g05805) was predicted for qSSL7-1 under BR-combined cold treatment in the current study. The expression level of OsBZR1 was upregulated both under cold stress and BR-combined cold treatment, and it showed higher transcript levels in the excellent BRILs than the non-excellent BRILs under cold stress and BR-combined cold treatment. Moreover, BZR1 regulated other genes uncoupled with CBFs, such as WKRY6, PYL6, and SOC1, etc., to improve cold tolerance in Arabidopsis [24,26]. As a negative modulator in BR signaling pathway, the protein kinase BIN2 interacted with and phosphorylated ICE1 under cold stress, which facilitated the interaction between ICE1 and the E3 ubiquitin ligase HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENE1 (HOS1), and thereby promoted ICE1 degradation. It was suggested that BIN2 mainly downregulated ICE1 abundance when the expression levels of CBF genes were attenuated [26]. In this study, OsWRKY77 (LOC_Os01g40260), a WRKY transcription factor, was also identified for qSSL1-2 under BR-combined cold treatment. It was shown that the expression level of OsWRKY77 was upregulated in the sensitive genotype, suggesting its negative role in cold tolerance at the germination stage in rice [52]. The expression levels of OsWRKY77 were upregulated both under cold stress and BR-combined cold treatment, and its expression levels in the excellent BRILs were lower than the non-excellent BRILs under cold stress and BR-combined cold treatment. In addition to the transcription factors and genes related to cold tolerance and BR signaling pathway mentioned above, some other genes related to cold tolerance were identified, such as zinc finger domain containing protein, auxin-induced protein, oxygen evolving enhancer protein, heat shock protein, etc. The precise function and regulatory mechanism of these genes still needs further investigation.
In conclusion, many QTLs and candidate genes related to cold tolerance or BR pathway were identified by mapping analysis and qRT-PCR. This study provided a basis for further mining the genes involved in low-temperature tolerance or BR signaling pathway and investigating the mechanism regulating low-temperature tolerance in rice. Further studies were needed to thoroughly investigate whether these genes are associated with both cold tolerance and BR signaling pathway in rice.

Plant Materials and Population Development
A super rice variety SN265 was used as the recipient parent to cross with the donor parent DXWR. F 1 plants were backcrossed with the recipient parent SN265 for 4 times to construct the BC 4 F 1 generation. Then, the BC 4 F 1 population was used to develop a BRIL population with 140 individuals in F 8 generation (BC 4 F 8 ).

Phenotypic Evaluation for Cold Tolerance
Evaluation of cold tolerance at seedling stage was performed according to the previous studies with minor changes [6,29,30]. For breaking dormancy, seeds were placed in a drying oven at 50 • C for 72 h. Then, the surface of the seeds, sterilized in 1% NaClO solution, was germinated and grew in a growth chamber with 75% humidity and 12 h light/12 h dark conditions. When the seedlings grew to the third leaf stage, the seedlings were watered and sprayed with sterilized distilled water (for cold treatment) and 0.1 µmol/L BR solution (for BR-combined cold treatment) for 24 h before the cold treatment, respectively. Furthermore, the seedlings were stressed at 4 • C for 7 days, and then moved to a greenhouse at 25 • C for 7 days to allow the seedlings to resume normal growth. Cold tolerance was evaluated based on the phenotypic changes of seedling shoot length (SSL), seedling root length (SRL), seedling dry weight (SDW), and seedling wet weight (SWW), respectively. SSL, SRL, SDW, and SWW of the abovementioned seedlings after recovering for 7 days were measured. The average of the three replicates of 10 seedlings for each treatment was analyzed. For the cold treatment, the data obtained under normal temperature (25 • C) were used as control. Cold tolerance scores based on the reduction rate of SSL, SRL, SDW, and SWW were calculated as follows: Reduction rate = ((data under normal temperature−data under cold treatment)/data under normal temperature) ×100%. For BR-combined cold treatment, both data obtained under normal temperature and cold treatment were used as controls, respectively. Cold tolerance scores were based on the reduction rate of SSL, SRL, SDW, and SWW = ((data under normal temperature−data under BR-combined cold treatment)/data under normal temperature) × 100% or ((data under BR-combined cold treatment−data under cold treatment)/data under BR-combined cold treatment) × 100%, respectively.
All of the experiments for evaluating cold tolerance were repeated three times under the same conditions, and the average cold tolerance scores from three replicates were used for the QTL mapping analysis. R software 3.5.1 (R Core Team. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing) was used to analyze the frequency distributions and correlations of all 4 phenotypic traits (www.R-project.org, accessed on 1 July 2022).

Genotyping-by-Sequencing and SNP Identification
Young leaves of 140 BRIL individuals and their parental lines were collected separately for genomic DNA extraction using a DNA extraction kit (Aidlab, Beijing, China). DNA libraries from BRILs were constructed using the GBS technology. The restriction enzymes PstI and MspI and T4 ligase were used for digestion and ligation, respectively. The DNA libraries were enriched and sequenced using Illumina HiSeq 4000 instrument (Huada Gene Technology Co., Ltd., Shenzhen, China). The Tassel software was used to analyze the raw data and high-quality SNPs between parents were termed by alignment with Nipponbare reference genome [53].

Construction of Linkage Map and QTL Analysis
The Illumina data of BRIL population was used to construct a linkage map using the QTL IciMapping sofware v4.163 (Meng et al. 2015, Beijing, China) [54]. The genotypic maps were aligned and split into recombination bins according to the recombination breakpoints, with the parameter of window size of 10 cm and walk speed of 1 cm.
QTLs were mapped using the inclusive composite interval mapping of the additive (ICIM-ADD) mapping method. QTLs were computed by a permutation test involving 1000 runs at a significance level of p = 0.05. The threshold for the logarithm of odds (LOD) scores was set to 3.0. QTL nomenclature was followed by the method of McCouch et al. (2008) [55].

Prediction and Validation of Candidate Genes
To predict potential candidate genes within QTL intervals, the physical positions of SNP markers flanking the QTLs were searched in the Rice Genome Annotation Project website (http://rice.plantbiology.msu.edu, accessed on 1 July 2022). The seedlings at the third leaf stage from 2 parental lines, 3 excellent cold tolerance BRIL individuals, and 3 non-excellent cold tolerance BRIL individuals at normal temperature (25 • C), for 12 h under cold stress (4 • C), and for 12 h under BR-combined cold treatment (pre-watering and pre-spraying with 0.1 µmol/L BR solution for 24 h, 4 • C) were collected to extract total RNA using the Trizol reagent. Quantitative real-time PCR (qRT-PCR) was used to investigate the expression patterns of 10 representative candidate genes. Moreover, 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 15 s, 60 • C for 30 s, and 72 • C for 1 min, and 1 cycle of 72 • C for 15 min. Each sample was analyzed in three biological and three technical replicates. The relative expression levels were calculated via the ∆Ct method [56]. The actin gene of rice (LOC_Os03g50885) was used as a reference gene. All of the primer sequences are listed in Table S1.