Skip to main content

Combining QTL-seq and linkage mapping to fine map a candidate gene in qCTS6 for cold tolerance at the seedling stage in rice



Cold stress caused by low temperatures is an important factor restricting rice production. Identification of cold-tolerance genes that can stably express in cold environments is crucial for molecular rice breeding.


In this study, we employed high-throughput quantitative trait locus sequencing (QTL-seq) analyses in a 460-individual F2:3 mapping population to identify major QTL genomic regions governing cold tolerance at the seedling stage in rice. A novel major QTL (qCTS6) controlling the survival rate (SR) under low-temperature conditions of 9°C/10 days was mapped on the 2.60-Mb interval on chromosome 6. Twenty-seven single-nucleotide polymorphism (SNP) markers were designed for the qCST6 region based on re-sequencing data, and local QTL mapping was conducted using traditional linkage analysis. Eventually, we mapped qCTS6 to a 96.6-kb region containing 13 annotated genes, of which seven predicted genes contained 13 non-synonymous SNP loci. Quantitative reverse transcription PCR analysis revealed that only Os06g0719500, an OsbZIP54 transcription factor, was strongly induced by cold stress. Haplotype analysis confirmed that +376 bp (T>A) in the OsbZIP54 coding region played a key role in regulating cold tolerance in rice.


We identified OsbZIP54 as a novel regulatory gene associated with rice cold-responsive traits, with its Dongfu-104 allele showing specific cold-induction expression serving as an important molecular variation for rice improvement. This result is expected to further exploration of the genetic mechanism of rice cold tolerance at the seedling stage and improve cold tolerance in rice varieties by marker-assisted selection.

Peer Review reports


Rice (Oryza sativa L.) is an important food crop that has adapted to a tropical climate. Therefore, cold damage of rice at high latitudes and high altitudes has become an important factor restricting its production [1]. Specifically during the early stage of vegetative growth, chilling directly causes physiological damage to seedlings, thereby hindering the vegetative formation of the plants and ultimately affecting yield [2]. Therefore, the key to solving this problem is to explore cold tolerance (CT) genes with strong advantages that contribute alleles to ultimately improve plant CT.

CT in rice is a quantitative trait controlled by multiple genetic and environmental factors. In the previous 20 years, numerous double haploid lines, back-cross inbred lines, recombinant inbred lines, and near-isogenic lines have been developed using CT varieties, such as Dongxiang wild rice [3], KunmingXiaobaigu [4], LijiangXintuanheigu [5], Koshihikari [6], M202 [7], and Norin-PL8 [8], as donors, and >250 CT-related quantitative trait loci (QTLs) have been identified on 12 chromosomes [9]. To date, >80 QTLs for CT at the seedling stage (CTS) in rice have been identified [2, 10,11,12]. Moreover, most QTLs that are linked with CTS are located within a range of 10 cM to 30 cM through primary populations, making it difficult to perform molecular-marker-assisted breeding. Thus, the key to solving this problem is to construct a fine-mapping population. To date, qCTS4 [13], qCtss11 [14], qSCT1 [15], qSCT11 [15], and qLOP2/qPSR2-1 [16] have been fine-mapped, with only seven genes responding to CTS isolated by map-based cloning (COLD1, qCTS-9, qPSR10, bZIP73, HAN1, OsLTPL159 and qCST10). COLD1 [17] encodes a G protein signal regulator, and its overexpression can significantly improve rice CT, whereas rice lines with no or low expression are cold-sensitive. The promoter region of qCTS-9 [18] has an indel marker significantly associated with low-temperature tolerance, and its overexpression increases rice tolerance to low temperatures. qPSR10 [19] was detected in a genome-wide association study, and its CT haplotype was found in japonica rice resources. bZIP73 [1] has a single-nucleotide polymorphism (SNP) on the genome of japonica and indica rice, with this SNP determining the tolerance of rice to low temperatures. HAN1 [20] regulates low-temperature response and CTS mediated by jasmonic acid and is a negative regulator of CTS. OsLTPL159 [21] encodes a nonspecific lipid-transfer protein, the coding sequence (CDS) for which is located in the mapped region of qCST10 [22]. These examples demonstrate that molecular regulation of CTS is extremely complex, and although these studies have laid the foundation for understanding the molecular mechanism of CT in early rice growth, dissecting the genetics behind CT in rice is still at its infancy. Identification of these QTLs by biparental cross-linkage mapping was a labor- and time-intensive process that included mapping the genotypes a large number of individuals in a segregated population; therefore, to obtain reliable QTLs, CTS improvement needs to be accomplished through fine mapping, gene cloning, functional verification, and finally application for molecular-marker-assisted breeding.

The rapid development of next-generation sequencing (NGS) technology has enabled the combination of bulked segregant analysis (BSA) with NGS to become an important method for mining QTLs and genes [23, 24]. Compared with traditional linkage mapping, QTL sequencing (QTL-seq) can improve work efficiency and provide high-density variation. QTL-seq has been successfully applied for many plants, including cucumber [25], soybean [26], rice [27], and tomato [28]. Although the technique enables rapid location of the major QTL, confidence interval (CI) resolution remains less than optimal. Therefore, researchers have had to combine QTL-seq with fine mapping [29] and RNA-seq [30] in order to identify candidate genes with a QTL. For example, CaqSW1.1 [31] was mapped from 1.37 Mb to 35 kb by QTL-seq and classical QTL mapping to allow identification of final candidate genes. Additionally, four candidate genes (SlCathB2, SlGST, SlUBC5, and SlARG1) associated with heat tolerance were detected by QTL-seq and RNA-seq in tomato [28]. Thus, utilization of a combined method for identifying a major QTL interval via QTL-seq and fine mapping to mine target genes provides a new perspective.

In this study, we employed two strategies (QTL-seq and fine mapping) to identify the genes for CTS in rice. An F2:3 population derived from a cross between Dongnong430 (a cold-sensitive variety) and Dongfu104 (a strong CT variety) were used for QTL-seq analysis, and qCTS6 was mapped in the 2.60-Mb genomic region on chromosome 6. Furthermore, we used a fine-mapping strategy to anchor qCTS6 in the 96.6-kb interval, where 13 non-synonymous SNPs (nSNPs) were identified within the seven-gene intervals, which were predicted to result in an amino acid change. Quantitative reverse transcription (qRT)-PCR and haplotype analysis revealed that OsbZIP54 was strongly induced by cold stress at the seedling stage, and that +376-bp (T>A)Dongfu104 can be considered a rare functional variation of japonica rice in northern China. These results revealed the mechanism by which cold sensitivity changes into CT in rice at the molecular level and improved the understanding of the mechanisms underlying CTS in rice.


Screening and evaluation of the CT score

Two parental varieties, DN430 and DF104 (Fig. 1A, B) along with their 460 F2:3 lines (Fig. 1C), were evaluated under normal and cold conditions of 9°C/10 days. The SR of the CT-variety DF104 was significantly higher than that of cold-sensitive DN430 (Fig. 1B), indicating that DF104 had a stronger CT than DN430. In the F2:3 population, the absolute values of skewness and kurtosis associated with the SR were close to 1 (Table S1), indicating the data as suitable for QTL analysis. Moreover, the F2:3 population showed different SRs from 0 to 100, with Fig. 1D showing that among the 460 lines, 30 CT and 30 cold-sensitive lines were identified. These were then selected to prepare the CT (T-pool) and sensitive pools (S-pool), respectively, which were then used for DNA re-sequencing.

Fig. 1
figure 1

Phenotype analysis of CTS in rice. A Performance of DN430 and DF104 under normal conditions. B Performance of DN420 and DF104 under cold treatment (9°C/10 days). C Identification of CT in the F2 population in an artificial climate chamber. D Frequency distribution of SR variation measured among 460 mapping individuals of F2:3 populations

Whole-genome re-sequencing and QTL-seq analysis

The average genome-coverage depth of the parents and two pools was 50×. When compared with the ‘Nipponbare’ reference genome, we obtained a total of 786,407 SNPs and 133,992 indels for four samples. After trimming and filtering, 521,630 SNPs and 85,824 indels were obtained (Table S2). From these SNPs, a total of 356,831 high-quality SNPs homozygous for each parent and showing polymorphism between the parents were developed for further QTL-seq analysis.

For all of the obtained SNPs, the ∆(SNP-Index) (Fig. 2A), Euclidean distance (ED) (Fig. 2B), G-value (Fig. 2C), and Fisher’s exact test (Fig. 2D) were used to intercept the association interval. We then identified the candidate regions associated with CTS in the rice genome (Table 1). For the ∆(SNP-Index), two major peaks on chromosomes 3 and 6 were identified for CTS (P < 0.01) spanning 2.58-Mb (1.50–4.08 Mb) and 3.24-Mb (28.00–31.24 Mb) intervals. The significant peak given by the ED algorithm completely covered the results returned by the ∆(SNP-Index) algorithm, and the two major intervals were subsequently named qCTS3 and qCTS6. The G algorithm similarly captured the results of the first two algorithms and identified another new interval (0.00–2.42 Mb) in chromosome 12 (named qCST12). Only qCTS6 was obtained when the four methods were applied. By comparing the association results of the four methods, their intersections of the genome regions were consistent with the association regions identified by Fisher’s exact test. A total of 2.60 Mb was found in the qCTS6 interval. Notably, in the three algorithms, the peak value obtained by qCTS6 was higher than that of qCTS3 (Table S3), indicating a significant difference in the allele ratio between the two mixed pools. Therefore, qCTS6 was considered the most significant CTS target for exploring candidate genes. According to statistical analysis, the qCTS6 coding region harbored 407 SNPs and 18 indels (Table S4).

Fig. 2
figure 2

QTL analysis of CTS using QTL-seq. A Manhattan plot showing the distribution of Δ(SNP-index) on chromosomes, B the distribution of the square of the ED on the chromosomes, C the distribution of G-value on the chromosomes, and D the distribution of −log10(P-value) on the chromosomes based on Fisher’s exact test. The blue and red lines represent 95 and 99% CIs, respectively. The black lines are average values of the four algorithms and were drawn by sliding window analysis. The number on the horizontal coordinate represents the chromosome number

Table 1 QTLs conferring cold tolerance in four method identified using QTL-seq

Narrowing of qCTS6 to a fine region

Because the 2.60-Mb region contained a large number of genes, we performed fine mapping of qCTS6 based on the base-variation information of the 2.60-Mb region to identify effective SNPs that might be associated with CTS according to the following criteria: Δ(SNP-Index) values >0.48, which would be significantly higher than the statistical upper bound of the CIs under the null hypothesis of “no QTLs at the P < 0.01 level”; and the selected SNP represents a non-synonymous mutation uniformly distributed within that region. We ultimately captured 27 excellent nSNPs within the qCTS6 interval, with these used for genotype scanning of the 460 lines and to obtain a qCTS6 linkage map (Fig. 3B). We used a logarithm of the odds (LOD) value >3.0 as a threshold for consecutive occurrence to determine the existence of the QTL. We then analyzed the seedling SR of the 460 lines under cold stress by linkage mapping. qCTS6 was simultaneously linked with SR and anchored to the 96.6-kb interval between 30,467,391 bp and 30,563,973 bp (Fig. 3B), thereby explaining 25.83% of the phenotypic variation in the SR (Table 2). The positive allele of qCTS6 was contributed by DF104, and CT6 was optimized to a physical interval of 96.6 kb. Other QTL intervals (29,739,644–29,771,312 bp) with minor effects (explaining 3.17% of the phenotypic variation) were not used in subsequent analysis due to their lower LOD value.

Fig. 3
figure 3

Further mapping of the qCTS6. A The Manhattan diagram shows the location of the qCTS6 on chromosome 6. The number on the horizontal coordinate represents the physical position (bp) on the chromosome. B qCTS6 detected by using the inclusive composite interval mapping (ICIM) module of QTL IciMapping 4.2. Linkage map based on SNP markers. The vertical lines on the blue bar represent the name of the SNP marker (all SNP marker information is given in Table S9). The LOD scores are shown with two distinct peaks corresponding to CTS. C Thirteen genes in the qCTS6 region were obtained from annotation information on the Nipponbare genome

Table 2 Identification of qCTS6 for cold tolerance at seedling stage by linkage analysis

Putative candidate genes for qCTS6

The 96.6-kb interval was intercepted from annotation information of the Nipponbare genome (, resulting in capture of 13 annotation genes (Fig. 3C). Except for FON1 [32] or the functional genes in the number of cloned flower organs, no other genes were evaluated. Previous studies reported important genes or gene families, such as those encoding pseudouridine synthase family protein [33], C3HC4 [34], hypothetical protein [18], protein phosphatase 2C [35], and the bZIP family transcription factor (TF) [36,37,38], related to CT in rice or other crops. Although similar genes or proteins were included in the list of 13 candidate genes, the 10 genes obtained from Gene Ontology (GO) analysis were not annotated with the phrase “response to cold”. The most enriched terms related to biological process, molecular function, and cellular component ontologies were metabolic processes (i.e., GO:0044267, GO:0008152, and GO:0034641), catalytic activity (i.e., GO:0003824 and GO:0140096), and cell (i.e., GO:0005623 and GO:0044464), respectively (Table S5 and Figure S1). Further, the 96.6-kb region on chromosome 6 harbors 24 SNPs located in the coding region between the parental lines and a Δ(SNP-Index), G-value, ED, and two-tailed Fisher’s exact test P-value higher than the statistical confidence at P < 0.01. Additionally, 89 SNPs and 18 insertions were located in the intergenic, intronic, untranslated, upstream, and downstream regions. Seven of 13 genes carried 13 nSNPs (Table S6), with four nSNPs identified using both parents as a reference and affecting four candidate genes encoding pseudouridine synthase (Os06g0717400), protein phosphatase 2C (Os06g0717800), a hypothetical protein (Os06g0718000), and a bZIP TF (Os06g0719500) (Table 3). As noted, these genes might play an important role in responding to cold stress in plants. Based on this finding, we identified the four final possible candidate genes using the following three approaches: the possible relationship between gene response and cold stress, non-synonymous mutations in the protein-coding region, and participation of the gene products in specific metabolic pathways that affect CT.

Table 3 Identification of SNPs in putative candidate genes for CTS

To identify strong candidate genes in qCTS6, we analyzed the expression levels of 13 genes under cold stress and normal conditions, revealing four genes showing significant changes in expression under cold stress. Among these, Os06g0761700, Os06g0717700, and Os06g0718400 levels differed significantly at a single time point, with only Os06g0719500 showing a strong induction of cold-stress response in parents with extremely strong CT (Fig. 4) and suggesting Os06g0719500 as a possible candidate gene for qCTS6. By contrast, under normal conditions, all 13 candidate genes showed no significant changes in expression (Figure S2). Notably, Os06g0719500 encodes the bZIP family TF annotated as OsbZIP54, with the SNP marker SNP-630563973 located in this coding region. Moreover, a strong correlation between mRNA levels and CT in these 60 lines indicated the regulatory function of OsbZIP54 on CT (Fig. 5). These findings suggested OsbZIP54 as a possible candidate CT gene on qCTS6 for further analysis.

Fig. 4
figure 4

Expression levels of the 13 candidate genes in DN430 and DF104 in response to cold stress and measured by qRT-PCR. **P < 0.01, Student’s t test

Fig. 5
figure 5

Correlation between OsbZIP54 expression and SR in 30 CT lines and 30 cold-sensitive lines following cold treatment

Haplotype analysis of candidate genes

To reveal the significance of variations in the OsbZIP54 coding region, we selected 60 lines from the T-pool and S-pool, respectively, and sequenced their OsbZIP54 genes by Sanger sequencing. We identified two SNPs (SNP-30563973 and SNP-30565007) on the OsbZIP54 exon and intron (Figure S3), respectively, with no polymorphism detected in the promoter region between the two parents (Figure S4). To confirm the association of haplotypes with the studied traits, violin plots were generated to assess trait distribution within haplotype groups (Fig. 6) We then found that 30 CT and 30 cold-sensitive lines showed consistent isolation in the +376-bp (T>A) coding region of OsbZIP54 (Table S7), suggesting that this might lead to phenotypic differences in CT. Moreover, the other SNP (SNP-30565007) was not associated with CT, as OsbZIP54 was isolated from the japonica rice variety DF-104. To clarify the distribution of +376-bp (T> A) in rice, we used data from the 3010 Rice Genomes Project and found that only 275 japonica rice genomes carried +376-bp (T>A), whereas indica rice genomes harbored >1600 instances. The proportion of OsbZIP54 functional variation in indica rice and japonica rice was 5.8:1, indicating that +376-bp (T>A) was a rare variation in japonica rice as compared with indica rice. We subsequently analyzed the variation types of OsbZIP54 in 295 northern japonica rice (Table S8), finding no instance of +376-bp (T>A) in the samples. This suggested that +376-bp (T>A) was a natural variation; therefore, +376-bp (T>A) can be used as a special rare functional variation in DF104. However, it cannot be ruled out that other japonica OsbZIP54 might harbor a similar single-base variation.

Fig. 6
figure 6

Haplotype–trait violin plots. Hap1 (TG), Hap2 (AA), Hap3 (AG), Hap4 (TA)


Rice is a staple food for ~50% of the world population [39]. Due to its origin in tropical and subtropical regions, rice is more sensitive to cold stress than other cereal crops. Notably, rice seedlings are sensitive to low temperatures in early Spring in temperate and subtropical zones, as well as in high-elevation areas, leading to delayed growth of seedlings, yellowing, withering, reduced tillers, and hindered growth [2]. To locate the major QTLs of CT, we used several approaches to obtain the most reliable QTL interval. First, we used a large population to identify their phenotypes under low-temperature conditions, thereby obtaining accurate CT phenotypes for each line. Second, two DNA pools with strict CT phenotypic differences were used to perform QTL-seq analysis. Moreover, we used four bioinformatics analysis approaches to map QTL regions at the 99% significance level. Finally, two highly significant peaks (qCTS3 and qCTS6) were detected on chromosomes 3 and 6, with the former mapped between 1.50 Mb and 4.08 Mb and the latter between 28.44 Mb and 31.04 Mb. QTLs and genes for CT on chromosomes 3 and 6 have been previously reported. The QTL qLVG3 was studied for CT at the bud stage of rice according to the ‘Nipponbare’ reference genome (chromosome 3 at a position 5.92-Mb away from qCTS3) [40]. qCTS3 was located in the 1,500,000- to 4,080,000-bp region, and CDPK13 [41] and OsCDPK13 [42] were detected in this region. In chromosome 6, the location of qCTS6 in the present study completely overlaps with the previously reported region of qLVG6 and is linked with CT during the budburst period [40]. Additionally, qRCT6b, a QTL controlling CT at the reproductive stage, was ~23.00-Mb distant from qCTS6. Moreover, we found no CT QTLs related to the seedling stage on chromosome 6 and no cold-response gene has been reported within the qCTS6 interval. Notably, implementation of three algorithms [∆(SNP-Index), ED, and G-value] resulted in a higher peak value for qCTS6 (Table S3) than that for qCTS3. Furthermore, compared with qCTS6, qCTS3 was not associated with Fisher’s exact test. This suggested the presence of a significant difference in the allele ratio between the two mixed pools. Therefore, we considered qCTS6 as the most significant target for CTS and exploration of candidate genes.

Using the strategy of traditional QTL mapping combined with QTL-seq, we mapped the Cf-10 gene to a 790-kb region in Cladosporium fulvum [43], and GaFzl was mapped to a 70-kb region containing seven annotated genes [44]. Therefore, we believe that strategies combining QTL-seq and traditional linkage mapping for uncovering very narrow candidate regions are well suited for faster targeting of a target gene. To screen potential candidate genes in the qCTS6 interval, we combined QTL-seq and fine mapping, which enabled reduction of the number of candidate genes within the qCTS6 interval defined by QTL-seq from 393 to 13 genes. Among these 13 genes, seven have 13 functional base variations. Furthermore, the expression level of OsbZIP54 in the two parents showed that it was strongly induced after cold treatment. Expression analysis of OsbZIP54 in 60 lines with differences in extreme CT revealed a correlation between CT and expression level (Fig. 5). Additionally, sequence analysis revealed two SNP differences in the coding-sequence region of OsbZIP54, in which +376-bp (T>A) is a functional variant, and that play a major role in determining the difference in CT (Table S7). However, as a newly identified CT regulator, the signaling-transduction system and the downstream pathways of OsbZIP54 remain unclear. DF104 carries a rare exonic SNP, and OsbZIP54 was highly induced under cold conditions, suggesting that the SNP might be involved in transcriptional regulation of OsbZIP54, which warrants further exploration. Interestingly, despite the strong expression of OsbZIP54 induced by cold stress, the promoters of DF-104 and DN-430 do not harbor polymorphisms. A previous study reported that bZIP73 is induced and upregulated by cold stress, and that there is a SNP difference at the 511-bp position of the coding region that determines CT of the rice [36]. In the present study, we identified no polymorphisms the bZIP54 promoter. However, another study showed that bZIP73–bZIP71 interactions might facilitate transcriptional activation of bZIP73, which subsequently confers a cold-signal switch in rice. This provides evidence that a possible reason for upregulated OsbZIP54 expression is cold stress. However, the expression level of CTB4a under cold stress is affected by promoter mutations and functional SNPs [45]; therefore, the genetic mechanism of rice in response to cold stress is complex, and significant alterations in the expression of genes under cold stress are not solely determined by their respective promoters. Next, we will further clone the alleles of the candidate gene, construct expression vectors for plants and yeast, and express these vectors in rice protoplasts, yeast, and Arabidopsis to compare the effects of this natural allele variation on gene expression and subcellular localization under cold stress. The biological functions of OsbZIP54 and its associated genomic variations need to be confirmed using gene editing and high-efficiency overexpression transformation systems. Although we found correlations between the nSNP in Os09g0444200 and the studied phenotypic traits, this is insufficient to assign a CT phenotype to a single SNP. By contrast, we also believe that a favorable allelic variant in one key gene is insufficient to provide CT, and that the final CT phenotype of the plant should be considered a result of a combination of favorable allelic variations from different key genes.

CT is a complex trait controlled by multiple genes [46], and those genes isolated thus far can be divided into two different types. One includes genes with specific functions and that participate in systemic metabolism to defend against cold stress [47]. For example, OsLTPL159 enhances the CT of rice at the early seedling stage by decreasing the toxic effect of reactive oxygen species, enhancing cellulose deposition in the cell wall, and promoting osmolyte accumulation [21]. Additionally, CTB4a, as a leucine-rich repeat receptor-like kinase that positively regulates ATP activity and content under cold stress by interacting with the β subunit of ATP synthase, thereby increasing pollen fertility [48]. The second type comprises genes that regulate gene expression during stress responses (i.e., signaling components and TFs) [47]. Together with their target genes, TFs can serve as regulators of signal transduction to activate or inhibit genes involved in cold-stress response. Therefore, TFs are excellent candidates for modifying complex traits of crops [49]. Conversely, there are 91 bZIP family member genes in rice [50], several of which are involved in the regulation of CT. Notably, as a transcriptional activator, the expression of OsbZIP52 was strongly induced at low temperature (4°C) [51]. Other bZIPs, such as bZIP73, interact with bZIP71 to regulate the level of abscisic acid and the balance of active oxygen, thereby improving the tolerance of rice to low temperatures [36]. However, in northern japonica rice, the potential of bZIP TFs in CT remains to be explored. In the present study, we identified a rare variation in OsbZIP54 carried in northern China japonica rice, and haplotype events did not occur in 295 northern cultivars (Table S8). Additionally, the 3010 Rice Genome Project database shows that OsbZIP54 has a very small occupancy rate in japonica rice as compared with indica rice. A rare thymine variation on the OsbZIP54 exon was found in DF104, which could serve as a genomic marker for improved CT. Moreover, the molecular mechanism of the interaction between OsbZIP54 and its targeted promoter elements, downstream genes, and cofactors warrant further investigation.


In conclusion, study of the genetic basis of CTS is important for breeding purposes. In this study, we used QTL-seq and linkage-map analysis to identify the qCTS6 QTL related to rice CT. Experimental verification by qRT-PCR and haplotype analysis indicated OsbZIP54 in the qCTS6 region as a candidate gene. Furthermore, the results showed that the Dongfu104 allele exhibiting a specific cold-induction expression pattern could serve as an important molecular marker for rice improvement.

Materials and methods

Plant materials

We used two japonica varieties [cold-sensitive female parent Dongnong 430 (DN430) and cold-tolerant male parent Dongfu104 (DF104)] as parental lines to develop the 460 F2:3 population. These materials were obtained from Northeast Agriculture University (Harbin, China).

CTS evaluation

The seeds were incubated at 40°C for ~36 h to break dormancy and then soaked in deionized water at 30°C for ~72 h for germination. The germinated seeds were transplanted into 96-well PCR plates for seedling growth in a greenhouse environment (26°C day/22°C night temperature, 12-h/12-h light/dark cycle, and 80% relative humidity). When the seedlings grew to trifoliate, they were treated to a low temperature of 9°C/10 days in an incubator, and all 460 F2:3 seedlings were allowed to resume growth at 28°C for 10 days to investigate the SR (measured as the percentage of total seedlings that survived relative to the total number tested) of 60 plants per replicate. Three biological replicates were performed.

Construction of segregating pools

All young leaves of 460 F2 individuals were collected separately for total genomic DNA extraction using a cetyltrimethylammonium bromide method [52], with minor modifications. The genomic DNA of 30 extremely CT (ECT) and 30 extremely cold-sensitive (ECS) individuals were selected as two bulked pools according to the SR of the F2:3 population (range: 0–100%). [To simplify the following description, we abbreviate DF104 as T, DN430 as S]. Isolated DNA was quantified using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). All DNA from the 30 ECT and 30 ECS plants were quantified on a Qubit 2.0 fluorimeter (Life Technologies, Carlsbad, CA, USA), and equal amounts of DNA from the ECT and ECS plants were mixed.

NGS sequencing and BSA-seq analysis

Total genomic DNA was extracted from bulked pools, and at least 3 μg genomic DNA was used to construct paired-end libraries with an insert size of 500 bp using the paired-end DNA sample prep kit (Illumina, San Diego, CA, USA). These libraries were sequenced using the HiSeq X10 (Illumina) NGS platform at Genedenovo (Guangzhou, China). Quality trimming is an essential step in generating high-confidence variant calling. Raw reads were processed to obtain high-quality clean reads according to three stringent filtering standards: 1) removing reads with ≥10% unidentified nucleotides; 2) removing reads with > 50% bases having phred quality scores ≤ 20; and 3) reads aligned to the barcode adapter.

To identify SNPs and indels, filtered reads were aligned to the Nipponbare reference genome sequence [53] using the Burrows–Wheeler Aligner (v.0.7.16a-r1181) with parameter ‘mem -M’; -M is an option used to mark shorter split-alignment hits as secondary alignments [54]. Variant calling was performed using the GATK UnifiedGenotyper (v.3.5; SNPs and indels were filtered using the GATK VariantFiltration function with proper standards (-Window 4, -filter "QD < 4.0 || FS > 60.0 || MQ < 40.0 ", -G_filter "GQ < 20"). All mutations were annotated for genes and functions, as well as genomic regions, using ANNOVAR [55]. Association analysis was performed using the SNP-Index [23], ∆(SNP-Index) [56], calculation of the G statistic [57, 58], ED [59], and two-tailed Fisher’s exact test [60] based on the SNPs. Average values were calculated using a 2000-kb sliding window with a step size of 20 kb and discarding of windows with less than 10 SNP/Indel. If the number of SNPs was insufficient, the results of this window were merged into the next window. The overlapping interval of the four methods was considered as the final QTL interval of CT.

Development of SNP markers and narrowing the candidate interval

To develop markers to validate the QTL-seq results and narrow the candidate region, 27 nSNPs with significant peaks were identified based on the CDS sequence of the two parents in the qCST6 region under the four association-analysis algorithms [∆(SNP-Index), G-value, ED, and Fisher’s exact test] and were considered as markers for fine mapping. Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA) was used to design primers (Table S9). All markers were selected for polymorphisms between parents, and genotyping of the 460 individuals was performed using the polymorphic markers, which were used to construct the linkage map and narrow the candidate region using the ICIM module of QTL IciMapping (v.4.2; The threshold of the LOD score for declaring the presence of a significant QTL was determined by the permutation test with 1000 repetitions at P < 0.01.

Prediction of candidate genes for the CTS

To predict possible candidate genes for the CTS, the following three strategies were comprehensively considered. First, candidate gene prediction was performed to compare the DNA sequences of genes within the QTL regions using the whole-genome DNA re-sequencing database of the parents. In the QTL regions, we focused on open reading frames with non-synonymous mutant SNPs between the two parents. Second, candidates with the highest probability of success were reselected according to their functional annotations and used to screen the rice genome database ( Third, the GO [61] and Kyoto Encyclopedia of Genes and Genomes [62] databases were used as references for the evaluation of candidate genes using BLAST software [63].

Verification of the expression of candidate genes

Seedling leaves of DN430 and DF104 were collected at 2 h, 6 h, 12 h, 24 h, and 48 h after cold treatment in three replicates, frozen in liquid nitrogen, and stored at −80°C until total RNA extraction. Control plants were also collected and stored similarly. qRT-PCR was used to quantify the expression levels of candidate genes under cold treatment. Total RNA was extracted from rice tissues using TRIzol reagent (Thermo Fisher Scientific) and treated with DNase I to eliminate any DNA contamination. RNA quality was assessed by electrophoresis and stored at −80°C until use. First-strand cDNA (10 μL) was synthesized according to the instructions for the PrimeScript RT Master Mix (Takara Biomedical Technology (Beijing) Co., Ltd., Beijing, China). Primers were designed with Primer Premier v. 5.0 (Premier Biosoft International) and are listed in Table S9. The housekeeping gene actin1 (Os05g36290) was used as an internal control [64], and qRT-PCR was performed using a LightCycler 2.10 system (Roche, Basel, Switzerland) with 2× SYBR Green I PCR Master Mix (Thermo Fisher Scientific). qRT-PCR analysis was performed as previously described [65].

Diversity analysis of candidate genes

Haplotypes of candidate genes were analyzed, as follows: 30 T-pool, 30 S-pool, and 295 China northern japonica rice. The 295 rice varieties were from the Heilongjiang, Jilin, and Liaoning provinces on China and other countries, including Japan, the Republic of Korea, the Democratic People’s Republic of Korea, and Russia [66]. Primer sequences were designed for candidate genes at SNP sites between parents, and the DNA fragments of the sequence in 60 lines were amplified by PCR to determine haplotype statistics of the distribution of target genes. All specific primer sequences are provided in Table S2. The PCR reaction mixture had a total volume of 20 μL and comprised 1.5 μL of forward primer (10 μm), 1.5 μL of reverse primer (10 μm), 2 μL of genomic DNA (50 ng/μL), 5 μL of ddH2O, and 10 μL of Pfu master mix (Cwbio, Beijing, China), which included Taq DNA polymerase, PCR buffer, Mg2+, and dNTPs. The PCR reaction was performed in an Eppendorf 5333 Mastercycler (Eppendorf, Hamburg, Germany) using the same protocol as that used for qRT-PCR [65]. The products were examined by 1% agarose gel electrophoresis. Direct sequencing of PCR products was performed by BGI Life Technology Co., Ltd. (Shenzhen, China). Haplotype trait data were plotted in the ggplot2 R package (ggplot2: elegant graphics for data analysis; Springer, NY, USA).

Statistical analysis

The difference between the two parents and progeny for CT was detected using SPSS18.0 software. Data represent means ± standard deviations. Correction analysis was performed to evaluate correlations between gene expression and seedling survival rates using SigmaPlot software (v.12.5; Systat Software, San Jose, CA, USA).

Availability of data and materials

The QTL-seq data for this study can be found in the National Center for Biotechnology Information Sequence Read Archive under the accession numbers SSR13319856, SSR13319857, SSR13319858, and SSR13319859 (



Bulked segregant analysis


Coding sequence


Confidence interval


Cold tolerance


Cold tolerance at the seedling stage


Euclidean distance


Inclusive composite interval mapping


Logarithm of odds


Non-synonymous SNP


Quantitative reverse transcription PCR


Quantitative trait loci


Single-nucleotide polymorphism


Survival rate


Transcription factor


  1. Liu CT, Schlappi MR, Mao BG, Wang W, Wang AJ, Chu CC. The bZIP73 transcription factor controls rice cold tolerance at the reproductive stage. Plant Biotechnol J. 2019;17(9):1834–49.

    CAS  PubMed  PubMed Central  Google Scholar 

  2. Andaya VC, Mackill DJ. Mapping of QTLs associated with cold tolerance during the vegetative stage in rice. J Exp Bot. 2003;54(392):2579–85.

    Article  CAS  PubMed  Google Scholar 

  3. Mao DH, Yu L, Chen DZ, Li LY, Zhu YX, Xiao YQ, et al. Multiple cold resistance loci confer the high cold tolerance adaptation of Dongxiang wild rice (Oryza rufipogon) to its high-latitude habitat. Theor Appl Genet. 2015;128(7):1359–71.

    Article  CAS  PubMed  Google Scholar 

  4. Dai LY, Lin XH, Ye CR, Ise K, Saito K, Kato A, et al. Identification of quantitative trait loci controlling cold tolerance at the reproductive stage in Yunnan landrace of rice, Kunmingxiaobaigu. Breed Sci. 2004;54(3):253–8.

    Article  CAS  Google Scholar 

  5. Shirasawa S, Endo T, Nakagomi K, Yamaguchi M, Nishio T. Delimitation of a QTL region controlling cold tolerance at booting stage of a cultivar,‘Lijiangxintuanheigu’, in rice, Oryza sativa L. Theor Appl Genet. 2012;124(5):937–46.

    Article  CAS  PubMed  Google Scholar 

  6. Endo T, Chiba B, Wagatsuma K, Saeki K, Ando T, Shomura A, et al. Detection of QTLs for cold tolerance of rice cultivar ‘Kuchum’and effect of QTL pyramiding. Theor Appl Genet. 2016;129(3):631–40.

    Article  CAS  PubMed  Google Scholar 

  7. Andaya V, Tai T. Fine mapping of the qCTS12 locus, a major QTL for seedling cold tolerance in rice. Theor Appl Genet. 2006;113(3):467–75.

    Article  CAS  PubMed  Google Scholar 

  8. Saito K, Hayano Saito Y, Kuroki M, Sato Y. Map-based cloning of the rice cold tolerance gene Ctb1. Plant Sci. 2010;179(1-2):97–102.

    Article  CAS  Google Scholar 

  9. Yang LM, Liu HL, Zhao HW, Wang JG, Sun J, Zheng HL, et al. Mapping quantitative trait loci and meta-analysis for cold tolerance in rice at booting stage. Euphytica. 2019;215(5):89.

    Article  CAS  Google Scholar 

  10. Han L, Qiao Y, Zhang S, Zhang Y, Cao G, Kim J, et al. Identification of quantitative trait loci for cold response of seedling vigor traits in rice. J Genet Genomics. 2007;34(3):239–46.

    Article  CAS  PubMed  Google Scholar 

  11. Suh JP, Lee CK, Lee JH, Kim JJ, Kim SM, Cho YC, et al. Identification of quantitative trait loci for seedling cold tolerance using RILs derived from a cross between japonica and tropical japonica rice cultivars. Euphytica. 2012;184(1):101–8.

    Article  Google Scholar 

  12. Lou Q, Chen L, Sun Z, Xing Y, Li J, Xu X, et al. A major QTL associated with cold tolerance at seedling stage in rice (Oryza sativa L.). Euphytica. 2007;158(1-2):87–94.

    Article  CAS  Google Scholar 

  13. Andaya VC, Tai TH. Fine mapping of theqCTS4locus associated with seedling cold tolerance in rice (Oryza sativaL.). Mol Breed. 2007;20(4):349–58.

    Article  CAS  Google Scholar 

  14. Koseki M, Kitazawa N, Yonebayashi S, Maehara Y, Wang Z-X, Minobe Y. Identification and fine mapping of a major quantitative trait locus originating from wild rice, controlling cold tolerance at the seedling stage. Mol Genet Genomics. 2010;284(1):45–54.

    Article  CAS  PubMed  Google Scholar 

  15. Kim S, Suh J, Lee C, Lee J, Kim Y, Jena KK. QTL mapping and development of candidate gene-derived DNA markers associated with seedling cold tolerance in rice (Oryza sativa L.). Mol Gen Genomics. 2014;289(3):333–43.

    Article  CAS  Google Scholar 

  16. Xiao N, Huang WN, Li AH, Gao Y, Li YH, Pan CH, et al. Fine mapping of the qLOP2 and qPSR2-1 loci associated with chilling stress tolerance of wild rice seedlings. Theor Appl Genet. 2015;128(1):173–85.

    Article  CAS  PubMed  Google Scholar 

  17. Yun M, Xiaoyan D, Yunyuan X, Wei L, Xiaoming Z, Dali Z, et al. COLD1 confers chilling tolerance in rice. Cell. 2015;160(6):1209–21.

    Article  CAS  Google Scholar 

  18. Zhao J, Zhang S, Dong J, Yang T, Mao X, Liu Q, et al. A novel functional gene associated with cold tolerance at the seedling stage in rice. Plant Biotechnol J. 2017;15(9):1141–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Xiao N, Gao Y, Qian HJ, Gao Q, Wu YY, Zhang DP, et al. Identification of genes related to cold tolerance and a functional allele that confers cold tolerance. Plant Physiol. 2018;177(3):1108–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Mao DH, Xin YY, Tan YJ, Hu XJ, Bai JJ, Liu ZY, et al. Natural variation in the HAN1 gene confers chilling tolerance in rice and allowed adaptation to a temperate climate. PNAS. 2019;116(9):3494–501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Zhao J, Wang S, Qin J, Sun C, Liu F. The lipid transfer protein OsLTPL159 is involved in cold tolerance at the early seedling stage in rice. Plant Biotechnol J. 2020;18(3):756–69.

    Article  CAS  PubMed  Google Scholar 

  22. Liu F, Xu W, Song Q, Tan L, Liu J, Zhu Z, et al. Microarray-assisted fine-mapping of quantitative trait loci for cold tolerance in rice. Mol Plant. 2013;6(3):757–67.

    Article  CAS  PubMed  Google Scholar 

  23. Abe A, Kosugi S, Yoshida K, Natsume S, Takagi H, Kanzaki H, et al. Genome sequencing reveals agronomically important loci in rice using MutMap. Nat Biotechnol. 2012;30(2):174–8.

    Article  CAS  PubMed  Google Scholar 

  24. Kadambari G, Vemireddy LR, Srividhya A, Nagireddy R, Jena SS, Gandikota M, et al. QTL-Seq-based genetic analysis identifies a major genomic region governing dwarfness in rice (Oryza sativa L.). Plant Cell Rep. 2018;37(4):677–87.

    Article  CAS  PubMed  Google Scholar 

  25. Lu HF, Lin T, Klein J, Wang SH, Qi JJ, Zhou Q, et al. QTL-seq identifies an early flowering QTL located near flowering locus T in cucumber. Theor Appl Genet. 2014;127:1491–9.

    Article  PubMed  Google Scholar 

  26. Song J, Li Z, Liu ZX, Guo Y, Qiu LJ. Next-generation sequencing from bulked-segregant analysis accelerates the simultaneous identification of two qualitative genes in soybean. Front Plant Sci. 2017;8:1–11.

  27. Wambugu P, Ndjiondjop MN, Furtado A, Henry R. Sequencing of bulks of segregants allows dissection of genetic control of amylose content in rice. Plant Biotechnol J. 2018;16(1):100–10.

    Article  CAS  PubMed  Google Scholar 

  28. Wen JQ, Jiang FL, Weng YQ, Sun MT, Shi XP, Zhou YZ, et al. Identification of heat-tolerance QTLs and high-temperature stress-responsive genes through conventional QTL mapping, QTL-seq and RNA-seq in tomato. BMC Plant Biol. 2019;19(1):1–17.

  29. Zhang XL, Wang WB, Guo N, Zhang YY, Bu YP, Zhao JM, et al. Combining QTL-seq and linkage mapping to fine map a wild soybean allele characteristic of greater plant height. BMC Genomics. 2018;19(1):1–12.

  30. Park M, Lee J-H, Han K, Jang S, Han J, Lim J-H, et al. A major QTL and candidate genes for capsaicinoid biosynthesis in the pericarp of Capsicum chinense revealed using QTL-seq and RNA-seq. Theor Appl Genet. 2019;132(2):515–29.

    Article  CAS  PubMed  Google Scholar 

  31. Das S, Upadhyaya HD, Bajaj D, Kujur A, Badoni S, Kumar V, et al. Deploying QTL-seq for rapid delineation of a potential candidate gene underlying major trait-associated QTL in chickpea. DNA Res. 2015;22(3):193–203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Suzaki T, Sato M, Ashikari M, Miyoshi M, Nagato Y, Hirano HY. The gene FLORAL ORGAN NUMBER1 regulates floral meristern size in rice and encodes a leucine-rich repeat receptor kinase orthologous to Arabidopsis CLAVATA1. Development. 2004;131(22):5649–57.

    Article  CAS  PubMed  Google Scholar 

  33. Lin D, Kong R, Chen L, Wang Y, Wu L, Xu J, et al. Chloroplast development at low temperature requires the pseudouridine synthase gene TCD3 in rice. Sci Rep. 2020;10(1):1–11.

    CAS  Google Scholar 

  34. Jung YJ, Lee IH, Nou IS, Lee KD, Rashotte AM, Kang KK. BrRZFP1 a Brassica rapa C3HC4-type RING zinc finger protein involved in cold, salt and dehydration stress. Plant Biol. 2013;15(2):274–83.

    Article  CAS  PubMed  Google Scholar 

  35. Tähtiharju S, Palva T. Antisense inhibition of protein phosphatase 2C accelerates cold acclimation in Arabidopsis thaliana. Plant J. 2001;26(4):461–70.

    Article  PubMed  Google Scholar 

  36. Liu C, Ou S, Mao B, Tang J, Wang W, Wang H, et al. Early selection of bZIP73 facilitated adaptation of japonica rice to cold climates. Nat Commun. 2018;9(1):1–12.

    CAS  Google Scholar 

  37. Liu C, Schläppi MR, Mao B, Wang W, Wang A, Chu C. The bZIP 73 transcription factor controls rice cold tolerance at the reproductive stage. Plant Biotechnol J. 2019;17(9):1834–49.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Liu C, Wu Y, Wang X. bZIP transcription factor RISBZ5: a potential negative regulator of cold and drought stress response in rice. Planta. 2012;235(6):1157–69.

    Article  CAS  PubMed  Google Scholar 

  39. Sasaki T, Burr B. International rice genome sequencing project: the effort to completely sequence the rice genome. Curr Opin Plant Biol. 2000;3(2):138–41.

    Article  CAS  PubMed  Google Scholar 

  40. Han LZ, Zhang YY, Qiao YL, Cao GL, Zhang SY, Jong Hwan K, et al. Genetic and QTL analysis for low-temperature vigor of germination in rice. Acta Genet Sin. 2006;33(11):998–1006.

    Article  CAS  PubMed  Google Scholar 

  41. Komatsu S, Yang GX, Khan M, Onodera H, Toki S, Yamaguchi M. Over-expression of calcium-dependent protein kinase 13 and calreticulin interacting protein 1 confers cold tolerance on rice plants. Mol Gen Genomics. 2007;277(6):713–23.

    Article  CAS  Google Scholar 

  42. Abbasi F, Onodera H, Toki S, Tanaka H, Komatsu S. OsCDPKI3, a calcium-dependent protein kinase gene from rice, is induced by cold and gibberellin in rice leaf sheath. Plant Mol Biol. 2004;55(4):541–52.

    Article  CAS  PubMed  Google Scholar 

  43. Liu G, Zhao T, You X, Jiang J, Li J, Xu X. Molecular mapping of the Cf-10 gene by combining SNP/InDel-index and linkage analysis in tomato (Solanum lycopersicum). BMC Plant Biol. 2019;19(1):1–11.

    Google Scholar 

  44. Feng X, Cheng H, Zuo D, Zhang Y, Wang Q, Liu K, et al. Fine mapping and identification of the fuzzless gene GaFzl in DPL972 (Gossypium arboreum). Theor Appl Genet. 2019;132(8):2169–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Zhang ZY, Li JJ, Pan YH, Li JL, Zhou L, Shi HL, et al. Natural variation in CTB4a enhances rice adaptation to cold habitats. Nat Commun. 2017;8(1):1–13.

    CAS  Google Scholar 

  46. Zhang Q, Chen QH, Wang SL, Hong YH, Wang ZL. Rice and cold stress: methods for its evaluation and summary of cold tolerance-related quantitative trait loci. Rice. 2014;7(1):1–12.

  47. Huang GT, Ma SL, Bai LP, Zhang L, Ma H, Jia P, et al. Signal transduction during cold, salt, and drought stresses in plants. Mol Biol Rep. 2012;39(2):969–87.

    Article  PubMed  CAS  Google Scholar 

  48. Zhang ZY, Li JJ, Pan YH, Li JL, Shi HL, Zeng YW, et al. Natural variation in CTB4a enhances rice adaptation to cold habitats. Nat Commun. 2017;8(1):1–13.

    CAS  Google Scholar 

  49. Yamaguchi-Shinozaki K, Shinozaki K. Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu Rev Plant Biol. 2006;57:781–803.

    Article  CAS  PubMed  Google Scholar 

  50. Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, et al. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013;6(1):1–10.

  51. Liu CT, Wu YB, Wang XP. bZIP transcription factor OsbZIP52/RISBZ5: a potential negative regulator of cold and drought stress response in rice. Planta. 2012;235(6):1157–69.

    Article  CAS  PubMed  Google Scholar 

  52. Murray M, Thompson WF. Rapid isolation of high molecular weight plant DNA. Nucleic Acids Res. 1980;8(19):4321–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Matsumoto T, Wu JZ, Kanamori H, Katayose Y, Fujisawa M, Namiki N, et al. The map-based sequence of the rice genome. Nature. 2005;436(7052):793–800.

    Article  CAS  Google Scholar 

  54. Kumar S, Agarwal S. Ranvijay: Fast and memory efficient approach for mapping NGS reads to a reference genome. J Bioinf Comput Biol. 2019;17(02):1950008.

  55. Wang K, Li MY, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164.

  56. Takagi H, Abe A, Yoshida K, Kosugi S, Natsume S, Mitsuoka C, et al. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J. 2013;74(1):174–83.

    Article  CAS  PubMed  Google Scholar 

  57. Magwene PM, Willis JH, Kelly JK. The statistics of bulk segregant analysis using next generation sequencing. PLoS Comp Biol. 2011;7(11):e1002255.

  58. Mansfeld BN, Grumet R. QTLseqr: an R package for bulk segregant analysis with next-generation sequencing. Plant Genome. 2018;11(2):1–5.

  59. Hill JT, Demarest BL, Bisgrove BW, Gorsi B, Su YC, Yost HJ. MMAPPR: mutation mapping analysis pipeline for pooled RNA-seq. Genome Res. 2013;23:687–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Fisher RA. On the interpretation of χ 2 from contingency tables, and the calculation of P. J R Stat Soc. 1922;85(1):87–94.

    Article  Google Scholar 

  61. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004;32(1):277–80.

    Article  CAS  Google Scholar 

  63. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Siahpoosh MR, Sanchez DH, Schlereth A, Scofield GN, Furbank RT, van Dongen JT, et al. Modification of OsSUT1 gene expression modulates the salt response of rice Oryza sativa cv. Taipei 309. Plant Sci. 2012;182:101–11.

    Article  CAS  PubMed  Google Scholar 

  65. Zhang Z, Li J, Tang Z, Sun X, Zhang H, Yu J, et al. Gnp4/LAX2, a RAWUL protein, interferes with the OsIAA3–OsARF25 interaction to regulate grain length via the auxin signaling pathway in rice. J Exp Bot. 2018;69(20):4723–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Li N, Zheng H, Cui J, Wang J, Liu H, Sun J, et al. Genome-wide association study and candidate gene analysis of alkalinity tolerance in japonica rice germplasm at the seedling stage. Rice. 2019;12(1):24.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


We thank Editage ( for English language editing.


We ensure that all rice seeds used in this study originated from Northeast Agricultural University in Heilongjiang Province, China. The legality of these seeds complies with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora. The rice seeds collected in the study are all cultivated rice in China rather than endangered and wild species. These varieties have passed the legal variety certification procedures in China and are licensed for production, planting, and market operations. Additionally, these varieties are members of the rice germplasm resources in our laboratory and are all conventional japonica rice types. The agronomic characteristics, cultivation characteristics, and varietal types of related varieties can be retrieved from the rice database (


This research was financially supported by the major science and technology project of Heilongjiang Province, China (grant no. 2020ZX16B010) and the China Postdoctoral Science Foundation (grant no. 2019M651249).

Author information

Authors and Affiliations



L.Y. and D.Z. conceived and designed the research; L.Y., J.W., and L.L. participated in data analysis; J.W., H.Z., W.X., and L.L. performed material development, sample preparation, and data analysis; L.Y. wrote the manuscript; H.Z and D.Z. corrected the manuscript; and the final manuscript was read and approved by all of the authors.

Corresponding author

Correspondence to Detang Zou.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1

. Clustering map of GO annotation of genes in 96.6- Kb regions. The abscissa is the content of GO categories, and the left of the ordinate is the number of genes. This figure shows the gene classification of GO secondary functions in the context of all genes in the associated region. Figure S2. Expression levels of the 13 candidate genes in DN430 and DF104 under normal condition measured by qRT-PCR. The results were statistically analyzed using Student’s t-test (**, P<0.01). Figure S3. The CDS region Sequence difference analysis of OsbZIP54. The gene structure of Os07g0569700 and sequence differences in OsbZIP54 between DF-104 and DN-430. Ref is the reference sequence of Nipponbare genome. Figure S4. The promoter Sequence difference analysis of OsbZIP54. The promoter structure of OsbZIP54 and sequence differences in OsbZIP54 between DF-104 and DN-430. Ref is the reference sequence of Nipponbare genome. Supplemental Table 1. Number of single nucleotide polymorphisms (SNPs) and InDels detected in samples. Supplemental Table 2. Survival rate of parents and F2:3 population under cold treatment. Supplemental Table 3. Significant peak statistics driven by four algorithms. Supplemental Table 4. Statistics of variation in the 2.60 Mb interval of qCTS6 based on resequencing data. Supplemental Table 5. GO annotation result of 96.6-Kb interval. Supplemental Table 6. Statistics of variation in the 96.6Kb interval based on resequencing data. Supplemental Table 7. The Haplotype analysis of Os06g0719500 in the T-Pool and S-Pool. Supplemental Table 8. Base Variation Statistics of 295 North China Japonica Rice in qCTS6 Interval. Supplemental Table 9. Primers used in this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yang, L., Wang, J., Han, Z. et al. Combining QTL-seq and linkage mapping to fine map a candidate gene in qCTS6 for cold tolerance at the seedling stage in rice. BMC Plant Biol 21, 278 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: